Driven Holographic CFTs

We study the dynamical evolution of strongly coupled field theories, initially in thermal equilibrium, under the influence of an external driving force. We model the field theory holographically using classical gravitational dynamics in an asymptotically AdS spacetime. The system is driven by a source for a (composite) scalar operator. We focus on a scenario where the external source is periodic in time and chart out the response of several observables. We find an interesting phase structure in the response as a function of the amplitude of the source and driving frequency. Specifically the system transitions from a dissipation dominated phase, via a dynamical crossover to a highly resonant amplification phase. The diagnostics of these phases include the response of the operator in question, entropy production, energy fluctuations, and the temporal change of entanglement entropy for small subsystems. We comment on evidence for a potential phase transition in the energy fluctuations of the system.


Introduction and Conclusions
The dynamics of quantum field theories driven far from equilibrium is a fascinating topic, owing to the complex interplay of quantum and statistical behaviours in the system. While a quantitative understanding of how field theories respond to nonlinear external sources remains in general an open problem, in recent years one has gained some insight into such phenomena.
On the one hand progress in this direction has been driven by experimental developments which allow for a detailed study. For instance the ability to simulate many-body dynamics in cold-atom systems has led to the opening of a new frontier in dynamical simulations, cf., [1] for a recent review. On the other hand, theoretical horizons have been broadened with the gauge/gravity duality providing an excellent arena to explore the dynamics of strongly interacting many-body systems using (classical) gravitational dynamics in a suitable limit (cf., [2] for a not so recent review). Coupled with the development of excellent numerical algorithms for studying dynamical problems in AdS gravity [3,4], the confluence of ideas and techniques provides an excellent opportunity to further our understanding of out-of-equilibrium dynamics.
A much studied protocol in this context is the quantum quench dynamics, wherein one takes a system initially in equilibrium, typically in the ground state, and subjects it to external sources which change the subsequent dynamics by modifying the Hamiltonian. The rate at which sources act on the system controls the features of the subsequent relaxation, assuming that the sources are non-vanishing for a finite amount of time. The analysis of such a quench protocol has benefited both from theoretical understanding using standard quantum field theory technology in low dimensions [5][6][7] and from a wide array of examples that have been studied holographically in the recent past [8][9][10][11][12][13][14][15][16][17][18][19][20][21]. In most cases the interest is in the approach to equilibrium at late times and the rate at which various observables thermalize [16,[22][23][24][25][26][27][28][29][30][31][32][33][34][35][36][37]. Note that since we inject energy in the process of the quench, even an initially pure state will appear to be well approximated by a thermal ensemble asymptotically (assuming that the field theory dynamics are sufficiently ergodic).
A slightly different but related scenario is one where we subject a system, again initially in an equilibrium configuration, to an external driving source which keeps doing work on it throughout the entire time period under study. More specifically, we will be interested in examining the behaviour when the initial state is chosen to be a thermal density matrix, so that one can simultaneously explore the response of a quantum dissipative system. For non-linear dynamical systems the response under such external driving can provide insight into the dynamics via the coherent build-up of the response.
Classical analogs of what we have in mind are situations where we drive a (damped) pendulum steadily or subject a viscous fluid to external forcing. The latter is particularly apposite, for the problem we study can be thought of as a hot deconfined plasma of a planar gauge theory disturbed by an external source, as studied in the hydrodynamic context in [38]. Rather than letting the driving grow without bound, we will subject our plasma to a periodic driving by turning on the source for a relevant operator. One therefore has two relevant dimensionful parameters characterizing the situation: (a) The amplitude of the external force whose scaling dimension is set by the conformal weight of the operator we exploit and (b) The frequency of the external driving. The third scale which is the temperature of the initial equilibrium state can be factored out, if we are interested in describ-ing the dynamics for a conformally invariant system, which is most natural in the gauge/gravity context. This scenario was explored in [20], who carried out a perturbative analysis for small amplitudes of the driving source. A related analysis of periodically driving a quantum system near a critical point was undertaken in [18].
Gravitationally the problem we study is the following: we have a Schwarzschild-AdS 4 black hole modeling our initial thermal density matrix of a three-dimensional CFT. At some instant of time on the boundary we turn on a periodic source for a relevant scalar operator, which we specifically choose to be of dimension 2 for simplicity. 1 The physics of the system is captured by examining the behaviour or various observables as we vary the amplitude A and the period P of the driving (measured e.g. in units of the initial temperature). We will in particular extend the perturbative analysis of [20] valid for A 1 to the non-perturbative regime A 1 for a wide range of driving frequencies. We find that the system naturally exhibits at least four different phases which are depicted in phase diagram Fig. 1; two of these (labeled IIb and III) are non-perturbative in A.
Before we describe the different phases, let us examine for a moment the physics of the gravitational system qualitatively. Initially we have a planar black hole in AdS 4 . When we turn on the scalar source, we are injecting energy into the bulk. This energy does work on the system and simultaneously heats it up. The latter is seen by the fact that some of the energy falls behind the horizon, which grows 2this is the gravitational response to the disturbance of the plasma. However, in this process we also induce an expectation value for the operator whose source we tweak. When we disturb the system 'slowly enough', the operative parameter measuring this being the product of the amplitude and the period, the system has time to catch-up. This is the dissipation dominated regime indicated by I in Fig. 1. In this regime the injected energy falls behind the horizon with little fanfare.
As we ramp up the disturbance, the plasma is driven more and more non-linear, with a dynamical cross-over visible as we move into phases IIa or IIb of Fig. 1. Note that the entire non-linear dynamics in the system is induced by the non-linearities of gravity, for we model the system simply by a free (massive) scalar field. In this phase the response gets more in-phase with the source. It is amusing to contrast this with non-linear scalar dynamics; we find that in this phase we can model the scalar 1PI effective action induced from the gravitational interactions to be well mimicked by a polynomial potential (see [17] for previous studies of self-interacting scalars in AdS). 1 This choice turns out to have several advantages as the dual scalar being conformally coupled to gravity in the bulk allows a certain level of technical simplification in various holographic renormalizations we need to do to extract physical data.
2 As we will be describing the dynamics of Einstein-scalar system with the scalar field satisfying the null energy condition, the areas of the event and apparent horizon (in the canonical foliation) have to grow monotonically -a consequence of the area theorem [39] (see [40] for an excellent overview). We will elaborate on this point in §2.3.
The "phase diagram" of the driven holographic plasma characterized by the period (P ) and amplitude (A) of the driving force, measured in units of the initial thermal scale T 0 . There are four distinct regimes marked on the diagram which are explained in the main text. σ in refers here to the imaginary (or in-phase) part of the conductivity defined in Eq. (3.3). As we move from southwest to northeast in the figure, the system is driven into a more non-linear regime; the crossing of the grey-dashed boundary is the turn on of the in-phase part of the conductivity σ in in regime II, and the crossing of the blue-dashed boundary signifies the entrance into the resonance phase of regime III i.e., |φ max 1 | → ∞. The character of the different regimes is further illustrated by displaying the phase portrait of the scalar operator (expectation value against source) used to drive the plasma.
In this regime there is less dissipation; the entropy production by the growth of the horizon area is slowed down relative to region I. The primary distinction between the two phases IIa and IIb themselves is the lag in the response seen as the period is increased (hence the tilt in the phase portrait).
For even larger disturbances, we enter region III, where the system response gets highly resonant and there is a steep growth in the response. As one might suspect this is the domain where the gravitational non-linearities are strongest and indeed one can check that such behaviour is not visible for a polynomially (self-) interacting scalar. In the course of our investigation we explore not just the phase portrait, but various other physical quantities of interest, such as the growth of entropy and dissipation in the system, the rate at which entanglement is produced, etc.. For instance, region IIb is characterized by enormous fluctuations in the energy of the system over a single period and continuous but non-differentiable behaviour in the entanglement entropy of a sub-system.
Let us contrast our results with the analysis in the perturbative regime of small amplitudes undertaken in [20]. 3 As one can see from phase diagram Fig. 1 for small amplitudes, A 1, one is largely in the dissipation dominated linear response regime. This is indeed consistent with the analysis of [20], who explore the dependence of observables on both the period of the driving as well as the dimension of the perturbing operator ∆. As for us the latter remains frozen and we are unable to check the detailed scaling relations they find, but in the common domain of overlap we do indeed have agreement. In particular, for perturbing operators of dimension ∆ = 2 in CFT 3 we expect to see that the energy dissipation as a function of the period scales as E diss ∼ P −1 (for A 1), independent of the initial temperature. Furthermore, we also expect that the work done in each cycle, measured by the entropy produced, to scale with the increased energy density. We find that in the slow driving regime this scaling closely tracks the prediction from local thermal equilibrium, but starts to grow more steeply as we transit into more interesting non-linear regimes.
While the various response functions provide us with a useful diagnostic of the phase structure of the dynamical evolution, we also attempt to gain insight into the non-equilibrium dynamics using entanglement entropy for small sub-systems. This non-local probe exhibits distinct characteristic features in the various regimes: for weak driving, the growth of entanglement is gradual (and appears to track the growth of thermal entropy), while for strong driving there are steep oscillations and glitches in its evolution. We should caution the reader that we have only examined entanglement entropy for relatively small sub-systems, owing to technical complications with numerical stability. Nevertheless these results suggest a rather rich structure in the temporal growth of entanglement with driving, which deserves further detailed exploration [41].
The outline of this paper is as follows. We begin in §2 by giving a quick overview of the basic set-up and the numerical solutions. Following this in §3, we set out the various observables we use to explore the behaviour of the system. In particular, we justify the rationale behind phase diagram Fig. 1 and how we should physically think of the different regimes. §4 is devoted to the study of entanglement entropy in this system where we focus on the region of an infinite strip and exploit the underlying homogeneity of the set-up. We conclude with a discussion in §5. Some technical results about holographic renormalization required for computing various observables is collected in the Appendices; Appendix A collects some useful information about holographic renormalization in our models while Appendix B provides details relevant for computing entanglement entropy.

Driven CFTs and their Holographic Duals
We first take the opportunity to set up the basic problem of a field theory driven out of equilibrium by turning on a source for a relevant operator. We then go on to describe how to model this in the holographic set-up and present the basic methodology and results from the numerical simulations.

Driving CFTs by Relevant Operators
We are interested in the dynamics of strongly coupled plasmas that are driven by an external source. The initial plasma is in equilibrium in some homogeneous thermal state at a temperature T 0 for t < 0. At t = 0 we introduce external sources with some specified spatial-temporal profile that we control. We focus exclusively on situations where the external sources are spatially homogeneous, but otherwise arbitrary and tunable at will.
To wit, the system under consideration can be modeled by an equilibrium density matrix, evolved under a time-dependent Hamiltonian, i.e., we take where O(x) is a single trace (gauge-invariant) relevant operator of conformal dimension ∆ < d. The source J (x) is chosen to have no spatial dependence and be temporally periodic and thus can be represented as Here Θ(t) is the Heaviside step function for turning on the periodic perturbation of amplitude A and driving frequency ω = 2π/P at t = 0; later in actual (numerical) implementations we will choose a suitable ramp factor to smoothly turn the perturbation on.
In the presence of the source, the Ward identities following from diffeomorphism and Weyl invariance get modified. A simple analysis shows that the boundary conservation equation now has an explicit source term indicative of the work done by the driving source on the CFT. Likewise the one-point function of the trace of the energy-momentum tensor no longer vanishes but satisfies Since the boundary theory is conformal, it does not have any intrinsic time scale. The time scales in the problem come from only the driving force, namely its amplitude and period. The situation of interest is thus characterized by three scales: • T 0 : the initial thermal scale for the homogeneous plasma.
• A: the amplitude of the source whose scaling dimension is d − ∆.
• ω: the driving frequency or the time-scale set by the period P = 2π/ω.

Holographic Driving
The gravity dual to this set-up is modeled by the dynamics of a scalar field φ with mass m 2 φ = −2, dual to a relevant perturbation of the boundary theory.
In our holographic implementation of this set-up we will work in d = 3 and consider a scalar operator with conformal dimension ∆ = 2. While this is rather specific, we will explore the phase structure of the driven system as a function of the ratio of scales outlined above. The qualitative features we believe are independent of these actual choices. 4 We have included a dimensionless gravity-scalar coupling α g which we can use to tune the amount of backreaction on the geometry; for the most part we will focus on α g = 0 or α g = 1, to model probe and interacting scalar fields respectively. We want to study gravitational dynamics driven by a scalar field whose nonnormalizable mode is turned out as dictated by the source J (x), i.e., take φ 0 (t) = A cos(ωt) and study the behaviour of the theory with varying amplitude A and frequency ω. The gravitational background is an asymptotically AdS 4 spacetime, which we write in ingoing Eddington-Finkelstein coordinates (sometimes called the Bondi-Sachs form) as: The coordinate dependences of the metric functions f , χ, ρ are explicitly indicated with homogeneity ensuring that ∂ x and ∂ y are Killing vector fields. Our initial state is a planar Schwarzschild-AdS 4 black hole with temperature T 0 = 3/π, corresponding to horizon size r + = 1. This bulk solution is given by (2.7) 4 We have also set AdS = 1 for simplicity.
For our choices of m 2 φ = −2 in d = 3, the amplitude A has mass dimension 1. Thus we have two interesting time scales associated with the external driving force: the period P and the inverse amplitude A −1 . To capture universal physics, we look at relatively late times of the non-thermalized system compared to both of these scales. Note also that in those late times the initial value of the temperature, T 0 , becomes irrelevant.
There has been much interest recently in holographic quenches, in which the system is initially driven to an excited state, and then is allowed to return to equilibrium, a process which exhibits some degree of universality. In contrast, we are interested in the dynamics of the steady state system while it is being driven. Hence, in our solutions we do not turn off the driving force at late times, and seek universal features associated with the driven steady state system. We will see that such dynamical features exist, and they strongly depend on the parameter ξ(P, A) ≡ P A , (2.8) the unique dimensionless parameter formed from the two time scales associated with the driving force. Below we refer to the regime ξ 1 as the weak driving regime, and ξ 1 as the strong driving regime (which is further divided into two separate dynamical regimes). We also measure time in units of the period P , thus we vary and discuss the dependence of observables on the two dimensionless parameters: the strength of the drive and time.

Bulk solutions
We solve the equations of motion resulting from the scalar-gravity Lagrangian (2.5) by direct numerical integration. The boundary conditions on the scalar are prescribed by the source and the metric is required to be asymptotically AdS 4 . The AdS boundary is attained as r → ∞ and the asymptotic behaviour of the fields is More specifically, we use the characteristic formulation of the resulting partial differential equations as explained in detail in [4] to numerically integrate for the solution.
The advantage of the method is that it allows us to use constrained evolution: at each time step we solve a nested set of ODEs to determine the time derivatives of all dynamical quantities, and then we use one of the standard time evolution schemes to march forward in time. While we follow the general logic of [4], in our implementation we found that some of elements described in [42] enabled for a more robust evolution.
To solve the radial ODEs we discretize the equations using a Chebyshev basis in the radial direction, typically taking a grid of 60 points. For time evolution we use an explicit Runge-Kutta method of order 4, with an adaptive step size. We filter at each time step by throwing out the top third of the Fourier modes for each dynamical variable to avoid artificial and unphysical growth in amplitudes of short wavelength modes associated with the UV cutoff.
In the regime of strong driving, we found it necessary to turn on the perturbation gradually from zero. Therefore we include a ramp-up time of 2 P , after which the amplitude reaches its intended value. Thus, the first few periods of each solution show behaviour sensitive to details of the ramp-up protocol. We look at observables only after this ramp-up time of 2 P .
In Fig. 2 we show one example of evolved bulk fields for a specific solution. As we perturb the system by a relevant operator, the scalar field grows towards the horizon. All fields are (at least approximately) modulated with the period of the source.
At this point it is worthwhile mentioning one important consistency check on the numerical scheme, which relies on the existence of a smooth horizon in the spacetime. Given a metric and a Cauchy slice in the bulk spacetime, one can find the outermost trapped surface on this slice. If we have a set of Cauchy slices that foliate the spacetime, then the future outermost trapping horizon, which we simply refer to as the apparent horizon by a common abuse of terminology, is typically defined by taking the union of the outermost trapped surface on all the slices. The apparent horizon thus defined is subject to an area law which was originally discussed by [39] we refer the reader to [40] for a concise modern summary and proof of the statement. It is however important to note that the statement relies on the existence of a sensible foliation of the spacetime by Cauchy slices. Indeed, it is possible as discussed in [43] to find exotic symmetry breaking foliations (which are however incomplete) in which even the Schwarzschild black hole solutions fails to have a trapped surface.
We mention this in passing, as [20] quotes the result of [43] to argue that apparent horizon areas need not be monotone generically. They however do not encounter such behaviour, for with the choice of ingoing coordinates in (2.6), there is a canonical choice of bulk Cauchy slices respecting the homogeneity of the disturbance. In this foliation the result quoted in [40] does apply and in fact simply follows from properties of null congruences using Raychaudhuri's equation. 5 Our results are indeed consistent with this expectation and we have checked that the area of the apparent horizon does grow monotonically in t (which labels the leaves of the foliation chosen), as we shall extensively see in the sequel. While initial results of [21] appeared to suggest otherwise, upon closer scrutiny, one finds that in numerical analyses so far the area of the apparent horizon does respect the second law as derived by [39]. 6

Driving Diagnostics
Having constructed the holographic duals we now turn to lessons that can be extracted from the geometry for the dynamics of strongly coupled field theories. Apriori there are a number of observables which are useful probes of the out-ofequilibrium situation and we will focus on those that offer most clear insight into the dynamics. Our primary goal is to quantify the behaviour of the system as a function of {P, A} and construct a phase diagram demarcating the various regimes in this phase space. Let us quickly enumerate the observables we will use and proceed to explain why they give us some insight into the dynamics: • The phase portrait of response φ 1 (t) as a function of the source φ 0 (t). Alternatively, this relation can be codified in a conductivity σ(t), as defined below in (3.3). We find 4 underlying phases regions that the system can fall into.
• The φ 1 -φ 0 phase portrait features for polynomial and non-polynomial potentials with the gravity-scalar coupling α g switched on and off.
• The cycle-averaged thermodynamics quantified by the energy density avg (t) and entropy density s avg (t), and the scaling relation s avg ∼ γ avg between them.
• The work done in each cycle, measured as the difference in average energy between two successive cycles, cycle = (n+1) avg − (n) avg . We typically take n to correspond to the penultimate cycle of our simulation.
• Fluctuations fluc (t) in the energy density around avg (t) and the maximal response |φ max 1 (t)|.
• Entanglement entropy and extremal surface evolution for fixed spatial strips A on the boundary.
When the system is driven by an external source, the most basic quantity is the response, which is characterized by the scalar one-point function in the presence of the source. In linear response theory, this can be obtained from the retarded Green's function of the operator O(x) evaluated in equilibrium. We are not just interested in the linear response regime, which would correspond in our set-up to A T 0 , but in the full non-linear response. To visualize the response of the strongly coupled plasma, especially in the non-linear regime, where its phase relative to the source is important, we will find it instructive to exhibit the phase portrait, the trajectory traced by the system in the φ 0 -φ 1 plane. We also codify the relation between scalar source and response by a complex conductivity, defined below.
In addition to the one-point function of the operator deforming the CFT, we are interested in the boundary energy-momentum tensor. This can be decomposed in to an energy density (t) and a pressure. In the holographic set-up one has The scale Ward identity (2.4) implies that pressure is not an independent observable since it can be obtained from knowledge of (t) and φ 1 (t), so we will not discuss the pressure separately. Additionally, to probe the local thermodynamics we will monitor the local entropy density s(t), obtained by computing the area of the apparent horizon at time t. 7 The dynamics of the bulk gravitational fields encode the heat production resulting from supplying external energy to the system. We monitor the explicit time dependence of the energy density (t) and the entropy density s(t) along with their values averaged over each driving cycle period P , and find for the most part that the averaged values are increasing with time. 8 These provide a useful diagnostic of the 7 Using the area of the apparent horizon (defined as the outermost trapped surface in the foliation respecting spatial homogeneity) results a causal boundary observable. One maps points on the apparent horizon to boundary points by Lie transport along radially ingoing null geodesics, which in the ansatz (2.6) are simply lines of constant {t, x, y}. On the other hand the teleological nature of the event horizon implies that its area would not provide a good measure for the boundary entropy density, cf., [3,44] for a discussion of this point. 8 Note that the averaging makes avg (t) and s avg (t) discrete in time.
departure from equilibrium, as one can monitor the scaling relation to infer the local thermodynamic equation of state. We define the thermodynamic scaling exponent γ when the system is in a steady state t > t s via Note that in thermal equilibrium, conformal invariance predicts γ 0 = 2 3 . We will encounter this and other scaling regimes in our driven system when conformal invariance is broken.
Note that one natural set of non-local observables we could use are the multipoint correlation function for gauge invariant local operators, perhaps for O itself. However, realistically this computation involves solving the wave-equation for the linearized scalar fluctuations on top of the background we have constructed, together with the imposition of suitable boundary conditions on the future horizon, to obtain sensible time-ordered correlation functions. These boundary conditions are somewhat tricky to implement (see however [45,46]) -we will therefore postpone a discussion of correlators to the future. 9 Below we describe the behaviour of the observables mentioned above in three distinct dynamical regimes, and comment on the bulk interpretation of those regimes. Once we have gained sufficient intuition from this exercise, we will then examine the entanglement entropy for a specified boundary region.

Dissipation Dominated Regime
The simplest situation occurs in the regime of weak driving ξ 1, which is best described as the dissipation-dominated regime (phase I). This includes the regime of small amplitudes, studied perturbatively in [20]. In this weak driving regime, the behaviour of all observables is dominated by dissipation, which we now demonstrate by looking at some specific observables.
As we drive the system by the scalar non-normalizable mode φ 0 it is instructive to divide the scalar response φ 1 to the part in-phase with the driving force, and the part completely out-of-phase with the perturbation. In analogy with an electromagnetic perturbation in linear response, we can complexify the time dependence of the scalar field 10 and define a complex conductivity 3) 9 We could following standard practice attempt to compute two-point correlation functions using the geodesic approximation [47]. However, as discussed in [48] and more recently in [49], this prescription doesn't generically reproduce correct time-ordered correlation functions (we really want in-in correlation functions in our set-up). As a result we will also refrain from computing geodesics in the numerical background. 10 That is, regard cos ωt and sin ωt as the real and imaginary parts of e iωt . φ 0 φ 1 Fig. 3: The phase portrait of the dimensionless responseφ 1 ≡ P A φ 1 versus the dimensionless sourceφ 0 ≡ 1 A φ 0 for ξ = 0.001 1 in the dissipation dominated regime (P = 0.001, A = 1) which we label as phase I. We evolve the solution for 10 periods with each colour segment representing one period. The early times t < 2 P show the effect of the perturbation ramp-up, and thus are numerical artefacts that we omit from the plot.
With this notation the out-of-phase and in-phase parts of the response correspond to the real and imaginary parts of the complex conductivity, σ out (t) and σ in (t), respectively. This is the usual convention for the more familiar conductivity, related to electromagnetic perturbations. As shown in Fig. 3 in the low driving regime the scalar response is precisely out of phase with the scalar source, σ in = 0, meaning all the energy is dissipated and none of it used to excite the internal energy associated with the scalar field i.e., no work is being done on the system. This is the quench limit and it matches with what we expect from the behaviour of the perturbation in linear response. The complex conductivity σ = σ out is purely real and has constant amplitude as a function of time at high frequencies. 11 This is manifested in the final steady state being reached almost immediately and consisting of closed untilted trajectories in phase space. As we shall see below, tilting of the trajectories in phase space is indicative of non-trivial response and work done onto the system. Fig. 4 shows what fraction of the complex conductivity σ out is present on each point on the (P, A) phase diagram, and for what we are concerned with currently, the system has the response being completely out-of-phase with the source when the period is low.
Both the energy and entropy density, averaged over each cycle, grow linearly with time in the dissipation-dominated regime . As the black hole grows, its entropy growth tracks its energy growth at a slightly higher rate than the equilibrium relation s avg ∼ and avg (t) versus time. their own evolution with time. Note that the expansion of the black hole horizon is not necessarily adiabatic (as measured e.g., by the rate of entropy increase 1 TṠ S ). In the low amplitude regime, one can also estimate in perturbation theory the amount of energy dissipated per cycle cycle which we define as the difference of the log 10 ω α Fig. 6: The dimensionless scaling parameter α(ω) from fitting cycle ∼ ω α for a small amplitude A = 1 in the linear response regime. It is expected for our choice of the scalar and dimension (∆ = 2 and d = 3) that α → 1 in both the small ( cycle ∼ ω) and large frequency ( cycle ∼ ω 2∆−d ) limits.
average energy avg between two successive cycles; for simplicity we take the result for the last two cycles of our evolution in quoting the results below. One expects the relation to take a scaling form cycle ∼ ω α . The scaling exponent α should be a nontrivial function of frequency itself; for low frequencies it is independent of the driving operator, but the high frequency limit cares about the spectral properties about the operator in question. Specifically, one finds that [20]: cycle ∼ ω for small frequencies and cycle ∼ ω 2 ∆−d for high frequencies. Since we are not scanning over different choices of the driving operator, we have a single shot at determining this result. As depicted in Fig. 6 we indeed find that the energy dissipated is linear both at low and high frequencies: α(ω) → 1 both for ω 1 and for ω 1 (a coincidence owing to our choice ∆ = 2 and d = 3). Interestingly there is some non-trivial intermediate frequency behaviour which appears to amplify the energy dissipated in a single cycle.
The bulk picture of the process is also very simple: as we send energy pulses, which are either weak or infrequent, they interact very rarely before falling into the black hole horizon. All injected energy from the boundary goes towards steadily increasing the black hole mass and the scalar field remains unexcited. The more diverse behaviour observed below can be attributed to gravitational interactions of those energy pulses before they fall into the black hole.

Dynamical Crossover Tilted Regime
We now discuss the qualitative changes in the system as we begin to move from the weak driving ξ 1 to the strong driving regime ξ 1 (from regime I to regime II through the grey-dashed line in phase diagram Fig. 1). Fig. 7 depicts a typical phase portrait of the system as we cross into the new dynamical regime. We see that this regime is characterized by an onset of excitations of the scalar field and breaking of discrete time translation symmetry. The left panel of Fig. 7 shows the transitioñ   In the right panel of Fig. 7 we see the effect of moving into the new dynamical regime at low amplitudes: there is a clear tilt in the phase portrait from the one in Fig. 3 with ξ 1 which indicates that the response is no longer completely out of phase with the source. The tilting of the trajectories at lower frequencies corresponds to the emergence of a finite in-phase contribution σ in > 0 in the conductivity; this sets the system somewhere between one with a purely out-of-phase conductivity (closed circular trajectories) and one with a purely in-phase conductivity (straight diagonal line trajectories). In other words not all of the injected energy is dissipated as was the case in regime I, but rather, work is actually being done on the system. As a result of having less dissipation in this regime, the energy and entropy of the black hole grow more slowly with time. Moreover, we find the scaling behaviour between the average energy and entropy, with a thermodynamic scaling exponent γ > 2 3 , for all values of (P, A), as shown in Fig. 8. In other words, while the work done in the system slows down the energy increase of the black hole, the entropy production is affected less. To understand this regime further, it is instructive to reproduce this type of phase portrait for a system without gravity. To that effect, we can study the special case of scalar field evolution in a fixed black hole background, with no backreaction on the geometry (i.e., α g = 0). To include non-linearity into the problem, we add self-coupling to the scalar field, to mimic the effect of the non-linearities due to gravitational interactions (see also [17]). Fig. 9 depicts the phase portrait of a selfcoupled scalar field with two types of polynomial potentials, which we took to be our original form (free massive scalar) and also one with quartic self-interactions: We can see that without non-linearity as in Fig. 9a, the phase portrait is tilted, but sharp features of the phase portrait are lost compared to the case with the same driving but also gravitational backreaction, depicted in Fig. 7b. Adding a polynomial non-linearity, as done in Fig. 9b, gives a phase portrait that starts to form slightly sharper features along with some amplification of the response. Thus, the simple system of self-interacting scalar field allows us sufficiently separate the two effects in regime II: we see that the tilt in the phase diagram is associated with decreased frequency, whereas the breaking of time-translation invariance is associated with increased amplitude. We note also that for this simple system, the third dynamical Fig. 9: The phase portrait of the dimensionless responseφ 1 ≡ P A φ 1 versus the dimensionless sourceφ 0 ≡ 1 A φ 0 for ξ(P = 10, A = 1) = 10 with α g = 0 and different polynomial potentials V (φ). The conventions are as described in Fig. 3.
regime of unbounded amplification discussed in the next subsection seems to be absent.
Thus, the bulk interpretation of this dynamical regime becomes clear: the pulses of energy injected at the boundary interact gravitationally before falling into the black hole. This results in additional physics to that of simple dissipation, modeled here by infalling the black hole. The gravitational interaction is due to perturbative exchange of gravitons, and can be mimicked by a polynomial self-interaction of the scalar field. In the next subsection we will see the effect of the gravitational interactions becoming strong when both A and P are large.

Unbounded Amplification Regime
As we increase the driving strength further in both A and P directions (from regime II to regime III through the blue-dashed line in phase diagram Fig. 1), we enter a dynamical regime no longer reproducible by polynomial self-interactions of the scalar field. We see the phase portrait of the scalar field in Fig. 10 for two instances of parameters in this regime. Moreover, we find this dynamical regime to be characterized by unbounded response and restoration of time translation symmetry.
As we increase the strength of the driving force ξ, the phase portrait becomes sharper and tilted, corresponding to an increased response and, again, less lag with the source as seen in Fig. 4. The 'slowness' of the energy injection from the boundary allows the scalar field to heat up as if the entire process were adiabatic, consequently allowing the scalar response to respond relatively quicker to the source. Note that although Fig. 4 shows |σ in /σ| ≈ 1 in this regime, the absolute value |σ| is actually very large in this unbounded amplification regime so that a small |σ out /σ| is still strong enough to keep the black hole perpetually growing in size.
The maximal response |φ max 1 | over our ten cycles of driving is plotted in Fig. 11 throughout the phase diagram. It is seen to increase rapidly with ξ past the dissipationdominated regime. This seems to indicate the presence of a non-linear resonance, which allows the scalar response to grow without bound. An interesting feature of Fig. 11 is that the maximal response does not grow in the high frequency regime regardless of how large ξ is by increasing A. It seems unlikely that unbounded behaviour is attainable even for amplitudes drastically higher than the bounds of numerical explorations reported in Fig. 11. Physically, this means that a rapid pulsing of small packets of energies can barely amplify the response of the system; the frequency of driving has to be below a certain bound for resonance to be possible -or in other words, a certain slowness in the sourcing is required. We conjecture that one should would see unbounded amplification only in the combined large P , large A regime which is slightly different from the traditional definition of resonance that depends only on frequency. An interesting curiousity is a slight dip in the response for moderate values of ξ preceding the rapid growth. This trough appears to demarcate the domains of bounded (regime II) and unbounded responses (regime III) empirically. It would be interesting to come up with a explanation for this phenomenon.
Finally, it is amusing to model the non-linear effects of gravity in terms of an effective scalar potential to see what is necessary to attain regime III. We find that while a scalar field with polynomial self-interaction does not seem to posses this log 10 P A This choice of scalar self-interaction is chosen to agree with our previous example (3.4) in the small field regime, but of course behaves differently for large field values. In Fig. 12 we see that indeed similar features of the phase diagram are reproduced: narrow closed trajectories and resonant response. We conclude therefore that the features of this dynamical regime are due to strong, non-perturbative gravitational effects occurring outside the black hole horizon. The fact that the non-linearities induced by gravity can be extremely strong, should perhaps be borne in mind while attempting to come up with simplified models of gravitational dynamics in AdS spacetime.

Energy Fluctuations
Another observable we monitor is the behaviour of energy fluctuations. More precisely, we consider the deviations from the average energy in a each cycle, fluc (t) = | (t) − avg (t)|. These cycle fluctuations are a crude proxy for genuine fluctuation information that can be extracted, for instance, by considering symmetrized twopoint functions of the boundary energy momentum tensor. Such ensemble-averaged fluctuations are known to exhibit phase transitions in periodically driven systems [50]. Some indication those transitions are possible in holographic systems is given in [20].
The results for our simulations in various regimes are plotted in Fig. 13. We observe a qualitative change in these cycle fluctuations between different regimes. While in the dissipation-dominated phase we do not see a lot of deviation from the mean, there is a steep growth in fluctuations as we enter the non-linear phases. The fluctuations are maximal in the unbounded amplification regime (regime III). We note that in contrast to the maximal scalar response, which also grows dramatically in that phase, the fluctuations do track the driving frequency, with there being more deviations in the large period limit.
It would be useful to confirm this behaviour directly with the computation of correlation functions, a task we leave for future investigation.

Entanglement entropy
Thus far we have discussed various local observables (response functions and thermodynamic data) which have served to help us chart the phase diagram of the driven system in Fig. 1. We now turn to other non-local field theory observables that are sensitive to the non-equilibrium dynamics. Since we are not going to examine the behaviour of higher point correlation functions, we will dive right into the dynamical behaviour of entanglement entropy.
In the boundary we have a density matrix ρ(t) which is time-evolving with respect to the perturbed Hamiltonian. At any given instant of (boundary) time, we pick a spatial region A and construct the matrix elements of the reduced density matrix ρ A (t) = Tr A c (ρ(t)) by tracing out the degrees of freedom in the complement (on the chosen Cauchy slice). The entanglement entropy is given by the von Neumann entropy of ρ A , i.e., S A (t) = −Tr A (ρ A log ρ A ) which we can monitor as a function of time.
Holographically computing the entanglement entropy for boundary regions in time dependent situations involves finding bulk codimension-2 extremal surfaces E A anchored on the said boundary region A [23]. We study the evolution of entanglement entropy focusing in particular on translationally invariant strip regions: (4.1) The bulk codimension-2 surface ends at x = ±a at some chosen instance of boundary time t A and is obtained by solving effectively a set of geodesic-like equations with our interpolated metric functions Σ, f , and χ (see Appendix B.1 for details). The covariant holographic entanglement entropy prescription [23] generalizing [51,52] states that Should there be multiple extremal surfaces, we choose the one with minimal area (homologous to A). The proper area of these surfaces diverges owing to the locality of the underlying QFT. In our case we encounter potential divergences not only from the surface reaching out to the asymptotic boundary, but also from the presence of the sources driving the system. The physical result we are after is the finite universal contribution S fin A , which will measure the entanglement created/destroyed as we drive the system away from thermal equilibrium. Fortuitously, for our choice of scalar operator, there are no contributions due to the source, and hence we can simply regulate by background subtraction. 12 As a result we will consider as our entanglement diagnostic, the following finite quantity where L y is the IR regulator in the non-compact translationally invariant direction.
Since we drive the system away from thermal equilibrium, S A (t = 0) is the corresponding value of the entanglement entropy computed in the Schwarzschild-AdS 4 geometry. In what follows we will simply quote the results of our numerical simulations both for the behaviour of the extremal surfaces themselves and ∆S A (t).

Extremal surfaces in the driven geometries
The extent to which the extremal surfaces penetrate into the bulk can for the most part be determined from the location of the cap-off point which we parameterize as (t * , u * = 1/r * , x = 0). 13 For very small regions we are reasonably close to the AdS boundary whence, the curves are approximately semi-circles u 2 + x 2 ≈ a 2 . As we increase to larger strip widths the extremal surfaces start to probe the interesting regions of the driven geometry and thus allows us to see qualitative differences between the four phases.
Generically we see that the following statements hold irrespective of the phases we consider: 1. The radial depth and the temporal extent spanned by the surface evolves nontrivially as a function of t A . One consequence of working with ingoing coordinates (2.6) is that the surfaces naturally dip back in time (see [33,53]).
2. The oscillatory driving of the system imprints itself in the profile of the extremal surfaces, with the scale of these oscillations set by the the driving parameters A and P . The periodic movement of the surface can be seen in pulsations of the turnaround point of the surface: u * and t * have oscillations of the same period superposed over some enveloping function.
horizon is at u + = 1. In these coordinates the proper size of the region A increases (due to ρ(t, r)) which means that the surfaces want to get closer to the horizon to extremize the area functional. The rate at which this happens depends on both the amplitude and the frequency of the driving. We also note that surfaces dip less temporally, i.e., t * − t A is increasing.
4. We also note that the location of the extremal surface appears to be consistent with causality of entanglement entropy [49]. While we have not explicitly checked that the surface lies in the casual shadow of the boundary region A, one simple consistency check visible from our results for t * is that t * < t A − a.
We remind the reader that in (2.6) lines of constant t and x are radially ingoing null geodesics. Causality at the very least requires that the cap-off point of the extremal surface lies below the ingoing null geodesic from the domain of dependence. Since for the strip region the boundary domain of dependence is a diamond anchored at (t A ± a, 0) and (0, ±a), we note that the ingoing light ray from the bottom tip of this diamond cannot signal to the cap-off point.
In the following discussion we will illustrate the behaviour of the extremal surfaces more explicitly in each of our phases. We have been reasonably conservative in our analysis and have chosen to work only with surfaces that do not get too close to the horizon (in fact u * < 0.2). This is to avoid both numerical issues as well as to avoid complications from the existence of multiple extremal surfaces. We follow a single branch of solutions as described at the end of Appendix B.1. The primary results of the extremal surfaces are shown in the plots Figs. 14, 15, 16, and 17, where we show the evolution of the extremal surface as well as u * (t A ) and t * (t A ).
Linear regime (small A): Although all phases display extremal surfaces that sink into the bulk with each driving cycle, the growth of u * in the linear regime of small amplitudes is most steady. We focus here on phases I (high frequency; dissipationdominated) illustrated in Fig. 14 and IIa (low frequency; tilted) illustrated in Fig. 15, which fall under this characterization. As the frequency is lowered and we pass from the dissipation-dominated phase to the tilted phase, there is drastic reduction in the growth of u * per cycle.
The evolution of t * in the two phases is also interesting; t * − t A is gradually increasing on average with time (recall that in the stationary geometry t * − t A would be constant). It turns out to be useful to look at a dimensionless parametert * ≡ (t * −t A )/P which measures the cap-off time relative to the boundary. In this context, there is more time-lag in phase I i.e.,t * I t * IIa 0, which hints at the cause for why the surfaces do not penetrate as far deep in the bulk in phase IIa as opposed to phase I. 14 In addition we see strong oscillatory patterns in phase II in spite of having only a steady increase in u * ; such a feature is absent in phase I.

Non-linear regime (large A):
We now turn to the phases III (low frequency; unbounded amplification) illustrated in Fig. 16 and IIb (high frequency; wobbly) illustrated in Fig. 17 in the non-linear regime of high amplitude. Some of the features seen in the linear regime continue to pertain: we see more pronounced oscillations iñ t * and a decreased tendency for the surfaces to lag behind in time at lower frequencies.
In the unbounded amplification regime (phase III), we see significant bursts of growth of the extremal surfaces. The oscillatory driving is felt rather acutely by the surfaces and the evolution is considerably violent. On average however, u * appears to advance more serenely despite having large amplitude oscillations per cycle.
In the dynamical crossover wobbly regime (phase IIb), there is a considerable amount of instability. We chose here to work with smaller strip widths a = 0.01 (instead of a = 0.05) to avoid complications of phase transitions between multiple competing extremal surfaces. The early part of the evolution is in line with what happens in the dissipation-dominated regime (phase I), but shortly after, there are discontinuities in thet * parameter with no noticeable effect in u * . Around t A /P ≈ 4.0 − 4.2 and t A /P ≈ 4.6, we see an exchange of dominance in the extremal surface, Fig. 15: Evolution of the extremal surfaces for a strip of width a = 0.05 with driving parameters ξ(P = 10, A = 1) = 10 (phase II; tilted). Conventions described in Fig. 14 apply. which starts out at a higher value oft * .
All in all, the extremal surfaces in the non-linear regime definitely has elements of intrigue owing to the large pulses of energy that affect the bulk geometry significantly. Although we do not delve into extremal surfaces that are positioned deeper into the bulk, we notice in the course of our analysis that the surfaces tend towards the horizon as expected. More curiously, we also find that for larger regions we cannot find extremal surfaces that stay outside the apparent horizon. This is not surprising since we expect based on earlier results that there will be surfaces that penetrate the apparent horizon of the black hole (cf., [24]). However, one of the disadvantages of our numerical scheme is that we are unable to explore this interesting regime due to the fact that the spacetime inside the apparent horizon has been excised. As explained in [4], this was to avoid complications with having caustics in the coordinate chart. Analysis of entanglement entropy however does require us to have the complete bulk geometry.

The evolution of entanglement
We now turn to the evolution of the entanglement entropy; the results are presented in Fig 18 for the regulated quantity ∆S A as introduced in (4.3).
In the dissipation-dominated regime (phase I), the entanglement entropy gradually increases, though in each cycle of forcing there is a time period for which the growth is negligible. We expect this feature is simply a consequence of the entanglement entropy tracking the thermal entropy. Even though we are not quite probing the full thermal contribution with the relatively small regions A, it bears to reason that the variation of the geometry is more or less equitable on all radial scales. This appears consistent with other probes of this phase. As we discussed in §3.1 the weak driving allows the system to efficiently dissipate the energy induced by the source and the conductivity σ(t) was purely imaginary. Basically the dominant effect here is the growth of the black hole horizon due to the driving and this in turn imprints itself into the growth of ∆S A seen in Fig. 18a.
On the other hand when we reach phase IIa (tilted regime) by way of small amplitudes, we start to see definite oscillatory evolution of ∆S A . In each oscillatory period we see a local reduction in ∆S A . On the other hand the temporal radial depth attained by the extremal surface as measured by u * is almost similar to that in phase I by juxtaposing the behaviour in Fig. 14 and Fig. 15. In phase IIa however, our extremal surfaces are closer to the boundary in contrast to phase I. We conjecture that the origin of the reduction in the ∆S A is associated with the sharp oscillations in t * or equivalentlyt * . These imprint themselves into the actual value of the area despite the surface not getting too far into the bulk (which is possible since even the asymptotics of the geometry is sensitive to the driving, cf., (A.17)). The onset of non-monotone growth of ∆S A in Fig. 18b characterizes the departure from the linear regime to the non-linear domain in line with the behavior of the phase portrait which in turns modifies the conductivity (which picks up a real part σ in > 0 in phase IIa).
The temporal change of ∆S A is much more pronounced in the non-linear regime. In the unbounded amplification phase III (see Fig. 18c) and the dynamical crossover wobbly phase IIb (see Fig. 18d), the ∆S A appears to track the time-coordinate of the cap-off pointt * quite efficiently. Indeed here we expect the non-linearities of the system to be the dominant effect. We know that the black hole grows quite rapidly in response to the energy injected into the system at the boundary from our discussion in §3.2 and §3.3. The behaviour in phase III is smooth with large amplitude oscillations, which qualitatively track quite well the behaviour oft * . The dynamical crossover wobbly phase (phase IIb) exhibits a lot more drastic behaviour. We encounter for the first time a jumps in the family of extremal surface that minimize the area (satisfying the boundary conditions and the homology constraint). These jumps translate into continuous but non-differentiable kinks in ∆S A visible in Fig. 18d. We again note that the radial position of the cap-off point of the extremal surface behaves much more smoothly and the glitches appear int * . Furthermore, the growth of the entanglement itself is rather steep as we see about an order of magnitude difference in ∆S A between the low amplitude and high amplitude regimes. It is interesting to contrast the change of entanglement entropy with the change in the thermal entropy to see how the two are correlated. As we have argued above, the fact that we have an ever increasing thermal entropy (the bulk black hole is constantly growing) implies that even for small sub-systems we will quickly see overwhelming thermal contribution. We display in Fig. 19 the functional dependence of ∆S A on the (normalized) instantaneous thermal entropy s(t)/s(t = 0). It is immediately apparently by eyeballing the plots that there appears to be nearperfect correlation in three phases with Fig. 19c corresponding to phase III being the only outlier. To get a quantitative feeling for the correlation we have also indicated the Pearson correlation coefficient ρ p as well as the Spearman rank coefficient ρ s . These are statistical markers for measuring correlations between two sets of data and are defined to take values in the interval [−1, 1]. The values ρ s , ρ p = 0, ±1 signify zero, perfect positive and perfect negative correlation respectively. While the Spearman coefficient indicates that the observables in question are monotonically related, the Pearson coefficient provides an accurate measure of linear correlation. Indeed from the results quoted in Fig. 19 we see that ∆S A (s) is a linear function to a very good approximation in phases I, IIa and IIb. It is curious that the linearity is respected even in the presence of the glitches in the growth of entanglement entropy (we do not see any drastic behaviour in the area of the apparent horizon). The unbounded amplification phase III clearly demonstrates the effects of non-linearities by decorrelating ∆S A and s(t).

Discussion
The non-equilibrium dynamics of strongly coupled field theories is amenable to detailed quantitative exploration using the AdS/CFT correspondence. We have exploited this set-up to study the behaviour when a homogeneous thermal plasma is driven away from equilibrium by a periodically sourcing a relevant (composite) scalar operator. The resulting dynamics exhibits a rather rich phase structure illustrated in Fig. 1.
We identified four distinct phases, characterizing them in terms of the frequency and amplitude of the external driving force. Of these the dissipation dominated phase I is perhaps most intuitive for here the weakness of the driving, allows the system to to catch up with the driving. This is clearly visible in the various observables we studied; the complex conductivity of the response is purely real owing to the phase lag between the source and response and the evolution of entanglement is pretty quiescent.
There is more structure when we ramp up either the period of driving, or the amplitude, for now the system departs quite rapidly away from equilibrium. The response therefore is more pronounced; we see more in phase response and greater temporal oscillations. In phases IIa to IIb there emerges a non-vanishing imaginary part to the conductivity, which in fact appears to capture the entire response for high values of the period and amplitude. We also notice that there are significant fluctuations in the energy density and the entanglement entropy and furthermore, the entropy density grows rather rapidly in this regime. Perhaps most intriguing is the unbounded amplification of phase III, where we see sharp fluctuations and a highly non-linear response. We argue that this response appears to be not captured by polynomial self-interactions of the composite operator; the intricate dynamics of gravity in AdS appears to induce effective non-polynomial couplings in the effective action for the operator O we use to perturb the system away from equilibrium. We believe this fact is significant and should be taken into account when attempting to construct effective models distilling the effects of gravitational interactions for strongly coupled systems .
While our focus has been on computing the simplest set of observables, essentially one-point functions and entanglement entropy for small sub-systems, the power of holography is that we can do much more. In time independent equilibrium scenarios it is straightforward to use the holographic map to compute correlation functions (at least two point functions). In the genuine non-equilibrium scenarios as those we have focused on the technology for computing such observables, whilst present [45,45] is still a bit cumbersome to work with (at least numerically). It would be interesting to develop these techniques further perhaps taking inspiration from the analytical models of [54,55]. This would allow us with a direct probe of fluctuations in the plasma, which can be contrasted with the dissipation in the system, the latter being measured by the entropy production through the growth of the horizon.
Likewise our exploration of the behaviour of entanglement entropy has been restricted to analysis of small sub-systems for pragmatic reasons. While the sub-system under consideration was chosen to have fixed size, the fact that we are continuously driving the system leads to an ever increasing thermal contribution to the entanglement. Geometrically this is easy to understand since the horizon for our bulk solution is ever growing (as we have indicated that both the event and apparent horizons are required to be monotonic in our set-up) and reaches out towards the boundary in the course of the evolution. As a result, the local thermal scale can overwhelm the relative smallness of the sub-region we choose. To have precise mapping of the entanglement structure we need to be able to ascertain the true minimum of the area functional in such scenarios bearing in mind that the extremal surface can (and often does) penetrate various horizons. A significant obstacle in ascertaining this is the fact that the characteristic method for solving Einstein's equations developed in [4] excises the region of the spacetime behind the apparent horizon. While this is a technical obstacle, overcoming it would not only enable us to probe the interior of a highly non-equilibrium black hole using holographic entanglement, but it could also allow us to explore other interesting scenarios such as the effect of perturbing the ground state of the system by external sources.

A Holographic Renormalization
We collect here some salient results for the computation of physical field theory quantities using standard holographic techniques.

A.1 Scalar deformations
The bulk action (2.5) should be supplemented by boundary counter-terms to ensure that (a) the bulk equations of motion follow from a consistent variational principle and (b) the on-shell action evaluated on the solutions is finite.
In standard Poincaré-AdS d+1 the scalar field behaves asymptotically as We will work with standard quantization (Dirichlet boundary conditions) for the scalar field, which involves treating the mode that fall-off as r ∆−d as the source for the scalar field.
In the presence of the source we let the metric to take the FG form, where g µν (z, x) = γ µν + O(z). If necessary we will denote by γ the induced metric on the surface z = z which differs from the boundary metric by a conformal transformation by z 2 . We will ignore this issue for most part and write the counter-terms in terms of γ µν below for simplicity.
With these conventions we find the following boundary counter-terms: We are using conventional AdS/CFT definitions: Our interest concerns conformally coupled scalar field which has a mass in AdS units given by To compute the boundary energy momentum tensor we vary where we should take care to include the appropriate radial dependence in the definition of γ µν . Lets split the contribution from the graviton and the scalar and write where the split is determined by the requirement that T µν φ ∝ φ. Then the two pieces can be computed efficiently as follows: which one can show evaluates to a nice covariant expression: where z is the location of the cut-off surface. The scalar contribution can be evaluated by using the fact that we are interested in the boundary variations to obtain: where · · · indicate the contribution from the higher order counter-terms and we have put back the powers of z now. The details now depend on the asymptotic expansion of φ. For general ∆ we have to worry about the fact that the Taylor series solution in the neighbourhood of z 0 looks like and we need to know the various intermediate pieces to complete the analysis. The case we are interested in is rather special, where there are no powers of z in the Taylor expansion between the source and the vev, so let us simply record the result for this case for now leaving a more general analysis for later. Before proceeding though, let us note that we can express (A.11) covariantly as follows (r = z −1 ): where n A is the unit normal perpendicular to the cut-off surface.

A.2 Specializing to
In this case the asymptotic expansion belongs to the special kind where where we are allowed to use the fact that φ 0 z ∆ − is the beginning of an independent Taylor series where the powers of z change by 2 units (use the fact that the Lagrangian has φ → −φ symmetry). This corresponds to the case we are interested where ∆ + = 2, ∆ − = 1 in d = 3.
In this circumstance we can simply use the terms written explicitly in (A.11) to obtain Then we find Now, we can get the final answer for the case of interest either by working with the Fefferman-Graham expansion in which case we need to know that The metric fall-offs allow us to compute the pieces in the boundary stress tensor directly since the z 3 term above is the correct answer. Using this or directly computing from the CY-ansatz (A.13) we claim to obtain (rescaled the result by a factor of 3/2).
We can check that this satisfies the Ward identities: where we used the boundary conservation law derived from the solutionȧ 3 = − 1 2 φ 0φ1 .

B Extremal surfaces and entanglement for strips
In this appendix we describe our methodology for finding extremal surfaces relevant for the computation of entanglement entropy. For simplicity we will focus on regions which exploit the symmetry of our set-up and consider A to be a strip extended along one of the translationally invariant directions, say y without any loss of generality, as in Eq. (4.1). We need a bulk codimension-2 surface that ends on the boundary of this region i.e., at x = ±a (at the chosen instant of boundary time t A ). We describe our strategy for finding this surface and computing its (regulated) area below.

B.1 Determining extremal surfaces
To find the extremal surface, we start by gauge fixing the reparameterization invariance on the surface. We take y to be one of the coordinates. Dimensionally reducing in this direction, we construct an effective action for a curve in the remaining directions and pick a proper-length parameter λ as the second coordinate. Thus, the extremal surface E A is embedded in the bulk as We choose the proper-length parameter to ensure that √ detγ ab = 1, which implies that the unregulated area of the extremal surface is given as in terms of parameter distance λ E A spanned by the curve and the IR regulator L y . In practical terms we work with the effective Lagrangian where the metric functions ρ, f , χ are obtained by interpolation of our numerical solutions. This is a geodesic problem, with some non-minimal coupling from the dimensional reduction along the translationally invariant direction of the strip. Instead if using the geodesic equations, we found it convenient to pass to a set of six first-order Hamilton-like equations by introducing P t = t , P x = x , and P + = r −ρ t which are related to the conjugate momenta. The equations we solve are the above three and We start from x = 0 in the bulk at some smooth cap-off point (x = 0, t * , r * ) where t = r = 0. 15 and propagate out to the boundary. We evolve until a with a fixed UV cut-off at r A and regulate the final answer for the entanglement entropy by background subtraction (see below).
In the main text we illustrate the temporal dependence of the extremal surfaces and S reg A for each of the four phases (I-IV) of Fig. 1 for fixed strip width a. Since we numerically control the data of the cap-off point we work iteratively: we start by fixing a suitable strip width a by tuning r * and t * , then we evolve the extremal surfaces by increasing t * and re-adjusting r * such that the strip width remains as a. We note that we assume that there are no discontinuities or multi-valuedness in the map from (r * , t * ) → (a, r A ), which we believe makes sense for small strip widths. 16 Finally, to work in a compact domain we choose u = 1/r ∈ [0, 1] which we will use to explain the properties of the extremal surfaces.

B.2 Regulated entanglement entropies
Since the extremal surfaces reach out all the way to the boundary, the proper area is divergent with the coefficient of the leading divergent term fixed by the area of the entangling surface ∂A. For a state of the CFT with vanishing sources for operators it is well known [52] that the entanglement entropy behaves as where S fin A is finite in the limit u A → 0. In normalizable states of the field theory S fin A is the universal contribution to entanglement which should be independent of the 15 This cap-off point is not necessarily the deepest point in the bulk; for the examples shown in this paper it however does turn out to coincide. 16 Such behaviour was noticed in extremal surface computation in global Vaidya-AdS by [53].
cutoff value u A . 17 One natural way for us to extract this quantity is to measure the entanglement relative to the t = 0 thermal Schwarzschild state ∆S A (t) = S A (t) − S A (t = 0), which can be extracted simply by vacuum subtraction. Usually, when we turn on sources for relevant operators, these can contribute additional divergences to the entanglement entropy [60]. In general in the presence of additional relevant scales one naively expects there to be logarithmically divergent terms polluting (B.5) and rendering vacuum subtraction meaningless. Fortuitously, this does no happen for the problem at hand. This can be extracted by examining the detailed discussion of [60], which we paraphrase below.
There is however a quick argument for the absence of logarithmic terms which we now describe. For scalar operators in CFT d with operator dimension d 2 < ∆ < d 2 + 1, as we have considered, it is well known in AdS/CFT that the corresponding bulk field has mass in the window where both asymptotic fall-offs are normalizable, i.e., m 2 ∈ (m 2 BF , m 2 BF + 1) with the Breitenlohner-Freedman bound mass m 2 BF = − d 2 4 as usual. 18 In this window note that ∆ − −∆ + < 2 and we have the Legendre transformed theory with an operator of dimension ∆ − by switching to alternate quantization [61].
Turning on a source for the faster-fall off mode ∆ + is equivalent, insofar as the leading back-reaction on the metric, to considering instead a state in the Legendre transformed theory where the alternate quantized operator with dimension ∆ − acquires a vacuum expectation value. However, since the divergence structure of entanglement is the same in all states of the field theory, and the conformal vacua of the two theories (standard and alternate quantization) coincide, it follows that the divergence structure of S A should be unchanged from (B.5), even with J (x) = 0. Our story is of course a special case with ∆ + = 2, ∆ − = 1 in d = 3. This observation is consistent with the results of [60] and the counter-terms used in [20].
To explicitly analyze the structure of the divergences in the entanglement entropy, let us consider the metric given in (A.17). Since the details of the divergences are blind to the boundary spatio-temporal behaviour of the sources we will examine the somewhat simplified setting where φ 0 = const to glean the relevant information.
With the time-translational symmetry restored by this choice, the Lagrangian for the extremal surface (which now is minimal) is simpler: 17 For the vacuum state of a CFT 3 with A being a circular disc S fin A would give the F-function [56,57] (the latter defined as the logarithm of the partition function of the theory a three-sphere). In fact, this can be used to define a UV finite quantity without recourse to background subtraction: following [58,59] we can just as well consider R d dR − 1 S A , with R being the disc radius, as the measure of entanglement growth. 18 Implicit in this statement is the fact that we are quantizing the scalar field with standard (Dirichlet) boundary conditions, so that the dimension of the dual operator is ∆ = ∆ + .
where g ii (z) is the spatial component of the metric in (A.17). This system has a conserved Hamiltonian, which we can exploit to write down an expression for the area directly. Introducing, z * which captures depth to which the minimal surface penetrates into the bulk, we have for the on-shell value of the area Using the explicit form of g ii , the second term is at least z 4 near the boundary so we can forget about it. The first term is all that matters, so lets look at √ g ii which has the z −1 divergence expected upon integration, but no further contribution of relevance in z → 0 limit. From the φ 2 0 term we get a contribution to the finite part of the entanglement, and this is indeed the physically relevant answer. It should be clear from this discussion is not specific to the choice m 2 = −2 in d = 3, but should hold for d 2 < ∆ + < d 2 + 1 as we argued abstractly above.