Conformal Floquet dynamics with a continuous drive protocol

We study the properties {of a conformal field theory} (CFT) driven periodically with a continuous protocol characterized by a frequency $\omega_D$. Such a drive, in contrast to its discrete counterparts (such as square pulses or periodic kicks), does not admit exact analytical solution for the evolution operator $U$. In this work, we develop a Floquet perturbation theory which provides an analytic, albeit perturbative, result for $U$ that matches exact numerics in the large drive amplitude limit. We find that the drive yields the well-known heating (hyperbolic) and non-heating (elliptic) phases separated by transition lines (parabolic phase boundary). Using this and starting from a primary state of the CFT, we compute the return probability ($P_n$), equal ($C_n$) and unequal ($G_n$) time two-point primary correlators, energy density($E_n$), and the $m^{\rm th}$ Renyi entropy ($S_n^m$) after $n$ drive cycles. Our results show that below a crossover stroboscopic time scale $n_c$, $P_n$, $E_n$ and $G_n$ exhibits universal power law behavior as the transition is approached either from the heating or the non-heating phase; this crossover scale diverges at the transition. We also study the emergent spatial structure of $C_n$, $G_n$ and $E_n$ for the continuous protocol and find emergence of spatial divergences of $C_n$ and $G_n$ in both the heating and non-heating phases. We express our results for $S_n^m$ and $C_n$ in terms of conformal blocks and provide analytic expressions for these quantities in several limiting cases. Finally we relate our results to those obtained from exact numerics of a driven lattice model.

We study the properties of a conformal field theory (CFT) driven periodically with a continuous protocol characterized by a frequency ωD. Such a drive, in contrast to its discrete counterparts (such as square pulses or periodic kicks), does not admit exact analytical solution for the evolution operator U . In this work, we develop a Floquet perturbation theory which provides an analytic, albeit perturbative, result for U that matches exact numerics in the large drive amplitude limit. We find that the drive yields the well-known heating (hyperbolic) and non-heating (elliptic) phases separated by transition lines (parabolic phase boundary). Using this and starting from a primary state of the CFT, we compute the return probability (Pn), equal (Cn) and unequal (Gn) time two-point primary correlators, energy density(En), and the m th Renyi entropy (S m n ) after n drive cycles. Our results show that below a crossover stroboscopic time scale nc, Pn, En and Gn exhibits universal power law behavior as the transition is approached either from the heating or the nonheating phase; this crossover scale diverges at the transition. We also study the emergent spatial structure of Cn, Gn and En for the continuous protocol and find emergence of spatial divergences of Cn and Gn in both the heating and non-heating phases. We express our results for S m n and Cn in terms of conformal blocks and provide analytic expressions for these quantities in several limiting cases. Finally we relate our results to those obtained from exact numerics of a driven lattice model.

I. INTRODUCTION
The study of driven quantum systems have attracted a lot of attention in recent years [1][2][3][4] . The theoretical interest in this area stemmed from the fact that such systems provide access to a gamut of phenomena that have no analog in their equilibrium counterparts. Some of these phenomena, for periodically driven systems, include dynamical phase transitions 5,6 , dynamical freezing 7-9 , realization of time crystals 10,11 , and the possibility of tuning ergodicity of the driven system with the drive frequency 12 . Moreover, periodically driven systems can lead to novel steady states which has no counterpart in equilibrium quantum systems 9,13 . The study of these phenomena has also received significant impetus from the possibility of realization of closed quantum systems using ultracold atom platforms; indeed, such platforms have recently been used to experimentally probe several nonequilibrium phenomena [14][15][16] .
A large set properties of such periodically driven systems can be inferred from studying its Floquet Hamiltonin H F which is related to the evolution operator U of the system via the relation U (T, 0) = exp[−iH F T / ] 3 . The computation of H F for a generic many-body system provides a significant challenge; indeed, an exact analytic computation of H F is generally possible for integrable models driven by discrete stepwise protocols (such a square pulse or kicks) 3 . This has led to development of several perturbative schemes for computation of H F ; some of these include Magnus expansion 3 , adiabaticimpulse approximation 17 , the Hamilton flow method 18 , and the Floquet perturbation theory (FPT) 19,20 . Out of these methods, the Floquet perturbation theory has the advantage of being easily applicable to a wide class of systems as well as being accurate over a large range of drive frequencies 21 .
More recently, several theoretical studies concentrated on the effect of Floquet dynamics on conformal field theories [23][24][25][26][27][28] . Usually driving a CFT is expected to generate an additional scale in the problem which drives the system away from its conformally invariant fixed point. However, it was recently realized that there is a class of models 29,30 where drive protocols need not break the conformal symmetry 23 . It is found that for CFTs with a sine square deformation (SSD) such dynamics can be initiated by a Hamiltonian whose holomorphic part is given by where L n for each integer n denotes a holomorphic generator of the Virasoro algebra, where c 0 denotes the central charge. These holomorphic generators are related to the stress tensor T µν of the CFT by dxT 00 (w) where w = τ + ix, x is the spatial coordinate and τ is the Euclidean time. The anti-holomorphic ones have a similar expression in terms of the complex conjugates.

arXiv:2101.04140v1 [hep-th] 11 Jan 2021
It is well-known that the class of Hamiltonians given by Eq. 1 are valued in su (1,1), giving rise to the evolution operator 23 valued in the group SU (1, 1) with ad − bc = 1. Using the isomorphism of SU (1, 1) and SL(2, R) we shall also have occasion to write the evolution operator in imaginary time in the generic form Its action on the complex plane C is given by the Möbius transformation whereã,b,c,d ∈ R andãd −bc=1.
The exact solution for U (T, 0) for Hamiltonians valued in su(1, 1) and driven by discrete protocols were discussed earlier for periodic kicks 22 and square pulse protocols [23][24][25][26][27][28] . It was found that such driven system could display two distinct phases depending on the drive parameters; these are termed as the heating (hyperbolic) and the non-heating (elliptic) phase and are found to be separated by a transition line (parabolic phase boundary) 22 . The presence of these phases can be shown to be a direct consequence of non-compact nature of the SU(1,1) group. However, such studies have not been extended to continuous drive protocols where exact analytic results are not available 31 ; in particular, the phase diagram of such a driven system for continuous drive protocols has not been studied so far.
It was noted in Ref. 23 that the action of the evolution operator U on a primary operator of the CFT can be understood, in the Heisenberg picture, in terms of a Möbius transformation of it's coordinates leading to where z is defined in Eq. 6 and the bar designates anti-holomorphic variables throughout. Using this relation, and a straightforward mapping from cylindrical or strip geometries to the complex plane following standard prescription 32 , several results were obtained on the energy density, correlation function and entanglement entropy of the driven system [23][24][25][26][27][28] . It was shown that in the heating phase, such a drive leads to emergent spatial structure of the energy density 27 . Moreover, it was found that the time evolution of the entanglement entropy shows linear growth with n (in the large n or longtime limit) in the heating phase. In contrast, it shows an oscillatory behavior in the non-heating phase and a logarithmic growth on the parabolic phase boundary. All of these studies focussed on evolution starting from the CFT vacuum on a strip geometry [23][24][25][26][27][28] ; the dynamics of the system starting from asymptotic states corresponding to primary operators of the theory, which necessitates computation of four-point correlation functions of the primary fields of the driven CFT, has not been studied so far.
In this work, we study dynamics of conformal field theories subjected to a continuous drive protocol. The protocol we use corresponds to the Hamiltonian given by Eq. 1 with f 1 (t) = 1 and where f 0 is the drive amplitude and δf is a constant parameter. In this study we shall focus on the large drive amplitude regime which corresponds to f 0 δf, 1. In this regime for ω D ≥ δf, 1, the FPT is expected to be accurate and we expect this to provide us with an analytic, albeit perturbative, understanding of the properties of the driven system.
The central results that we obtain from such a study are as follows. First, we chart out the phase diagram of the driven system as a function of δf and ω D using exact numerics. Our results show re-entrant heating and nonheating phases separated by a parabolic phase boundaries as a function of δf and T . Second, we provide analytic expressions of the evolution operator U (T, 0) using FPT. The phase diagram obtained from the perturbative result provides a near-exact match with its exact numerical counterpart over a wide range of frequencies. Our perturbative results indicate the existence of a parameter α, given within first order FPT by (where we have put = 1) which indicates proximity of the system to the phase boundary. The system stays in the non-heating phase for α 2 < 1 and heating phase for α 2 > 1; the phase boundary between these two phases is given by α = ±1.
We also use the first order FPT results to explain the re-entrant transitions between the heating and the nonheating phases in the phase diagram. Third, we use the obtained expressions for U within FPT, to obtain analytic expressions for energy density, equal and unequal time two-point correlation functions, entanglement entropy, and return probability of generic sine-square deformed CFT starting from a primary state. Our results expresses these quantities in terms of α and allows us to characterize their behavior near the transition from the heating phase. For example, we find that the return probability of any primary state after n drive cycles shows an universal behavior below a crossover time n c both in the heating and the non-heating phases near the transition line; for n > n c , the probability decays exponentially for the heating phase and remains an oscillatory function in the non-heating phase. We provide analytic estimate of n c as a function of α. Fourth, we discuss the nature of the emergent spatial structure of the energy density and the correlation function of primary operators for the drive protocol. We analytically show that the energy density of any primary state in the heating phase displays peaks which shifts from L/4 and 3L/4 to L/2 as one moves from deep inside the heating phase to the transition line. For unequal-time correlation function starting from the CFT vacuum, we find a line of such peaks both in the heating and the non-heating phases; the position of these peaks can be analytically found within first order FPT. Finally, we provide expressions of the equal-time correlation function C n and half-chain m th Renyi entropy S m n of the driven CFT after n drive cycles starting from an initial primary state in terms of conformal blocks V p . In general, these blocks do not have analytical expression for arbitrary CFTs; here, we provide their analytical forms in several asymptotic limits. We discuss the applicability of these limits to the driven CFT and discuss the properties of C n and S m n in these asymptotic regimes. Finally, we relate some of our results to those obtained by exact numerical study of the SSD model on a 1D lattice 29,30 .
The plan of the rest of the paper is as follows. In Sec. II, we derive expressions of U and provide the phase diagram for our drive protocol. This is followed by Sec. III where we compute energy density, correlation functions and return probabilities starting from a primary state. Next, in Sec. IV, we relate our results to those obtained from numerical study of driven SSD Hamiltonian on a 1D lattice. Finally, we discuss our main results and conclude in Sec. V. Some details of the perturbative FPT calculations and representation independent derivation of the Mobius transformation are presented in the Appendices.

II. PHASE DIAGRAM
To find U (T, 0) corresponding to the Hamiltonian given by Eq. 1 with f (t) given by Eq. 8, we first note that the holomorphic generators (Eq. 2) for n = −1, 0, 1 form an su(1, 1) subalgebra (10) A representation of this algebra is furnished in terms of the Pauli matrices as where σ α for α = x, y, z are standard Pauli matrices. In this representation (Eq. 11), the holomorphic part of the Hamiltonian H(t) becomes which corresponds to a Zeeman Hamiltonian of a single spin-half particle with a time dependent magnetic field alongẑ and an imaginary constant magnetic field along y. The latter feature is a consequence of the SU(1,1) group structure and is crucial in realization of the heating and the non-heating phases 22 . The expression for U (T, 0) corresponding to Eq. 12 can be found exactly for periodic kicks 22 and square-pulse 23 protocols. In contrast, for the continuous protocol given by Eq. 8, an exact analytic expression for U (T, 0) does not exist. However, the numerical result for U can be obtained in a straightforward manner. Such a numerics is carried out by dividing the time period T into N Trotter steps with width δt i ∼ T /N . The maximal allowed width of each of these steps depend on system energy scale and the drive frequency; they are chosen so that H does not vary appreciably within any step. This allows one to write This procedure leads to U (T, 0) = exp[−iH F T / ]; we find numerically that the Floquet Hamiltonian is We note that σ x does not appear in H F . The condition for the different phases are thus obtained by 22 |TrU (T, 0)| > 2 (heating phase), < 2 (non-heating phase), and = 2 (phase boundary). The exact numerical phase diagram obtained from this procedure is shown in the left panel of Fig. 1. We find several re-entrant transitions between the heating and the non-heating phases leading to multiple lobes whose boundary correspond to the transition lines. We note that the transition between these phases can be induced by tuning either δf or ω D .
In particular, we note that large ω D δf, 1, the phase transition between the heating and non-heating phase occurs at δf = 1.
The absence of σ x in H F can be shown to be the consequence of an emergent dynamic symmetry of the evolution operator at stroboscopic times which follows from the periodicity of H(t). To see this we note that since H(t) = H(T − t), one has U j = U N −j+1 for all j in the Trotter product. Further, we note for any of these Trotter steps, one has σ x U j σ x = U −1 j . Thus one can write where we have used the relation U −1 This forbids presence of terms ∼ σ x in H F . It is to be noted that this symmetry is absent for t = nT where n is an integer.
Next, we develop an perturbative analytic expression of U in the limit when f 0 is the largest scale in the problem. Here the diagonal term in H(t) is treated exactly and the effect of the off-diagonal term is taken into account within standard time-dependent perturbation theory 20 . In this scheme, the first term for U and H F is given by where s = arccos[cos(δf T π/L)]. We note that H F retain the periodic structure of U 0 .
To obtain the first order correction to U (T, 0), we use the standard result 20 where I denotes the 2 × 2 identity matrix and U 1 denotes the first-order correction to U in the interaction picture given by where H 1 = − π L iσ y is the perturbative term of the Hamiltonian. A simple calculation detailed in App. A yields where α is given by Eq. 9. We note that U 1 (T, 0) is not unitary; this is a well-known issue with perturbation theory for U . In what follows we unitarize U 1 (T, 0) as follows. We note that this can be done by first writing (I − U 1 (T, 0)) exp[−iH F T ], where the matrix H F = (i/T )U 1 . The evolution operator is then given by U 0 (T, 0) exp[−iH F T ] which is unitary. However, in the present case, we find that evolution operator obtained by this method retains terms ∼ σ x . This is clearly inconsistent with the dynamic symmetry discussed in Eq. 14. Thus we chose to unitarize U 1 (T, 0) directly since it does not have any term σ x . These two alternative routes to obtaining H F coincides with each other in the large ω D limit, where s → 0.
The unitarization of U 1 (T, 0) given by Eq. 18 can be achieved by writing Comparing Eqs. 13 and 19, we find that the first order FPT result provides approximate analytic expressions for p(T ) and q(T ) given by In what follows we shall use Eq. 20 to compare between exact numeric and approximate analytic results.
The evolution operator obtained in Eq. 19 indicates several interesting features. First, we find that n y is imaginary; this is a direct consequence of the SU(1,1) group structure. Second, we note that Tr[U 1 (T, 0)] = 2 cos(θ); thus the condition |Tr[U 1 (T, 0)]| > (<)2 for realization of heating (non-heating) phases translates to θ being imaginary(real). This happens when α 2 > (<)1 which leads to our identification of α as the parameter whose value determines the phase the system. Third, the parabolic phase boundary between these two phases is given by |Tr[U 1 (T, 0)]| = 2 which leads to θ = 0. This is realized for α = ±1 on the boundary between the heating and non-heating phases. From the expression of α in Eq. 9, we find that for ω D f 0 , the condition α = 1 is satisfied for δf = 1. This can be verified directly from the exact phase diagram in the left panel of Fig. 1. Fourth, the re-entrance of the heating and non-heating phases as a function of δf shown in the right panel of Fig. 1 displaying the phase diagram obtained from first order FPT, can be easily explained by periodic nature of | cos(θ)| since the phases must repeat for θ → θ + π. We note the two phase boundaries in each of the lobes of this phase diagram correspond to α = 1 (lower branch of the lobes) and α = −1 (upper branch of the lobes). The center of these lobes correspond to α → ∞ which occurs along the lines δf /ω D = n 0 /2 for n 0 ∈ Z. This can be directly seen from Eq. 9 where the term in the sum corresponding to n = n 0 diverges in this limit. Finally, a comparison between the left and the right panels of Fig.  1 shows that the two phase diagrams match qualitatively for ω D /(π/L) > 1. This allows us to justify the use of the analytical method for a qualitative understanding of the dynamics within first order FPT. A computation of the second-order results of FPT is carried out in App. A; we find that it retains all features of the first order theory and provides a near-identical phase diagram. In App. B, we chart out a representation independent derivation of the first order perturbation results in the high frequency limit which matches the results obtained here using the SU(1,1) representation of the Virasoro generators.
Before ending this section , we obtain the elements of U n = U (nT, 0) (where n ∈ Z is an integer). Using Eqs. 4 and 13, we find the elements of U j for any integer j obtained using exact numerics in the non-heating phase (where p > q) to be where d j = a * j and c j = b * j . For the heating phase, for which p < q, one obtains For the transition line where p = q, a careful evaluation of the limit shows The corresponding analytic, first order FPT, expressions of elements of U (jT, 0) in terms of α and s can be directly read off from these equations by using the correspondence expressed in Eq. 20. For example for the transition line, the first order FPT result is a" j = 1 − ijs and b" j = −ijs. Also, the corresponding elementsã j ,b j ,c j and d j of the evolution operator in Euclidean time can be obtained for each of these phase from Eqs. 21, 22, and 23 via standard analytic continuation T → −iτ . We shall use these expressions in the next section computation of return probability, energy density, correlation functions and the entanglement entropy.

III. RESULTS FOR DRIVEN CFT
It was pointed out in Ref. 23 that the operation of the evolution operator U on any primary operator of the CFT can be understood in terms of a Möbius transformation of its coordinates and is therefore given by Eq. 7. This is further discussed in App. B. In this section, we shall use this result to study the properties of the driven CFT with a cylindrical geometry which corresponds to periodic boundary condition of a driven 1D chain of length L. In our setup, the coordinate of the cylinder is given by w = τ + ix, where x is the spatial coordinate and τ is the Euclidean time; we use w = (L/(2π)) ln z for mapping such a cylinder to the complex plane. Thus for any primary operator O(w,w) on the cylinder, we can write 23,32 where z n is the transformed coordinates given by Eq. 6 andã n ,b n ,c n , andd n are the elements of U n in Euclidean time which can be obtained from analytic continuation of T → −iτ in Eqs. 21, 22 and 23. In all computations, we shall work in Euclidean time and carry out the analytic continuation to real time using τ = iT at the end of the calculation. Apart from the primary operators we shall also use transformation properties of the stress tensor. The holomorphic part of the stress tensor denoted by T (w) transforms, under a general coordinate transformation, as Here we note that S zn,z = 0 since S vanishes for any Möbius transformation and S z,w = c 0 /(24z 2 ) for transformation from a cylinder to the complex plane. The transformation for the antiholomorphic part of the stress tensor can be obtained in the similar manner by replacing (w , w) → (w ,w).
In what follows, we shall use these general results to obtain the stroboscopic time (n) dependence of the return probability, energy density, equal and unequal time correlation functions, and the entanglement entropy of the driven CFT. The initial state for this purpose is chosen to be an asymptotic in-state of the CFT denoted by where |0 denotes the CFT vacuum, φ is a primary field with dimension (h,h) and we note w,w → −∞ corresponds to τ → −∞ which in turn implies z,z → 0 on the complex plane. This defines for us a primary state. The corresponding out-states are given by which corresponds to τ → ∞ and hence z,z → ∞ on the complex plane. The reason for choosing these primary states are as follows. First, we note that since L 0 and L ±1 annihilates the CFT vacuum, the vacuum state does not evolve. This renders all equal time correlation functions of primary operators, entanglement entropy and energy density to be fixed at their equilibrium values under action of U ; only the unequal-time correlation function shows nontrivial dynamics. This feature of driven CFT usually compels one to work with the strip geometry 23,26 or use different Möbius transformation 27 where one can obtain non-trivial time evolution starting from the ground state.
Here we use an alternative approach by using the cylindrical geometry but starting from the primary states of the CFT (|h,h ) as initial states; these are the simplest states of the model which have non-trivial dynamics under action of U .

A. Return probability and energy density
The return probability amplitude A n after n cycles of the drive is given by where the denominator corresponds to the normalization of the asymptotic states. To this end, we first look at the contribution of the holomorphic part of the return probability amplitude given by To compute A hol n , we take w i = τ i + ix i for i = 1, 2 and consider the limits τ 1 → −∞ and τ 2 → ∞ at the end of the calculation. Using the fact that 0|U † = 0|, we can rewrite as N = 0|U n † φ(w 2 ,w 2 )U n φ(w 1 ,w 1 )|0 . We then use Eq. 24 to obtain where z n2 = (ã n z 2 +b n )/(c n z 2 +d n ) is the transformed coordinate,ã n ,b n ,c n andd n are obtained from Eqs. 21, 22 ,23 after analytic continuation T = −iτ and we have used ∂w/∂z = L/(2πz) and ∂z n2 /∂z = (c n z 2 +d n ) −2 . After a careful evaluation of the limits, we finally obtain A hol n =ã −2h n . A similar computation shows the contribution from the anti-holomorphic part to be A ant−hol n = a −2h n . Thus one obtains the return probability P n , after analytic continuation to real time, to be The exact numerical value of P n can thus be obtained by using Eqs. 21, 22, 23. Here and for all quantities in the rest of this section, we shall provide expressions for the prediction of FPT; the numerical result in terms of p and q can be directly read from these expressions by using Eq. 20(i.e. with the mapping α → p/q and s → pT ). This procedure yields , non − heating The behavior of ln P n is shown in the left panel of Fig.  2 as a function of drive frequency after n = 50 cycles of the drive and for δf = 0.5. We find that ln P n shows dips at large n in the heating phases; in contrast, it shows oscillatory behavior near the transition in the non-heating phase as can be seen from the inset of the left panel of Fig. 2. We also note that the return probability shows an exponential decay for large n in the heating phase and an oscillatory behavior in the non-heating phase according to standard expectation. On the transition line, it shows a power law decay ∼ n −2(h+h) for large n. These qualitatively different behaviors of P n can be clearly seen in the right panel of Fig. 2 where ln P n is plotted as a function of n for three values of drive frequencies corresponding to three different phases. Here we also note that near the phase transition line, in either phase, the behavior of P n is identical to that on the critical line below a crossover timescale n ≤ n c . This is most easily seen by expanding the cosh(nθ) term in powers of nθ in Eq. 32. This procedure yields, in the heating phase, where the ellipsis denote terms higher order in n. We note that n c here is estimated as n for which the contribution of the term ∼ n 4 becomes equal to the leading term ∼ n 2 . We find that n c diverges at the transition as 1/ √ α 2 − 1 as one approach the transition from the heating phase; it can also be made large by tuning either the drive frequency or amplitude since s ∼ δf T . For n n c , the behavior P n is identical to that on the transition line. A similar result can be obtained if the transition line is approached from the non-heating phase. Thus we find that characteristic behavior the return probability on the transition line can be observed upon approaching the line below a crossover timescale n c which diverges at the transition; P n has a universal dependence on s for n n c . This behavior is shown in Fig. 3 in the high drive frequency regime ( ω D )/(π/L) = 100 where δf = 1 gives the position of the transition line. We clearly see that ln P n remains indistinguishable for n ≤ n c 50. The divergence of n c at the transition line is shown in the right panel of Fig. 3.
Next, we study the behavior of the energy density of the driven system given by To evaluate this, we first consider the holomorphic contribution to the energy density. To this end, we first move from the cylinder to the complex plane so that one can write where the constant term c 1 = (2π/L) 2 c 0 /24 comes from the Schwarzian given by Eq. 25 and c 0 denotes the central charge. In what follows we shall ignore this term since it does not change upon driving the system. Using this and Eqs. 24 we find Eq. 36 reduces the computation of the energy density to that of computation of φT φ correlator in the CFT vacuum state. This can be done in a straightforward manner using standard identity 33 and yields We substitute Eq. 37 in Eq. 36 and find that in the limit z 2 → ∞ only the second term in the right side of Eq. 37 contribute. Taking the z 1 → 0 limit, we finally obtain the holomorphic contribution to the energy density to be The antiholomorphic contribution can be simply read off from this to be E hol n (x) = E ant−hol n (−x). Thus the total energy density which is the sum of the holomorphic and anti-holomorphic parts are given by We note that the expressions of the energy density of the asymptotic states is identical to that obtained for the vacuum state in the strip geometry in Ref. 26. The only difference appears in the prefactor; the asymptotic states have a prefactor h while the vacuum in the strip geometry has a prefactor of c 0 /32. Thus we expect similar emergent spatial structure for E n (x) as discussed in Ref.
26. Below, we analyze this phenomenon for the continuous drive protocol. To this end, we analytically continue to real time, substitute Eq. 21 in Eq. 40 and, using Eq. 20, obtain for the heating phase We note from Eq. 41 that for a generic x, E n (x) decays exponentially with n with the large n limit: E n (x) ∼ e −4nθ . This behavior can be in Fig. 4 where E n (x) is plotted for several n as a function of x. The plot shows that E n (x), at large n, decays at all positions except at two places for which the leading terms of Q 2 1 (x) and Q 2 2 (x) in the large n limit cancels each other. A straightforward calculation shows that within first order FPT. Thus the position of the peaks of E n (x) in the large n limit moves from L/2 to L/4 and 3L/4 as one moves from the phase boundary to deep inside the heating phase. This behavior is supported for exact numerics as can be seen in the right panel of Fig. 4 where we plot the peak positions (±δx/L) as a function of δf for ω D = 100(π/L).
In contrast for the non-heating phase we find using Eq. 22, E n (x), in real time, is given by Eq. 41 with Thus E n (x) is an oscillatory function of n as shown in the left panel of Fig. 5 where E n (x) is plotted as a function of x of several n. Finally on the transition line, using Eq. 23, we find   where we have analytically continued to real time. Thus we find the peak of E n (x) in the large n limit occurs at x = L/2 for which Q" 1 (x) = 1 and Q" 2 (x) = 0. In contrast, for x = 0, L, we have Q" 1 (x) = 1 + 4n 2 s 2 and Q" 2 (x) = 0, so that E n (x) ∼ 1/(1 + 4n 2 s 2 ) 2 for all n. This behavior is consistent with the fact that n c diverges at the transition. A plot of E n (x) as a function of x for several n, shown in the right panel of Fig. 5, confirms this behavior.
We also note that E n (x) takes particulary simple forms at the center and end of the chains where it is given by From Eq. 45, we once again find the existence of a crossover timescale n c = n c (1 + α)/(1 + 4α) below which E n (0) decays as 1/(1 + µn 2 ) 2 where µ = 2s 2 α(1 + α). This behavior is verified in the left panel of Fig. 6 which shows the behavior of E n (0) (with E n (0) = E n (L)) as a function of n; we find that the curves for the heating and non-heating phase become identical to the transition line for n = n c . For n n c , the decay becomes exponential in the heating phase. In contrast, for E n (L/2), plotted in the right panel of Fig. 6, there is a clear distinction between behavior of E n in the heating phase and on the transition line. This can be attributed to the fact that x = L/2 coincides with the position of the peak on the transition line, while E n (L/2) decays in the heating phase.

B. Correlation functions
In this subsection, we present results for both equaltime and unequal-time correlation functions of primary operators and the stress tensor.

Equal time correlation function
We begin with the analysis of the equal-time correlation function starting from the asymptotic state |h,h given by where w j (w j ) = +(−)ix j for j = 1, 2.
We first compute the holomorphic part of C n (x 1 , x 2 ) which is given by where z jn = (ã n z j +b n )/(c n z j +d n ). Here we choose the dimension of the primary fields φ to be same that of the state |h . The anti-holomorphic part an be written similarly in terms ofz andz jn . To compute the four-point correlators of the primary fields, we use the standard result which expresses these in terms of the cross ratio of the complex coordinates z i . To this end, we define z ij = z i − z j with z k = z kn for k = 1, 2. Using this notation, we define the cross ratio We note that for z 3 → ∞ and z 4 → 0, we have = (ã n e iπx2/L +b n e −iπx2/L )(c n e iπx1/L +d n e −iπx1/L ) (ã n e iπx1/L +b n e −iπx1/L )(c n e iπx2/L +d n e −iπx2/L ) A similar result forȳ n can be obtained for the antiholomorphic part by replacing To compute the four-point correlator, we first define the quantity To evaluate this, we first note that in the non-heating phase and on the transition line, |1 − z 2n /z 1n |, |1 − z 2n /z 1n | → 0 for large n in Euclidean time. To take advantage of this limit, we make the conformal transformation (where z 1 ≡ z 1n and z 2 ≡ z 2n ) so that in the new coordinate z 2n = 0. A similar transformation is done forz i . Then using the standard transformation rule of operators φ(z,z) → (∂z i /∂z i ) −h (∂z i /∂z i ) −h φ(z i ), one can write, after some straightforward algebra 34 where y n = z 2n /z 1n andȳ n =z 2n /z 1n . The advantage of using this form for F which admits conformal block decomposition is that one can write down a perturbative expansion of the blocks around |1 − y n |, |1 −ȳ n | = 0. This allows us to obtain an analytic, albeit perturbative expression of C n (x 1 , x 2 ) for arbitrary h,h. For the driven problem, |1 − y n |, |1 −ȳ n | 0 (in Euclidean time) in the non-heating phase and on the transition line for all x 1 , x 2 in the large n limit. Also, for all phases, this limit holds for all n only for |x 1 − x 2 |/L 1. The perturbative results that we chart out next is expected to be accurate in these limits. For the present case, one obtains 35 where p denotes sum over primaries and V p denotes the p-th conformal block with dimensions h p ,h p and the coefficients C hhp that depend on the details of the CFT. We note here that the identity block corresponds to h p = 0 and for this block C hhI = 1.
Substituting Eq. 53 in Eq. 47 (and its corresponding anti-holomorphic part), one obtains where only the leading term of the k F k x k 1 is retained.
Next, we discuss the contribution of the identity block which is universal since C hhI = 1 = Chh I . For this block, h p =h p = 0. We note that if we retain only the first order term, we find from Eq. 54 that C n (x 1 , x 2 ) becomes independent of n. The first non-trivial contribution of the identity block arises from the x 2 term in the expansion of k F k x k . This universal contribution is given by Thus in this limit, the identity block contribution to the deviation of C n from its equilibrium value provides a measure of the central charge. A plot of J n (x 1 , x 2 ) = (C n (x 1 , x 2 ) univ /C 0 (x 1 , x 2 ) − 1)/(2h 2 /c 0 ) (after analytic continuation to real time) is shown in the top panels of Fig. 7 as a function x 1 /L = x/L for x 2 /L = x 1 /L + 0.03 after n = 100 cycles of the drive. For these plots we have chosen ω D = 100π/L. The top left panel shows the behavior of C n (x 1 , x 2 ) univ in the non-heating phase (δf = 1.5) displaying a broad oscillatory structure. In contrast, in the top right panel, for the heating phase (δf = 0.9) it displays two sharp peaks consistent with the behavior of E n (x). The position of these peaks shift to L/2 on the transition line as can be seen from the inset of top right panel of Fig. 7. The bottom panels show the behavior of C n (x 1 , x 2 ) univ for large n = 100 in the non-heating phase (left panel) and on the transition line (right panel) where the perturbative expansion of F is expected to be accurate. We find clear emergence of spatial pattern in the non-heating phase in contrast to the behavior of E n . We shall discuss this behavior in more details in the context of unequal-time correlation function.
In the large central charge c 0 limit, with the conformal dimensions held fixed, the Virasoro conformal blocks reduce to global conformal blocks which are given in terms of hypergeometric functions, by 36 , (56) Hence the correlator is given by where C 0 (x 1 , x 2 ) is defined in Eq. 54. Note that since the Virasoro blocks reduce to global block in this limit, there is no identity block, i.e., h p = 0.
Finally, we note for CFTs with large central charge c 0 1, when the asymptotic states have large conformal dimension H h, with h/c 0 and H/c 0 both held fixed, a closed form answer is available via the monodromy methods.
Therefore we compute the equal-time correlation function C n (x 1 , x 2 ) = H,H|φ(w 1 ,w 1 )φ(w 2 ,w 2 )|H,H / H,H|H,H . This can be done exactly in the same way as charted out above; the only difference is that one needs to keep track of two operator dimension H and h. A straightforward calculation shows that in this case one has 38 with a 0 = 1 − 24H/c 0 . A similar expression can be obtained forVp by substituting z jn →z jn and h, h p → h,h p . One can then write C n in terms of V p andV p as A motivation for studying such large c 0 CFTs comes from the AdS/CFT correspondence. These are CFTs which are expected to have semiclassical gravity duals. The global block answer for light correlators is reproduced in bulk AdS by the geodesic Witten diagrams 37 . The Virasoro vacuum block is non-trivial in two dimensions unlike its higher dimensional versions as it contains the stress tensor and its descendants. In the large c 0 limit, the dynamics of the Virasoro vacuum block matches with results from semiclassical gravity 38 . In CFT 2 , since one has Virasoro, one can use the geodesic Witten diagrams to interpolate between the global and the semiclassical monodromy answer (when pair of operator conformal dimensions scale with the central charge) by taking into account backreaction due to the heavy geodesics 39 . The time-dependent drive of an inhomogeneous metric will have implications for the physics of black holes in the dual gravitational theory. We are not going to explore this issue further here.

Unequal time correlation function
In this subsection, we compute the unequal-time correlation function of the primary fields G n (x 1 , x 2 ). We note that unequal-time correlation functions of the vacuum state, unlike-their equal-time counterparts, display non-trivial dynamics 25 . In what follows, we chart out the result for these correlation functions for the CFT vacuum given by (60) The holomorphic part of this correlator yields after mapping to the complex plane, The antiholomorphic part can be computed in a similar manner with z i →z i . Thus one finally gets for h =h Defining the dimensionless center of mass and relative coordinates as x cm = π(x 1 + x 2 )/L and x rel = π(x 1 − x 2 )/L, analytically continuing to real time, and substituting Eqs. 21, 22 and 23 in Eq. 62, one obtains A straightforward analysis of G n (x cm , x rel ) in the heating phase reveals that it will display peaks, in the large For x rel = 0, we find that the peak position coincide with those of E n (x). In general, the position of the peaks trace a curve in the (x 1 , x 2 ) space and can be tuned by changing α; for α → ∞, the peaks only occurs when x cm = x rel + π/2. Thus this phenomenon constitutes another example of emergence of spatial structure in driven CFTs which was initially found by analysis of E n (x) 26 . For all pairs of values (x 1 , x 2 ) which do not satisfy Eq. 64, G n decays exponentially with n in the large n limit: ]. This behavior is clearly seen in the left panel of Fig. 8 where G n (x 1 , x 2 ) is plotted as a function of x 1 /L and x 2 /L. The position of the divergence of G n (x 1 , x 2 ) coincides with the solution of Eq. 64 as can be seen from the right panel of Fig. 8.
In contrast, the peaks of G n (x 1 , x 2 ) in the non-heating phase do not occur at an fixed positions independent of n in the large n limit. Here the divergences for occur when nθ = mπ (for integer n and m) if x 1 = x 2 ; for all x 1 = x 2 , they occur when n, x 1 , and x 2 satisfies

C. Entanglement
In this subsection, we consider the evolution of the m th Renyi entropy S m n after n drive cycles starting from |h,h states. We note that the ground state has S m n = S m 0 since the state does not change under the drive in the cylindrical geometry. This entanglement can be computed most simply by considering the correlation of the twist operators T (w,w) 40 . Here we shall concentrate on the m th Renyi entropy which is given, after n cycles of the drive, in terms of the twist operator T m (w,w) as where ρ n ( ) is the reduced density matrix of state after n cycles of the drive corresponding to an initial state |h,h , is the spatial dimension of the subsystem, w i (w i ) = +(−)ix i , we choose x 1 = 0 and x 2 = , and the twist operator T m represents a primary field with dimension h m = c(m−1/m)/24 40 . In what follows, we shall focus on the half-chain entanglement entropy which corresponds to = L/2 for the sake of simplicity; we note however, that the method can yield results for arbitrary .
As before we evaluate the holomorphic part of the α m n ( = L/2). This is given by The computation of C m 4n = C m hol 4n C m anti−hol 4n and hence S m n thus reduce to a problem similar to that worked out for the correlation function. However, here the operator dimension of T m are different from h, and hence all the expressions obtained in Sec. III B can not be directly used. Nevertheless, the computation procedure is similar, and we present the main results here. We find that C m 4n can once again be written in terms of sum of contribution over conformal blocks. In the perturbative limit, where |1 − y n |, |1 −ȳ n | 1, one has We note that this perturbative result is expected to be accurate for large n and in the non-heating phases or on the transition line as discussed earlier. Also in these cases since |1 − y n |, |1 −ȳ n | 1, only the leading term in the sum may be retained. Substituting Eq. 68 in Eq. 67 and after some straightforward algebra one obtains, assuming h =h and using h m =h m If we retain only the identity contribution which is universal, then using C aaI = 1 we can simplify, This universal contribution, after analytic continuation to real time, yields a simple expression for the evolution of δS m n = S m n − S m n=0 , given by In the non-heating phase and on the transition line, this yields, after analytic continuation to real time, Thus the oscillatory behavior of δS n in the non-heating phase and its logarithmic growth on the transition line are expected qualitative features that are reproduced in this perturbative approach. We note that the linear growth of the entanglement in the hyperbolic phase is also reproduced by this procedure; however, here we can not ascertain the accuracy of this result since the perturbation expansion of F over conformal blocks are uncontrolled. In the large central charge limit, for states with fixed conformal dimensions, we can use the global block approximation to express (assuming h m =h m and h =h) Defining α m0 n (L/2) = (2π/L) 4hm 4 −(hm+h) and δS m n = ln[α m n (L/2)/α m0 n (L/2)]/(1 − m), we find, using Eq. 73, This provides the drive induced contribution to the halfchain m th Renyi entropy after n drive cycles. Finally we consider the case for large c 0 CFTs such that h/c 0 h m /c 0 with h/c 0 and h m /c 0 held fixed. Here once again, it is possible to obtain analytic expression of α m n ( L 2 ) using the monodromy block Eq. 58. The universal contribution (from the identity block) non-perturbative in n is given by, where we have analytically continued to real time. We note here that the exact conformal block contribution can also be determined numerically using the Zamolodchikov recursion relations 41 ; however, we are not going to address this in this work.

