{'fraction_non_alphanumeric': 0.048510616928497896, 'fraction_numerical': 0.018363059617171218, 'mean_word_length': 3.963855421686747, 'pattern_counts': {'":': 0, '<': 29, '<?xml version=': 0, '>': 23, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 27, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Periodically, Quasi-periodically, and Randomly Driven Conformal Field Theories: Part I

In this paper and its sequel, we study non-equilibrium dynamics in driven 1+1D conformal field theories (CFTs) with periodic, quasi-periodic, and random driving. We study a soluble family of drives in which the Hamiltonian only involves the energy-momentum density spatially modulated at a single wavelength. The resulting time evolution is then captured by a M\"obius coordinate transformation. In this Part I, we establish the general framework and focus on the first two classes. In periodically driven CFTs, we generalize earlier work and study the generic features of entanglement/energy evolution in different phases, i.e. the heating, non-heating phases and the phase transition between them. In quasi-periodically driven CFTs, we mainly focus on the case of driving with a Fibonacci sequence. We find that (i) the non-heating phases form a Cantor set of measure zero; (ii) in the heating phase, the Lyapunov exponents (which characterize the growth rate of the entanglement entropy and energy) exhibit self-similarity, and can be arbitrarily small; (iii) the heating phase exhibits periodicity in the location of spatial structures at the Fibonacci times; (iv) one can find exactly the non-heating fixed point, where the entanglement entropy/energy oscillate at the Fibonacci numbers, but grow logarithmically/polynomially at the non-Fibonacci numbers; (v) for certain choices of driving Hamiltonians, the non-heating phases of the Fibonacci driving CFT can be mapped to the energy spectrum of electrons propagating in a Fibonacci quasi-crystal. In addition, another quasi-periodically driven CFT with an Aubry-Andr\'e like sequence is also studied. We compare the CFT results to lattice calculations and find remarkable agreement.


Introduction
Non-equilibrium dynamics in time-dependent driven quantum many-body systems has received extensive recent attention. A time-dependent drive, such as a periodic drive, creates a new stage in the search for novel systems that may not have an equilibrium analog, e.g., Floquet topological phases [1][2][3][4][5][6][7][8][9][10][11][12][13] and time crystals [14][15][16][17][18][19][20][21][22]. It is also one of the basic protocols to study non-equilibrium phenomena, such as localization-thermalization transitions, prethermalization, dynamical Casimir effect, etc [23][24][25][26][27][28][29][30][31]. Despite the rich phenomena and applications in the time-dependent driving physics, exactly solvable setups are, in general, very rare. Usually, we have to resort to numerical methods limited to small system size. In this work, we are interested in a quantum (1 + 1) dimensional conformal field theory (CFT), which may be viewed as the low energy effective field theory of a many-body system at the critical point. The property of conformal invariance at the critical point can be exploited to constrain the operator content of the critical theory [32,33]. In particular, for (1 + 1)D CFTs, the conformal symmetry is enlarged to the full Virasoro symmetry, which makes tractable the study of non-equilibrium dynamics, such as the quantum quench problems [34,35]. For a time-dependent driven CFT, however, relatively little is known.
Most recently, an analytically solvable setup on the periodically driven CFT was proposed in Ref. [36]. The authors implement the periodic driving with two non-commuting Hamiltonians H 0 and H 1 for time durations T 0 and T 1 respectively, where H 0 = L 0 h(x) dx is the uniform CFT Hamiltonian on a line of length L, and H 1 is obtained from H 0 by deforming the Hamiltonian density h(x) as H 1 (x) = L 0 2 sin 2 πx L h(x)dx, which is also called sine-square deformation (SSD) in literature [37][38][39][40][41][42][43][44][45][46][47][48][49][50]. Interestingly, it was found that different phases can emerge during the driving, depending on duration of the two time evolutions. As depicted in Fig. 1, there exits a heating phase with the entanglement entropy growing linearly in time, and a non-heating phase with the entanglement entropy simply oscillating in time. At the phase transition, the entanglement entropy grows logarithmically in time. Later in Ref. [51], these emergent phases and the phase diagram were further confirmed by studying how the system absorbs energy. More explicitly, the total energy of the system grows exponentially in time in the heating  phase, oscillates in the non-heating phase, and grows polynomially at the phase transition. Furthermore, the system develops interesting spatial structures in the heating phase. The energy density forms an array of peaks 1 with simple patterns of entanglement as shown in Fig. 2.
In this work and its sequel, we introduce and study a general class of soluble models of driven CFTs with a variety of driving protocols. We determine their dynamical phase diagrams of heating versus non-heating behavior, particularly when the periodicity of the drive is absent. We extend the previous study on periodic driving to quasi-periodic 2 and random drivings, and make a connection to the familiar concepts of crystal, quasi-crystal, and disordered systems. The connection is based on the coincidence of the group structures underlying the two problems: 1. The driving protocol we considered for the CFTs involves SL 2 deformed Hamiltonians.
These are generalizations of the SSD Hamiltonian protocols, where the deformed Hamiltonians H q are chosen as Here h(x) and p(x) are the energy and momentum densities, and f q (x) (g q (x)) are real functions of the form a + b cos 2πqx L + c sin 2πqx L , with q ∈ Z. The remarkable aspect of these protocols, which is the key to their solubility, is that the time evolution of many physical quantities after a prescribed time is captured simply by a 2 × 2 matrix transformation, i.e. a SL 2 or Möbius transformation. This simplification occurs despite the fact that we are discussing Figure 3: A local view of the operator evolution on a Riemann surface. By choosing a suitable coordinates, each step of the driving can be characterized by a Möbius transformation that is determined by the SL 2 deformed Hamiltonian. a spatially extended system. In this case, the operator evolution can be recast into a sequence of Möbius transformations on a suitable Riemann surface (See Fig. 3), 2. The hopping problem in tight binding model can be solved using transfer matrix method, namely reformulating the discrete Schördinger equation Eψ j = [Hψ] j = ψ j+1 +ψ j−1 +V j ψ j by product of the transfer matrices where Ψ n = (ψ n+1 , ψ n ) T represents the corresponding wave-function.
Both problems are now solved by analyzing products of SU(1, 1) SL(2, R) matrices, creating intriguing analogies. In fact, the main part of the paper is to dive into the analogies and exam whether the rich phenomenons in solids can reassemble in the time domain.

Outline and main results of this paper
• In Sec. 2 we explain the details of the general setup of our study, which is a time-dependent driven (1 + 1) CFT with arbitrary SL 2 deformations. As mentioned in the introduction, the physical consequence of such driving is encoded in the product of a sequence of SU(1, 1) matrices that correspond to the driving steps.
• In Sec. 3, we introduce the main diagnostics of our driven CFT: the Lyapunov exponent and group walking. The former is a useful characterization to quantify the growth of Π n w.r.t. the number of driving step n, i.e.
where · is a matrix norm. Applying to our driven system, the Lyapunov exponent has the meaning of the heating rate and serves as a good "order parameter" in carving the phase diagram. For example, λ L > 0 represents a heating phase, and we show that total energy of the system grows exponentially as E(n) ∝ e 2λ L ·n and the entanglement entropy of the subsystem that includes the energy-momentum density peaks grows linearly in time as S(n) ∝ λ L · n. One interesting universal phenomenon here is that the total energy and the entanglement are not distributed evenly in the system, instead the driven state will develop an array of peaks of energy-momentum density in the real space. This phenomenon has been reported in Ref. [51] for special setups, now we verify the universality in a larger class.
For λ L = 0, the system is either in the non-heating phase where total energy and the entanglement entropy oscillate or at the phase transition where the total energy grows polynomially and the entanglement grows logarithmically.
The second diagnostic we introduce is the notion of group walking, which is particularly useful in analyzing and visualizing the details of the spatial structures. This tool is necessary in the cases such as quasi-periodic and random driving when we need to resort to the numerics to identify the universal features.
• In Sec. 4, we study the properties of the periodic driving, providing criteria of the heating phase, non-heating phase, and the phase transition. We discuss the generic features in each phase. This section generalizes the minimal setup in Ref. [36,51], and also provides the necessary tools for the discussions in quasi-periodic driving where technically we approach the quasi-periodic limit via a family of periodic driving.
• In Sec. 5, we consider the quasi-periodic driving using two examples: Fibonacci type and Aubry-André type. The Fibonacci driving is the main focus. In the Fibonacci driving, we use the Fibonacci bitstring/word (see Appendix. B) X j=1,2,3... = 10110101 . . . and two distinct unitary operators U A = e −iH A T A , U B = e −iH B T B to generate a quasi-periodic driving sequence U j = X j U A + (1 − X j )U B . The simplest way to generate the Fibonacci bitstring is through the following substitution rule: Begin with a single bit 1, and apply the substitution rule 1 → 10, 0 → 1 at each step, then we will generate the following sequence 1 → 10 → 101 → 10110 → 10110101 → ..., which approaches the Fibonacci bitstring in the infinite step limit. Denoting the n-th Fibonacci number as F n , namely F n = F n−1 + F n−2 with F 0 = F 1 = 1, the Fibonacci bitstring/word satisfies: X j+Fn = X j , where n 2 and 1 j < F n . In the Fibonacci driving, we find the following features: 1. In the heating phase, the Lyapunov exponents (heating rates) exhibit self-similarity. This also implies there exist heating phases with arbitrarily small positive Lyapunov exponents. At these points, the growth of entanglement entropy/energy can be arbitrarily slow. In addition, there are very rich patterns in the time evolution of entanglement/energy in the heating phase. In particular, the locations of the energymomentum density peaks exhibit even/odd effects at the Fibonacci numbers.
2. Exact non-heating fixed points. We find that there always exist exact non-heating fixed points in the phase diagram, as long as both of the two driving Hamiltonians are elliptic. At the non-heating fixed point, the time evolution of the entanglement entropy and the total energy can be analytically obtained at the Fibonacci numbers F n . They exhibit an oscillating feature of period 6, i.e., S A (F n ) = S A (F n+6 ) and E(F n ) = E(F n+6 ). At the driving steps that are not Fibonacci numbers, the envelope of the entanglement entropy grows logarithmically in time, and the total energy grows in a power law.
3. We find an exact mapping between the phase diagram of a Fibonacci driving CFT and the energy spectrum of a Fibonacci quasi-crystal. More precisely, the nonheating phase in the parameter space of a Fibonacci driving CFT corresponds to the energy spectrum of a Fibonacci quasi-crystal. Both form a Cantor set of measure zero.
As a complement, we also investigate the quasi-periodic driving with an Aubry-André like sequence, where the phase diagram has a nested structure that resembles the famous Hofstadter butterfly found in the Landau level problem [56]. We also exam the measure of the non-heating phase and show it vanishes similar to the Fibonacci driving.
• In Sec. 6 we conclude with discussions. We also provide several appendices with details of calculations and examples.
2 Time-dependent driven CFT with SL 2 deformations In this section, we introduce the general setup and basic properties of a time-dependent driven CFT with SL 2 deformations. The formalism in this section is general, i.e. suitable for arbitrary driving sequence. In the end of this section, we will explain the three classes of driving that we will focus on in this paper and its sequel [57]: the periodically, quasi-periodically, and randomly driven CFTs as advertised in the introduction. More technical details can be found in Appendix A (See also Refs. [43,36,51]). We are mainly interested in the time-dependent driven CFT with discrete time steps. That is, we drive the CFT with H 1 for a time interval T 1 , then with H 2 for a time interval T 2 , and so on, where H 1,2,... are SL 2 deformed CFT Hamiltonians that we will explain momentarily. Starting from an initial state |Ψ 0 , the wavefunction after n steps of driving has the form: The initial state here is not limited to a ground state. For instance, it can be chosen as a highly excited pure state or a thermal ensemble at finite temperature, as will be studied in detail in Ref. [58]. It is found that the emergent phase diagram of the time-dependent driven CFT is independent of the choices of the initial state, and only depends on the concrete protocols of driving, namely the driving sequences {U j } here. For simplicity, throughout this work we will choose the initial state |Ψ 0 as the ground state of a "uniform CFT", i.e. with uniform Hamiltonian density where T (x) (T (x)) are the chiral (anti-chrial) energy-momentum tensor with translation symmetry, L is the total length of the system. Now let us specify the choices of the Hamiltonians {H j } in Eq. (5), we require them to be generated by a deformation {(f j , g j )} as follows: where f j (x) and g j (x) are two independent real functions with periodic boundary conditions. 3 That is to say, in general we can deform the chiral and anti-chiral modes independently in a system with periodic boundary conditions. An alternative way to view the deformation in (7) is to rewrite (7) using energy density T 00 (x) = 1 2π (T (x) + T (x)) and the momentum density T 01 (x) = 1 2π (T (x) − T (x)) as follows Although the formulas and results we obtain in the this work hold for the general case, in many places of this paper we will choose f j (x) = g j (x) such that the deformed Hamiltonian takes the following simple form The study of the energy spectrum of such Hamiltonian can be found in [44]. In particular, the so-called sine-square deformation (SSD) with f j (x) = sin 2 ( πx L ) in Eq. (9) has received extensive study in both condensed matter physics and string theory recently [37][38][39][40][41][42][43][44][45][46][47][48][49][50]. In fact, the initial study of the Floquet CFT in Refs. [36,51] is also based on SSD.

SL 2 deformation
A convenient parametrization of the deformed Hamiltonian H j in Eq. (7) is to use the Fourier components of T (x) and T (x) denoted as L n and L n L n := c 24 The operators L n (L n ) form a Virasoro algebra with c being the central charge of the underlying CFT. For example, the uniform Hamiltonian H 0 defined in (6) can be expressed as and what we will call a "SL 2 deformed" Hamiltonian corresponds to the following enveloping function and similarly for g j (x). In this case, the corresponding H j is a linear superposition of {L 0 , L ±q } and {L 0 , L ±q }, which are the generators of the SL (q) (2, R) subgroup. 4 To be concrete, we have where we have defined L q,+ := 1 2 (L q +L −q ) and L q,− := 1 2i (L q −L −q ). Note the SSD deformation mentioned above corresponds to the special case when Therefore, the SL 2 deformation can be thought as a generalization of the SSD deformation, while retaining the analytic tractability. 5 In general, by defining the quadratic Casimir the (chiral and anti-chiral) SL 2 deformed Hamiltonians can be classified into three types as follows [40,41,60,61] Different types of Hamiltonians will determine the operator evolution in different ways (See Appendix.A.1). With the SL 2 deformation, many physical properties of the driven system including the phase diagram, the time dependence of the entanglement entropy [43,36] and the energy-momentum density [51,52] have been obtained in a periodically driven CFT system. Here we generalize the driving to an arbitrary sequence {U j } as shown in Eq. (5). In the following, we will derive general formulas based on the SL 2 deformation sequence, and later apply to the periodically, quasi-periodically, and randomly driven CFTs.

Operator evolution
For SL 2 driven quantum states, it is convenient to compute observables via Heisenberg picture, namely the correlation functions are given by Ψ 0 |O 1 (x 1 , t 1 ) · · · O n (x n , t n )|Ψ 0 , where the Heisenberg operators O(x, t) are defined by discrete time evolution from the w-cylinder (where w = τ + ix, and x = L and x = 0 are identified) to the q-sheet Riemann surface z.
For each step, U j = e −iH j T j is generated by the SL 2 deformed Hamiltonian and t = m j=1 T j is only defined for a discrete set of times in our setting.
The virtue of the driving Hamiltonian in (7) is that the operator evolution can be represented by a conformal mapping (z, z) → (z , z ), under which the primary operator O(z, z) transform as where h (h) are conformal dimensions of O. Then the full unitary U (t) is a composition of a sequence of conformal mappings. For the special type of enveloping function (13), a convenient coordinate is given as follows (see Fig. 4 for an illustration.) under which the evolution generated by U j can be expressed as a Möbius transformation: 6 The explicit form of M j is determined by the Hamiltonian H j and the time interval T j . An important observation is that the driving protocol (13) we use in fact generates M j matrix in the following specific form which is a SU(1, 1) matrix. Note SU(1, 1) ∼ = SL(2, R), both are subgroups of SL(2, C). The isomorphism is expected since we start from a SL(2, R) action on the states. Thus, the net effect of the full evolution U = U n · · · U 2 · U 1 is given by the product of n SU(1, 1) matrices 6 More details can be found in Appendix. A. See also Refs. [43,36,51].
Note that the later matrix acts on the right since we are using the Heisenberg picture of evolution. To summarize, the operator evolution under a sequence of driving {U j } is given by the following formula where z n is related to z through the Möbius transformation in (21) with the matrix Π n ∈ SU(1, 1).

Time evolution of entanglement and energy-momentum density
To characterize the possible emergent phases, we study the time evolution of the entanglement entropy and the energy-momentum density of the system. In terms of correlation functions, the former is determined by the two point function of twist operator, while the latter is determined by the one point function of energy-momentum tensor.
For example, the m-th Renyi entropy of the subsystem A = [x 1 , x 2 ] can be obtained by the following formula S (m) where |Ψ n is the time-dependent wavefunction in Eq. (5), and T m (T m ) are twist (anti-twist) operators that are primary, with conformal dimensions h = h = c 24 (m − 1 m ). For initial state |Ψ 0 being the ground state of H 0 with periodic boundary conditions, the time evolution of the entanglement entropy for the subsystem A = [(k − 1/2)l, (k + 1/2)l] where k ∈ Z and l = L/q is given as 7 S Here α n and β n are the matrix elements appearing in the operator evolution in Eq. (25). α n and β n are the corresponding matrix elements for the anti-chiral part. One can also study the time evolution of energy-momentum tensors based on the operator evolution as discussed in the previous subsection. However, since T (x) is not a primary field, the operator evolution in Eq. (19) should be modified as where the last term represents the Schwarzian derivative. The expectation value of the chiral energy-momentum tensor density is [51] 7 For a general choice of subsystem A = [x 1 , x 2 ], the exact expression of the entanglement entropy under a time-dependent driving will be quite involved. See, e.g., the appendix of Ref. [36]. However, if the CFT is in a heating phase, one can obtain an approximated expression of the entanglement entropy of A = [x 1 , x 2 ] by keeping the leading order [51]. 8 Hereafter, for convenience of writing, we write T (x, t = n j=1 T j ) as T (x, n).
For the anti-chiral component 1 2π T (x, n) , the expression is the same as above by replacing α n (β n ) → α n (β n ) and e 2πix l → e − 2πix l . The total energy and momentum of the system are We would like to make a few remarks here: 1. For the periodic boundary conditions we considered here, the time evolution with q = 1 deformations are trivial as H j also annihilates the ground state of H 0 . 9 In contrast, if we consider an open boundary condition, the ground state of H 0 will no longer be the eigenstate of the deformed Hamiltonian H i . Then one can have a non-trivial time evolution, as studied in Refs. [36,51].
3. If we only deform the Hamiltonian density in Eq. (8), i.e. let f j (x) = g j (x), one can find that α n = α n and β n = β n [36], and therefore P (n) = 0, i.e., the total momentum stays zero. In this case, both the left movers and right movers are excited, but they carry opposite momentum and the total momentum are canceled to be zero. 4. In later sections, we will compare the CFT calculations with the lattice calculations. An efficient way to perform numerical calculations on the lattice is to consider q = 1 with an open boundary condition, since for larger q, the length of the wavelength of deformation l := L/q is effectively suppressed for a fixed L. In this case, by deforming the Hamiltonian in Eq. (9), where only the Hamiltonian density is deformed, one can find the time evolution of the entanglement entropy as follows: [36] The expectation value of the chiral energy-momentum density is: The anti-chiral part has the same expression as above with the replacing e 2πix L → e −2πix L . Then one can find 1 , based on which one can obtain the total energy as One can find the similarity and difference by comparing with the case with periodic boundary conditions. For example, in the ground state with α n = 1, β n = 0 in (33), one can obtain the Casimir energy E = − πc 24L , which is different from that in periodic boundary conditions. Nevertheless, the dependence of the entanglement entropy/energymomentum denstiy on the matrix elements α n (β n ) in the Π n in Eq. (25) are similar.
There is rich information contained in the formula discussed above. As will be seen later, if the CFT is in a heating phase, there will be energy-momentum density peaks emerging in the real space. The locations of these peaks are determined by β n /α n . It turns out that both the quantities S A (n) and E(n) can serve as 'order parameters' to distinguish different emergent phases in the time-dependent driven CFTs. For example, for the periodically driven CFT as studied in Ref. [36], it is found there are two different phases with a heating phase and a non-heating phase, where the time evolution of entanglement entropy exhibits qualitatively different features as shown in Fig. 1. Also, it is found n Ref. [51] that the total energy grows exponentially fast as a function of driving cycles n in the heating phase and simply oscillates in the non-heating phase.
As a short summary, once we know the operator evolution in Eq. (24) or equivalently the matrix form in Eq. (25), one can study the entanglement/energy-momentum evolution based on Eqs. (27) and (30) or Eqs. (31) and (33).

Periodic, quasi-periodic, and random driving CFTs
In general, the sequence of unitary operators {U j } in Eq. (5) can be chosen in an arbitrary form. In this work, we are interested in three classes: periodic, quasi-periodic, and random drivings.
1. Periodical driving: the sequence of unitary operators {U j } in (5) are chosen with a 'period' p (p ∈ Z + ) such that U j = U j+p , ∀j ∈ Z. Then the time evolution of wavefunction in (5) can be written as To obtain the physical properties of the system under periodic driving, we only need to analyze the corresponding transformation matrix M 1 M 2 · · · M p ∈ SU(1, 1) within a period.
2. Quasi-periodic driving: {U j } form a quasi-periodic sequence. Quasiperiodicity is the property of a system that displays irregular periodicity, where the sequence exhibits recurrence with a component of unpredictability (For example see the review [62] for a more rigorous mathematical definition of quasi-periodic sequence). In this paper, we will focus on the following two protocols of quasi-periodical driving: (a) Fibonacci type. This is the type of quasi-periodic driving we will study in detail in Sec. 5.1. We use the Fibonacci bitstring/word (see Appendix. B) X j=1,2,3... = 10110101 . . .
and two distinct unitaries U A , U B to generate a quasi-periodic driving sequence U j = X j U A + (1 − X j )U B , i.e. we apply U A (U B ) if the bit is 1 (0).
(b) Aubry-André type. In this case, we generate the quasi-periodic driving sequence as follows That is to say, we consider two Hamiltonians H 0 and H 1 and fix the driving period T 1 for H 1 while let the driving period of H 0 increase with driving cycle T 0 = nωL where ω is an irrational number and L is the total length of the system. Note in terms of the unitary U = exp(−iH 0 T 0 ), its action on the operator only depends on T 0 mod L.

3.
Random driving: {U j } form a random sequence. More concretely, each U j is drawn independently from the ensemble {(u k , p k )} k=1...m , where u k = e −iH k T k is the unitary matrix and p k is the corresponding probability, with the normalization k p k = 1.
In brief, for all the three kinds of time-dependent drivings, our goal is to describe the behavior of the physical properties of the CFT in the long time driving limit n → ∞, where n is the number of driving cycles.
As a remark, one can find that the types of driving sequence are similar to those in the potentials in crystals, quasi-crystals, and disordered systems. One can find interesting relations between different phases of time-dependent driven CFTs and different types of wavefunctions in a lattice, as briefly discussed in the introduction. Furthermore, both types of quasi-periodicities we mentioned above have been discussed in the quasi-crystal literature, e.g. see Refs. [63][64][65].

Diagnostics
The previous section explains how the physical properties of an SL 2 driven CFT state can be extracted from the conformal mapping generated by the driving sequence. The mapping is further encoded in an SU(1, 1) matrix, denoted as Π n in (25), which itself is a product of n SU(1, 1) matrices.
Mathematically, the long time asymptotics of the driven state now can be understood by the n-dependence of Π n . In this section, we will introduce two useful diagnostics to characterize such dependence: (1) Lyapunov exponent; (2) Group walking. The former is a simple scalar quantifying the growth of Π n , while the latter is more refined and uses two points on the unit disk to track the trajectory of Π n . Although not independent, both of them will be useful and used in the later sections.

Lyapunov exponent and heating phase
For all the three classes of drivings we introduced in the previous subsection, the problem is reduced to the study of the product Π n (defined in (23)) of a sequence of SU(1, 1) matrices that encode the conformal mappings. One useful and simple characterization for the growth rate of this matrix product is the so-called Lyapunov exponent (for a review of the subject, see e.g. [66]). Generally, we can consider a product of n matrices Π n = M 1 · M 2 · · · M n , where M j ∈ SL(d, R). Then the (upper) Lyapunov exponent is defined as where · is a matrix norm. We would like to make a few comments about the definition here: 1. Here, the specific choice of norm · · · is not essential. To be explicit, we will choose the Frobenius norm in this paper, i.e.
2. The definition also applies to SL(d, C), as one can always embed SL(d, C) in SL(2d, R).
Following Eqs. (27) and (30), we find the asymptotics in long time limit where we have neglected the contribution from the anti-chiral mode for the moment. 10 From this perspective, we may interpret the Lyapunov exponent λ L as the heating rate in the heating phase. If λ L > 0, then the time-dependent driven CFT must be in a heating phase, with the total energy exponentially growing in time. We will also explain in detail momentarily that since the norm of the ratio β n /α n approaches 1 in the long time limit when λ L > 0, an array of energy-momentum peaks will emerge in real space whose exact locations will be determined by the phase of the ratio β n /α n . On the other hand, if λ L = 0, the system is either in a non-heating phase or at the phase transition. We emphasize here that the vanishing of Lyapunov exponential allows a sub-exponential growth of the matrix norm ||Π n || as a function of n, e.g. this could happen at the phase transition/boundary.
To summarize, using the Lyapunov exponent λ L , we can classify the phases as follows: λ L > 0 : heating phase (with exponentially growing energy), λ L = 0 : non-heating phase or phase transition.
To further identify the detailed properties of entanglement/energy evolution in the non-heating phase and at the phase transition, one needs to study the finer structure of the matrix Π n , which we will pursue in the next subsection. The Lyapunov exponent works for general matrices. When specialized to SU(1, 1) or SL(2, C), another commonly used classifier is the trace of the matrix. Namely, | Tr M | > 2, | Tr M | = 2 and | Tr M | < 2 correspond to the hyperbolic, parabolic and elliptic types of matrix, respectively. This criterion was used to identity different phases for the Floquet driving CFT studied in previous works [36,51].
For periodic driving, it follows from the definition of matrix norm that this trace classifier is equivalent to the Lyapunov exponent. We have | Tr M | > 2 if and only if λ L > 0, | Tr M | < 2 if and only if λ L = 0. One can also extend it to the quasi-periodic driving as follows. As will be detailed discussed later, any quasi-periodic driving corresponding to an irrational number w can be considered as the limit of a sequence of periodic driving, which is generated by the continued fractions of ω. For each element in the sequence, we can apply the trace classifier to obtain a sequence of phase diagram, with its limit being the true phase diagram for the quasi-periodic driving system.
For the random driving, the Lyapunov exponent will become a more appropriate definition, which we use exclusively in the corresponding discussion.

Group walking: fine structures of the time-dependent driving
The Lyapunov exponent defined in the previous subsection is a single number. To view the 'internal structure' in the matrix product Π n in Eq. (25), it is helpful to study how the matrix elements evolve in time, which determines the time evolution of the entanglement entropy and the energy-momentum density.
A convenient parametrization of the SU(1, 1) matrices such as M j and Π n in (23) is given as follows and N ρ = 1 − |ρ| 2 is the normalization factor. 11 The unit disk D := {z ∈ C, |z| < 1}, the boundary(or edge) of the disk ∂D := {z ∈ C, |z| = 1}, and the complex numbers ρ and ζ are depicted as follows: Thus, the evolution of matrix Π n as a function of step n can be captured by the evolution of a pair of points (ρ n , ζ n ) on the unit disk. We will call this process 'group walking' for brevity. An equivalent but more convenient parameterization of the trajectory (ρ n , ζ n ) is to use (ρ n , ρ n ζ n ). For example, the total energy (30), locations of the energy-momentum density peaks (29) and entanglement entropy (27) are expressible using (ρ n , ρ n ζ n ): 1. The energy formula in Eq. (30) only depends on ρ n and increase monotonically w.r.t. |ρ n |.
In the heating phase, the exponential growth of E(n) as a function of n is tied to the phenomenon that |ρ n | approaches exponentially close to the boundary ∂D. On the other hand, if the total energy simply oscillates in n, e.g., in the non-heating phase of a periodically driven CFT, then |ρ n | should follow the same oscillation pattern. In between, as a we approach the phase transition, the orbit of ρ should be closer and closer to the boundary. As a summary, the behaviors of ρ described above can be visualized by the following cartoon: Non-heating phase near phase transition (46) Here we sketch the rough features of the group walking of ρ n , this cartoon are not meant to be exact. As we will see in Sec. 5 on the quasi-periodically driven CFT, in general there are many rich fine structures on the orbit of ρ.
2. The locations of the energy-momentum density peaks are determined by (ρ n ζ n ) only.
Recall that the poles of the (chiral) energy-momentum density (29) the former is indeed isomorphic to the Z 2 quotient of the latter, namely SU(1, 1)/Z 2 and agrees with our parametrization.
3. The entanglement entropy in Eq. (27) depends on both ρ n and ρ n ζ n where we have neglected the contribution of the anti-chiral mode.
In summary, the group walking of ρ n and (ρ n ζ n ) on a unit disk D determine the behaviors of the energy-momentum/entanglement evolution as follows: 1. |ρ n | determines the growth of total energy.
2. In the heating phase, (ρ n ζ n ) in the long driving limit (n 1) determines the location of peaks of the energy-momentum density.
3. |ρ n | and (ρ n ζ n ) together determine the time evolution of the entanglement entropy S A (n).

Periodic driving
This section is a generalization of previous works in Refs. [36,51] by considering a more general setup of periodic drivings. Apart from its own interesting features, this generalized setup can be used to analyze the quasi-periodically driven CFTs in Sec. 5.
In Refs. [36,51], a minimal setup of a periodic driving with two driving steps within one driving period was considered. The two different driving Hamiltonians are chosen as H 0 and H 1 = H SSD being the sine-square deformed Hamiltonian. Here we generalize this minimal setup in two aspects: one is to consider arbitrary SL 2 deformed Hamiltonians, and the other is to consider more general periodic sequences. In general, as the number of driving steps within a driving period increases, the phase diagram will become quite rich. 12

General protocol for periodic driving
For a periodical driving with period p ∈ Z + , we have U j = U j+p for all j ∈ Z + . Then the time evolution of wavefunction after np driving steps is determined by the unitary operators (U p · · · U 2 · U 1 ) as follows, Figure 5: A general protocol for a periodically driven CFT. There are p steps of driving within each driving period. In the i-th step of driving, we consider the driving with (H i , T i ), where H i is a SL 2 deformed Hamiltonian in Eq. (7) and T i is the corresponding time interval of driving.
In terms of conformal mapping, the operator evolution after np driving steps only depends on the the matrix product Let us denote the matrix elements of Π p and (Π p ) n as follows Next, we will determine the phase diagrams and relevant physical quantities based on the operator evolution given by Π p , or equivalently the following Möbius transformation and similarly for z.

Phase diagram and Lyapunov exponents
The matrix Π np = (Π p ) n has three distinct asymptotics depending on the trace of Π p , for convenience, let us classify the types of SU(1, 1) matrices in parallel to the classification of Möbius transformation we used in Ref. [51], also see Fig Figure 6: Illustration for the locations of fixed points of Möbius transformation in the three phases. In the non-heating phase, the two fixed points are inside and outside the unit circle respectively. They will merge at the same point on the unit circle at the phase transition. Then the two fixed points will split but still sit on the unit circle in the heating phase.
As a reminder, the fixed points of the Möbius transformation are convenient way to characterize the transformation when we repeat it multiple times. 13 We rewrite the Möbius transformation into the following form where γ 1,2 are the fixed points we mentioned, and η is the multiplier. For Π p parametrized in (52), we have the following explicit formulas Note the sign of the discriminant depends on the trace of Π p , which can be used to categorize the Möbius transform (54) (see Fig. 6 for an illustration. ) For parabolic class when | Tr(Π p )| = 2, we have γ 1 = γ 2 , the transformation (54) becomes trivial and we need to invoke When repeating n times, we only need to modify η → η n for | Tr(Π p )| = 2 case and β * p → nβ * p for | Tr(Π p )| = 2. And therefore, we have a simple expression for the matrix elements of (Π p ) n defined in (52): 13 Therefore, this is the main tool we used in the previous study [51] to visualize the effects of periodic driving.
Another advantage of the representation using fixed points γ 1,2 and multiplier η is that the Lyapunov exponent now only depends on η as follows That is to say, the hyperbolic Π p with | Tr(Π p )| > 2 implies a positive Lyapunov exponent and therefore heating phase; while the elliptic and parabolic classes both have λ L = 0. By analyzing the corresponding group walk in the next subsection we will confirm that | Tr(Π p )| < 2 corresponds to non-heating phase while | Tr(Π p )| = 2 is the phase transition as expected.
2. For | Tr(Π p )| > 2, the multiplier |η| = 1 and we have the following limit at n → ∞, Recall that both γ 1 and γ 2 live on ∂D as shown in Fig. 6. Therefore, in this case, both ρ n and (ρ np ζ np ) will approach exponentially close to the boundary of the unit disk ∂D.
3. For | Tr(Π p )| = 2, we have the following limit where γ * ∈ D as shown in Fig. 6. From Eq. (62), we notice that in this case, both ρ np and (ρ np ζ np ) will approach ∂D polynomially (in n) close.

Entanglement/energy evolution
Given the explicit expressions of the matrix elements of Π np = (Π p ) n in (58) and (59), we can further obtain the time evolution of the entanglement entropy S A (N = np) and the total energy E(N = np) based on formulas (27) and (30), respectively. For the total energy, it grows exponentially in the heating phase and the exponent is exactly twice the Lyapunov exponent, the latter is given in (60). In the non-heating phase and the phase transition, the total energy oscillates and grows polynomially.
Here we only consider the contribution of the chiral modes, the anti-chiral modes follow parallel discussions.
As noted in Ref. [51], the energy-momentum density has interesting spatial structures. In fact, as mentioned in Sec. 3.1, a positive Lyapunov exponent λ L indicates there is an array of peaks in the energy-momentum density T (x, n) in real space. The same spatial structure, namely the array of peaks, is also present at the phase transition with λ L = 0, although the growth is polynomial in n, significantly slower than the heating phase.
Following Eqs. (29), (58) and (59), we find the locations of the (chiral) energy-momentum peaks are given as follows where we have assumed 0 < η < 1 in the above formula, for η > 1, we need to replace γ 2 by γ 1 . Here γ 2 corresponds to the unstable fixed point in the Möbius transformation in the heating phase, and γ is the unique fixed point at the phase transition. A cartoon plot of the energy-momentum density distribution in real space is shown as follows: where different colors represent different chiralities. For simplicity, let us keep the anti-chiral part (red) undeformed, and only deform the chiral part (blue). Then the entanglement entropy in the heating phase depends on the choice of subsystem A as follows [51]: If one also deforms the anti-chiral part and let it live in heating phase with Layapunov exponent λ L , then we need to add up two contributions when A also includes any anti-chiral peaks. Note in general λ L = λ L as they can be deformed independently in the CFT with periodic  boundary condition. One can further check the entanglement pattern by looking into the mutual information as studied in [51], and find each peak is mainly entangled with the two peaks of its nearest neighbor with the same chirality, as schematically shown in Fig. 2. At the phase transition, similar to the energy, the spatial structure of the entanglement persists, while the growth is slower One final remark is that in the above discussions, the entanglement cuts are chosen to avoid the centers of the energy-momentum density peaks. In Appendix. A.2.2, we also consider the cases when the entanglement cuts are located at the center(s) of the energy-momentum density peaks. Then some interesting features in the entanglement entropy could arise.
To summarize, we put the phase diagrams and related quantities in the periodically driven CFT in Table. 1.

A minimal setup
Now we consider a minimal setup of the periodically driven CFT to demonstrate the main features in the previous discussions. In this setup, we consider only p = 2 driving steps within one period That is, we drive the CFT with (H 0 , T 0 ) and (H 1 , T 1 ), where T 0 and T 1 are the time intervals. We consider a SL 2 deformed Hamiltonian with q = 1 with open boundary conditions 14 : We choose H 0 and H 1 as H θ=0 and H θ =0 . Note that H θ=0 corresponds to the uniform Hamiltonian, and H θ=∞ corresponds to the SSD Hamlitonian in Eq. (15) up to an overall factor 2. The regime with λ L > 0 corresponds to the heating phase, while the dark blue regime with λ L = 0 corresponds to non-heating phase, and boundary between them is the phase transition.
Denoting the time interval of driving as T θ , then the corresponding Möbius transformation M (H θ , T θ ) has the following form Here L eff = L cosh(2θ) denotes the effective length of the total system. Physically, it characterizes the effective distance that the quasiparticle needs to travel to return to its original location [43].

Phase diagram and Lyapunov exponent
The Lyapunov exponent λ L is determined by the trace of the transformation matrix Π p as shown in (60). In our setting p = 2 and where x 0 = πT 0 L and x 1 = πT 1 L eff , with L eff = L cosh(2θ). Therefore, inserting into we obtain the result shown in Fig. 7. From the figure, we can also read out the phase diagram straightforwardly, namely the regime with λ L > 0 corresponds to the heating phase, while the dark blue regime with λ L = 0 corresponds to non-heating phase, and the boundary between them is the phase transition. We also show in Fig. 8 the group walking pictures for θ = 0.2 with different choices of T 1 /L. To gain some analytical understanding of the formula for Lyapunov exponent, let us consider a simple example when θ = ∞, namely H 1 = H SSD . Along the line T 0 = L/2, (73) simplifies to Tr(Π p=2 ) = 2πT 1 L . And therefore, the Lyapunov exponent is a function of T 1 given as follows In particular, in the limit T 1 L, we have That is to say, along the line T 0 = L/2, the Lyapunov exponent grows logarithmically with T 1 /L in the large T 1 limit (T 1 /L 1). In fact, for θ = ∞, the result in (76) holds for arbitrary T 0 /L = nπ (n ∈ Z) in the large T 1 /L limit. Now we would like to make a few comments on Fig. 7 1. The area of regime with larger Lyapunov exponent grows when we increase the parameter θ in Hamiltonian. The heuristic argument is that when the 'difference' between H 0 and H θ is greater, the driving protocol is easier to heat the system.
2. When we push θ → 0, the area of heating phase, namely the regime with λ L > 0 decrease to zero. However, there is always at least a point (T 0 /L, T 1 /L eff ) = (1/2, 1/2) staying in the heating phase for arbitrary θ = 0. This point corresponds to x 0 = x 1 = π 2 in Eq. (73) where we find Tr(M 0 M 1 ) = 2 cosh(2θ) > 2 for arbitrary θ = 0. Now, let us take a closer look at the special point (T 0 /L, T 1 /L eff ) = (1/2, 1/2) in the heating phase. In terms of the transformation matrix M 0 and M 1 given in (72), we have These two matrices are actually special examples of a larger class of SU(1, 1) matrix defined below Reflection. Let M be an SU(1, 1) matrix parametrized as (43), namely is a reflection. This condition is equivalent to demanding matrix M is traceless Tr(M (ρ, ζ)) = 0. One important property of the reflection matrix is that it squares to −1, i.e.
Apparently, a reflection matrix is elliptic since its trace is smaller than 2. But the product of two distinct reflection matrices are hyperbolic, i.e, | Tr(M ( The reflection matrix will also play an important role in later discussions on both the quasiperiodically driven and the randomly driven CFTs [57]. Applying to our case where we have two distinct reflection matrix M 0 and M 1 , we conclude that Π 2 = M 0 M 1 is hyperbolic and therefore induces a heating phase for (T 0 , T 1 ) = (L/2, L eff (θ)/2). 15 Indeed, we can explicitly check that the Lyapunov exponent has a simple form (assuming θ > 0) In addition, we have which further fixes the location of the energy-momentum peak to be at x = 0 (x = L) for the chiral (anti-chiral) mode (c.f. Eq. (32) for the open boundary condition discussed here). In fact, the chiral and anti-chiral peaks switch positions after each driving period. The total energy and entanglement entropy are also expressible using θ (c.f. (31) and (33)) where p = 2 and we consider A = [0, L/2] here. We will compare this CFT result with the numerical calculation on a lattice model in the next subsection.

Numerical simulation on lattice
In Ref. [36], the authors compare the CFT and lattice calculations on the entanglement entropy evolution in a periodically driven CFT. It was found that the comparison agrees very well in the non-heating phases, but deviates in the heating phase. The heuristic reason is that the two driving Hamiltonians in [36] are chosen as H 0 and H θ=∞ , which result in a large Layapunov exponent in the heating phase (See Fig. 7). Then the system can be easily heated up with only a few driving steps. Now, by considering the general H θ , we can tune the system to have a small heating rate by choosing a small θ. The lattice model we consider is a free fermion lattice, which has finite sites L with open boundary conditions. We prepare the initial state as the ground state of with half filling. The SL 2 deformed Hamiltonian has the form where f (j) = 1 − tanh(2θ) · cos 2πj L , c j are fermionic operator satisfying the anticommutation One can refer to the appendix in Ref. [36] for the details of calculation of the entanglement entropy and correlation functions. The comparison of the numerical and CFT calculations on both the entanglement and energy time evolution can be found in Fig. 9. The agreement is remarkable.
One can also refer to Appendix. A.2.3 for the interesting case that when the entanglement cut lies at the center of both the chiral and anti-chiral energy-momentum density peaks, the entanglement entropy can decrease linearly in time.

Quasi-periodic driving
In this section, we will study the non-equilibrium dynamics in a quasi-periodically driven CFT with SL 2 deformed Hamiltonians. We would like to understand the following two questions in this section: Fibonacci quasi-periodical driving 2. What is the generic feature of the entanglement/energy evolution in the quasi-periodically driven CFT?
As an initial effort to answer these questions, we will mainly focus on the case of quasiperiodical driving with a Fibonacci sequence, which is simpler to handle compared to a more general quasi-periodic sequence. Our setup is closely related to the Fibonacci quasi-crystal, which was proposed in the early 1980's by Kohmoto, Kadanoff, and Tang [63], and Ostlund, Pandit, Rand, Schellnhuber, and Siggia [64]. It was observed and later proved that the spectrum of Fibonacci Hamiltonian is a Cantor set of zero Lebesgue measure [63][64][65]. Since then, the Fibonacci dynamics has been extensively studied in both physics and mathematics. See, e.g., Ref. [67] for a recent review. For simplicity, in the following we may call the quasi-periodically driven CFT with a Fibonacci sequence as a Fibonacci driven CFT.
In the end of this section, we also discuss another kind of quasi-periodic driving with Aubry-André-like sequence by focusing on the properties of its phase diagram.

Fibonacci driving and relation to quasi-crystal
We start with an introduction to the setup and tools we use to analyze the Fibonacci driving, many of which are borrowed from the rich literature of Fibonacci quasi-crystals.

Setup and trace map
A Fibonacci driving in this paper is generated by two unitaries The Hamiltonians H A , H B are chosen to be the SL 2 deformed Hamiltonian same as the ones used in the previous sections. Therefore, each unitary U A(B) corresponds to a conformal map M A(B) ∈ SU(1, 1) and the final conformal map Π n that determines the operator evolution is given as a product For example, the first few matrices An useful property of the Fibonacci driving Π n is that for n being a Fibonacci number 16 F k with k 3 there is a recurrence relation for its trace This relation was used in quasi-crystal literature, e.g. see Ref. [63]. Also see Appendix B for a derivation following the substitution rule of the Fibonacci bitstring and the property that det Π n = 1.
The initial conditions for this recurrence relation can be taken as It is sometime convenient to define an auxiliary The auxiliary element x F 0 is defined such that the recurrence relation (90) also holds for k = 2.
For SU(1, 1) matrices, we have x F k ∈ R. To visualize the trace map, let us introduce a three dimensional vector ( , then the trace map in (90) can be expressed as the following mapping between points in three dimensional space with the initial condition l = (x F 3 , x F 2 , x F 1 ) given in Eq. (91), or alternatively we can use (x F 2 , x F 1 , x F 0 ) with the auxiliary element x F 0 . Remarkably, the trace map has a constant of motion [63] see Appendix B for an explicit check that I is independent of k.

Example with H θ and fixed point
Let us now take explicit example of SL 2 deformed driving Hamiltonians.
where H 0 is taken as the CFT Hamiltonian with a uniform Hamiltonian density, and H θ is taken as the SL 2 deformed one in Eq. (71) The corresponding conformal transformation M A and M B has been computed in (72) and copied here (95) 16 Our convention for the Fibonacci number is that And L eff = L cosh(2θ) denotes the effective length of the system under H θ . Therefore, the initial conditon for the trace map is given as follows And the invariant I defined in (93) is Generally speaking, the invariant I constrains the motion of (x F k , y F k , z F k ) on a two dimensional manifold M. For I ∈ R, there are three topologically distinct scenario as shown in Fig. 11 1. I = 0: The manifold M can be decomposed into five parts. The central part is the curvilinear tetrahedral ('island'), with the vertices/singularities at A (1, 1, 1), B (1, −1, −1), C (−1, 1, −1), and D (−1, −1, 1). The tetrahedral is parameterized by θ 1 and θ 2 with x = cos(θ 1 + θ 2 ), y = cos θ 1 , z = cos θ 2 . The left four parts are funnels. The first funnel is parameterized by x = cosh(θ 1 + θ 2 ), y = cosh θ 1 , and z = cosh θ 2 , with its vertex at the point A. The other three funnels are similar defined with the vertices at B, C, and D. In the Fibonacci driven CFT, this case corresponds to T 0 /L ∈ Z or T 1 /L eff ∈ Z. Physically, this corresponds to a single quantum quench which is not our focus here.
2. I > 0: The four vertices A, B, C and D are replaced with four necks, which connect the central part ('island') of the manifold to the four funnels. The whole manifold is therefore non-compact. This case corresponds to all the nontrivial choices of (T 0 , T 1 ) in our setting (97). It turns out that for almost all the initial points on the manifold, they will flow to infinity under the trace map in (90). [62] 3. I < 0: The central part ('island') becomes disconnected to the outside funnels and therefore compacted. This case is absent in our setting for the Fibonacci driving.Nevertheless, this case may be related to some non-Hermitian Hamiltonian or non-unitary time evolution and deserves a careful study in future. For a fixed I > 0, one can tune two of the three parameters (θ, T 0 , T 1 ) to move the initial point (x F 3 , x F 2 , x F 1 ) on the surface M then the orbit under the trace map is completely determined. As we will show in the following sections, most of the orbits will escape to the infinity and resulting an heating phase. However, there still exists returning orbit, e.g. when we have two zeros in the initial condition (x F 3 , x F 2 , x F 1 ), we will end up with a period 6 orbits with a = (1+I) 1/2 , see Fig. 12 for an illustration. We will call such initial points that correspond to the non-heating point as "fixed point", in the sense that those points are fixed under T 6 action.

Phase diagram: From periodic to quasi-periodical driving
In this section, we show the shape of the phase diagram of a Fibonacci driven CFT via numerically approaching the Fibonacci bitstring by its finite truncation. This strategy has been proven useful in the analysis of the energy spectrum of a Fibonacci quasi-crystal [63]. In the quasi-crystal case, the energy spectrum forms a Cantor set of zero Lebesque measure. In this section, we will show numerical evidence of such "fractal" structure, while in the next section we will map our phase diagram to the energy spectrum of quasi-crystal and establish the claim. Recall that (in Appendix. B) we generate the Fibonacci driving using the Fibonacci bitstring where χ(t) = χ(t + 1) is a period-1 characteristic function and ω = is an irrational number with a simple continued fraction representation Now to approach the Fibonacci bitstring from a periodic string, we can truncate the continued fraction of ω at finite order n and obtain a rational number (principal convergent) ω n = F n−1 /F n , namely the ratio of two nearby Fibonacci number. The corresponding bitstring {X j } now has periodicity F n and therefore produce a periodic driving. We can now use the tools introduced in Sec. 4 to obtain a phase diagram for each ω n .
In Fig. 13, we show the evolution of phase diagrams of periodically driven CFTs with protocol (H A , T A ) = (H θ=∞ , T 1 ) and (H B , T B ) = (H 0 , T 0 ). The phase diagram is periodic in T 0 direction with period L, we only show the phase diagram within one unit cell 0 T 0 L. As we increase n, there are two notable features Figure 14: Phase diagrams in a periodically driven CFT with the sequence generated by finitely truncated Fibonacci bitstring, i.e. {X j } with ω n = F n−1 /F n . Here we choose n = 2, 4, 5, 6, 10, 20, 100, and 1000 for 8 plots respectively. The two Hamiltonians we use are H 0 (θ = 0) and H 1 (θ = 0.5) in (71). The phase diagram is periodic in T 0 direction with period L and in T 1 direction with period L cosh(2θ) 1.54L. The blue (yellow) regions correspond to the heating (non-heating) phases.
1. The number of regions of the non-heating phases increases with n, and tends to infinity as n → ∞.
2. The measure of the non-heating phases decreases with n, and tends to zero as n → ∞.
These two features suggest that the non-heating phases in the quasi-periodical driving limit may form a Cantor set of measure zero, analogous to the feature of the energy spectrum in a Fibonacci quasi-crystal. In fact, this is indeed the case, as we will discuss in detail in the next subsection.
We also present the evolution of phase diagrams by the Hamiltonians H 0 and H θ with finite θ. See Fig. 14 for θ = 0.5, and Fig. 39 in the appendix for θ = 0.2. The two features mentioned above are also observed in these cases. It is noted that for a finite θ in H 1 (θ), the phase diagram is also periodic in T 1 direction, with the period L cosh(2θ).

Exact mapping from a Fibonacci driven CFT to a Fibonacci quasi-crystal
The features in the phase diagrams in Fig. 13 and Fig. 14 suggest that the non-heating phases in the quasi-periodical driving limit may form a Cantor set of measure zero. In this subsection, we verify this by performing an exact mapping between the phase diagram of a Fibonacci driven CFT and the energy spectrum of a Fibonacci quasi-crystal. The latter has been proved mathematically that the energy spectrum is indeed a Cantor set [65] (See also Ref. [62] for a review). Figure 15: Phase diagrams in a periodically driven CFT with the sequence chosen in (101) where ω n = F n−1 /F n .The parameters are the same as those in Fig. 13 except that now we change the variables to V CF T and E CF T (See Eq. (112)). The blue (yellow) regions correspond to the heating (non-heating) phases.
Before introducing the mapping, let us first briefly review the background of the Fibonacci quasi-crystal. We consider the discrete Schrödinger operators of the form where ψ j := j|ψ is the position-space wavefunction, with j labeling the j-th site, and V j is the onsite potential. For eigenvalue problem Hψ = Eψ, it is useful to consider the transfer matrix Denoting Ψ j = (ψ j+1 , ψ j ) T , we have 17 In the Fibonacci quasi-crystal, the potential V j can also be generated by the Fibonacci bitstring {X j } (191) as follows The allowed energy spectrum E is determined by requiring that Defining T Fn := T Fn T Fn−1 · · · T 1 , and x Fn := 1 2 Tr( T Fn ), it turns out the traces {x Fn } satisfy the same recurrence relation in Eq. (90). The only difference between the Fibonacci driving CFTs and the Fibonacci quasi-crystals is the initial conditions, which we will specify now. By taking V B = −V A = V , one can find the initial conditions for the Fibonacci quasi-crystal are 18 The invariant I in the constant of motion in Eq. (93) becomes I = V 2 . In a quasi-crystal, the potential V is fixed, and therefore each E specifies an initial condition, which may flow to infinity by iterating the trace map (E is in the gap), or is bounded (E is in the spectrum).
Next, let us compare the initial conditions in the Fibonacci driven CFT. We consider the phase diagrams in Fig. 13, which correspond to H 0 (θ = 0) and H 1 (θ = ∞). By taking the limit θ → ∞, the initial conditions in Eq. (96) become 19 and the invariant I in Eq. (93) is To compare with the initial conditions of Fibonacci quasi-crystal, here we choose (x F 3 , x F 2 , x F 1 ) instead of (x F 2 , x F 1 , x F 0 ) as the initial condition. Based on Eq. (90), one can obtain Now by defining E CFT := 2 cos πT 0 L + 3V CFT , V CFT := − πT 1 L · sin πT 0 L , then the initial condition line can be written as: with the invariant I in Eq. (110) expressed as That is, by redefining variables, we can find a map between the initial conditions in Eq. (108) and Eq. (112). With this map, the allowed energy E in the spectrum of a Fibonacci quasi-crystal  Fig. 15 where the non-heating phases vanishes for larger V CFT .
is mapped to the non-heating phase in a Fibonacci driven CFT specified by E CFT (T 0 , T 1 ), and vice versa. One should note that, however, on the CFT side, (E CFT − 3V CFT ) ∈ [−2, 2] always lives in a window of finite width. On the quasi-crystal side, we have (E − 3V ) ∈ (−∞, +∞). This means the non-heating phases in a quasi-periodically driven CFT are only mapped to part of the energy spectrum in the Fibonacci quasi-crystal. This can be intuitively seen by considering the initial condition lines in Eqs.(108) and (112) on the two dimensional manifold M determined by Eq. (93). As shown in Fig. 16, the overlap of l (E CFT ;I) and l (E;I) is always a straight line of finite length 2. For smaller V or V CFT , l (E CFT ;I) overlaps with l (E;I) mainly in the region with |x|, |y| 1, |z| = 1 in the middle 'island'. For the initial conditions in this region, they are much easier to be bounded as we iterate the trace map [68]. As V or V CFT increases, the overlap of l (E CFT ;I) and l (E;I) moves gradually from the 'island' in the middle to the 'funnel' outside. Then it becomes more difficult for the initial conditions to stay bounded as we iterate the trace map. This analysis agrees with the fact that in the phase diagrams in Fig. 15, there are no non-heating phases observed for large V CFT .
This "inclusion map" for small V CFT is totally fine for our goal: Since the energy spectrum of a Fibonacci quasi-crystal forms a Cantor set of measure zero, then part of the energy spectrum 18 Here we use the recurrence relation to infer the value of x F0 and x F−1 from (x F3 , x F2 , x F1 ), the reason we choose to start with x F−1 = 1 for quasi-crystal is that we need a convenient base point to map to the CFT initial point, whose x F1 happens to be 1 as well. It should be clear later when we construct the mapping. 19 Note that by taking the limit θ → ∞, we always consider finite T 0 and T 1 such that T 0 , T 1 L cosh(2θ) when θ → ∞. In this case, the initial conditions (x F2 , x F1 , x F0 ) CFT form a straight line with y = 1, rather than a closed loop in Fig. 12. See Fig. 16 for the initial conditions (x F3 , x F2 , x F1 ) CFT with θ = ∞. (which is a connected and finite region in the parameter space) is also a Cantor set of measure zero. Then with the exact mapping discussed above, we conclude that the non-heating phases in the quasi-periodically driven CFT form a Cantor set of measure zero.
Furthermore, in Fig. 17, we also check explicitly the measure of the non-heating phases in the phase diagrams in Fig. 15 as we approach the quasi-periodic limit. The procedure of obtaining the measure is as follows: Fixing a V CFT (or equivalently the invariant I) in Fig. 15, for each ω n = F n−1 /F n , there are many 'energy bands' of non-heating phases. Denoting the band width of the j-th band as d j (n), this band width depends on both I and θ (which is ∞ here). Then the measure of non-heating phases with ω n = F n−1 /F n is defined as where E max (n, I, θ)−E min (n, I, θ) is the total width of the energy window, which is 4 for θ = ∞.
That is, the measure of the non-heating phases decreases exponentially as a function of n, and tends to become 0 in the limit n → ∞. This agrees with the fact that the non-heating phases in the quasi-periodical driving limit form a Cantor set of measure zero. The decaying rate λ(I, θ) may also be interpreted as the escape rate, since it describes the rate of initial conditions in Fig. 12 escaping into the infinity. Also, we remind here that the real driving steps are F n Figure 18: Phase diagrams in a periodically driven CFT with the sequence generated by finitely truncated Fibonacci bitstring, i.e. {X j } with ω n = F n−1 /F n . The parameters are the same as those in Fig. 14 except that now we change the variables to V CF T and E CF T (See Eq. (112)). The blue (yellow) regions correspond to the heating (non-heating) phases.
rather than n. And F n ∼ ω −n at n → ∞, the measure of non-heating phases depends on F n as σ n (I, θ) ∝ F λ(I,θ)/ log ω n for large n. That is, σ n (I, θ) decays polynomially as a function of the driving steps F n . In addition, we check how the decaying rate λ(I, θ) depends on the invariant I. As shown in Fig. 17 (right), it is found that λ(I, θ) depends on I as λ(I, θ = ∞) = a · I b , with a 0.285 and b 0.552. This monotonic dependence is reasonable in the sense that a smaller I corresponds to a narrower neck connecting the 'island' and 'funnel' (See Fig. 11 and Fig. 16), which may suppress the escape rate from the island to the funnel.

Cases that cannot be mapped to Fibonacci quasi-crystal
The exact mapping studied in the previous subsection applies for the case of θ → ∞ in H 1 (θ). For a finite θ, we do not have such an exact mapping. Here we hope to study the common features among the phase diagrams with different θ (See, e.g., the phase diagrams in Fig. 13,  Fig. 14, and Fig. 39 in the appendix).
As analyzed in the previous subsections, to study the measure of the non-heating phases or the escape rate of initial conditions to infinity on the manifold M(See Fig. 12), it is more appropriate to fix the invariant I in Eq. (93). This is because the trace map in Eqs.(90) or (92) holds for a fixed invariant I. In other words, the points (x Fn , y Fn , z Fn ) in Eq. (92) move on the manifold M with a fixed geometry. For this reason, we can replot the phase diagram in Fig. 14 by changing variables in the initial conditions in Eq. (96) as follows: where L eff = L cosh(2θ). with the invariant With the above procedure, now we map the phase diagram in the region {(T 0 , T 1 )|0 T 0 /L 1, 0 T 1 /L eff < 1/2} in Fig. 14 to Fig. 18. The merit of this mapping is that for each V CFT in Fig. 18, the invariant I is fixed. Then we study the measure of the non-heating phases as defined in Eq. (114), with the result shown in Fig. 17. There are two interesting features: 1. Similar to the case of θ = ∞, the measure of the non-heating phases depends on n as σ n (I, θ) ∝ e −λ(I,θ)·n . That is, the measure of the non-heating phases decays exponentially (power-law) as a function of n (F n ), indicating that the measure will become zero in the quasi-periodical driving limit n = ∞.
2. Interestingly, the decay rate (or escape rate) λ(I, θ) for θ = 0.5 and θ = ∞ collapse to the same curve with λ I = a · I b , where b 0.552 (See the right plot in Fig. 17). This means λ(I, θ) is only a function of I, and is independent of θ.
In addition, in Fig. 17, we also present the results for the measure of non-heating phases for the case of θ = 0.2 (See Fig. 39 and Fig. 40 for the corresponding phase diagrams). The decaying rates λ(I, θ = 0.2) as a function of I again fall on the same curves as that of θ = 0.5 and θ = ∞, as seen in Fig. 17 (right plot). This means the decay rate λ(I, θ) is only a function of the invariant I, but is independent of θ which characterizes the concrete deformation of Hamiltonians.

Lyapunov exponents in the quasi-periodical driving limit
The phase diagrams as studied in the previous subsections simply tell us whether the CFT is in the heating or non-heating phases. As n → ∞, the measure of the non-heating phases becomes zero, and one can only "see" the heating phase in the phase diagram. In this subsection, we will use Lyapunov exponents to further characterize the fine structures in the heating phases in the limit n → ∞.
Let us first consider a periodical driving with ω = ω n = F n−1 /F n in Eq. (101), where the period of driving is F n . The Lyapunov exponent in the heating phase can be obtained via Eq. (60) as: where x Fn = 1 2 Tr(Π Fn ) can be efficiently computed using the recurrence relation. Now we consider the distribution of Lyapunov exponents in Fig. 15 in the quasi-periodical driving limit. To be concrete, we fix V CFT (or equivalently I) in Fig. 15, and scan the Lyapunov exponents along E CFT . As shown in Fig. 19, it is found that the Lyapunov exponents exhibit self-similarity structures. That is, by zooming in the distribution of Lyapunov exponents, one can find the same distributions (in different scales). One can zoom in the distribution all the way and see the self-similarity structure, as long as a large enough n is taken. The self-similarity structure of Lyapunov exponents also indicates that the Lyapunov exponents can be arbitrarily small. In other words, in the heating phases of a Fibonacci quasi-periodic driving CFT, there exist some regions with arbitrary small heating rates for the entanglement/energy growth.
We also study the distribution of Lypunov exponents in the T 0 /T 1 parameter space in Fig At last, in Fig. 21, we give a color plot of the distribution of Lyapunov exponents in the parameter space (T 0 /L, T 1 /L). One can find the patterns inherit some features of the periodically driven CFTs (See Fig. 7). It is also helpful to compare these three plots with the phase diagrams in Fig. 13, Fig. 14, and Fig. 39, respectively. We emphasis that although there are large areas of regions with almost zero Lyapunov exponents, they are actually in the heating phase (See Fig. 13, Fig. 14, and Fig. 39). If we zoom in these regions, one can observe the self-similarity structure(See, e.g., Figs.19 and 20).

Fixed point in the non-heating phase: Entanglement and energy dynamics
From the previous discussions, we conclude that the measure of the non-heating phases shrinks to zero when we approach the quasi-periodical driving limit, i.e. without special guide it is hard 3. The envelopes of the entanglement entropy and total energy will grow logarithmically and in a power-law as a function of the driving steps n (not the Fibonacci index), respectively.
We will first illustrate the above statements with simple examples in the following discussions, and then prove them in Sec. 5.2.4.

Entanglement and energy evolution at Fibonacci numbers
In the following, we will study the exact non-heating fixed point and its properties with the simple choice of H A = H θ and H B = H θ=0 , where the form of H θ is given in Eq. (71). The the initial condition has the following form where a := cosh(2θ). In fact, for this fixed point, not only the traces (recall x F k = 1 2 Tr(Π F k )) have periodicity x F k = x F k+6 , the corresponding matrices themselves are also returning periodically total energy at the Fibonacci numbers are where S A (0) denotes the entanglement entropy in the initial state (which is the ground state of H θ=0 here), and E F 6k = − πc 24L corresponds to the Casimir energy. Also see Fig. 22 and Fig. 23 for examples with c = 1 and θ = 0.5 and 0.2.
An interesting remark is that the special initial condition we choose that forms the fixed point are conformal maps that correspond to the "reflection matrix" (since we require the initial traces to vanish, see discussions near (77)). Note, the product of two distinct reflection is hyperbolic, i.e. if we drive the system with periodic driving we will end up heating the system. However, what we present just now is that if we drive it in a Fibonacci pattern, they happen to return and form a non-heating point in the quasi-periodic driving phase diagram.

Entanglement and energy evolution at non-Fibonacci numbers
As shown in Fig. 22 and Fig. 23, between the Fibonacci numbers, the numerical simulation shows that entanglement and energy still increase for initial condition at fixed point. In the content of Fibonacci quasi-crystal, it was found that the wavefunction amplitude in the energy spectrum and some physical observables (e.g., the resistence) have a power-law growth as a function of the lattice site n [69][70][71]38]. In our setup, at the driving steps that are non-Fibonacci numbers, we expect the entanglement entropy or total energy also grows in a certain sub-exponential way. Now we provide analytic understanding using the property of the Fibonacci driving protocol. The idea is that, for any integer n which can be written as a sum of distinct Fibonacci numbers, the corresponding conformal transformation matrix Π n can be written as where each Π F k with Fibonacci number has been obtained in Eq. (122). In particular, we can find some simple sequence 1. Let be an integer that increase with m. Correspondingly Recall from (122), the product Π 6 · Π 3 and Π n can be evaluated explicitly The corresponding entropy for the half system A = [0, L/2] and the total energy are where we have constrained θ > 0. The plot of Eq. (129) can be found in Fig. 22 and Fig. 23. Note that S A (n) and E(n) grow linearly and exponentially respectively as a function of m for large m. However, m is not the actual driving step number. We need to convert to the actual step number n ∼ e 6m log ϕ , which grows exponential with m for large m, and ϕ = √ 5+1 2 is the golden ratio. Therefore That is to say, at the non-Fibonacci numbers n in Eq. (126), the entanglement entropy S A (n) grows logarithmically in time, and the total energy grows in a power-law in time. This corresponds to the feature of phase transition (or critical phase) in the periodically driven CFT (See Table.1).
2. Now we choose a difference sequence, to demonstrate that the growth rate at this fixed point depends on the sequence we pick when approaching the long time limit. Let And correspondingly The entropy and energy formulas are In the two examples above, the entanglement entropies all grow with n. One can observe in Fig. 22 that at certain points the entanglement entropy may decrease. We will investigate these points using the following examples Note the last element is important. The corresponding matrix is given as The corresponding entropy and energy formulas are where θ > 0. The results are similar to Eq. (133), with a minus sign difference in the entropy formula. In other words, we have a logarithmic decrease in the entanglement evolution and a power-law growth in the total energy evolution (See Fig. 22 and Fig. 23).
The entanglement decrease might look bizarre, but this could happen in a system with infinite entropy to start with, e.g. in the continuous field theories where a UV regulator is required in the entropy calculation, which itself is a manifestation of the large entanglement in the vacuum state.
Technically, we may explain the origin of the decreasing entropy as follows: the form of the conformal transformation in (136) indicates that the energy-momentum density (see (32)) locates exactly at x = L/2, which coincides with the entanglement cut we choose. As discussed in detail in Appendix A.2.2, in this case, the entanglement entropy will decrease in time. Physically, it is because the degrees of freedom that carry the entanglement between two regions are accumulated at the entanglement cut. We emphasize that the points with decreasing entanglement entropy are due to the coincidence of the energydensity peak and the entanglement cut, and therefore are not generic. In general, at the non-heating fixed points, the envelopes of the entanglement entropy and total energy will grow logarithmically and in a power law in time, respectively.
2. Another example with decreasing entropy we present in the Fig. 22 and Fig. 23 is that and the corresponding entropy and energy Using the same procedure above, one can find many other series of discrete points with different growing (and decreasing) rates in the entanglement/energy evolution in Fig. 22 and Fig. 23, these series together form the fan structure in the Figures.

Comparison of CFT and lattice calculations
In this part, we compare the CFT and lattice calculations for the time evolution of entanglement entropy/energy at the non-heating fixed point as discussed in the previous subsections. The lattice model we use is the same as that studied in Sec. 4.4.2. The two lattice Hamiltonians under consideration are where L is the total length of the lattice and f (j) = 1 − tanh(2θ) · cos 2πj L , with the initial state chosen as the ground state of H 0 . The corresponding driving time intervals are T 0 = L/2 and T 1 (θ) = L eff /2, where L eff = L cosh(2θ). Then we drive the system with the Fibonacci sequence as introduced in Sec. 5.1.1. Fig. 24 presents the comparison of the lattice and CFT results on the entanglement/energy evolution. We find that the agreement is remarkable. A comparison on the entanglement entropy evolution with larger driving steps n can be found in Fig. 25. In general, the agreement will break down for a large enough n, since more higher-energy modes will be involved as n increases (Recall that the envelope of the total energy growth is power-law in time at the nonheating fixed point). On the lattice model, the high-energy modes are no longer well described by a CFT, and therefore there must be a breakdown at certain n. 21 21 More precisely, let us denote n * as the driving step at which the agreement between CFT and lattice calculations break down. From Fig. 24, one can observe that n * is a monotonically increasing function of L. This dependence can be understood as follows: One may consider the wavefunction in the 'Fock space' (which is Verma module here) of a CFT of finite length L. The initial state is the ground state |0 . As we drive the system, higher energy modes |N (N > 1) will be involved. It is noted that N is independent of the length L of the CFT. Since the energy spacing is proportional to 1/L. one can find the energy E N (L) corresponding to |N

Exact nonheating fixed points in more general cases
With the concrete examples illustrated in the previous discussions, now we are ready to prove the statements as mentioned in the beginning of Sec. 5.2, which we rewrite here: 1. If both the driving Hamiltonians are chosen as elliptic types, one can always find exact fixed points in the non-heating phases.
2. At these (non-heating) fixed points, both the entanglement entropy and the energy evolution are of period 6, i.e., S A (F k ) = S A (F k+6 ) and E(F n ) = E(F n+6 ).
3. The envelopes of the entanglement entropy and total energy will grow logarithmically and in a power-law as a function of the driving steps n, respectively.
Proof of claim 1. The non-heating fixed point has initial condition ( is higher than E N (L ) if L < L . For a small L, E N (L) may be in the high-energy region which are no longer described by a CFT. However, by increasing L to a large enough L , we can push E N (L ) into the low-energy region which are well described by a CFT. This is why we have a better agreement in Fig. 24 for a larger L.
In other words, to prove claim 1, we only need to find suitable initial conditions T A and T B for two elliptic Hamiltonian H A and H B such that the SU(1, 1) matrices M A and M B are traceless, i.e. are reflection matrices (see (79) is the wavelength of deformation (See, e.g., Eq. (13)). One can obtain the reflection matrix by choosing T = 1 2 L eff in Eq. (143), where L eff := l/C is the effective length. Then Eq. (143) becomes which is traceless obviously. That is to say, to arrive the non-heating fixed point, we need to set T A(B) = l/2C A(B) and the corresponding reflection matrices M A and M B take the form of (144) with subscripts A and B.
Proof of claim 2. Having two reflection matrix M A and M B , we immediately have the following useful property Now let us use this to prove the claim 2.
To verify the periodicity of entanglement entropy S A (F n ) and energy E(F n ), it is sufficient to show the periodicity of the conformal transformation matrix Π Fn . Let us first example the case for n = 0 . . . 6. Following the definition, the first 3 are Note the Π F 0 = Π F 1 is introduced as a convenient notation which satisfies the recurrence relation. For n = 3, 4, we can use the recurrence relation and find where in the last equation we have used the property that M 2 A = M 2 B = −1. The result two n = 5, 6 are also obtained analogously Note we have already verified the periodicity for And the recurrence relation for Π Fn only depends on the previous two elements, therefore we conclude The claim 2 follows immediately.
Proof of claim 3. First Ref. [71,72] generally prove that the norm of Π n is polynomially bounded at λ L = 0. Therefore, we only need to show there is at least a sequence of point that approaches infinity polynomially similar to what we showed in last section for a specific example. Note, for any two distinct reflection matrix M A and M B , their product is hyperbolic. Then we can use their product to generate a sequence of points, e.g. For Therefore, the total energy will grow exponentially w.r.t. m as if in the heating phase. However, remember the step index n is also exponential in m and therefore in terms of physical steps n, the energy grows in power law. Similarly, we have the entanglement entropy grow logarithmically. 22

Entanglement and energy dynamics in the heating phases
In this subsection, we present the rich features in the heating phase, including the time evolution of total energy, energy-momentum density and entanglement entropy.

Group walking and entanglement/energy evolution
As discussed in Sec. 3.2, the group walking of ρ and (ρ ζ) of the SU(1, 1) matrix in (43) reflect the time evolution properties of the total energy, energy-momentum density, and the entanglement entropy. In the following, we will study the group walking of ρ and ρζ, as well as the energy/entanglement evolution. The detailed discussion on energy-momentum density is left to the next subsection. First, we consider the group walking of ρ in Π n Π n = 1 As shown in Fig. 26 is a sample plot of the group walking of ρ and ρζ nearby a certain point (which is T 0 /L = 1/2 and T 1 /L eff = 0.04 here, with L eff = L cosh(2θ)) in the parameter space. See also Fig. 27 for another sample plot. The two driving Hamiltonians are H 0 and H θ=0.5 , with H θ given in Eq. (71). As the driving steps n increase, one can find that ρ will walk to a certain point ρ n=∞ on ∂D. This is quite an interesting feature, since ρ does not move around on ∂D in the limit n → ∞. 23   The parameters are the same as those in Fig. 26. Note that the total energy and entanglement entropy (in the left plot) are plotted in a small window to see the detailed oscillating structure.
As we have discussed in Sec. 3.2, as ρ approaches the boundary ∂D, both the entanglement entropy and the total energy of the system will grow accordingly, as seen from Eqs. (45) and (49). In Fig. 26 and Fig. 27, one can find that for different (T 0 , T 1 ), ρ approaches the boundary ∂D with different rates, which correspond to different Lyapunov exponents and different growth rates in the entanglement entropy/total energy. This can be intuitively seen by looking at the entanglement/energy growth in Fig. 28 and Fig. 29, respectively. Also, the patterns in |ρ n | in Fig. 26 and Fig. 27 result in the oscillating features in the entanglement/energy growth in Fig. 28 and 29. 24 As a remark, it is noted that for the case of T 1 /L eff = 0.0401 in Figs. 26 and 28, ρ does not walk to the boundary ∂D, and the total energy seems to simply oscillate in n. This is because the Lyapunov exponent at this point is too small. One needs to take a longer time (larger n) to observe the growth behavior.
Second, Let us consider the group walking of (ρζ). As shown in Figs.26 and 27, different from the behavior of ρ, which walks to a fixed ρ ∞ in the n → ∞ limit, (ρζ) will walk around on the boundary ∂D even in the n → ∞ limit. This will result in two interesting features: 1. The locations of energy-momentum density peaks will move around during the quasiperiodic driving, as seen from Eq.(48).
2. More oscillating structures will be introduced in the time evolution of the entanglement entropy. As discussed in Sec. 3.2, the total energy only depends on |ρ|, but the entanglement entropy depends on both |ρ| and ρζ, as seen in Eq. (49). This extra oscillating structure can be found in Figs.28 and 29, in particular for the case with large Lypunov exponents.  In a short summary, there are rich patterns in the group walking of ρ and (ρζ) in the heating phase of a Fibonacci quasi-periodically driven CFT. The velocity of ρ walking towards ∂D determines the value of Lyapunov exponents, and therefore the growth rate of the entanglement entropy and total energy. The behavior of |ρ| determines the concrete features of the energy growth E(n) through Eq. (45). In general, one can observe various oscillating features in the growth of E(n). Both |ρ| and (ρζ) determine the features of S A (n) through Eq. (49), where more patterns of oscillations can be observed comparing to E(n).

Locations of energy-momentum density peaks
There are several features on the distribution of energy-momentum density we hope to point out in the heating phase of a Fibonacci driven CFT: 1. In the heating phase, there is an array of peaks of energy-momentum density distributed in the real space. The locations of these peaks are determined by Π n through Eq. (48) In fact, as we have shown in Sec. 3.1, a positive Lyapunov exponent λ L always indicates an array of peaks in the energy-momentum density T (x, n) in real space.
2. Different from the periodically driven CFT, where the peaks are located at the same positions after each period of driving (See Eq. (66)), in the quasi-periodical driving, since there is not a regular driving period, the locations of the energy-momentum density peaks will in general move around after each driving step. This can be seen based on the group walking of (ρζ), which determines the locations of energy-momentum density peaks [see Eq. (48)], in Figs.26, 27, and 30.
3. Although the locations of the energy-momentum density peaks will move around, in the long time driving limit (n 1), we can still find some regular patterns. In particular,   Fig. 26 (Fig. 27). The total number of driving steps are F 25 = 121393. This means there are in total 25 steps of moving (purple lines) at the Fibonacci numbers. (Bottom) Arg(ρζ) ∈ (−π, π] as a function of n, where n denotes the n-th Fibonacci numbers F n . The parameters are the same as the top panel, but with a larger driving number F 32 .
we find there is an even/odd effect in the distribution of these energy-momentum density peaks when the driving steps are Fibonacci numbers. More concretely, let us denote x peak as the locations of peaks of the energy-momentum density. Then we observe that as shown in Fig. 30.
We find that the even/odd effect as mentioned above is closely related to the group walking of ρ as discussed in the previous subsection, i.e., in the long time driving limit n → ∞, ρ will flow to a certain point ρ ∞ on ∂D. By taking this as an assumption, i.e, lim n→∞ ρ n = ρ ∞ ∈ ∂D, one can prove that there are indeed even/odd effects in the locations of energy-momentum density peaks, as follows.
Let us denote the matrix elements of Π F n−1 and Π Fn as follows Then applying the recurrence relation, we have The first formula in Eq. (157) indicates that According to Eq. (48), this implies that in the long time driving limit there is an even-odd effect in the locations of peaks of the energy momentum density at the Fibonacci numbers (See Eq. (154)).

Other quasi-periodic driving: Aubry-André like
The other well-known model for the quasi-crystal is the Aubry-André model [73,74], which describes free electron hopping on a one-dimensional lattice with the following Hamiltonian where ω is an irrational number that is incommensurate to the lattice periodicity and thus characterizes a quasi-periodic on-site potential. It has been used for the study of localization [75] and also appears in the two-dimensional integer quantum Hall effect [56]. Given an irrational ω, this model has a localization transition at λ = 2 and the two phases are related by a duality transformation. This is rigorously proved in the mathematics literature [62,75,76] by studying the so-called almost Mathieu operator, which is equivalent to the Aubry-André model at the single particle level. Motivated by the Aubry-André model, here we introduce another type of quasi-periodic driving, which we call an Aubry-André (quasi-periodically) driven CFT.

Setup
Let us first recall our minimal setup of the periodically driven CFT in Sec.4.4. Our protocol for the quasi-periodic driving will be a modification to that. In the minimal setup, the system takes the open boundary condition with q = 1. As depicted in Eq.
Note that T 0 and T 1 appear in M 0 and M 1 through cos and sin, which implies that T 0 → T 0 + L and T 1 → T 1 + L eff yield the same set of SU(1, 1) matrices 25 . This observation leads to the following way of introducing the quasi-periodicity. We keep T 1 fixed and let T 0 have a dependence on the number of cycles. In the n-th cycle, H 0 is applied for time T 0 = nωL, where L is the system size. The whole sequence is depicted in Fig. 31. If ω = p/q is a rational number, the unitary evolution generated by e −iH 0 T 0 will repeat after every q cycles, and protocol is reduced to a periodic driving with each period consisting of q cycles.
If ω is irrational, such periodicity disappears which gives rise to a quasi-periodic driving. One can also fix T 0 and vary T 1 instead, which does not change the result qualitatively and thus will not be discussed.

Phase structure for a single irrational number
It is subtle to directly access the dynamics with an irrational ω, especially when there is no analytical tool. For example, to extract the Lyapunov exponent by numerically computing the matrix product can be unstable due to numerical errors. Therefore, we will track a sequence of rational numbers to approach the physics at the irrational numbers. Each rational number corresponds to a periodic driving system and we can apply the trace classifier to identify the phase structure. This is exactly the same as our discussion for the Fibonacci driven CFT in Sec.5.1.3.
To illustrate the general features, let us consider ω = ( √ 5 − 1)/2 being the inverse golden ratio as a concrete example, which can be approximated by the sequence ω n = F n−1 /F n , n = 1, 2, 3 · · · . For a given n, the driving repeats after every F n cycles. We can compute the matrix that corresponds to the evolution for one period and determine the dynamics from the trace.
The phase structures for different values of n and θ are shown in Fig. 32, which shares a lot similarity with the Fibonacci case but also has some difference. Given n, the region for the non-heating phase gradually decreases as θ increases. Given θ, the non-heating region splits into many 'bands'. Notice that T 1 = 0 and T 1 = L eff are actually identified, the number of  Figure 32: Phase structure and Lyapunov exponents for the sequence of periodically driven CFT determined by ω n = F n−1 /F n . In (a) and (b), the blue (white) regions are the region where the system is non-heating (heating). In (c), for a given n, the blue lines (blank region) correspond to the non-heating (heating) phase.
bands is exactly equal to F n when F n is even and F n − 1 when F n is odd. If one compute the Lyapunov exponent for a fixed n, it shows peaks with equal spacing, as shown in Fig. 32(d). It implies that as ω → ( √ 5 − 1)/2, the self-similarity structure will not appear as in the Fibonacci case. The reason why the peaks have equal spacing will explained later.
We can also check explicitly the measure of the non-heating phase as approaching the quasi-periodic limit. Following the same prescription as (114), we call σ n (θ) the measure of the non-heating phase, which is a function of n and θ. The result is shown in Fig. 33(a). For a given θ, the measure is exponentially decaying at large n, similar to what we have found for the Fibonacci case. We introduce the decaying rate λ(θ) as σ n (θ) ∝ e −λ(θ)n , and find that it does not has a strong dependence on θ.
These two features also hold for generic rational and irrational numbers and can be summarized by the following statements.
1. If ω = p/q is a rational number with p, q being co-prime, then the region for the nonheating phase splits into multiple bands. The number of bands depends on q as Note the periodic boundary condition in T 1 when count the number of bands, namely the first band and last band (in vertical order) are considered to be the same band.

2.
When ω = a n /b n approaches an irrational number as n → ∞, the measure for the nonheating phase decreases exponentially with n. The second statement is only empirical and based on the numerical observation. The first statement can be understood by a perturbative argument as follows.
Perturbative proof of the statement 1. Let us assume ω = p/q and analyze the matrix of the unitary evolution for one period, which consists of q cycles. Notice that L and L eff are merely the units for T 0 and T 1 , we can set both of them to be 1 without changing the physics.
We first consider the limit θ = 0, so that M 0 and M 1 are both diagonal matrices of pure phases. For the k-th cycle of the period, T 0 = k p q and we have The lower index denotes the cycle and the upper index means it is the zero-th order term in the small θ expansion. The matrix for the whole period Π (2) · · · M (0) (q) and its trace can be easily computed, and we have where we have used the condition that p and q are co-prime. A typical curve for | Tr Π (0) q | is shown in Fig. 34(a). Therefore, even without adding θ, the trace can touch the critical value | Tr Π q | = 2 at the following positions q ∈ 2Z : q ∈ 2Z + 1 : T 1 = r q , r = 0, 1, · · · q − 1 . (162) Note that T 1 has a period 1, which fixes the range of r. Then, we turn on a tiny θ and consider its contribution to Tr Π q perturbatively. In particular, if | Tr Π q | exceeds 2 for a certain range of T 1 , a heating phase appears there and will continue to exist for larger θ. At small θ, such thing is more likely to happen at those special positions where | Tr Π that such thing does always happen. One typical example is shown in Fig. 34(b). For large q, we consider the following perturbative calculation. Each M (k) can be expanded to the second order in θ The first non-vanishing contribution comes from the second order of θ, which has two terms. One is the cross term of M (k) , the other only contains M (k) . In the limit of large q, second term dominates and we have Tr Π q ≈ 2 cos πqT 1 + π p(q + 1) 2 − 4qθ 2 sin πT 1 sin π(q − 1)T 1 + π p(q + 1) 2 .
Notably, this is a good approximation even for small q, the q = 5 case is shown in Fig. 34(b) as an illustration. If q is an odd number, one can check that at those special positions T 1 = r/q, we have from which we can see that the value of | Tr Π q | indeed exceeds 2 except at r = 0 (T 1 = 0). The proof for q being even is similar and one can find that | Tr Π q | exceeds 2 for all T 1 = (r + 1 2 )/q. So far, we have perturbatively shown that | Tr Π q | will exceed 2 at some equally spaced special positions. Therefore, the heating phase will appear there as long as one turns on θ, which accordingly opens 'gaps' in the phase structure. This gives rise to the multiple bands of the non-heating phase and also explains the equally spaced peaks observed for the Lyapunov exponent. Furthermore, when q is even, gaps will appear at all those special positions which leads to q disconnected bands. When q is odd, the gap cannot appear at T 1 = 0 which yields q − 1 bands (after identifying T 1 /L eff = 0 with T 1 /L eff = 1). This completes the proof of our first statement. We shift the origin of T 1 /L eff by 1/2 when presenting the data.

Phase diagram and nested structure
We can also study the phase diagram for generic ω including both rational and irrational numbers. The result is shown in Fig. 35, with the colored region being the non-heating phase and blank region being the heating-phase. Here are some comments on the features of the phase diagram: 1. The whole diagram is symmetric with respect to ω = 1/2, which is a direct consequence of our analysis above. Namely, for a rational number ω = p/q, the gap opening is independent of the numerator, which implies that the phase diagram should be invariant under p/q → 1 − p/q. Our following discussion will focus on the part ω 1/2.
2. The diagram has (infinitely) many empty regions, whose center is at simple rational numbers ω = 1/n, n ≥ 2. Some representatives are drew explicitly by the gray lines in Fig. 35. The reason why they all sit in the relative empty region is that one has to use a rational number with large denominator to approach one of them, the non-heating bands of which are too fragmented to read by eyes.
3. In between every two neighboring simple rationals, the structure of the subregion resembles the original phase diagram. This implies that the phase diagram has a nested structure, which is similar to the famous Hofstadter butterfly [56] and can be understood in the following way.
As reviewed in Appendix. B, every rational and irrational number ω ∈ [0, 1] can be uniquely represented by a continued fraction. Here we adopt the idea and consider a generalization We will see that it provides a useful guidance to resolve the diagram layer by layer. Those with only N 1 being nonzero are dubbed as the principal series, those with nonzero N 1 , N 2 as the first descendants and so on.  The nested structure of the phase diagram exactly follows such an organization: 1. First, the principal series {1/N 1 }, N 1 = 2, 3, · · · form the skeleton of the phase diagram, which is shown in Fig. 36(a). They also sit in the relatively empty region in the full phase diagram.
2. The first descendants { 1 } fill in the blank region between the principal series and serves as the skeleton for the next descendants. For example, 1 2+1/N 2 and 1 3−1/N 2 fills the region between 1/3 and 1/2, as shown in Fig. 36(b). 3. Notice that all the descendants of 1/2 that are smaller than 1/2 are 1 2+p/q , 1 < p < q. Similarly all the descendants of 1/3 that are larger than 1/3 are 1 3−p/q , 1 < p < q, which can also be written as 1 2+(1−p/q) . These two series, filling the subregion between ω = 1/3 and ω = 1/2, can be considered as the 'mirror reflection' to each other with ω = 2/5 being the reflection center. This is similar to what is observed for the full diagram, with the difference that each pair { 1 2+p/q , 1 2+(1−p/q) } do not have the same number of bands. This explains the similarity and also difference between Fig. 35(a) and (b). The nested structure follows from continuing such kind of game.

Features of the group walking
So far, the results resemble what have been found for the well known Aubry-André model or almost Mathieu operator. The unique advantage of our CFT setup is that it brings physical meaning to the group walking. As discussed in Sec. 3.2, the group walking of ρ and ρζ are related to the energy, energy-momentum density and the entanglement. We will give a corresponding discussion in this section.
To illustrate the general feature, let us choose the inverse golden ratio ω = In the numerics, we choose ω = F 29 /F 30 as an approximation to plot (a) and (b). In (b), the scattered dots in the circle represent (ρζ) n for n ≤ 500. The inset of (b) is Arg(ρζ)/π as a function n, with n denoting the n-th Fibonacci number F n . exponential growth respectively, which is consistent with our general claim for the heating phase.
The behavior of ρζ is quite different. In the short time, ρζ seems to have a random distribution in the disk. In the late time, it flows onto ∂D, which implies the formation of energy-momentum peaks. Although (ρζ) n and (ρζ) n+1 does not show a strong correlation, as shown by the background scattered dots in Fig. 37(b), the sub-sequence {(ρζ) k+Fn }, n ≥ 1 for any fixed k does have a definite limit as n → ∞ and the detailed value of the limit depends on k. Physically it means that the peaks observed at the time k + F n with fixed k will appear at the same position.
Such behavior of ρζ is generic as long as the irrational number is ω = ( √ r 2 + 4 − r)/2, r = 1, 2, · · · . For example, one will observe the same feature if choosing ω = √ 2 − 1, ( √ 13 − 4)/2 and so on. An intuitive reason is that the n-th principal convergent for this type of irrational numbers can be written as b n−1 /b n . (For the inverse golden ratio, b n is the n-th Fibonacci number.) Such a sequence b n provides us with a natural choice of the observation time. For a generic irrational number, its n-th principal convergent is p n /q n with {p n }, {q n } being two different sequence and a "natural choice" of the observation time becomes less clear.

Discussion
In this paper, we propose a general framework to study the non-equilibrium dynamics of (1+1)D CFTs with SL 2 deformation. We exam the details of the dynamical phases that emerge in the periodic and quasi-periodic driving using the tools we propose. In the sequel of this paper [57], we will apply the framework to the random driving sequence, where the use of the Lyapunov exponent and group walking becomes a necessity rather than a convenient option.
Here we highlight some of the unexpected features that we have found: 1. For the driving protocol that uses the SL 2 deformed Hamiltonian, when the total energy and entanglement start to grow, they always grow in a pattern where energy is concentrated in discrete points and form peaks; while the entanglement is shared within nearest neighbours. This phenomenon was first observed in Ref. [51] and is found to persist in the more general setting here.
2. Introducing irregularity in the driving protocol usually enlarges the heating phase, which is what we have observed in the quasi-periodic driving where the non-heating regime shrinks to a set with measure zero. 26 However, for the Fibonacci driving, the phase diagram is found to have a special non-heating fixed point, where the total energy and entanglement will return in a pattern following the Fibonacci sequence. More explicitly, if we observe the system only at the steps coinciding with the Fibonacci numbers, what we see is a state returning to itself with period 6.
It is worth to mention that if we pick out the two unitaries U A , U B that underlies the aforementioned Fibonacci driving, and apply them in a periodic fashion . ., the system will end up with a heating state. This is a surprise since it implies that at this special point, the "irregularity" actually converts the heating protocol to a non-heating one. The reason is that the pattern in the Fibonacci driving sequence manages to conspire with the special unitaries in a way that they happen to cancel each other and result in a return, which is explained in Sec. 5.2.4.
3. Besides the last point, the "order" of the quasi-periodic driving manifests itself in another way. In the heating phase of the quasi-periodic driving, ordinary stroboscopic observation of the energy peaks is featureless. However, if we consider the Fibonacci driving and observe the state at the steps coinciding with the Fibonacci numbers, the energy peaks oscillate between two fixed positions rather than randomly distribute (which is what we expect for a general irregular driving). Similar feature is also observed in the Aubry-André driving, where the energy peaks return to the same positions.
In the rest of this Discussion, we would like to comment on the special setting we use and some future directions. SL 2 deformation 27 is kind of a "shortcut" in analyzing the driven systems, because its effect in a single driving period can be characterized by a conformal transformation without introducing external sources to the system. The simplicity of the single driving allows us to pursue the "complexity" in the pattern of the driving sequence as we do in this paper. For future directions 1. Within the SL 2 deformation framework: so far we have been only probing the driven state by simple observables such as one-point function of energy-momentum tensor for energy distribution or two-point function of twist operator for the entanglement entropy, both of which only depend on the central charge. To explore more CFT data such as the operator content and the OPE coefficients, we need to consider more complicated observables. For example, we may consider measuring multi-point functions during the driving and ask how could a carefully designed driving protocol help us extract more CFT data.
2. Beyond the SL 2 deformation framework: the features we observe in this paper and in Ref. [51] is robust against perturbation in the initial state. However, it is not obvious what will happen when we perturb the driving Hamiltonian away from simple modulations, e.g. introducing operators other than energy-momentum in the driving Hamiltonian and ask their effects.
Another related question is that if we treat the CFT we have as a low energy effective theory, then in the heating phase we will finally drive the system to an energy scale where we need to consider its UV completion, i.e. we need to include some irrelevant operators in the driving Hamiltonian. We could ask what will happen at that point? For example, how does the energy peak and entanglement pattern get modified? These questions are relevant in explaining the data from lattice simulation beyond the conformal regime.
Note added: During the preparation of this manuscript, we learnt that the Fibonacci quasiperiodically driven CFT is also studied in [77], which will appear on arXiv on the same day. We thank the authors for sending us their manuscript before posting.

Acknowledgement
We thank for helpful discussions with Bo Han, Daniel Jafferis, Eslam Khalaf, Ivar Martin, Shinsei Ryu, Hassan Shapourian, Tsukasa Tada, Michael Widom, Jie-Qiang Wu, Yahui Zhang and Di Zhou. In particular, we thank Yahui Zhang for suggestions on the study of the quasiperiodic driving CFTs, and thank Michael Widom for pointing out Ref. [63], which stimulated our interest in considering the Fibonacci sequence in quasi-periodically driven CFTs. XW is supported by the Gordon and Betty Moore Foundations EPiQS initiative through Grant No.GBMF4303 at MIT. Y.G. is supported by the Gordon and Betty Moore Foundation EPiQS Initiative through Grant (GBMF-4306) and DOE grant, de-sc0019030. AV and RF are supported by the DARPA DRINQS program (award D18AC00033) and by a Simons Investigator Award.

A More on time-dependent driven CFTs
In this appendix, we give more details on some formulas/results as used in the main text.

A.1 Operator evolution with arbitrary SL 2 deformations
In this appendix, we introduce the procedures of obtaining the concrete form of Möbius transformation in Eq. (21) in the main text. Some related details can be found in Refs. [43,36,51,60,61]. Let us illustrate the calculation with a simple example, and then give results for an arbitrary SL 2 deformation.
The illustrative example we consider has the following deformed Hamiltonian: where T 00 (x) is the Hamiltonian density with T 00 (x) = 1 2π T (x) + T (x) .
As a remark, for θ = 0, H deform corresponds to a uniform one without any deformation; for θ = ∞, H deform corresponds to a SSD Hamiltonian, whose energy spectrum has been recently studied in detail in Refs. [40][41][42]44].
To evaluate the correlation function, such as the simplest one Ψ(t)|O(x)|Ψ(t) , where |Ψ(t) = e −iH deform t |Ψ 0 , one can study the operator evolution O(x, t) = e iH deform t O(x)e −iH deform t , as follows.
The correlation function Ψ(t)|O(x)|Ψ(t) can be considered as the path integral on a wcylinder with the operator O inserted, as depicted in Fig. 4, where w = τ + ix. This cylinder can be mapped to a q-sheet Riemann surface with a conformal map z = e where α = cos πt l eff + i cosh(2θ) · sin πt l eff and β = −i sinh(2θ) · sin πt l eff . When q = 1 such that l = L, we get the result as presented in Eq. (72) in the main text.

A.2 Entanglement entropy evolution
In this appendix, we give some details on the time evolution of the entanglement entropy in a time-dependent driven CFT.
With the same procedure, if one chooses the subsystem as A = [kl, (k + j)l] where k, j ∈ Z and j < q, then one can obtain S (m) log |α n + β n | + log |α n + β n | . The difference between this result and Eq.(180) reflects the non-uniform property of the driven CFT.

A.2.2 Linear decrease of the entanglement entropy
In this appendix, we show that if the subsystem is chosen in such a way that the two entanglement cuts lie on the centers of two chiral (anti-chiral) energy-momentum density peaks, and at the same time we keep the anti-chiral (chiral) part undriven, then the entanglement entropy may decreae in time. The choice of subsystem A can be understood as follows: Let us consider the setup of periodically driven CFT in Sec.4, such that the locations of peaks are fixed in the long time driving limit. Without loss of generality, we drive the chiral modes in time, but keep the anti-chiral modes undriven. As seen from Eq.(66), the distance between two chiral enerngy-momentum density peaks are x peak,2 − x peak,1 = jl, where j ∈ Z, j < q and l = L/q. In addition, since the locations of peaks correspond to the fixed point γ 2 (i.e., z peak = γ 2 , where we have assumed 0 < η < 1) in Eq.(54), this means x peak,1 and x peak,2 will not move around in the stroboscopic sense. Then based on Eqs.(177) and (178), one can find that Now by considering the expressions of α np and β np in Eq. (58), one can find that and therefore Ψ n |T m (w 1 , w 1 )T m (w 2 , w 2 )|Ψ n = 2π 4h , based on which one can find that where 0 < η < 1. Recalling that the total driving step number is N = np, we can write the entanglement entropy evolution in terms of the Lyapunov exponent in Eq.(60) as follows: It is interesting to compare this formula with the result in Eq. (68). Here the linearly decreasing entanglement entropy is due to the coincidence of the entanglement cuts with the centers of the chiral energy-momentum density peaks (while keeping the anti-chiral modes undriven). During the driving, the degrees of freedom that entangle A and A will flow and accumulate at the energy-momentum density peaks, which are located at the entanglement cut. Intuitively, the Bell pairs that are non-local in space now become local, which results in a decrease in the entanglement entropy evolution. It is emphasized that although the entanglement entropy decreases in time, which is due to the choice of entanglement cuts, the total energy of the system still grows in time (in the heating phase) in a periodically driven CFT.

A.2.3 Comparison of CFT and lattice calculations
To confirm the linearly decreasing feature of the entanglement entropy evolution in the previous subsections, we compare the CFT and lattice calculations in this appendix. The results here are also related to the entanglement evolution at the non-heating fixed point in a Fibonacci driven CFT in Sec. 5.2.2 (See Fig. 22).
As an example, let us consider the minimal setup of periodically driven CFTs in Sec.4.4. We drive the CFT periodically with H θ=0 and H θ in Eq. (71), with open boundary conditions. Now we consider two driving protocols: Protocol I: Driving with H 0 first, and then H 1 = H θ , Protocol II: Driving with H 1 = H θ first, and then H 0 , The driving time intervals are chosen as T 0 = L/2 for H 0 , and T 1 = L eff /2 for H 1 , respectively. The resulting Möbius transformations in one driving period correspond to the following SU (1, 1) matrices: One can find |Tr(M )| = |Tr(M )| = |2 cosh(2θ)| > 2, i.e., both M and M are hyperbolic. Therefore, in both protocols, the CFT is in a heating phase. The difference is that, the chiral and anti-chiral energy-momentum density peaks are located separately at x = 0 and x = L in protocol I, but are located at the same position x = L/2 in protocol II, as pictorially shown as follows: x 1 Protocol I x 1 Protocol II (188) By choosing the subsystem A = [0, x 1 ] with x 1 = L/2, the entanglement cut will not cut any peaks in protocol I, but will cut both the chiral and anti-chiral peaks in protocol II. Intuitively, one can unfold the CFT with open boudary conditions to a single copy of chiral CFT with a periodic boundary condition. After the unfolding, the entanglement cut in protocol II lies on the centers of two chiral energy-moemntum density peaks. Based on our study in the previous subsection, the entanglement entropy will decrease linearly in this case. Next, we show this is indeed the case rigorously .
Based on Eqs. (31) and (33), one can obtain the entanglement entropy/total energy evolution as follows: where A = [0, L/2]. One can find that the entanglement entropy grows (decreases) linearly as a function of n in Protocol I (II). On the other hand, the total energy of the system grows in both protocols: · cosh(4nθ), for both Protocols I and II. Figure 39: Phase diagrams in a periodically driven CFT with the sequence generated by finitely truncated Fibonacci bitstring, i.e. {X j } with ω n = F n−1 /F n . Here we choose n = 2, 4, 5, 6, 10, 20, 100, and 1000, respectively. The two Hamiltonians we use are H 0 (θ = 0) and H 1 (θ = 0.2) in (71). The phase diagram is periodic in T 0 direction with period L and in T 1 direction with period L cosh(2θ) 1.08L. The blue (yellow) regions correspond to the heating (non-heating) phases.
Now we compare the CFT and lattice calculations on the entanglement entropy/total energy evolution. The lattice model we consider is the same as that in Sec.4.4.2. That is, the two lattice Hamiltonians under consideration are H 0 = 1 where L is the total length of the lattice and f (i) = 1 − tanh(2θ) · cos 2πi L , with the initial state chosen as the ground state of H 0 . The comparison is shown in Fig.38. The agreement between CFT and lattice calculations are remarkable.

A.3 Phase diagrams from periodic to quasi-periodical driving
In this appendix, we present one more group of results on the evolution of phase diagrams as we use a periodical driving to approach the Fibonacci quasi-periodical driving limit, as shown in Fig. 39 and Fig. 40.
In Fig.39, we consider a periodic driving with the sequence generated by finitely truncated Fibonacci bitstring, i.e. {X j } with ω n = F n−1 /F n . The two driving Hamiltonians are H θ=0 and H θ=0.2 (See Sec.5.1.2 for more details).
In Fig.40, we replot the phase diagrams in Fig.39 with the new variables as defined in Eq.(116). Based on Fig.40, we obtain the measure of non-heating phases in Fig.17. Figure 40: Phase diagrams in a periodically driven CFT with the sequence generated by finitely truncated Fibonacci bitstring, i.e. {X j } with ω n = F n−1 /F n . Here we choose n = 2, 4, 5, 6, 10, 20, 100, and 1000, respectively. The two Hamiltonians we use are H 0 (θ = 0) and H 1 (θ = 0.2) in (71). The phase diagram is periodic in T 0 direction with period L and in T 1 direction with period L cosh(2θ) 1.08L. The blue (yellow) regions correspond to the heating (non-heating) phases.

B Fibonacci bitstring/word and recurrence relation
In this appendix, we review some basics of the Fibonacci bitstring/word (in this paper, we will call it Fibonacci bitstring instead of Fibonacci word 28 ) for the readers' convenience. In particular, we explain the equivalence of two ways to generate the Fibonacci bitstring: (1) the quasi-periodic potential and (2) the substitution rule. Based on the substitution rule, we explain the recurrence relation of the traces and the constant of motion that are used in the main text.
This definition is straightforward but not useful in our application. Instead, we will follow the presentation in Ref. [79] to show that the above bitstring can be generated by a substitution rule. The equivalence is based on the following two properties of {X j } 1. Let us use the notation " x := max{m ∈ Z|m x}" for the floor function, and we have To prove this statement, let us first use the relation ω 2 + ω − 1 = 0 to rewrite ω 2 = 1 − ω and −ω 3 = 1 − 2ω. Then according to the rule (192) we have X j = 1 ⇔ ∃m ∈ Z : m − 2ω (j − 1)ω < m − ω ⇔ ∃m ∈ Z : jω < m (j + 1)ω 2. Let us use F n to denote the n-th Fibonacci number, namely F n = F n−1 + F n−2 with F 0 = F 1 = 1, then we have X j+Fn = X j , for n 2 , 1 j < F n .
To prove this statement, it is sufficient to show that for n 2, the difference of two parts in X j = (j + 1)ω − jω is unchanged while shifting the argument of the floor function · by F n ω for 1 j < F n : we can write F n ω = m + r where m is the integer that is closest to F n ω and |r| = dist(F n ω, Z) denotes the distance (with sign) between F n ω and the nearest integer m. Obviously, shifting by an integer will not affect the difference, so we only need to check the effect of shifting by r. From the fact that the convergent 29 F n−1 /F n is a best Diophantine approximation 30 of the irrational number ω, we have the and rename it as Fibonacci bitstring to avoid confusion. 29 For irrational real number x, we always have an infinite continued fraction representation x = a 0 + 1 a 1 + The n-th principal convergent is the rational number p n /q n obtained by a truncation at a n . The irrational number ω = has a particularly simple continued fraction representation a 0 = 0, a 1 = a 2 = . . . = 1, and its finite truncation is the ratio of two nearby Fibonacci numbers F n−1 /F n as one can easily check. 30 For a real number x, a rational number p/q is a best approximation (of second kind) if |qx − p| < |q x − p | , ∀q > q > 0 .
For the proof that the best approximations are given by the convergent, see e.g. Ref. [80].
following inequality |r| = dist(F n ω, Z) < dist(jω, Z) for 1 j < F n (199) in other words, for jω with 1 j < F n , adding r will not be able to fill the gap between jω and a nearby integer. This statement further holds for j = F n , since ω is irrational, F n ω + r can not be an integer. To summarize, neither (j + 1)ω nor jω will change its value after a shift of r for 1 j < F n and therefore we have proved (196).
The second property (196) provides an efficient algorithm to generator the 0, 1 bitstring for {X j }. For instance, let us denote the first F n bits as string A n B n where A n stands for the first F n−1 bits and B n stands for the next F n−2 , and they together has length F n = F n−1 + F n−2 as required, here are examples for first few n Then, we group the string A n B n , rename it as A n+1 = A n B n and according to (196), the B n+1 is obtained by copying the first F n−1 bits of A n+1 which is exactly A n , i.e. we have the following recurrence relation A n+1 = A n B n , B n+1 = A n .
One may concern about that (196) actually produces a longer bitstring than the above recurrence relation since it also generates an additional segment B n after A n+1 B n+1 which will overlap the first F n−2 bits of B n+2 . One can check that the overlapping part is consistent with the rule here, because the first F n−2 bits of B n+2 is indeed A n−1 = B n . A final comment is that the above recurrence relation can be recast into a "local" substitution rule that is closer to the rabbit populations problem Fibonacci originally considered. Let us start with a single bit 1, and apply the following substitution rule 1 → 10 , 0 → 1 (202) at each step, then we will generate the following sequence 1 → 10 → 101 → 10110 → 10110101 → . . .
which approaches to the Fibonacci bitstring after infinite steps.

B.2 Recurrence relation and constant of motion
In the main text, we are interested in the product Π n of n SU(1, 1) matrices where matrix M j depends on X j in the Fibonacci bitstring. The substitution rule (201) directly leads to the following recurrence relation for Π n The relation can be extended to k = 2 by defining an auxiliary Π F 0 = M 2 that is distinct from Π F 1 = M 1 , although strictly speaking F 0 = F 1 = 1. A key observation made in Ref. [81] is that their traces obey the following recurrence relation To derive this relation, we start with (205) and find Then we insert identity Π F k + Π −1 F k = Tr(Π F k )I for the uni-determinant 2 × 2 matrice Π F k , and obtain whose trace gives (206).
Using the trace relation, Ref. [81] further notes a constant of motion Indeed, one can check that the change of r.h.s. under the shifting k → k + 1 is zero: