On periodically driven AdS/CFT

We use the AdS/CFT correspondence to study a thermally isolated conformal field theory in four dimensions which undergoes a repeated deformation by an external periodic time-dependent source coupled to an operator of dimension Delta. The initial state of the theory is taken to be at a finite temperature. We compute the energy dissipated in the system as a function of the frequency and of the dimension Delta of the perturbing operator. This is done in the linear response regime. In order to study the details of thermalization in the dual field theory, the leading-order backreaction on the AdS black brane metric is computed. The evolution of the event and the apparent horizons is monitored; the increase of area in each cycle coincides with the increase in the equilibrium entropy corresponding to the amount of energy dissipated. The time evolution of the entanglement entropy of a spherical region and that of the two-points function of a probe operator with a large dimension are also inspected; we find a delay in the thermalization of these quantities which is proportional to the size of the region which is being probed. Thus, the delay is more pronounced in the infrared. We comment on a possible transition in the time evolution of the energy fluctuations.


Introduction
The framework of gauge/gravity duality is a particularly useful setup to study how classical gravitational singularities are resolved in a quantum theory of gravity. Decay of a false vacuum in AdS [1] generically results in a big-crunch classical geometry. Studies have shown that some of such crunches are unresolved and lead to unbounded ill-defined boundary theories [2][3][4][5][6]. Others seem unresolved with impunity and are nevertheless expressed by well-defined boundary conformal field theories (CFT) [7,8]. The class of periodic timedependent Hamiltonians which are unbounded for half of the time and bounded the other half, provide an interesting situation that is "midway" between the case of bounded and unbounded potentials. Depending on the frequency of the time-dependent part of the Hamiltonian, it can happen that the crunch is shielded by a black hole (BH) horizon [9]. This process is dual to thermalization in the boundary field theory. Finite temperature can act as a stabilization mechanism by providing a positive thermal mass on top of possible tachyonic direction(s) in the field space. Moreover, understanding the dynamics of thermally isolated systems that are driven or quenched by an external force, modeled by a time-dependent Hamiltonian, is an active area of research in condensed matter physics (see e.g. [10] for a review). Theoretical research is motivated by many recent advances in experimental studies of cold atom systems. Near quantum phase transitions, long-distance physics of many-body quantum systems is described by conformal field theories (CFT). Many universal properties of driven quantum systems can then be modeled by quantum field theories with time dependent couplings, see e.g. [11][12][13][14][15][16][17]. It remains in general challenging to study such non-equilibrium processes at strong coupling in higher dimensions; this motivates us to investigate these issues in the framework of the AdS/CFT correspondence (see e.g. [18][19][20][21][22][23][24][25][26][27][28][29][30] for previous work on this topic).
In this paper we consider a thermally isolated field theory which is driven by a periodic external perturbation ξ(t) in the Hamiltonian: H = H 0 + ξ(t)δH , ξ(t) = ξ(t + 2π/ω) . (1.1) One can think of the periodic source ξ(t) also as the vacuum expectation value (VEV) of some field interacting with the CFT. We analyze the case in which the undeformed Hamiltonian H 0 belongs to a conformal field theory in d dimensions with an AdS d+1 dual, and δH is a relevant deformation of the form: where O ∆ is a generic relevant scalar operator with dimension ∆ > d−2 2 , i.e. above the unitarity bound; the form of the time dependence is chosen here for simplicity as ξ(t) = ξ 0 cos ωt. The initial state of the theory is taken to be at finite temperature T , and we work in the limit in which the total amount of work done on the system is much smaller than the total initial internal energy. The physical response of the CFT to the deformation depends significantly on the dimension ∆ of the perturbing operator. From dimensional analysis, in the limit ω ≫ T , the energy dissipated per unit of volume in each cycle W c is expected to scale as ω ζ where ζ = 2∆ − d; this is an increasing (decreasing) function of ω for ∆ > d/2 (for ∆ < d/2). For ω ≪ T instead W c generically scales as ωT ζ−1 ; this is an increasing (decreasing) function of T for ∆ > (d + 1)/2 (for ∆ < (d + 1)/2). Both these expectations are confirmed in strongly-coupled theories with AdS duals.
We use the AdS/CFT correspondence to study how the CFT reacts to the deformation (1.1). The initial thermal state is described by an AdS d+1 black brane geometry, and the operator O ∆ is dual to a scalar field φ with mass squared m 2 = ∆(∆ − d). The periodic deformation (1.1) is then implemented as a time-dependent boundary condition for φ. In the linear approximation, W c can be determined by computing scalar two-point functions (see e.g. [31][32][33][34]; for a review of linear responses in AdS/CFT, see e.g. [35]). In practice, this amounts to solving the wave equation for φ in the black brane geometry. We specialize to the case of d = 4 and compute W c numerically as a function of ω and ∆. In both the limits ω ≫ T and ω ≪ T , the results are consistent with the ones obtained by dimensional analysis and by the Kramers-Kronig relations. In the limit where the dimension of the driving operator approaches the unitarity bound ∆ → 1, W c becomes a sharply peaked function at a small frequency ω m ; numerically we find that ω m scales as (∆ − 1) 1/2 . While the energy W c injected into the system by the external source during one cycle tends to drive the system away from its equilibrium thermal state, the self interactions of the system lead to a restoration of thermal equilibrium by entropy production, which on the gravity side is dual to an increase of the area of the horizon 1 . In order to study the details of thermalization in the dual field theory, we compute the leading order gravitational backreaction in the expansion parameter ξ 0 on the bulk AdS 5 metric. The number of cycles is taken to be large but finite, in such a way that the total work performed on the system is small compared to the initial temperature; this can be achieved by taking the parameter ξ 0 to be sufficiently small. We indeed find that the thermalization takes place in a very efficient way.
Equipped with the calculation of the leading correction to the metric, we monitor the time evolution of a few physical observables. We follow the time evolution of the event horizon and we check and confirm that it increases with time. It is also possible to follow another quantity, the apparent horizon, which was advocated to play an important role in the CFT [45,48]. The area of the apparent horizon is not always a monotonic function of the time t for intermediate values of t [49]. However in this case we found that also the area of the apparent horizon increases. The increase of area in each cycle is the same for both definitions of horizon, and it coincides with the increase in the equilibrium entropy corresponding to a change in the internal energy given by the work done in a cycle; this is true for all the values of 1 < ∆ ≤ 4.
We next investigate the time evolution of the two-point function of a probe operator with a large dimension ∆ p , which in the holographic dual is related to the geodesic length, and of the entanglement entropy of a spherical region, which is probed by the volume of minimal surfaces. For both the quantities there is a delay in the thermalization which is proportional to the size L of the region which is probed, in the limit L ≫ 1/T . The thermalization time is longer for infrared observables; this kind of behavior has been observed before in different systems in [37][38][39][40]. The two-point function and the entanglement entropy have also an oscillatory term in time. The amplitude of this oscillation is not extensive in the spatial region which is probed, and so it is negligible in the limit of large region.
It was pointed out [50] that the energy fluctuations of driven systems are significantly affected in a universal way by the protocol by which the periodic driving force is applied. In our theoretical setting we show that these universal features are determined by ∆ and ω. The results in [50] raise the possibility of a transition in the time dependence of energy fluctuations. We discuss these features in the context of quantum field theory and we find that the operator spectrum of N = 4 Super Yang-Mills (SYM) allows for such a transition.
In section 2 we make some general considerations based on dimensional analysis which apply to a generic driven CFT; we also comment on the relation to the case of quenches. In section 3 the theoretical setting is introduced, including the equations used to compute the backreaction of the metric. Linear response functions are computed numerically in section 4. In section 5 we discuss the apparent and event horizons and we compute the leading correction to the two-point function and the entanglement entropy. In section 6 we initiate the study of energy fluctuations. Section 7 concludes with a discussion of the results. Finally, appendix A contains some technical details for the special cases of integer and half-integer ∆s and appendix B deals with the effective Hamiltonian of a free field theory for large frequencies.

Periodically driven CFT
One of the main tools in studying the reaction of a medium to a weak perturbation which tends to drive it out of equilibrium is the linear response theory. In this section we review some basic definitions and we make some universal considerations based on dimensional analysis. In the linear approximation, the response to the deformation (1.2) is encoded in the retarded Green function, which can be defined in terms of a general spacetimedependent deformation: (2.1) The source J ∆ (x) has dimension d − ∆. In the linear response regime, where J ∆ (x) is sufficiently small, the VEV of the operator is given by where G R is the retarded Green function: and θ H is the Heaviside step function. Here we will focus on a time-dependent sourcewhich is constant in space -of the form J ∆ (t) = ξ 0 Re(e −iωt ); for small ξ 0 , we will be in the linear-response regime where The energy per unit volume, which has dissipated in a period 2π/ω, denoted by W c , is proportional to ξ 2 0 Im G R (ω, 0).
The functions Re G R (ω, 0) and Im G R (ω, 0) are respectively even and odd (this is a consequence of the fact that G R (ω, 0) is the Fourier transform of a real function). For frequencies ω ≫ T , the behavior W c ∝ ω ζ where ζ = 2∆ − d is expected from dimensional analysis. This expectation has been confirmed in CFTs with an AdS dual, see e.g. Appendix A of Ref. [31]; this will be reviewed in sections 3 and 4, see eqs. (3.41,4.13,4.14). In the same class of theories, in the case of integer ∆s in even dimensions and of half-integer ∆s in odd dimensions, Re G R (ω, 0) instead scales as ω ζ log(ω/T ), while it scales as ω ζ for generic dimensions, see eqs. (4.13,4.14). For d−2 2 < ∆ < d 2 the energy dissipated per unit volume in a cycle, W c , is a decreasing function of ω for large ω ≫ T . In the opposite limit ω ≪ T one expects that Re G R (ω, 0) approaches a constant (the equilibrium value), while Im G R (ω, 0), being an odd function of ω, is expected to generically start off as ω T ζ−1 . This behavior will be confirmed in the numerical analysis of section 4. Examples of the behavior of W c are shown in figure 1. In the case of a free field theory, the retarded propagator is: where P denotes the principal value. In the limit m → 0, the two δ functions collide at (ω, k) → 0. We will find a similar behavior in AdS/CFT in the limit ∆ → 1, with the difference being that the δs are smeared into a continuous peak. Previous studies on driven systems have focused on the case of quenches, where the time dependence of the coupling interpolates between two different asymptotic Hamiltonians, i.e. ξ(t) = ξ i + (ξ f − ξ i )θ H (t). The effective Hamiltonian approach [51], [9] suggests that for large frequencies, ω, the effect of introducing a periodic deformation can be described as a quench. A simple case where ∆ = 2 is discussed in appendix B, where we show that the averaged long-time behavior is that of the quench to leading order in 1/ω 2 ; in this approximation the energy dissipation is neglected.
The authors of [22,26,30] recently studied quenches in AdS/CFT as a function of the operator dimension ∆, in d = 4. In order to approach the limit of θ H (t) in a smooth way, they considered the class of functions: where t 0 is a time scale, as a regularization of θ H . In the limit t 0 → 0, they found by numerical analysis that the energy dissipated in the quench is a divergent quantity for ∆ ≥ 2, which scales as 1/t 2∆−4 0 (as log 1/t 0 for ∆ = 2). We point out that this divergence is a direct consequence of the fact that W c scales as ω ζ in the large ω limit of the periodically driven system; the quantity 1/t 0 acts as a UV cutoff of the Fourier modes of ξ(t). The Fourier transform of the source f (t) is: In the case of ζ > 0, the energy dissipated at small t 0 , corresponding to a very fast quench, scales as . This is in agreement with the numerical result obtained in [22,26]. In the case of ζ < 0 there is no divergence, so the energy dissipated in the process should approach a constant in the limit of a sudden quench. For ζ = 0, which in d = 4 corresponds to ∆ = 2, the divergence is logarithmic. We will now turn to the study of periodically driven strongly coupled systems with the aid of the AdS/CFT correspondence.

Theoretical setting
In this section we review work on the gravity dual of a conformal field theory perturbed by a relevant scalar operator of dimension ∆ and then adopt it to the purpose of studying the time dependent periodic driving force dealt with in this work. We use as an expansion parameter the magnitude of the amplitude of the time-dependent perturbation. The gravitational theory is Einstein gravity coupled to a massive scalar field: We consider generic values of the mass m 2 ≤ 0, which on the boundary corresponds to different operator dimensions ∆(∆ − d) = m 2 , giving the usual two solutions ∆ ± = d/2 ± m 2 + d 2 /4. In the range m BF = −d 2 /4 ≤ m 2 ≤ −d 2 /4 + 1 we can either set ∆ = ∆ + or ∆ = ∆ − ; for m 2 > −d 2 /4 + 1 we are forced to set ∆ = ∆ + , in order to respect the unitarity bound for the dimension of a scalar operator. We ignore higher order terms in the potential because they only affect dynamics beyond the leading order in our expansion parameter.
For example, in the case of the AdS 5 dual of the N = 4 theory, we may consider an operator O 2 which gives mass to the scalar fields, or an operator O 3 which gives mass to the fermions. The operator O 2 corresponds to a field φ in AdS 5 with m 2 = −4, while O 3 corresponds to a field with m 2 = −3. The value of the Newton constant G N is related to the number of colors N c of the N = 4 theory: G N = π/(2N 2 c ). In the following we will consider the four dimensional case in detail.
We choose to study the CFT in flat Minkowski space and so we work in Poincaré patch. In order to model a translational invariant state in the boundary CFT, we consider the following metric: For convenience, we work in Eddington-Finkelstein (EF) coordinates. The AdS boundary is at r → ∞, where A(v, t) → r 2 , Σ → r. The variable v, for large r, corresponds to the time of the boundary theory. The Einstein equations read [22,45]: where h ⊙ ≡ḣ + Ah ′ /2 and ′ and˙denote derivatives with respect to r and v. The Klein-Gordon equation reads: We will solve these equations in a perturbative expansion, taking as a zeroth order solution a black brane with φ = 0. We will use as expansion parameter the amplitude of the field φ.
The initial thermal state of the boundary CFT is dual to a black brane geometry, which corresponds to: where T is the temperature. The expression in Schwarzschild coordinates is: which can be obtained from EF by changing variables from v to t, with where the integration constant is chosen in such a way that v = t for r → ∞. The location of the horizon, r 0 , is given by the value of r for which A vanishes; the temperature is then T = A ′ (r 0 )/(4π).