IV. RELATION TO LATTICE MODELS
In this section, we relate our results obtain using conformal field theory in the last section to those obtained by exact numerics on a specific lattice model. The model chosen is the sine-square deformed (SSD) fermionic model whose Hamiltonian is given by 23,29 where c j is the fermion annihilation operator on site j, δ = 2π/L, L is the chain length, the lattice spacing is set to unity, we have assumed that the system to be at halffilling, and J is the hopping strength of the fermions.
In what follows, we shall use periodic boundary condition for this Hamiltonian. We note that due to the local phase factor exp[±iδj], H ± are not Hermitian operators. However, their sum is still Hermitian and leads to where J is the hopping strength of the fermions. In what follows, we shall implement the drive via a timedependent hopping J → J(t) = J 1 + J 0 cos(ω D t)) + δJ. It is well known 23,29 that the low-energy sector of Hamiltonian can be expressed as so that one can identify δf = δJ + J 1 and J 0 = f 0 .
In what follows, we shall scale all energies by J 1 and compute the time evolution of the instantaneous energy density of the system, E n (x) after n cycles of the drive as follows.
To this end, we first compute U (T, 0) The procedure for this identical to the one carried out in Sec. II and involves decomposition of U into N time steps of width δt = T /N : U (T, 0) = j=0..N −1 U j , where U j = U (t j + δt, t j ). The width δt is chosen such that H SSD does not vary appreciably within this interval. One then diagonalizes H j and expresses U j in terms of its eigenvalues and eigenvectors. The matrix U is then constructed by taking product over all U j s. Finally one diagonalizes U to obtain the Floquet eigenvalues F SSD m and eigenvectors |m SSD . In terms of these after n drive cycles, the wavefunctions of the driven chain can be written as n E(L) where c m = m SSD |ψ init denotes the overlap of the Floquet eigenstates with the initial state. The initial state is chosen to be one of the primary CFT states; for the lattice model studied here, these states are tabulated in Ref. 42. We note since f 0 1, the initial primary states corresponding to H SSD is expected to be accurately described by those charted in Ref. 42. Here we choose the state corresponding to h =h = 1/2 which for the lattice correspond to the state where |FS is the half-filled Fermi sea. Thus |ψ init corresponds to two particles populating the lowest available energy states over the half-filled Fermi sea 42 . One then computes the instantaneous energy density at any given site as E n (j) = ψ init |h j (n)|ψ init , where h j (n) = U † (nT, 0)h j (t = 0)U (nT, 0). In what follows we shall study the behavior of E n (j) and relate it to the corresponding CFT results obtained in Sec. III A. The results obtained from this procedure is shown in Fig. 11 for a chain of length L = 200. The top left panel shows a plot of E(L) as a function of n for ω D = 40J 1 / for several representative values of δf . At this frequency, the transition line is at δf c 1. One therefore finds that below a threshold number of drive cycles, n c , E(L) as obtained from the lattice shows universal behavior for both the phases and mimics the behavior of the system on the transition line. Moreover, these results show an excellent match with the corresponding CFT results outlined in Eqs. 41, 43, and 44 (with L in Eq. 1 identified to the chain length of the lattice Hamiltonian). The long time behavior of E(L) as a function of n in all the phases is shown in the top right panel of Fig. 11. We find that the lattice model reflects a clear distinction, as predicted by CFT, between the behavior of E(L) in the heating and non-heating phases at long times (large n). The bottom left panel shows the energy density E(x) as a function of x/L in the heating phase corresponding to ω D = 40J 1 / and δf = 5π/L after n = 1500 cycles of the drive. We find that the peak positions predicted by CFT are correctly captured by the lattice model; however, the peak amplitudes are lower which is due to lattice effects that are expected to cause deviation of lattice dynamics from the CFT results in the heating phase at large n 24 . Finally, in the bottom right panel of Fig. 11, we find that n c diverges at the transition; this divergence seems to coincide with the predicted 1/ |1 − α 2 | behavior since α ∼ 1/δf for ω D π/L 1.

V. DISCUSSION
In this work, we have studied the dynamics of driven CFTs using a continuous protocol. Our analysis shows that such a drive protocol, characterized by amplitudes f 0 , δf and frequency ω D yield heating and non-heating phases separated by transition lines. Such phases were obtained for discrete protocols earlier in Refs. 22 and 23 where the evolution operator U (T, 0) admits an exact analytic solution. In contrast, for continuous drive protocol, there is no exact analytic result for U . We therefore first present the phase diagram in the limit of large drive amplitude as a function of δf and ω D showing several re-entrant transitions between the heating and the nonheating phases. We also develop an analytic, albeit perturbative, approach to these driven systems using FPT. We find that for ω D ≥ δf, 1 (in units of π/L), the first order FPT results provide an excellent match with the numerical results. Our analysis allows us to identify a parameter α as a function of f 0 , δf , ω D whose values determine the phase of the system; for |α| < (>)1, the system is in the non-heating (heating) phase. The transition lines correspond to α = ±1.
We have also investigated the return probability, energy density, correlation function and the Renyi entropies of the driven CFT and have provided perturbative analytic expressions for several quantities using FPT. These expressions provide reasonable match with exact numerics. We show that the return probability P n of a primary state displays decaying (oscillatory) behavior in the heating (non-heating) phase as expected. On the transition line P n shows a power-law decay with n. Our analysis identifies a crossover stroboscopic timescale n c ; for n ≤ n c , the return probability shows a universal behavior analogous to that of a system on the transition line.
We show that n c ∼ 1/ |1 − α 2 | and can thus be tuned by changing both δ and ω D .
For energy density of a primary state, we find emergence of spatial structure as identified earlier in Ref. 25. The peaks of the energy density for large n occurs at L/4 and 3L/4 when the system is deep inside the heating phase (|α| 1); they move towards L/2 as one approaches the transition line (α = ±1). In contrast, there are no such peaks in the non-heating phase which are independent of n in the large n limit; here E n (x) shows an oscillatory behavior as a function of x for all n. For small n, we find that E n (x) obeys universal behavior similar to that when the system is on the transition line and identify a crossover time scale till which this behavior persists. However this phenomenon does not occur if x corresponds to the position of the peaks in the heating phase or on the transition line.
We have also computed the equal-time correlation functions of primary fields, C n (x 1 , x 2 ) of the driven CFT starting from a primary state. These correlation function requires evaluation of the four-point function of the CFT and are therefore expressed in terms of F which admits decomposition into Virasoro conformal blocks V p . The analytical expression of V p for arbitrary (h, h p ) does not exist. Here we have identified several limiting case where analytical results may be presented. The first of these is the case when the cross ratios that appear in the argument of F (|1 − y n | and |1 −ȳ n | in our case) are small. This limit is applicable for large n if |x 1 − x 2 | L for all phases; it is also applicable for all x 1 and x 2 in the non-heating phase and on the transition line provided n is large. Our analysis in this limits shows emergent peaks in the hyperbolic phase analogous to the energy density. In addition, we also find emergent spatial structure indicating a line of divergence in the non-heating phase and on the transition line. We also provide analytic expression using large c 0 limit for the block with fixed h,h where the Virasoro blocks can be replaced by the global conformal block. Finally, we note that using monodromy methods, it is possible to find analytic expressions of the correlation functions in the large c 0 when the dimension of the primary state H h with H/c 0 and h/c 0 held fixed. We point out that this regime may be relevant for large c 0 CFTs used in AdS/CFT correspondence.
The structure of the unequal time correlator G n (x 1 , x 2 ) in the presence of the drive provides another example of emergent spatial structure. We note that here one can study the dynamics starting from the cylinder vacuum state since the unequal-time correlation function for such an initial state, in contrast to its equal-time counterpart, shows non-trivial evolution. Our analysis shows that G n , in the heating phase, diverges along a curve in the x 1 , x 2 plane; we provide an analytic expression for this curve within first order FPT which shows reasonable match with exact numerics. We also find that in contrast to E n (x), G n also shows divergences along a curve in the non-heating phase; the shape of this curve depends on n through Eq. 65. This constitutes an ex-ample of emergent spatial structure in the non-heating phase which does not exists for E n (x).
Next, we provide a computation of the half-chain entanglement entropy (m th Renyi entropy) for the driven CFT staring from a primary state. The computation of S m n is similar to that of the equal-time correlation function since it can indeed be viewed as equal time correlation function of the twist operator T m with conformal dimension h m . We express S m n in terms of the conformal blocks and discuss limits in which their analytic expressions are available. Such a limit constitutes the case of long-time (n 1) limit of S m n in the non-heating phase and on the transition line. Here we show that universal contribution to F (and hence S m n ) from the identity block has oscillatory dependence of n in the non-heating phase and a logarithmic growth on the transition line. We also provide analytic expression for S m n in the large c 0 limit where the Virasoro blocks can be replaced by global conformal blocks. Finally for h h m , and c 0 1 with h/c 0 and h m /c 0 held fixed, we find analytic expression for S m n (L/2) using monodromy methods; our results here may be relevant to CFTs used in AdS/CFT correspondence.
Finally, we point out that our results show excellent match with exact numerics of a 1D lattice model of fermions on a finite chain of length L with Hamiltonian H SSD (Eq. 77). In particular, we find that the emergent spatial structure of the energy density in the heating phase at long times and its universal behavior below a crossover scale n c is accurately reflected in such lattice dynamics. We note that we study the system in the presence of a global drive. In a typical lattice system which obeys Galilean invariance, such a drive does not usually lead to emergent spatial structure of correlation functions or energy densities. The fact that we find such an emergent structure here clearly shows the necessity of a CFT based interpretation of such a dynamics where space and time are intertwined 24 .
In conclusion we have studied driven CFTs using a continuous periodic protocol and have provided a phase diagram showing re-entrant transitions between heating and non-heating phases. We have also studied the return probability, energy density, correlation functions and Renyi entropies of such a driven CFT starting from primary states. Our results indicate several features of these quantities such as the universal behavior of the return probability and the energy density below a crossover stroboscopic timescale and emergence of spatial structure in both heating and non-heating phases as found in the correlation function of primary fields. We discuss relations of these results to a recently studied lattice model and find excellent match between exact numerical lattice model based results with analytic prediction of the CFT.

