Emergent Spatial Structure and Entanglement Localization in Floquet Conformal Field Theory

We study the energy and entanglement dynamics of $(1+1)$D conformal field theories (CFTs) under a Floquet drive with the sine-square deformed (SSD) Hamiltonian. Previous work has shown this model supports both a non-heating and a heating phase. Here we analytically establish several robust and `super-universal' features of the heating phase which rely on conformal invariance but not on the details of the CFT involved. First, we show the energy density is concentrated in two peaks in real space, a chiral and anti-chiral peak, which leads to an exponential growth in the total energy. The peak locations are set by fixed points of the M\"obius transformation. Second, all of the quantum entanglement is shared between these two peaks. In each driving period, a number of Bell pairs are generated, with one member pumped to the chiral peak, and the other member pumped to the anti-chiral peak. These Bell pairs are localized and accumulate at these two peaks, and can serve as a source of quantum entanglement. Third, in both the heating and non-heating phases we find that the total energy is related to the half system entanglement entropy by a simple relation $E(t)\propto c \exp \left( \frac{6}{c}S(t) \right)$ with $c$ being the central charge. In addition, we show that the non-heating phase, in which the energy and entanglement oscillate in time, is unstable to small fluctuations of the driving frequency in contrast to the heating phase. Finally, we point out an analogy to the periodically driven harmonic oscillator which allows us to understand global features of the phases, and introduce a quasiparticle picture to explain the spatial structure, which can be generalized to setups beyond the SSD construction.

However, a CFT as a gapless many-body system is expected to be vulnerable to a generic driving. If we start from the ground state of the original Hamiltonian, then Floquet driving might lead it to an infinite temperature state easily. To avoid that happening, we need to choose special protocols. In this paper, we are going to use the Virasoro symmetry generators as our driving Hamiltonian so that we can take maximal advantage of the conformal symmetry to constrain the system. As one of the most canonical choices, we will use the sl(2, R) subalgebra, the exact meaning of which will be discussed later. Although this choice may look a bit special, it is powerful enough to reveal some universal features of the problem that apply more generally.
We will follow the setup used in [34], where the authors consider an open chain and implement the driving with the sine-square deformed Hamiltonian [37][38][39][40][41][42][43][44][45][46][47]. It was shown that if we start from the ground state and turn on the Floquet drive, we can identify a non-heating phase in the high-frequency driving regime and a heating phase in the low-frequency driving regime by looking at the entanglement entropy growth. The fact that we have these two phases has an algebraic reason which can be understood by using a quantum mechanical model, as we will discuss later.
The main part of this paper will present a more detailed study on what happens in the Floquet dynamics, paying special attention to the spatial structure that emerges, that has not previously been discussed.
Let us summarize the main phenomena. In the heating phase, although the total energy and entanglement keep growing, the system does not evolve into a featureless state. We find that in this phase, the system heats up in a very non-uniform way. The energy pumped in concentrates at two points, one of which has purely chiral excitations and the other one only has anti-chiral excitations. The entanglement entropy also comes from the entanglement between the excitations at these two points. Furthermore, the energy and entanglement entropy are related by a simple formula. All these features above are universal and only depend on the central charge of the CFT. In the non-heating phase, if we do a stroboscopic measurement, we can find that energy excitation will move back and forth in the system with the total energy and entanglement entropy oscillating in time.
Since a real experiment will inevitably have noise, we are also interested in the question that how stable those phenomena are to noise. For example, we could start from an excited state or have local perturbation during evolution. Furthermore, the driving frequency could have a small fluctuation. By combining analytical and numerical analysis, we will argue that the non-heating phase is delicate but all the reported features in the heating phase are quite robust to these perturbations. For the non-heating phase, an arbitrarily tiny noise in the driving frequency will eventually heat the system. The dimensionless heating rate is proportional to α 2 , where α characterizes the magnitude of randomness.
The paper is organized as follows. In Sec.2, we will briefly review the set-up in [34], summarize the method of studying the evolution of operators and see how to interpret the Hamiltonian by sl(2, R) algebra. In particular, in Sec.2.3, we discuss the phase diagram from a different angle using the mapping to a driven harmonic oscillator. In Sec.3, we will present our main result of this paper on various features of the heating phase. We will focus on how the energy is absorbed, how the entanglement is generated and their relation. We will also draw intuition from this special setup and make a few comments on what to expect for the case of more general initial conditions, boundary conditions and Floquet drives. In Sec.4, we will analyze the stability of these phenomena against driving with random periods. Finally in Sec.5, we give some conclusions and outlook.

Setup for a Floquet CFT
In this section, in the interest of completeness, we review the set-up and some results of prior work in [34] that are relevant to our discussion.

Floquet driving
We start with a (1 + 1)D CFT with an open boundary condition. Let us denote its total length by L, and its central charge by c. We will consider the following time-dependent Hamiltonian where H 0 is the ordinary Hamiltonian that can be written as an integral of energy density T 00 (x) along the real space as follows, H 1 is the so-called sine-square deformed (SSD) Hamiltonian [37][38][39][40][41][42][43][44][45][46][47][48] For simplicity, the initial state |ψ 0 is chosen to be the ground state of H 0 , i.e. |ψ 0 = |GS . It is also useful to introduce the Floquet operator F = F 0 F 1 = e −iH 0 T 0 e −iH 1 T 1 to characterize the unitary evolution for a single cycle. In the stroboscopic measurement, the Floquet dynamics is determined by the state |ψ(nT ) = F n |GS . For example, the two point function of local operators O 1 (x 1 ) and O 2 (x 2 ) after n cycles is given by In the "Heisenberg" picture, the calculation amounts to determining the operator evolution O(x, nT ) = F −n O(x)F n . For general Floquet drives, this is a difficult problem. However, for the SSD Hamiltonian defined in (3), the operator evolution has a simple expression in terms of the Möbius transformation.

Operator evolution and Mobius transformation
In this section, we will derive the explicit expressions for the operator evolution. It is convenient to work in Euclidean coordinates, and the Lorentzian correlator can be obtained by analytic continuation. We will use three coordinates in this paper, denoted by They correspond to the stripe geometry, complex plane and upper half-plane respectively. See Fig. 1 for an illustration. In the imaginary time, the Floquet operator is given by F = e −τ 0 H 0 e −τ 1 H 1 . Let us first check how the operator evolves after one cycle, namely Here we assume O(w, w) to be a primary operator with the conformal dimension (h, h). On the strip w, the algebraic relations between H 0,1 and O are complicated. It will be easier to work in z = e 2πw/L coordinate instead, where H 0,1 are expressible as contour integrals of the stress tensor. More explicitly, we have The term −cπ/6L comes from the Schwarzian derivative and will not affect the operator evolution. The contour C is shown in Fig. 2 (a). The subtlety is that C is not closed due to the branch cut arising from the open boundary condition. The branch cut can be treated as follows. First, we use the Baker-Campbell-Hausdorff formula to expand Eq. (5) and write it as commutators. Each term can be depicted as a double contour integral shown in Fig. 2 (b). The conformal boundary condition requires T (z) = T (z) right above and below the branch cut, respectively. Thus, we can attach two horizontal lines along with the branch cut (as the red horizontal lines in Fig. 2(b)) for free since the contributions exactly cancel. After the above manipulations, the new contour can be deformed to enclose operator O as shown in Fig. 2(c). Therefore, on the z-plane, the Floquet operator acts on the operators as if there is no branch cut.
As a consequence, the operator evolution that is driven by the stress tensor will be determined by a two-step conformal transformation where (∂z/∂w) corresponds to the transformation from the strip (w) to the complex plane (z) and (∂z 1 /∂z) is the transformation generated by the Floquet dynamics F . To determine the map z 1 (z), we notice that without the branch cut, H 0 and H 1 in Eq. (6) can be written as Virasoro generators L 0,±1 and their anti-holomorphic patterns as follows, we use tilde to emphasize that the identification only works for operator evolution. The generators L 0,±1 form an sl(2, R) algebra. Therefore, the corresponding Floquet operator F generates a Möbius transformation on z, namely The coefficients a, b, c, d are determined by the dimensionless driving periods τ 0 /L and τ 1 /L as follows, More explicitly, the evolution induced by H 0 acts as a dilation on the z-plane, namely z goes to z = e 2πτ 0 /L z, which explains the e πτ 0 /L factors. The evolution by H 1 is also a dilation but in a different coordinate χ. 1 For the Floquet problem, we would like to study the operator evolution for n repeated cycles of Möbius transformations, namely z n = f (f . . . f (z)) and we will denote it as, 1 Since the H 1 acts on z coordinate in a complicated way, we can instead look for a new coordinate χ, on which H 1 acts as a simple dilation. Namely we assume a coordinate change χ(z) and accordingly T (z) = χ 2 T (χ), Requiring H 1 generates a dilation amounts to the following condition, Under the evolution of H 1 , χ goes to χe 2πτ1/L and correspondingly the z transforms as, Inserting z = e 2πτ0/L z, we get Eq. (10).  The successive application of Möbius transformation is better described using the fixed points f (γ) = γ and the "rotations" η relative to the fixed points With these new variables, Eq. (14) can be rearranged into the following form.
For our physical application, c = πτ 1 L e πτ 0 /L is non-zero and there are in general three possible scenarios depending on the position of the fixed point: 1. Elliptic class: the quadratic equation f (γ) = γ has two distinct roots that are conjugate to each other γ 1 = γ * 2 , the rotation parameter is determined by the following formula In this case, η is a pure phase, namely |η| = 1. See Fig. 3 (a) for an illustration. The corresponding Mobius transformation z n = f n (z) can be represented as an SL(2, R) matrix as follows, 2. Hyperbolic class: the two distinct roots are purely real and the parameter η defined above is also a real number. See Fig. 3 (c) for an illustration. The parameter 0 < η < 1 represents the rescaling near the fixed points. The Mobius transformation matrix is in the same form as (18).
3. Parabolic class: two roots are merged together γ 1 = γ 2 = γ. Therefore (16) does not apply. For this case, we introduce a new parameter β = a−d 2c such that 1 The corresponding transformation matrix is which can not be diagonalized. The parabolic class may be thought as the marginal case of either elliptic or hyperbolic class, see Fig. 3 (b) for an illustration.
We remark here that the Möbius transformation also applies to quasi-primaries such as the stress tensor. In that case, although we will obtain a Schwarzian derivative term when transforming between different geometries, the operator evolution driven by Möbious transformation on the complex plane is still determined by Eq. (7). More explicitly, the stress tensor on the strip after n-cycle driving becomes Finally, for the operator evolution in real (Lorentzian) time, we perform the analytic continuation τ 0 → iT 0 , τ 1 → iT 1 . In real time, a space-time position (x, t) on the strip maps to z = e i2π(x+t)/L on the z-plane, which is always on the unit circle. Therefore, the operator evolution is geometrically related to the automorphism of a unit circle under the conformal mapping. Although the Möbius transformations after the analytic continuation generally belong to SL(2, C), the basic structure remains the same. (Naively we may expect an additional class known as loxodromic class shows up where η is a general complex number, not necessarily a phase or purely real. However, the physical parameters that appear in the Floquet setting do not fall into such class.)

Parametric oscillator (swing) analogy
The last section explained the relation between the operator evolution and the Mobius transformation, which is further classified into three classes: elliptic, hyperbolic and parabolic. The corresponding Floquet dynamics are also classified into the non-heating, heating, and the critical classes respectively, and the phase diagram was first obtained in [34]. For reader's convenience, we reproduce the phase diagram in Appendix. A.
It is instructive and amusing to gain intuition into this classification in a more elementary setting with the same SL(2, R) structure. The example we would like to use is the parametric oscillator with the following Hamiltonian [49], where f (t) = f (t + T ) and g(t) = g(t + T ) are periodic functions. One familiar example is the Mathieu oscillator, which corresponds to f (t) = 1, g(t) = g 0 − 2g 1 cos(2t). Classically, they are useful in explaining the motion of a playground swing, see Fig. 4(a) for a classical picture. Furthermore, the recognition of the SL(2, R) structure in the problem also has interesting consequence in the ultra-cold quantum gases, e.g. see Ref. [50,51].  For the quadratic Hamiltonian, the Heisenberg operators ( Therefore the stroboscopic evolution of (x, p) is represented by a SL(2, R) transformation F (x,p) , whose classification determines the stroboscopic trajectory of (x(nT ), p(nT )). More explicitly, to compare with the Möbius transformation used in Eq. (9) we treat (x, p) as a point on the complex projective plane CP 1 which can be more conveniently parametrized by z = x/p. Then the SL(2, R) action on the point (x, p) shown in Eq. (23) is equivalent to the Möbius transformation Eq. (9) and also have three classes. The fixed points γ 1,2 and the rotation angle η of the Möbius transformation can be translated to the eigenvectors v 1,2 and the ratio of eigenvalues λ 1,2 of the SL(2, R) matrix, respectively. Their correspondence is given explicitly below.

Möbius transformation
This explains the different dynamics. In the elliptic class, (x, p) as a real vector only keeps rotating on the x − p plane. The energy measured by p 2 /2 + x 2 /2 just oscillates E(nT ) ∼ cos(nθ + ϕ) with a period controlled by the angle θ of the eigenvalue. In the hyperbolic class, F (x,p) has two real right eigenvectors: v 1 with an eigenvalue λ > 1 and v 2 with an eigenvalue 1/λ < 1. Therefore unless the initial condition (x 0 , p 0 ) is along v 2 , (x(nT ), p(nT )) will flow to infinity along v 1 exponentially fast, which causes the energy to grow exponentially in the long time limit In the parabolic class, F (x,p) only has one right eigenvector v with eigenvalue 1 thus becomes singular. To determine the dynamics, we can look at the Jordan normal form of F (x,p) . 3 Unless the initial condition (x 0 , p 0 ) is along v, (x(nT ), p(nT )) will flow to infinity linearly, which causes the energy to grow quadratically As a concrete example, the Mathieu oscillator introduced at the beginning of this section can support all of the three different dynamics, and its phase diagram is presented in Fig. 4. The Floquet CFT studied in the current paper is richer than its oscillator analog. In particular, the (1 + 1)D CFT has locality in space, which will lead to features in the energy density and entanglement that are the focus of the following sections.

Energy and Entanglement
Energy and entanglement are the most straightforward and fundamental diagnostics of states evolving under Floquet driving. Fortunately, both can be studied in (1 + 1)D CFT analytically using the operator evolution method we have discussed. In this section, we will present the stroboscopic measurement of the energy and entanglement under Floquet driving. We will also provide a semi-classical picture of the phenomenon and point out an interesting relation between energy and entanglement.

Energy density and total energy
The energy of a state under the Floquet evolution can be measured by the expectation value of the stress tensor T 00 = T + T , whose time evolution can be obtained by Eq. (21). The n dependence of the energy arises from the first term. To compute G|T (z n )|G , we need to perform another conformal transformation to the upper half-plane via ξ = √ z. This mapping generates a Schwarzian derivative 3 Since F (x,p) is singular, its Jordan normal form has a nonzero off-diagonal element The off-diagonal element will increase linearly with the driving cycles, i.e. F n (x,p) = P (c) Heating phase Figure 5: The evolution of energy density profile in different phases. In the plot, we choose the system size L = 2π and the central charge c = 1. Different colors means different times. The green, blue, brown and red curve corresponds to t = T, 2T, 3T, 4T respectively. (a) T 0 = 0.5L, T 1 = 0.1L which is in the non-heating phase. The energy density only oscillates (b) T 0 = 0.9L, and T 1 is tuned to make the system right at the phase boundary. The position of the peaks in the plot are not given by γ. This is because we are not at the late time regime. One can check that as we increase n, the peaks will move towards log γ. (c) T 0 = 0.9L, T 1 = 0.1L which is in the heating phase. We can clearly see the formation and growth of two peaks. and leaves a second term G|T (ξ)|G . On the one hand, Ward identity and scale invariance constrains G|T (ξ)|G ∝ 1/ξ 2 . On the other hand, it is invariant under the horizontal translation. Therefore this term has to vanish and all the contribution comes from the Schwarzian term, namely, where w = τ + ix is the complex coordinate for the stress tensor T on the strip, z n is the coordinate on the z-plane after n-cycle driving and c is the central charge. After analytic continuation τ 0 → iT 0 , τ 1 → iT 1 , the expectation value of T has the following form with A, B, C, D depending on T 0 /L, T 1 /L and n through the prescription described in section 2.2. Replacing z with z gives us the expectation value of T . These formulae allow us to look at the evolution of energy density directly, which is found to have different behaviors in different phases, as shown in Fig. 5. In the non-heating phase, the energy density just fluctuates without a definite period. In the heating phase, the energy density quickly develops two peaks, which grows with time vary fast. The positions of the energy peaks are determined by the unstable fixed points of the Möbius transformation, i.e. e 2πix peak /L = γ 2 or γ * 2 . At the phase boundary, there are also two energy peaks but growing much slower.
These phenomena can be understood from the perspective of the fixed points of Möbius transformation. The n dependence enters Eq.  rescaling factor (∂z n /∂z) 2 , whose different behaviors in three phases explain the feature shown in Fig. 5.
1. Non-heating phase: The two fixed points sit on different sides of the unit circle. In our convention, γ 1 is inside the unit circle and γ 2 is outside the unit circle, as depicted in Fig. 6 (a). Since z n is constrained on the unit circle, it cannot flow to either of them but just keeps rotating around them. That is the reason that energy density fluctuates in this phase. Since the rotation angle η, defined by Eq. (17), is not a rational phase, these fluctuations do not have a definite period.
2. Heating phase: Both two fixed points are now on the unit circle. In our convention, γ 1 is a stable fixed point and γ 2 is an unstable fixed point. For the chiral stress tensor, when z = γ 2 , although z n doesn't change the rescaling factor (∂z n /∂z) 2 = η −2n will grow exponentially with n. For the anti-chiral stress tensor, the same thing happens at z = γ 2 . Therefore we observe two energy peaks at two symmetric positions. On the other hand, for a generic position z, z = γ 2 , z n will flow to the stable fixed point and the rescaling factor (∂z n /∂z) 2 will decrease exponentially with n. Therefore the stress tensor shrinks, making the two peaks sharper and sharper.
3. Critical line: The two fixed points merge to γ on the unit circle, which is a marginal case. One can show that the maximal value of the rescaling factor keeps growing but in a power-law fashion, which explains the slowly growing peaks. The position for the maximum gradually moves to the position corresponding to γ.
Note the phenomena here do not rely on the initial state, as long as it is not a common eigenstate of H 0 and H 1 . For example we may consider a generic initial state |φ with the expectation value of stress tensor φ|T (z)|φ = E φ (z), then the Eq. (27) generalizes to and the discussions above still hold. In this scenario, the boundary condition is also irrelevant since the operator evolution discussed in Sec. 2.2 is independent of the choice of boundary conditions. 4 The idea of relating the fixed points to the heating/non-heating phenomena also applies to more general setups. For example, we can use H 0 and H 2 = 2 L 0 dx sin 2 2πx L T tt (x) to generate the Floquet dynamics. H 2 is related to L ±2 and thus the operator evolution is still a conformal transformation but with four fixed points. When none of the fixed points are on the unit circle, the system is in the non-heating phase without energy peaks. In certain parameter regime, there are two unstable fixed points locating on the unit circle, which implies heating dynamics and correspondingly four growing energy peaks. Furthermore, we can define a Hamiltonian by a generic deformation H = L 0 dxf (x)T tt (x). As long as f (x) is a smooth real function and has a Fourier decomposition, H can be represented as a linear combination of Virasoro generators and the operator evolution can be written as a conformal transformation. 5 If f (x) = sin 2 kπx L , k ≥ 1, the conformal mapping is essentially the same as what we discussed here, which supports a nonheating and heating phase. However, for a generic f (x), determination of fixed points and the corresponding dynamics is a hard problem, which we leave for a future study.
Besides the energy density, we can also look at the total energy For stroboscopic measurement, we can plug in the Eq. (27) and have In either non-heating or heating phase, the Möbius transformation has two fixed points and we need to use Eq. (18) to get, For the non-heating phase, since η is a pure phase the total energy will oscillate with time. Generally, η is a non-rational phase factor, thus we do not expect any periodicity. Since the energy is oscillating, we cannot talk about the long-time behavior itself but the average, which is a finite number. For the heating phase, 0 < η < 1 and η n becomes exponentially small at the late time regime. Therefore, to the leading order we can drop the η n and η 2n terms in  the numerator and find the energy grows exponentially, When the system is at the boundary between the non-heating and heating phase, the two fixed point merges together and we need to use Eq. (20). Noting that βγ is pure imaginary, the total energy can be written as The total energy grows quadratically in cycle number n. This long-time asymptotics of the total energy provides a direct diagnostic of the different phases. The oscillation, exponential and quadratic growth behavior matches the simple picture obtained in the driven harmonic oscillators, as shown in Sec. 2.3. In particular, noticing that η = λ −2 , the heating rate in Eq. (24) and Eq. (34) are exactly the same. This is because the growth behavior only depends on the underlying algebra and not its detailed realization.

Entanglement pattern in the heating phase
Besides the total energy, the entanglement entropy of the left/right half system also has different behaviors in different phases, which was shown in [34]. Here, we will focus on the spatial structure of entanglement that has not previously been discussed. in particular we examine the heating phase and discuss the relation between the energy peaks observed above and the entanglelment.
Given a pure state |Ψ , the reduced density matrix of a subsystem A is defined by the partial trace ρ A = Tr A |Ψ Ψ| and its von Neumann entanglement entropy is given as S A = − Tr ρ log ρ. In our setting, we consider a time dependent state |ψ(nT ) = F n |GS that evolves under the Floquet driving and study the corresponding entanglement entropy S A (nT ) as a function of driving cycle n. Figure 8: Entanglement cuts for different cases. In (a) and (c), the subsystem includes only one energy peak. In (b), the subsystem doesn't include any energy peak.
For a subsystem A = [0, x] starting from the left end and end at position x ∈ (0, L), we plot the results in Fig. 7 and keep the details of the calculation in Appendix. B. The entanglement entropy has a background value from the initial state. As time increases, the curve quickly develops two kinks, the positions of which exactly coincide with the energy peaks. Only the curve between the two kinks grows with time while the curve outside does not. This implies only when the subsystem includes one of the energy peaks, does the entanglement grow with time. If the subsystem includes either none or both peaks, the entanglement remains at its background value and does not grow at all.
This statement can be further verified by studying the entanglement of the subsystem A with ending points x 1 , x 2 ∈ (0, L). We fix x 2 to sit between the two energy peaks and study how the long-time behavior of entanglement growth depends on the choice of x 1 . Without loss of generality, we assume the chiral energy peak is on the left and the anti-chiral peak is on the right in the following discussion. In general, there are three different choices of x 1 : 1. 0 < x 1 < x C . In this case, the subsystem A only includes the chiral peak, as depicted in Fig. 8 (a). The entanglement entropy is, where the first term grows linearly with time (i.e. the driving cycle n), which is consistent with the result in [34]. As long as x 1 < x C , the slope only depends on the central charge and the characteristic constant η but not on the positions of entanglement cuts. This behavior is universal and does not depend on the operator content. The non-universal terms are sub-leading in the n 1 limit.
2. x C < x 1 , x 2 < x A . In this case, the subsystem is between the chiral and anti-chiral peak, as depicted in Fig. 8 (b). To the leading order, one can show that it saturates to an O(1) value, which depends on the operator content and position of insertion. The exact value is not relevant but the most important is that the entanglement entropy does not have interesting time dependence in the long time limit.
3. x A < x 1 . In this case, the subsystem A only includes the anti-chiral peak, as depicted in Fig. 8 (c). The entanglement entropy grows linearly as in Eq. (36).
Using the results above, we can also infer the bipartite mutual information between the two energy peaks. Let us choose two disjoint regions X and Y . X only includes the left peak and Y only includes the right peak. We call their complement as Z, which is composed of three disjoint regions Z 1 , Z 2 and Z 3 , located on the left of the chiral peak, between the two peaks and on the right of the anti-chiral peak, respectively. Since we are studying a pure state, the mutual information between X and Y is Following the prescriptions in Appendix. B, the calculation of S XY requires a four-point function of twist operators in the boundary CFT which is in general unknown. Instead, we can use the subadditivity of entanglement to bound S Z where in the last step we use the fact that none of S Z j , j = 1, 2, 3 grows with time and saturates to an O(1) value. Therefore S Z itself can only be O(1) value and does not make an important contribution to the time dependence of the mutual information. On the other hand, since each of X and Y includes one peak, S X and S Y grows linearly with time as shown in Eq. (36). Thus the mutual information also linearly grows with time, where the non-universal terms are sub-leading in the n 1 limit. All of the results above provide strong evidence that the state prepared by this Floquet driving only contains bipartite entanglement. We can think of the entanglement pattern as being described by many EPR pairs accumulating at the two peaks, i.e. one member of the pair is at one peak and the other member of the pair is at the other peak. In each Floquet cycle, there are c 3 log 2 1 η pairs created. This suggests a quasi-particle picture which is developed in the next section and will help us understand the phenomena outlined by the calculations.

The quasi-particle picture
In this section, we provide a quasi-particle picture to understand the formation of the peaks and the entanglement pattern similar to the discussions in Calabrese and Cardy [52,53]. It is not surprising that such a quasi-particle picture exists since our analysis above should apply to any (1 + 1)D CFT, including the one realized by (1 + 1)D massless free fermion. What is interesting is that the predictions from the quasi-particle picture agree quantitatively with the CFT calculations. From the quasi-particle picture, in each Floquet driving cycle, when we suddenly change the Hamiltonian, we expect that there will be quasi-particle excitations emitting from different points. The pairs of particles moving to the left and right from a given point are highly entangled. For example, at the beginning of each cycle, we change H 0 to H 1 , which creates EPR pairs in the system. Then they move together with all other quasi-particles that have been created in previous cycles with velocity v(x) = 2 sin 2 (πx/L). The velocity is determined by the sin-square envelope we defined in Eq. (3). In the second part of each driving cycle, the  quasi-particles will be governed by H 0 and the velocity will now change to v(x) = 1. Therefore, we can determine the distance that a quasi-particle travels in one cycle by the following formula: This distance depends on T 0 , T 1 , and the initial position x i of the quasi-particle, thus uniquely determining the final position x f of the quasi-particle. For the plot in (40), we assume the quasi-particle bounces back once from the boundary 6 . In general, the number of bounces is n if (n − 1)L < T 0 < nL, where n is a positive integer. One can find that T 1 does not come in, because under H 1 the quasi-particles will never reach the boundary due to the vanishing of velocity at the boundary 7 . With this quasi-particle picture, we can determine the positions of the two energy peaks as observed in Fig. 5. Recall that in the long time limit n 1, after each driving cycle, the positions of the two energy peaks stay the same. There are only two possibilities as follows (without loss of generality let us focus on the chiral energy peak and track its position x C ), 1. If the number of bounces n is odd, then the chiral energy peak will become an anti-chiral energy peak due to the bounces at the boundary. To keep the positions of the chiral/antichiral energy peaks the same, we have to do the switching: 2. If the number of bounces n is even, then the chiral energy peak is still chiral after each driving cycle. Then one has The tracking of the anti-chiral energy peak can be analyzed in the same way. Noting that the two energy peaks are symmetric about x = L/2, i.e., x A = L − x C , we can determine the positions of the two energy peaks by the "quantization" condition 6 The conformal boundary condition ensures the magnitude of the velocity remains the same after the bouncing. 7 More exactly, this statement is only true for a vanishing function at least faster than linear, because t ∼ 0 dx x a diverges when a 1. Here for the SSD, we have v(x) = 2 sin 2 (πx/L)2 ∼ x 2 for x → 0.
Evaluating the above equation explicitly, one can find that x C and x A are determined by the following equation: with x C = x * and x A = L − x * . One can check that the solution x * in Eq. (42) matches the repulsive fixed point γ 2 = e 2πix * /L of the Möbius transformation in Eq. (15) exactly. This quasi-particle picture turns out to be useful. On the one hand, it gives us a semiclassical explanation for the formation of the peaks. In the heating phase, the system keeps absorbing energy by creating many EPR pairs. Due to the fixed point solution in the semiclassical motion, those EPR pairs will accumulate at the fixed point, with chiral part staying at x C = x * and anti-chiral part staying at x A = L − x * . If we keep track of what happens within one cycle, we will find that the particles at the two peaks will switch their position after one cycle if the number of bounces n is odd, but will stay the same if n is even. In the non-heating phase, the system does not absorb much energy and there is no fixed point in the equation of motion. Therefore we only observe energy oscillation.
On the other hand, it also provides us insight into the growth of entanglement entropy. The EPR pairs generated, not only carry energy but also share entanglement. Therefore, as the system absorb energy, the entanglement also grows. Based on this semi-classical picture, it is not hard to conjecture that all the entanglement is shared by the two peaks. Since the energy and entanglement are carried by the same objects, it is also natural to expect some relationship between them. In the following section, we will derive such a relation.

Energy-entanglement relation
As we said before, since we have this interpretation that the energy and entanglement are all carried by those EPR pairs at the two peaks, it will be natural to ask whether there is any relation between them. The result for entanglement entropy always contains a divergent non-universal piece due to the absence of a UV cutoff in a field theory, which is absent on the energy side. Therefore, in the following, we only compare their universal time dependence and dispense with the non-universal part.
First, let us look at the results for the heating phase. By comparing Eq. (34) with Eq. (36), we can find the following equation, which relates the total energy growth to the entanglement growth for the chiral or anti-chiral peaks. We use proportion instead of equality because we only keep the universal information and drop all the other non-universal details. Then, let us see whether this relation also holds in the non-heating phase and the critical case. Because we already know that the entanglement comes from the two peaks, it is sufficient and technically easier if we just use the result for the left half and right half entanglement, as has been computed in [34]. Since they used a different notation, we reproduce their results in Appendix. B using our notations so that readers can compare with the total energy more conveniently. For the non-heating phase, the entanglement entropy keeps oscillating as #η n around a non-zero average value, which matches this result. For the critical case, the entanglement entropy grows logarithmically in time S ≈ c 3 log n, which also matches this result. Several remarks follow below. This equation only contains the central charge and thus is true for any CFT. Generically, we would not expect such a universal relation. It only holds here because the states prepared are related by a conformal transformation that makes the universal relation possible [54]. It is also, to some extent, reminiscent of the Cardy formula in an equilibrium CFT, which says energy is proportional to the square of the entropy. However, what we found here is that the energy is the exponential of the entropy, thus much larger than the entropy, which suggests the state we prepare is far from the equilibrium state.

Effect of randomness in non-heating phase
In a real experiment, local perturbations and imperfections of the pulse sequences are inevitable. In the section, we will discuss the effect of having some small randomness on the driving period T 0 and T 1 , i.e. in each cycle T 0 and T 1 are independently drawn from a certain distribution, the final results are obtained after doing a "disorder" average. Here are some detailed explanations of the protocol, 1. The driving time for each cycle is given by where δT 0 and δT 1 denote the deviation from the (constant) mean values T 0 and T 1 . For simplicity, we consider the case that δT 0 and δT 1 are uniformly distributed in the following domain: where α characterizes the magnitude of randomness. L is the total length of the system, which is the fundamental time scale of the system.
2. Given a sequence of randomized driving time, the operator O on the z-plane under n-cycle imaginary time evolution is given by the familiar formula The derivative terms are calculated using the chain rule, where (c.f. Eqs. (9) and (10)) The difference comparing to the previous sections is that the parameters a, b, c, d here are determined by the randomized driving cycles, in particular, the terms in the chain rule formula are independent. We comment that in principle for the real-time evolution, we need to analytically continue for each cycle in order to keep track of the trajectories of z n and z n to determine the branch cut crossing. The disorder average is done after the analytic continuation. Here we only consider the stress tensor and the entanglement entropy of the left (right) half system. Both quantities are free of the branch cut issue and we can safely perform the analytic continuation.

Energy and entanglement growth
In this section, we present numerical results for the energy and entanglement growth affected by the randomness. We will focus on the non-heating phase and briefly comment on the heating phase.
The first result is presented in Fig. 10, where we choose parameters T 0 , T 1 and α 1 such that each individual sample (T 0 , T 1 ) belongs to the non-heating phase, namely corresponds to an elliptic Möbius transformation. With the randomized protocol, we find that the total energy of the system grows with time (even for α 1). Asymptotically, it grows exponentially with n as numerically verified in Fig. 10 (a). At the same time, the entanglement entropy grows linearly with n, as shown in Fig. 10 (b). That is to say, the non-heating phase will disappear immediately for arbitrarily weak randomness and we are only left with the heating phase. The total energy and entanglement entropy grow in the same way as what we find for the heating phase without any randomness. As a side note, if we implement the random driving set-up Each κ is extracted from the CFT calculation of S A (n) by averaging over N sample = 1000.
for a Mathieu oscillator, small randomness also leads the energy to grow exponentially (see Appendix. C for details), which suggests that the phenomenon here might be a generic feature of the SL(2, R) algebra and not special to our CFT setting. Then, it is desirable to ask how the heating rate is related to the magnitude of randomness α. We will study this problem based on the entanglement entropy, as follows. With randomness, there are four dimensionless parameters in the calculation of entanglement, l A /L, T 0 /L, T 1 /L, and α. Recalling that S(n) also depends on the total length L of the system through its initial value, we consider the quantity S A (n) − S A (0) only depends on the dimensionless ratios we introduced above. 8 Let us introduce the heating rate κ and write the entropy growth as For simplicity, we consider the choice of T 0 and T 1 in Eq. (44) with T 0 = T 1 = T . Then for As shown in Fig. 11(a), we study κ as a function of T /L with different α. There are several interesting features: (i) Fixing α, as T increases from T /L = 0, κ will increase accordingly. In particular, κ grows the fastest near the phase transition T * /L 0.416 (note that the phase transition is defined for the case with no randomness). This indicates that the system is more sensitive to the randomness near the phase transition. (ii) Fixing T < T * , one can find that κ will increase with the randomness α. That is, with larger fluctuations in the driving periods, the system will be heated up more easily. (iii) For T > T * , κ collapse to the same curve, indicating that the heating phase (defined before adding noise) is robust under the effect of fluctuations in the driving periods.
With the analysis above, now we are interested in how κ depends on the magnitude of randomness α for a fixed T with T < T * , which corresponds to the non-heating phase (before adding randomness). As shown in Fig. 11 (b), it is found that κ depends on α in the following way: where the coefficient κ/α 2 depends on T , as can be seen in Fig. 11 (b). With this observation, one can alternatively plot κ/α 2 as a function of T . It is found that κ/α 2 for T < T * only depends on T , with the concrete value 0.1 ∼ 1.

Energy density distribution
After determining how the randomness in the driving periods T 0 and T 1 affects the phase diagram, let us look at how it changes the energy density distribution. If we start from the non-heating phase regime and add randomness, we find that the energy density peaks at the two ends of the system, as shown in Fig. 12. This phenomenon can be understood from the semi-classical quasi-particle picture. Without randomness, the quasi- The energy density distribution is the same as that without any randomness, except that the peaks are smeared out a little.
particles are created and moved in the system coherently. After adding randomness, those motions become irregular. However, since the group velocity of the quasi-particles is smaller near the boundary, accordingly we have a higher probability to see more quasi-particles there. However, if we start from the heating phase regime, the total energy keeps growing exponentially with time and the energy peaks will not disappear for moderate randomness, as shown in Fig. 13. We can first drive the system with a fixed period and let the energy peaks form. Then we turn on sufficiently weak randomness. Now the energy peaks will not be moved back perfectly but with a small discrepancy. As a result, the energy peaks will be smeared a little bit but still there. Therefore, just as the time crystal with MBL [20][21][22], all the features for the heating phase we find here are also robust even though we slightly perturb it away from the fine-tuning (randomness-free) point.
Before we close this section, we want to emphasize that in our CFT calculation of the energy and entanglement evolution in the presence of randomness (see Fig. 10), the average is performed numerically. It is desirable to derive an analytic result, for example, for the heating rate κ in Fig. 11. We leave this problem for future study.

Summary
In this paper, we presented a detailed study of the Floquet CFT introduced in [34]. The phase diagram obtained in that paper can be understood by mapping the problem to a Floquet harmonic oscillator. The reason for such a mapping arises from the fact that these two problems share a sl(2, R) algebra and the classification of dynamics becomes the classification of the linear combination of SL(2, R) generators. Because of that, although we use the entanglement entropy and total energy to explicitly determine the phase diagram, the calculation of which depends on the choice of the initial state, the result is actually a property of the driving Hamiltonian and doesn't depend on the initial state choice.
In the non-heating phase, the energy profile and total energy keep oscillating. In the heating phase, although the energy increases exponentially fast, the system is heated in an extremely non-uniform way, i.e. only two points absorb the heat. What is more, all the entanglement entropy is also shared by these two peaks. These peaks are determined by the fixed points of the relevant Möbius transformation that is defined by the dynamics. On the phase boundary between the heating and non-heating phases, we still observe two peaks but the total energy only increases quadratically with time.
Although the question of whether the system absorbs energy relies on the detailed calculation, the energy density and entanglement structure can be understood by a quasi-particle picture. The questions of how energy distributes can be mapped to solving a pure classical motion. Inspired by this picture, we find a relation between the total energy and the entanglement between the two peaks, E(t) ∝ c exp 6 c S(t) , which says the quasi-particle carries much more energy than entanglement. Such a relation, as contrary to the classic Cardy formula, is a clear manifestation of a non-equilibrium state. It will be interesting to understand whether this relation is special to this set-up that only involves SL(2, R) or is true for more general cases.
To make some connection to the real experiment, we examine the robustness of all these features against random driving. Even if we add tiny randomness to the driving period, the non-heating phase completely disappears and we only have the heating phase, where the total energy grows exponentially with time. After we know whether the total energy grows, the energy density can be analyzed perturbatively. If we start from a (T 0 , T 1 ) that is deep inside the non-heating phase and turn on the randomness, the energy density will peak near the boundary. This is because the quasi-particles move incoherently and have smaller velocity near the boundary. On the other hand if we start from a (T 0 , T 1 ) that is deep inside the heating phase and turn on moderate randomness, we expect the energy peaks will remain although they are smeared out a little.
Most of the phenomena above, in particular including the existence of heating and nonheating phase and the features about the energy profile, not only hold for this special set-up but should also occur for any generic Floquet driving that only uses Virasoro generators as the Hamiltonian. The reason is that this type of Floquet driving can always be thought of as a conformal mapping on the complex plane. The energy density calculation to a large extent can be reduced to the problem of finding fixed points of the conformal mapping. If none of the fixed points is on the unit circle, the system must be at the non-heating phase and energy density just oscillates. Once it has a repulsive fixed point on the unit circle, we will see two energy peaks, one of which is purely chiral and the other is anti-chiral. The system is generically heated up. If there is also an attractive fixed point on the unit circle, the energy density will decrease to make the energy peaks sharper and sharper. For this case, one has to do a more detailed calculation to determine whether the system is heating or not. Consequently, the problem of classifying dynamics is equivalent to the problem of classifying conformal mappings. This Floquet CFT using sine-square deformed Hamiltonian is the first and simplest example that explicitly realizes this. It will be interesting to generalize this special set-up to more general protocols and give a more thorough discussion on the connection between dynamics and geometry.
Another interesting problem is to consider the Floquet CFT from a thermal initial state at finite temperature β −1 , which is closely related with experiments. Since there are now three length scales, i.e., the total length L of the system, the driving periods T , and the finite temperature β, then the time evolution of entanglement and energy density may exhibit more rich features in particular in the early time of driving. In the long time limit, we expect there are still two phases, i.e., the heating and non-heating phases. One intuition is based on the quasi-particle picture as presented in Sec.3.3. One can find that the existence of fixed point or not in the solution of equation of motion, which determines the system is in heating or non-heating phases, is independent of the introduction of finite temperature β −1 . We expect these two phases will persist even if the system is prepared at a thermal initial state. We leave the detailed study in a future work.
The heating phase discussed here realizes a highly nonequlibrium state where entangled EPR pairs are continuously produced and localized at specific locations. Given the utility of entanglement as a resource for quantum information processing, experimental realization of the protocols discussed here may be desirable. Indeed given the high tunability of ultracold atomic systems in optical lattices [55], and the ability to measure both energy and entanglement entropy [56], an important future direction will be to find routes to implement these protocols in the lab.

A Phase diagram of the Floquet CFT
The Floquet CFT defined in Sec. 2 was known to have two different phases, which was first shown in [34]. Here, we reproduce the phase diagram in Fig. 14 for the sake of being self-content. The high frequency regime is a non-heating phase, where the entanglement entropy and energy oscillate. The low frequency regime is a heating phase, where the entanglement entropy linearly grows and the energy exponentially grows. On the phase boundary, the entanglement entropy grows logarithmically and the energy grows quadratically. Fig. 14 only shows one domain, and the phase diagram repeats itself when we increase T 0 /L with a period 1.

B Entanglement growth for a subsystem and branch cut crossing
In this section, we will present the details of the entanglement entropy growth calculation for a subsystem in the heating phase. As the early-time regime contains non-universal information, our analytical analysis will focus on the late-time regime (i.e. n 1).
The one point correlation function T m (z n , z n ) in a boundary CFT can be mapped to a two point function through the mirror trick on the whole plane, i.e.
T m (z n , z n ) ∝ Note the derivative term in Eq. (53) decreases exponentially as a function of n in the long time limit 9 ∂z n ∂z hm ∂z n ∂z which is related to the linear growth of entanglement entropy. While the behavior of T m (z n , z n ) depends on an intersting branch cut structure that will lead to the spatial feature (the kink) of the entanglement entropy plotted in Fig. 7.
More explicitly, the branch cut arises from the factor √ z n − √ z n . The subtlety is that although both z n and z n flow to the same attractive fixed point γ 1 at long time limit, their

B.3 Two entanglement cuts
In this section, we present the details of computing the entanglement entropy for a subsystem that does not end at the boundary, i.e. with ending points x 1 , x 2 ∈ (0, L). In other words, we need to insert two twist operators Tr ρ m A (t) = C m (x 1 , x 2 , t) = ψ(t)|T m (x 1 )T m (x 2 )|ψ(t) = G|T m (x 1 , t)T m (x 2 , t)|G .
We follow the strategy in Sec. 2 to do the calculation first in the imaginary time and analytic continue to real time in the end. The n-dependent part of the above formula is given as follows, j=1,2 ∂z j,n ∂z j hm ∂z j,n ∂z j hm T m (z 1,n , z 1,n )T m (z 2,n , z 2,n ) .
The mirror trick maps the two-point function T m (z 1,n , z 1,n )T m (z 2,n , z 2,n ) in a boundary CFT to a four-point function on the whole plane without boundary. The important n-dependent part is given by the following formula,

C Random driving Mathieu oscillator
In this section, we discuss the random driving Mathieu oscillator. The classical Newton's equation for a Mathieu oscillator is h controls the amplitude of the driving force and T is driving period. The intrinsic period of the harmonic oscillator is 2π. For a weak driving force h 1, the first smallest unstable driving period is T unstable = π. The system is stable(non-heating) for any T < T unstable . These are shown in Fig. 16 (a).
For a random driving Mathieu oscillator, we let the driving period T uniformly distribute in an interval T = T + δT, δT = [−α, α], where α controls the strength of the randomness. In each cycle, we randomly choose a T from the distribution and evolve the system accordingly. The final results will be averaged over "disorder". If we choose T < T unstable and α 1, we find that the averaged energy grows exponentially with time, as is demonstrated in Fig. 16 (b).