Klein-Gordon equation
We will expand around the black brane metric, taking the amplitude of the scalar φ as an expansion parameter. It is useful to parametrize the scalar field as follows: where λ is an expansion counting parameter, which is proportional to the amplitude of the driving force ξ 0 on the field theory side. The functions A, Σ receive corrections only at order λ 2 . It is convenient to introduce the variables: The equation of motion for the scalar φ on the black brane background, to leading order in λ, is: Note that this equation contains only a first order time derivative, while in Schwarzschild coordinates it contains also second order time derivatives 2 . In practice eq. (3.10) will be solved numerically, see section 4. In order to set the correct boundary conditions and also to compute one-point functions, it will be useful to expandφ as a power series in ρ, around the boundary at ρ = 0. If m 2 is such that ∆ ± is neither integer nor half-integer, the scalar φ can be written as: The case of integer and half-integer values of ∆ ± is treated in appendix A. The following iterative relations [26], which are valid for both ∆ = ∆ ± , can be derived from the Klein-Gordon equation: where˙denotes a derivative with respect to τ . There is no a priori algebraic relation between a ∆ − +j and a ∆ + +j ; the relative coefficient is fixed by the requirement that the field φ has ingoing-wave boundary conditions at the horizon. 2 In Schwarzschild coordinates, (rescaling ρ = µ/r, τ = µt), the equation of motion is:

Backreaction to leading order
Next we compute the backreaction of the scalar φ on the metric at leading order in the expansion parameter λ. The metric functions A, Σ receive corrections at order λ 2 : 13) and the corresponding metric at this order reads: (3.14) From (3.3), one obtains a set of linearized equations forÃ 2 ,Σ 2 [26]: Here ′ and˙denote derivatives with respect to ρ and τ . These equations have source terms, which are due to insertions ofφ 1 = Re(φ 1 e −iω T τ ). The solutions forΣ 2 to the homogeneous version of eq. (3.15) are a constant and 1/ρ; the term proportional to 1/ρ corresponds to a linear rescaling in the ρ coordinate, while the constant inΣ 2 can be gauged away using the residual diffeomorphism invariance ρ → ρ + f (τ ) and thus we can set both these terms to zero. Eq. (3.15) gives all the information that we need to solve forΣ 2 . We can then insert Σ 2 into eq. (3.16) and solve forÃ 2 ; the solution to the homogeneous version isÃ 2 ∝ ρ 2 , so this equation fixesÃ 2 only up to a function of time which is proportional to ρ 2 in space; this is fixed by eq. (3.18). We show that this term is actually a linear function of time τ . After reviewing the general equations for the leading order backreaction [26], we specialize to the case of a time-dependent periodic source. It is not consistent within the approximation to consider a source that is switched on for an infinite number of cycles, because at some point the total backreaction will grow large, no matter how small the amplitude λ is taken. We can hence imagine switching on the source at time τ * 0 and switching it off after a large number of cycles at time τ * f . This while keeping λ as small as needed in order to trust the perturbative expansion. We neglect also the effect of the transient period, which should be negligible in the limit of a large number of cycles. In order to make the time dependence of the problem explicit, we can use the following Ansatz, valid in the time window τ * 0 < τ < τ * f : It turns out thatÃ 2 , as a function of τ , contains a constant term A 2,c (ρ), a linear one A 2,l (ρ), and a periodic part A 2,p (ρ);Σ 2 contains a constant term Σ 2,c (ρ) and a periodic one Σ 2,p (ρ). The metric (3.14) resembles the AdS Vaidya metric: which is a solution of the Einstein equations in the presence of the energy-momentum tensor with the other components being zero; this describes the absorption of null dust by a black brane. If we set Σ 2,c,p = A 2,c,p = 0 in our metric, then we would recover a Vaidya metric with M (τ ) ∝ τ . Several studies of non-equilibrium physics in the AdS/CFT have focused on the Vaidya metric, see e.g. [36][37][38][39][40][41][42][43][44].
The function A 2,l can be obtained by solving eq. (3.18): This expression is proportional to ρ 2 , which can be shown by forming a linear combination of the real and imaginary part of eq. (3.10). The functions Σ 2,c,p can be obtained by solving the following equations: with the boundary condition that the constant and the 1/ρ terms in the expansion at small ρ of Σ 2,p,c are zero. The periodic part A 2,p can be expressed in terms of φ 1 and Σ 2,p using eq. (3.18); the result is The constant part A 2,c is obtained by solving the following differential equation: Let us first focus on the case of a generic m 2 (except for the cases of integer and halfinteger ∆s, which are treated in appendix A); the following form of the series expansion in ρ can be used [26]: where the αs and σs are functions of time. Among all these coefficients, α 0 is special; note that eq. (3.16, 3.17) contain only spatial derivatives of A, and in a combination which is blind to the coefficient α 0 . The only equation that fixes it is eq. (3.18). The coefficient α 0 is the only one with a secular term which increases linearly in time; all the other coefficients in eqs. (3.27,3.28), are superpositions of a constant and a periodic term in time, with period 2ω T : From eq. (3.15) one obtains (j ≥ 0 is an integer): Replacing some of the coefficients of eqs. (3.30) in eq. (3.18), a differential equation which gives the time evolution of α 0 can be found [22,26]: The coefficient α 0 has an important physical significance: its variation in time is proportional to the total work performed on the system by the external force. Both Σ 2,c,p are divergent at small ρ for ∆ + > 7/2; in this regime it is useful to use the small ρ expansion to set up the right boundary conditions for solving eq. (3.24): For ∆ + > 3, both A 2,c and A 2,p diverge at small ρ; it is useful to use the small ρ expansion to set up the boundary conditions for the differential equation (3.26): In this subsection we have studied the leading-order response of the metric to the scalar perturbation; the general form of the O(λ 2 ) correction to the metric is parametrized by eqs. (3.19,3.20) in terms of five functions A 2,l,c,p (ρ), Σ 2,c,p (ρ); we have shown that the function A 2,l (ρ) is quadratic in ρ. In order to fix the boundary conditions and also to compute the one-point correlators, we have expanded the metric functionsÃ 2 ,Σ 2 as power series in ρ, see eqs. (3.27,3.28), and we have written recursion relations for the coefficients of these expansions.

One-point correlators
Once the background metric and the scalar field solutions have been determined, one can derive the expectation value of the boundary scalar operator O ∆ and of the energymomentum tensor. The parameter ξ 0 which parametrizes the source J ∆ = ξ 0 cos ωt in the boundary field theory is proportional to the expansion parameter λ on the gravity side; by dimensional analysis ξ 0 = λµ 4−∆ (in our conventions we set the proportionality coefficient equal to one) . In the case in which the dimension of the source, ∆ = ∆ + , we use the following parametrization: where the Klein-Gordon equation for the field φ is solved using an ingoing boundary condition at the horizon [31]. Holographic renormalization [53] can be used to obtain an expression for the expectation values of the energy-momentum tensor. Explicit expressions for generic ∆ = ∆ + were found in [26]: In our case, the time dependence of E and P has a periodic part with frequency 2ω T and a linear part in time, which corresponds to the work done on the system. These expressions are consistent with the Ward identity In the case where the dimension of the source, ∆ = ∆ − , we use the parametrization: The expressions (3.36,3.37) for O ∆ and E are valid also for ∆ = ∆ − ; this can be checked by inspecting the consistency of eq. (3.32) with eq. (3.39). The cases with integer and half-integer dimensions are treated in detail in appendix A. In our conventions, the retarded Green function and the function A 2,l (both in the case where ∆ = ∆ + and ∆ = ∆ − , as well as for ∆ integer and half-integer, with the exception of ∆ = 2 case) are: For ∆ = 2 these quantities are given in eq. (A.12). The work done on the system per unit volume and time is: In our conventions, this expression is valid for any dimension ∆. Specializing to different values of the operator dimension: (3.43)

Linear response
We are now prepared to study the physical properties of the system, by numerically computing the response function χ ∆ (ω T ) for any values of ω T and for operator dimension 1 < ∆ ≤ 4. We check the numerical calculations using the Kramers-Kronig relations and comparing with the analytical results in both the limits ω T ≫ 1 and ω T = 0.

Numerical study
In order to study the time-dependent case numerically, let us write the field φ 1 in the following form: where the integer M is chosen in such a way that the sum includes (for a given m 2 ) at least all the exponents which are smaller or equal to ∆ + . The first coefficients a ∆ − +k read: By substituting the Ansatz of eq. (4.1) into eq. (3.10), one finds an equation of the kind: where J 0 is a source term, whose explicit form depends on the number of terms M in eq. (4.1). For example, if M = 3 The equation (4.4) can be solved by a shooting method, imposing that the solution is regular at the horizon and that the expansion at order ρ ∆ − ofφ 1 vanishes at ρ = 0. The value of χ(ω T ) can then be extracted from the coefficient of the term proportional to ρ ∆ + inφ 1 , using eq. (3.35). In the case of ∆ = ∆ − , we can use the same numerical solution that we used for ∆ = ∆ + , using instead eq. (3.40) to compute χ(ω T ). Numerical results with ∆ = ∆−. In the limit ∆ → 1 the function Im χ becomes sharply peaked at a low frequency.     proceeds as follows: In order to find χ numerically, it is convenient to use: Numerical results are shown in fig. 6. the real part is proportional to log ωT at large ωT , the imaginary part goes to a constant.
• m 2 = −3, where ∆ − = 1, ∆ + = 3 As a tool to enforce the desired boundary conditions, it is then convenient to introduceφ 1 : Numerical results are shown in fig. 5. There is a discontinuity in the real part of χ as a function of ∆, when it approaches the value ∆ = 3.
• m 2 = 0, where ∆ − = 0, ∆ + = 4, which is the case of deformation by a marginal operator. For the numerical calculation we set: and we proceed as in the previous cases. The numerical results for this case are shown in fig. 5.

Checks
A cross check on the numerics comes from the Kramers-Kronig (KK) relations: where P denotes the principal value. In the case of ζ < 0, we can apply these relations directly, as χ(ω, 0) vanishes for large ω. In the case of ζ ≤ 0, we should first subtract the divergent part at large ω T ; this is more difficult, because we need to compute χ with a rather good accuracy. Practically, we managed to check KK relations only up to ∆ < 3.5.
For ω T = 0, the Klein-Gordon equation (3.10) can be solved in terms of hypergeometric functions in the static case; the solution with an ingoing-wave boundary condition at the horizon is [54]: From this exact solution the following relation can be extracted [26]: The equilibrium value of the response function χ ∆ (0) can be found using eqs. (3.35,3.40).
The large ω T = ω/(πT ) regime is obtained in the small temperature limit, which corresponds to the vacuum AdS metric at T = 0. Calculations are done for example in [31,52]: The case of integer ∆s is special and contains logarithms: 14) We find agreement at the per-mille level between the numerical and analytical results both in the ω T ≫ 1 and in the ω T = 0 regimes. In the limit ∆ → 1, the theory is just above the unitarity bound; this is reflected in a singular behavior of Im χ(ω T ), which becomes a function which is sharply peaked at a very low value of the frequency ω T,0 . The position of the maximum ω T,0 tends to zero for ∆ → 1; by using a numerical fit we find that ω T,0 ≈ 0.966(∆ − 1) 1/2 .
The function Im χ ∆ is generically expected to start off with a linear term in ω T for ω T ≪ 1, because it is an odd function of ω T . The numerical calculation confirms this behavior.

Apparent horizon
We will next use the calculation of the metric backreaction to study the time evolution of the area of the horizon, which is related to the entropy in the boundary CFT. It is still an open question whether a valid notion of local entropy density exists in out-ofequilibrium situations. It has been argued [55] that the area of the horizon, projected onto the boundary along an in-falling null geodesic (which corresponds to a line with constant v in EF coordinates) could be identified with the entropy density in the dual field theory. However, it has been suggested in [45,48] that apparent horizons provide a better definition of entropy density than global horizons. We will compute and compare the area of both the apparent and the event horizon on solutions taking the leading order backreaction into account.
Let us start with the apparent horizon; it is associated with the presence of a trapped surface. Considering the following light-like vectors: written in the coordinates (v, r, x), we compute the expansion of the null geodesics: where w = l, n. It follows that where ′ ,˙are derivatives with respect to r, v. The position of the apparent horizon is detected by the vanishing of θ n . Using the variables ρ, τ , the location ρ a of the apparent horizon to order λ 2 is determined by: ( At this order one can write ρ a = 1 + λ 2 δρ a + O(λ 4 ) and the variation of the position of the apparent horizon ρ a is given by: where the functions A 2,c , A 2,l , A 2,p , Σ 2,c , Σ 2,p are evaluated at ρ = 1. The area of an element of the apparent horizon A a at order λ 2 is then: Using the expression (3.25), we find: The area of the apparent horizon always increases as a function of time, independently of the value of φ 1 (1) and ω T ; the relative coefficient of the periodic and of the linear part in time is such that in every cycle there is an instant at which the time derivative of A a vanishes. The increase of area in a cycle is proportional to ω T |φ 1 (1)| 2 . Some plots of |φ(1)| as a function of ω T for various values of ∆ are shown in figs. 7,8. By comparing eq. (3.23) and eq. (3.41), we find that the value of |φ 1 (1)| 2 is directly related to Im χ: (in our conventions, this expression is valid for ∆ = 2; for ∆ = 2, it gets modified to |φ 1 (1)| 2 = Im χ 2 (ω T )/ω T ). The increase in area and in energy in a cycle are given by: Taking the area of apparent horizon [45,48] as a representative of entropy density (S a = A a /(4G N )), we see that the variation of entropy in a cycle is: δS a /S a = 3 4 δE/E. This is the same variation in entropy as one would expect from the equation of state S ∝ E 3/4 of the undeformed CFT at equilibrium for a variation in internal energy equal to the work done on the system. This is consistent with thermalization. It should be stressed that this result can be made valid for arbitrary values of ω T by considering the limit of a sufficiently small amplitude ξ 0 . This is done in a way that keeps the total work per unit of volume performed on the system much smaller than the initial energy density. For ω T ≫ 1, this does not overlap with the regime of validity of the hydrodynamic approximation, which instead is an expansion around ω T = 0 and is valid for arbitrary amplitude ξ 0 . The hydrodynamic description of a driven system in the case of ∆ = 4 was studied, at second order in the boundary derivative expansion, in [18]; this corresponds to keeping just the linear and the quadratic order in χ 4 (ω T ), see eq. (3.19) in [18].
The quantity Arg(φ 1 (1)) gives the difference between the phase associated with the source and that given by the moment at which the time derivative of A a vanishes. In the limit ω T → 0, the evolution of A a is in phase with the source and then Arg(φ 1 (1)) → 0; this is the limit in which the hydrodynamic approximation is also valid. Some plots as functions of frequency, for several values of ∆, are shown in figs. 9,10; from the numerical calculations, we find that this phase difference approaches a constant value (which can have either sign depending on ∆) at large ω T .

Event horizon
The event horizon corresponds to a null geodesic which separates trapped null geodesics from the ones that can escape to the boundary. In the unperturbed geometry, it sits at ρ(τ ) = 1. At order λ 2 , the geodesic equation receives the correction: Let us expand the position of the event horizon as ρ e = 1 + λ 2 δρ e (τ ) + O(λ 4 ); the leading correction ρ e satisfies: whose solution is: Note that one could add a function of the form Ce 2τ , with C an arbitrary constant, to the solution δρ e in eq. (5.14), and still get a solution of eq. (5.13). The curves with C > 0 correspond to trapped null geodesics, while the ones with C < 0 correspond to null geodesics which can escape to the boundary. This is different than the position of the apparent horizon; the only term which is the same is the one linear in time. The area of the event horizon A e is (5.15) which, using (3.25), can be written as: and it increases monotonically in time, as expected from the area theorem of general relativity. The time derivativeȦ e is always strictly positive; this is different from the area of the apparent horizon A a , for whichȦ a = 0 once for every cycle. The difference between the two areas is: which is bigger than zero. This is consistent with the expectation that the area of event horizon is bigger than that of the apparent one at a given moment in time. A comparison in an illustrative example is shown in fig. 11. Since the two definitions of horizon area show the same linear growth in time, the corresponding increase of entropy density δS = δA/(4G N ) in a cycle is in both cases the same as the one given by the equation of state of the undeformed CFT: δS/S = 3 4 δE/E. It should be emphasized that the calculation of the position of the event horizon in eq. (5.14) is valid only in the regime where the total work per unit of volume done on the system remains much smaller than the initial energy density. It is also important that a large number of cycles is performed by the source in such a way that the initial transient period can be neglected. Both these conditions can be achieved by keeping the amplitude ξ 0 of the source sufficiently small.

Geodesics and equal-time two-point functions
We will probe the back-reacted geometry by computing the renormalized equal-time Wightman function of an irrelevant scalar operator with large dimension ∆ p . This quantity is directly related to the space-like geodesics of the geometry which join points of the same time at the boundary. For calculations in the case of the Vaidya metric see [40,41,44].
Let us first discuss the λ = 0 case, and parametrize the geodesics with the coordinate X and the functions τ g (X), ρ g (X). We recall that all the distance and time scales are measured in units of the inverse temperature T −1 . The geodesic length at zeroth order reads: Translation invariance of τ and X gives the conserved quantities: Geodesics which join points of the same time at the boundary satisfy the relation τ ′ g (ρ * g ) = 0, where ρ * g is the value of ρ g at the midpoint; for this class of geodesics K = 0 and This family of geodesics can be parametrized by the boundary length in the X direction 2L or, alternatively, by the maximum of ρ g , which is ρ * g ; this class of geodesics does not penetrate the horizon (ρ g < 1, and ρ * g → 1 for L → ∞). We consider geodesics stretching from X = −L and X = L at the boundary; the integration constant for τ g is chosen in such a way that τ g (±L) = 0. Some examples of these geodesics are shown in fig. 12.
Near the boundary, ρ → 0, the asymptotics of the geodesics are: The geodesic length is then a divergent quantity; it is convenient to define a renormalized length L 0,R = L 0 + 2 log(ρ cut /2) , (5.22) where ρ cut is a cut-off. The renormalized equal-time Wightman function of an operator O ∆p with large dimension, ∆ p ≫ 1, is proportional to: where x = (L, 0, 0) and x ′ = (−L, 0, 0). The numerical result as a function of X for the unperturbed AdS 5 black hole is shown in fig. 12. In the limit L ≫ 1, the quantity L 0,R is dominated by a configuration which is situated mostly nearby the horizon, that is at ρ g ≈ 1; expanding ρ g (X) = 1 +ρ g (X), one obtains the following geodesic equation to leading order inρ g : which can be solved exactly: Then one can extract an approximation for τ g : Substituting back in eq. (5.18), one obtains that L 0,R scales like 2L. In order to compute the leading order correction in λ 2 to the geodesic length L R , it is not necessary to compute also the leading order correction O(λ 2 ) to the geodesic. It suffices to evaluate the change in length of the unperturbed geodesic; the contribution of the correction to the geodesic will be of higher order, due to the fact that the leading order variation of the action vanishes on solutions of the Euler-Lagrange equation.
We can then expand the renormalized geodesic length as follows: where τ * is the time at the boundary. From eq. (5.21) one can check that there are no further divergences in (5.27) at order λ 2 . It is important to remember that the functions A 2 (τ, ρ) andΣ 2 (τ, ρ) are given by eq. (3.19) in the time window τ * 0 < τ < τ * f ; for time τ < τ * 0 the system is in equilibrium andÃ 2 =Σ 2 = 0 (near τ ≈ τ * 0 there is a short transient period which we ignore).
The geodesic probes a time in the past which is proportional to the length L, see eq. (5.26) and fig. 12. For boundary times τ * such that so that we can use the expressions (3.19) along all the geodesic extension, we can decompose the quantity L R (τ * ) into a constant, a linear and a periodic part in τ * − τ * 0 : where A 2,c , A 2,l , A 2,p , Σ 2,c , Σ 2,p are functions of ρ g . It is convenient to split: The quantity L R,eq corresponds to the geodesic length in the conformal field theory at equilibrium, with an energy density given by the initial energy plus the total work done on the system: The deviations from the equilibrium are parametrized by δL R ; this quantity has a constant part L C and a periodic part L P . The ratio τ d,g = L C /L L can be interpreted as a time delay in the thermalization of the two-point function. In the limit L ≫ 1, using the approximations (5.25,5.26), one finds In the regime L ≫ 1, the delay τ d,g is linear in L: From this result one can deduce that thermalization time becomes longer with the scale L. This kind of behavior in strongly coupled systems has been observed before in [37][38][39][40]. An analogous situation takes place in quenches [11][12][13] where, due to causality, L R is linear in time after the quench up to times of order L/v, where v is the maximal speed of propagating signals. This is true also for other quantities, like the entanglement entropy (see the next section). For realization of the same situation in Vaydia geometry, see e.g. [37,40,41]. The periodic part of the geodesic length L P can be interpreted as a periodic oscillation around the equilibrium value; in the limit L ≫ 1, L P is negligible compared to the equilibrium value: L P /L R,eq → 0. This can be checked by inserting the approximations (5.25, 5.26) into eq. (5.27). This gives where A 2,p , Σ 2,p are evaluated at ρ = 1. This contribution does not scale linearly with L for L ≫ 1; for this reason the present contribution is of the same order as the boundary   effects. This quantity L P was computed numerically as a function of ω T and L, see figs. 13,14,15. It approaches to a constant at large L, with some wiggles obeying δω T × δL ≈ 1. The analytic expression eq. (5.36) gives a qualitative explanation for the wiggles in the figs. 13, 14, 15. In this section we studied the time evolution of the equal-time two-point function of a probe operator with large dimension, see eq. (5.23). The delay in thermalization is proportional to the size L of the region which is probed, see eq. (5.35); this is true in the limit L ≫ 1 in units of T −1 . This result does not depend on the dimension of the driving operator ∆ and on the frequency ω T . The two-point function has also an oscillatory term with frequency ω T ; in the regime L ≫ 1 this term is negligible compared to the equilibrium value. The dependence of the amplitude of this oscillatory term on ω T , ∆ and L is shown in figs. 13, 14, 15.

Entanglement entropy
After discussing geodesics, we skip the two-dimensional probes, i.e. Wilson loops, and we study the time evolution of the entanglement entropy of a spherical region on the boundary. The results will turn out to be similar to those found using geodesics as a probe. A proposal to compute the entanglement entropy in strongly coupled theories with an AdS dual was introduced in [56] and generalized to the time-dependent case in [36]. For the derivation of the time-independent case, see [57,58]. For a review of entanglement entropy in AdS/CFT see [59].
In the gravity dual, the entanglement entropy is computed as follows. Consider first the boundary ∂B of the region of space B whose entanglement entropy we wish to compute and then construct the extremal surface (for AdS 5 it is a 3-dimensional surface) in AdS space which ends on ∂B; the entanglement entropy of B is then given by the volume of this surface. In this section we will choose a region B which is a sphere with radius L. The evolution of the entanglement entropy in the case of a Vaidya metric was studied in [36][37][38][39][40].
Let us first consider the extremal 3-surface parametrized in terms of ρ v (X), τ v (X) and X = X(sin θ cos ϕ, sin θ sin ϕ, cos θ). The minimal volume for λ 2 = 0 is: Translation invariance of τ gives a conserved quantity: For the class of minimal surfaces that we are considering, K = 0; this gives τ ′ v = ρ ′ v /(ρ 4 v −1) as for the geodesics considered in the previous section. We choose the integration constant for τ v in a such a way that τ v (L) = 0. Examples of minimal surfaces are shown in fig. 16. Near the boundary one can use an approximate solution: .
By insertion of this expansion, one can check that V 0 has divergences, which we can regularize by subtracting the divergent part: where ρ cut is some UV cutoff. This quantity is shown in fig. 16. In the limit L → ∞ (where lengths are measured in units of inverse temperature), V 0,R should reduce to the thermal entropy of the region, which is an extensive quantity. Let us check this limit, which is dominated by a minimal surface sitting at ρ v ≈ 1 for most of its extension. We denote by ρ * v the maximum of ρ v ; in the limit that we are considering 1 − ρ * v = ǫ is a very small positive value. Let us introduceρ v = ρ v − 1, then the leading order equation for smallρ is: This differential equation can be solved analytically; the solution gives: where the value of ǫ can be extracted from the requirement that ρ 0 (computed using the approximation (5.42)) vanishes: One can find τ 0 (X), by solving the equation ; it turns out that in the limit L ≫ 1 in units of inverse temperature, it is a good approximation to use: By insertion into eq. (5.37), one can indeed check that the extensive part of V 0,R reduces to the thermal entropy in the L → ∞ limit. As before, it is not necessary to compute the change in the minimal volume as a function of λ 2 ; it suffices to evaluate the O(λ 2 ) action on the unperturbed minimal surface: whereÃ 2 ,Σ 2 are taken as functions of ρ 0 and τ 0 +τ * , where τ * is the time on the boundary.
By insertion of eq. (5.39), one can check that there is an extra divergence for ∆ > 3, proportional to ρ 2(3−∆) cut , (which is a log ρ cut divergence for ∆ = 3): 12(2∆−7) and α ∆ − = 2σ ∆ − . The coefficients (a ∆ , α ∆ , σ ∆ ) are defined in eqs. (3.11,3.27,3.28). For ∆ = 7/2, the functionΣ 2 has also a term proportional to log ρ, so an extra log divergence is expected; for ∆ > 7/2 there is another power-like divergence in addition to the one canceled by V ct,1 ; the following extra counter-term is thus needed: where we have used α ∆ − +1 = ∆−3 3(7−2∆) a ∆ −ȧ ∆ − , and σ ∆ − +1 = 5−∆ 6(2∆−9) a ∆ −ȧ ∆ − . We have to subtract these time-dependent divergent pieces in order to obtain a finite result. In the static case, this kind of divergent contributions to the entanglement entropy were studied in [60]; here we find new divergences which are flagged by their dependence on ω T . These divergent contributions do not depend on the state of theory, but only on the deformation in the Lagrangian.
In a similar way to L R (τ * ) in eq. (5.29), for boundary times τ * such that where −τ v (0) is the maximal time in the past which is reached by the surface, we can decompose the quantity V R (τ * ) such that it has a constant part, a linear part in τ * and a periodic part: where In these expressions A 2,c , A 2,l , A 2,p , Σ 2,c , Σ 2,p are functions of ρ v and the counterterms in eqs. (5.46,5.47) are included.
In analogy to the case of the two-point function, it is convenient to split: The quantity V R,eq corresponds to the equilibrium value of the entanglement entropy in the undeformed conformal field theory at equilibrium, with an energy density given by eq. (5.33). The deviations from the equilibrium are parametrized by δV R . The ratio τ d,e = V C /V L can be interpreted as a time delay in the thermalization of the entanglement entropy. For large L, using the approximations (5.44): The delay τ d,e is linear in L, in the regime L ≫ 1:   As in the case of the two-point functions, the periodic part V P is negligible compared to the equilibrium value in the limit L ≫ 1. This can be checked substituting the approximations for ρ ≈ 1 of eq. (5.44) into eq. (5.52): where A 2,p ,Σ 2,p are evaluated at ρ = 1. This expression scales as L 2 for L ≫ 1, and so it is of the same order as the boundary effects. The quantity V P was computed numerically, see figs. 17, 18, 19. In this section the time evolution of the entanglement entropy was studied. As for the two-point functions, for L ≫ 1, the delay in thermalization is proportional to the size L of the region which is probed, see eq. (5.55). This result does not depend on ∆ and ω T . The entanglement entropy has also an oscillatory term, which is negligible compared to the equilibrium value in the regime L ≫ 1. The dependence of this quantity as a function of ω T , ∆ and L is shown in figs. 17, 18, 19.

Energy fluctuations
The external periodic perturbation can introduce new dynamical features in the behavior of the system; in particular the energy fluctuations may exhibit a phase transition in their behavior as a function of dimension of the driving operator ∆. In [50] the case of a thermally isolated driven system was analyzed. It turns out that the resulting energy fluctuations have a universal distribution different from the one given by the Gibbs measure 3 . Two qualitatively different regimes are predicted, depending on the details of the protocol by which the time-dependent Hamiltonian is employed.
This behavior depends on two universal parameters: 1. the power γ which appears in the entropy S = E γ , several systems have this behavior over a significant range of energies.
2. the power s which appears in the energy dependence of the work made per cycle by the driving force, W = E s .
As shown in section 2, the protocol resulting from a periodic perturbation of a CFT in flat space is such that the work done in each cycle W c is proportional to E s in the ω ≪ T limit, where: The two regimes then depend on η = 2(s − 1)+ γ: for η < 0 the variance of the energy σ E is such that σ 2 E /σ 2 E,eq approaches a constant after a large number of cycles, while, for η > 0, σ 2 E /σ 2 E,eq ∝ E η . The transition between these two regimes occurs when η = 0; for this value σ 2 E /σ 2 E,eq ∝ log E. Close to this transition, a divergent time scale is required in order to reach the asymptotic regime. This is somehow similar to a second-order phase transition, with the diverging time scale being analogous to the divergent correlation length. It would be interesting to study this transition in a strongly-coupled CFT and to explore how the deviation from the Gibbs measure manifests itself in the bulk gravity dual.
In the case of a conformal field theory in flat d-dimensional space, γ = (d − 1)/d; our choice of protocol is determined by the dimension of the operator ∆ and by the frequency ω. In the limit ω ≪ T , s is given by the following expression: The parameter η becomes positive for ∆ > 3(d + 1)/4; i.e. ∆ > 15/4 = 3.75 in four dimensions (corresponding to a scalar with mass m 2 = −15/16 in AdS 5 ) and ∆ > 3 in three dimensions (corresponding to a massless scalar in AdS 4 ). N = 4 SYM on a sphere is a typical example of a system in which S(E) ∝ E γ with different exponents γ for different energy ranges. In the micro canonical ensemble there are four different regimes [61]: free gravitons (γ = 9/10), free strings (γ = 1), small black holes (BH) (γ = 8/7) and large BH (γ = 3/4). In each of these phases, depending on the detailed features of the protocol, it should be possible to achieve both η > 0 and η < 0, and this should correspond to a different kind of quantum gravitational behavior in the bulk. The result eq. (6.1) has been derived in the Poincaré patch, so it is actually applicable only to the case of a large BH. To get a sense of what is possible, we can apply the result of in thermal equilibrium: eq. (6.1) which is valid for γ = 3/4, also to the other three regimes. This yields the critical ∆ for which the transition occurs: 3) The spectrum of N = 4 SYM (which contains operators of dimension ∆ = 2, 3, 4) then allows both for the η > 0 and η < 0 regimes.

Conclusions
In this paper we have used the AdS/CFT correspondence to explore how a strongly-coupled d = 4 CFT in a thermal state reacts to a deformation by a relevant operator with dimension ∆ which is periodic in time; we studied the problem to leading order in the deformation parameter ξ 0 . We computed the amount of energy which is dissipated in the system as a function of ∆ and ω in the linear response regime. The leading order backreaction on the metric is studied; it turns out that, as a function of time, it can be decomposed into a constant, a linear and a periodic part, see eqs. (3.19,3.20).
The entropy can be monitored in terms of the area of the horizon. The entropy increase in a cycle is the same as the one that is found from the equation of state of the undeformed CFT, given the amount of energy dissipated in a cycle. This is true both for the entropy defined in terms of the event and of the apparent horizons, which turn out to be different within short-time behavior.
We have studied the time evolution of the two-point functions of a probe operator with large dimension and of the entanglement entropy of a spherical region. Both the quantities increase in time proportionally to the energy dissipated in the system. There is a delay in achieving the equilibration, which is linear in the size L of the region which is probed. In the limit L >> 1 in units of the inverse temperature, the delay in the thermalization of the two-point function and of the entanglement entropy is respectively: This result does not depend on the dimension ∆ of the driving operator and on the frequency ω T . Both the quantities have also a term which is periodic in time, which is subleading as a function of the size of region L compared to the equilibrium value. The thermalization time becomes longer the further one ventures into the IR domain. This feature was observed in other studies of different out of equilibrium systems, see e.g. [37][38][39][40][41][42][43][44]. All this body of work lends extra credence to the thermodynamical features of black holes.
From the results of [50], we expect a transition in the time evolution of energy fluctuation to occur as a function of the protocol parameters ω T and ∆. In the regime of small ω T , we expect that for ∆ > 15/4 the variance of the energy has a run-away behavior of the type σ 2 E /σ 2 E,eq ∝ E η , where σ E,eq is the Gibbs-like variance and η > 0. It would be interesting to understand the bulk dual of such a run-away behavior; it could be that this is related to Hawking radiation leaking out of the boundary due to coupling to the external source, as in [62,63]. We leave this as a topic for future investigation.
A Appendix: some special values of m 2 In this appendix we will consider separately the value of m 2 which gives integer or halfinteger dimension ∆. Generically we can expand eq. (3.10) in powers of k, replacing φ of Here k is not necessarily integer, it is just a discrete label. The expansion is (note that it is exact, not just leading order): It turns out that the b k are all zero except for integer dimension ∆. For half-integer ∆s there are still some differences with respect to the generic case, concerning the metric backreaction. We will consider each of the five cases separately next.
A.1 m 2 = −4 For m 2 = −4, all the coefficients are determined in terms of a 2 , b 2 : The functionsÃ 2 ,Σ 2 can be expanded as follows: α j ρ 2+j +ᾱ j ρ 2+j log ρ +ᾱ j ρ 2+j log 2 ρ , (A.4) σ j ρ 3+j +σ j ρ 3+j log ρ +σ j ρ 3+j log 2 ρ , where the first of these coefficients are given by: The coefficient α 0 is determined by the differential equation: Holographic renormalization gives [22]: There is an ambiguity in the counterterm in the renormalization procedure; here we use the one which corresponds to δ 1 = 0 in the conventions of [22]. The source and the response function are then parametrized as follows: which gives: For m 2 = −15/4, both a ∆ − = a 3/2 and a ∆ + = a 5/2 are free parameters; the remaining a ∆ + +k are given by eq. (3.12). No logarithms appear in the expansion; some aspects of the source and VEV are still peculiar, so we discuss it as a separate case. The functionsÃ 2 , Σ 2 can be expanded as follows: The first terms are determined as: There are also renormalization ambiguities; here we use the conventions corresponding to δ 1 = δ 3 = 0 and δ 2 = −1/4 in the notation of [22]. The source and the response functions are: a 1 = Re e −iω T τ , a 3 = Re χ 3 (ω T )e −iω T τ . (A.29) Let us first see what we expect to happen if ω ≫ √ 2κ or equivalently if |q| ≪ 1. The effective Hamiltonian after time t ≥ 0 can be calculated by the method of separation of time scales [51], [9] and is to order 1/ω 2 : Hence, we expect the system to quench at time t = 0 with the harmonic oscillator frequency jumping from κ to 1 2 κ |q|, where |q| ≪ 1 for the approximation to hold. Let us check our expectations. The equation of motion for the Hamiltonian H + reads ψ + κ 2 cos(ωt)ψ = 0 , (B.5) which has the solution The Mathieu cosine C and sine S functions depend also on the standard parameters q and a [64]. Here a = 0 and q is given in eq. (B.3). The stability bands for the Mathieu functions are symmetric in q → −q. We have here adhered to the convention of taking C(0) = 1 and S ′ (0) = 1 and we leave q negative. The large-ω limit of eq. (B.6) gives: where we used the fact that (for a = 0 which yields the Mathieu exponent ν = q/ √ 2, see [64], [9]) C(τ ) = 2 − q cos 2τ 2 − q cos ντ + O(q 2 ) S(τ ) = 2 − q cos 2τ ν(2 − q) sin ντ + O(q 2 ) (B.8) ≃ 1 + q sin 2 τ cos ντ + O(q 2 ) , ≃ 1 ν 1 + q sin 2 τ sin ντ + O(q 2 ) . (B.9) This shows that the amplitude to leading order in the 1/ω 2 expansion is given by waves, completely analogously to the situation of a quench [15]. This is also what one would expect by comparing the Hamiltonian H − with the effective Hamiltonian H eff + to leading order in 1/ω 2 . Subleading terms will however reveal the difference between a quench and the suddenly turned-on oscillating mass term.
Let us now calculate the propagator for the system at hand without making any assumptions about ω. We have Ψ 0 |ψ(t 1 )ψ(t 2 )|Ψ 0 = Ψ 0 |ψ 2 (0)|Ψ 0 C ωt 1 2 C ωt 2 2 + Ψ 0 |π 2 (0)|Ψ 0 4 ω 2 S ωt 1 2 S ωt 2 2 + Ψ 0 |{ψ(0), π(0)}|Ψ 0 1 ω C ωt 1 2 S ωt 2 2 + S ωt 1 2 C ωt 2 2 where we used the bracket [ψ(0), π(0)] = i. Since at time t = 0 we have the harmonic oscillator, we find Ψ 0 |ψ 2 (0)|Ψ 0 = 1 2κ , Ψ 0 |π 2 (0)|Ψ 0 = κ 2 , Ψ 0 |{ψ(0), π(0)}|Ψ 0 = 0 , (B.11) from which we obtain the propagator where the sign ǫ = +1 for t 2 > t 1 while ǫ = −1 for t 1 > t 2 . This step implements the time-ordering of the operators in the definition of the propagator. If we expand in |q| ≪ 1, with the characteristic exponent ν = q/ √ 2 + O(q 3 ), we obtain the propagator of the form where A(t) ≡ 1 + q sin 2 ωt 2 . (B.14) In the limit |q| ≪ 1, the function A(t) → 1; this limit reproduces the propagator of a quench in a harmonic oscillator (see e.g. eq. (8) of [15]), in which the frequency is suddenly changed from κ to κ f . The final frequency of the effective quench is: in agreement with the expectation from the effective Hamiltonian, see eq. (B.4). Notice that the q corrections on top of the quench propagator are a) small due to the smallness of q and b) are very fast ω 2 ≫ √ 2κ 2 which is equivalent to |q| ≪ √ 2. Consider now the case of a free field theory. The Hamiltonian density is: The propagator for a given k-vector thus reads ∆ k (t 1 , t 2 ) = 1 2 κ A(t 1 , ν)A(t 2 , ν) cos |k|t 1 cos |k|t 2 + In the limit q → 0, the above propagator (B.21) again coincides with the one of a quenched free field theory, see [15]. For each momentum mode, the initial frequency of the effective quench is κ, while the final frequency is: The order q correction to the propagator has a more complex structure compared to the quantum mechanical case, where it is contained in the prefactor A(t 1 )A(t 2 ) of the propagator (B.13). As discussed in [15], it is possible to define a momentum-dependent effective temper- At large times the system tends to a state with thermal-like correlation functions, having an effective temperature which is a function of the absolute value of the momentum | k|.