VI. ACKNOWLEDGEMENT
The authors thank Shouvik Datta and especially, Koushik Ray for several stimulating discussions. RG acknowledges CSIR SPM fellowship for support. DD acknowledges supports provided by SERB, MATRICS and Max Planck Partner Group grant, MAX-PLA/PHY/2018577.

Appendix A: Floquet perturbation theory
In this appendix, we provide details of the Floquet perturbation theory used in the main text. We begin our analysis starting from Eq. 1 of the main text with f (t) given by Eq. 8. We shall put π/L = = 1 in this section.
In the limit of large f 0 , U 0 , the zeroth order term in the perturbative expansion of U is given by , Since sin(ω D T ) = 0, where T = 2π ω D is the time period of the drive, we find (restoring π/L and ) U 0 (T, 0) = e −iδf T σz , H F = sσ z /T (A2) where s = arccos(cos(δf T )) (where π/L is set to unity) is defined in the main text. We note that H F reflects the periodicity of U .
The first order term for U in the perturbation expansion is given by Since H 1 ∼ iσ y and U 0 only depends on σ z (Eq. A2), a straightforward calculation yields U 1 (T, 0) = = −i(I 1 σ + + I * 1 σ − ) where J m (x) are m th Bessel functions. This leads to Eq. 18 and then, following the unitarization procedure discussed in the main text, to Eq. 19 for H Once again using the Pauli matrix dependence of H 1 and U 0 we find From Eq. A7 we note that U 2 − U 2 1 /2 depends only on Im[I 2 ]. Thus to second order in perturbation theory, U 2 is given by U 2 = U 0 (I + U 1 + (U 2 − U 2 1 /2)) = e −iH (2) F T = e −is (1 + iβ) −iα sin y iα sin y e iy (1 − iβ) where α is defined in Eq. 9 in the main text. We unitarize U 2 following the same procedure as H (2) F = θ (2) n (2) z σ z + in (2) y σ y /T sin(θ (2) T ) = (sin s − β cos s) 2 − α 2 sin 2 s n (2) z sin(θ (2) T ) = sin s − β cos s n (2) y sin(θ (2) T ) = α sin s (A9) The phase diagram is obtained from Eq. A9 by imposing condition on Tr[exp[−iH (2) F T ]] as discussed in the main text. This translates to the conditions | cos(θ (2) )| > (<)2 for the heating (non-heating) phases and | cos(θ (2) )| = 2 on the transition line. The phase diagram, shown in Fig.  12 as a function of δf and ω D turns out to be qualitatively similar to that obtained using first order FPT (Fig. 1 in the main text).
To obtain the next order correction to U , we write To evaluate this, we use the Baker-Campbell-Hausdorff relations for L 0 and L ±1 which states e iaL0 (L 1 + L −1 )e −iaL0 = e ia L 1 + e −ia L −1 (B7) Using Eq. B7 and identifying a = −2π/L(f 0 sin(ω D t)/ω D + δf t) (Eq. B5), one can evaluate U 1 (T, 0) in a straightforward manner. The integrals involved are similar to those in App. A and one obtains where α is given by Eq. 9 of the main text. Thus the expression of the evolution operator, up to first order in perturbation theory, is given by We note that this expression matches with Eq. 18 in the main text for L 0 = σ z /2 and L ±1 = ∓σ±.
Next we briefly comment on the unitarization procedure to be followed here. Since U 1 is non-unitary within first order FPT, we would like to unitarize it. This seems difficult without using the SU(1,1) representation of the Virasoro operators. Here we therefore concentrate on the high frequency regime, where the two routes to unitarization of U 1 , discussed in the main text, leads to identical result. In the high frequency regime, s → 0 and sin s/s → 1. Furthermore, in this regime once can write U 1 1 − iH F T / . Thus expanding Eq. B9 in powers of s and retaining only first order terms, we find where n z , n y , and θ are given by Eq. 19 of the main text. We note that H F coincides with that obtained within first order FPT obtained in the main text with the identification L 0 = σ z /2 and L 1 + L −1 = iσ y and is also consistent with the dynamic emergent symmetry discussed in the main text.