Dynamics of Holographic Entanglement Entropy Following a Local Quench

We discuss the behaviour of holographic entanglement entropy following a local quench in 2+1 dimensional strongly coupled CFTs. The entanglement generated by the quench propagates along an emergent light-cone, reminiscent of the Lieb-Robinson light-cone propagation of correlations in non-relativistic systems. We find the speed of propagation is bounded from below by the entanglement tsunami velocity obtained earlier for global quenches in holographic systems, and from above by the speed of light. The former is realized for sufficiently broad quenches, while the latter pertains for well localized quenches. The non-universal behavior in the intermediate regime appears to stem from finite-size effects. We also note that the entanglement entropy of subsystems reverts to the equilibrium value exponentially fast, in contrast to a much slower equilibration seen in certain spin models.


Introduction
In recent years we have seen enormous progress in qualitative and quantitative understanding of out-of-equilibrium quantum dynamics. Theoretical and numerical methods have been very effective to unearth the generic behaviour of a variety of observables in such systems. Coupled with the rapid growth of experimental techniques in cold atom and many-body systems to probe such dynamics, one can furthermore ratify our theoretical understanding. Motivated by these considerations we continue our explorations of dynamics of strongly coupled non-equilibrium quantum systems using holographic methods.
-1 -One simple scenario of interest in many circumstances is a situation where we start with a QFT in global equilibrium and deform it by turning on external sources for relevant operators. The sources provide external dials which can serve to do work on the system and drive it out of equilibrium. We could consider sources that act homogeneously in space (but localized in time), which is often referred to as global quench, or have it act locally in spacetime, which corresponds to a local quench. Both types of protocols are well studied in literature in the past decade or thereabouts. In either case we are considering deformations of the form where O(x) is a (composite) operator of the QFT and J the classical source we dial. The distinction at this level between local and global quenches is simply in the spacetime support of the source J (x). Much of the analytic progress in this front has been in 1+1 dimensional CFTs, where the quench protocols of the form (1.1) can be incorporated into a Euclidean path integral, and studied efficiently by computing correlation functions of the deforming operator O(x) in the unperturbed state of the CFT, cf., [1,2] for the original discussion and [3] for a review.
Our primary interest is in exploring the dynamics of strongly coupled QFTs subject to such protocols in higher dimensions. A natural framework to explore this question is provided by the holographic AdS/CFT duality which maps the QFT problem onto the dynamics of a gravitational system in asymptotically AdS spacetime. For concreteness we will focus on 2+1 CFTs which are originally in global thermal equilibrium and subject them to a quench by a local scalar operator O of dimension ∆. The gravitational problem then comprises of Einstein-Hilbert gravity coupled to a massive scalar, whose mass m is related to the conformal dimension by the standard formula, viz., ∆ = 3 2 + 9 4 + m 2 2 AdS . 1 The initial global equilibrium state maps onto a planar Schwarzschild-AdS 4 black hole and the problem at hand involves analyzing the deformation of this said black hole consequent to turning on a boundary source for the scalar field. This then amounts to a gravitational infall problem. The pulse of scalar on the boundary propagates into the bulk and dissipates through the black hole horizon. Of interest to us are the observables in the interim process.
While there are many quantities that could be, and indeed have been , studied in this context, we will for definiteness focus our attention on entanglement entropy. While strictly not an observable, the entanglement entropy for a particularly chosen spatial region of the QFT captures important aspects of the field theory dynamics. Not only does it provide a measure of how correlations in the system evolve following the quench, but it furthermore is also a simple quantity to compute in the holographic context. The holographic entanglement entropy proposals of [44,45] and their covariant generalization [4] provide an extremely simple route to its computation. All we are required to do is solve a classical problem of finding areas of extremal surfaces anchored on the said region of interest.
In what follows we will explore how holographic entanglement entropy evolves following a local quench. We will restrict our attention to a very specific scenario, wherein we quench a CFT 3 with a ∆ = 2 operator. The disturbance will be taken to be localized in space and time -we pick exponential damping in space and an inverse Pöschl-Teller switch on/off in time, cf., (3.1). We retain translational invariance in one spatial direction, breaking homogeneity in the other. We study entanglement entropy for strip-like spatial regions that are aligned with the symmetry we retain, so that the problem of finding extremal surfaces can be mapped to effectively finding geodesics in an auxiliary three dimensional spacetime. Of interest to us are how the entanglement entropy growth is correlated with the position and size of the strip relative to the quench location.
To appreciate the question, let us recall some well known facts. The classic analysis of [1] of entanglement entropy growth following a global quench in CFT 2 has spurred lots of activity on the subject. While the two dimensional case can effectively be described by a quasiparticle picture, since the entanglement growth is linear due to left and right movers decoupling (following an initial quadratic ramp up [22,28,30]), the holographic models present a much different picture in higher dimensions. 2 The results of various analyses of global quenches have been beautifully encapsulated in the 'entanglement tsunami' picture developed by Liu-Suh in [28,30] and further explored recently in [47]. Following an initial quadratic growth in time, the entanglement entropy for any region grows linearly at a rate dictated by the tsunami velocity v E . To define this quantity unambiguously the authors chose to normalize the local value of entanglement entropy relative to the final thermal entropy expected for the same region once equilibration is complete. This does leave a single parameter which 2 We note here that oftentimes global quenches are holographically modeled by considering a Vaidya-AdS geometry (see [4,46] for early discussions) that corresponds to infalling null matter in the bulk, which does not accord a clean CFT interpretation. A more cleaner perspective is offered by either solving the non-linear dynamics of gravity coupled to realistic matter like a scalar field, or more simply by implementing an end of the world brane boundary state [25] explicitly in holography. The results for the growth of entanglement entropy are however independent of the particularities of the modeling.
-3 -is the aforementioned velocity. It was found not only v E ≤ 1 as required by causality with equality in d = 2 consistent with the CFT 2 analysis, but one could further bound it by a universal dimensional dependent constant v * E (d). 3 This upper bound on velocity was attained holographically for matter that collapsed into a Schwarzschild-AdS d black hole at late times.
Given this rather clear situation for global quenches, we are interested in ascertaining the behaviour when we localize the quench protocol to a finite spatial domain. We in principle could focus on deformations by sources delta-function supported at point. This is natural when studying this problem in QFT as one can map the computation to that of computing correlation functions on some background, however for our purposes of carrying out numerical investigations we choose to smear out the source. We expect firstly that the underlying locality of the QFT forces entanglement entropy to behave causally; as explained in [48,49] this means that the source makes its presence felt only when it acts in the causal past of the entangling surface (the boundary of the region of interest). This is indeed what one sees in explicit computations in CFT 2 . The entanglement entropy only starts changing after a time lag set by the time it takes for the quench disturbance to propagate between the region of interest and its complement. As long as the quench front is localized either in the region or in the complement, we only have the initial state entanglement.
Previous analyses of holographic local quenches by [23] involved modeling the system by the infall of a massive particle -this is effectively an eikonal approximation wherein one is assuming that the wavepackets of the quench are tightly collimated. Moreover, the authors chose to work with very heavy operators ∆ 1 which could then be approximated in terms of worldlines of a small black holes. The relevant geometry can be obtained by applying a suitable symmetry transformation to the global Schwarzschild-AdS black hole and with it in hand properties of holographic entanglement entropy were explored. This picture was further supported by field theory analysis of such deformations at large central charge [50,51]. Our aim to tackle this problem from a different perspective by studying the entanglement evolution in a quenched gravitational background as explained above. We will recover most of the results mentioned above in our analysis.
We can moreover explore quantitative features of the entanglement evolution. We see that the propagation of entanglement is confined to a an effective light-cone. We extract an entanglement velocity v E from this emergent causal structure. Unlike the case of the global quench, the velocity depends on the details of the quench. It appears to grow monotonically with increase in the amplitude of the quench source as well as with the increase of the initial temperature. For a certain range of parameters is appears to track the tsunami velocity bound v * E (3) of [28], while for others it reaches close to the speed of light.
One can understand this behaviour qualitatively as follows: for well localized quenches in large regions A one is in the eikonal limit. Here the growth of the entanglement entropy is linear with an 'entanglement velocity' that is close to the speed of light. On the other hand for regions which are confined within a broad quenching pulse, the situation is nominally similar to a global quench experiment. As we are collapsing scalar matter, we should expect that the behaviour in this domain is isomorphic that seen by [28], and indeed we recover the tsunami velocity. 4 Away from these limiting cases we see contamination from edge effects both from finite size of A and the finite width of the quench source. We have not examined the detailed non-linear effects that cause the velocity to grow from the tsunami bound towards the speed of light, but display some examples which illustrate the pattern.
While our numerical results are constrained to probing small spatial regions relative to thermal scale, 5 we nevertheless are able to extract both this entanglement velocity as well as examine the return to the equilibrium. In contrast to studies in lattice models in low dimensions which display a logarithmic return of entanglement entropy to its equilibrium value after the quench, we find that the holographic systems prefer to equilibrate exponentially.
The outline of the paper is as follows. In §2 we describe the basic set-up for holographic local quenches, describing the general methodology and the determination of entanglement entropy from the gravitational background. In §3 we give the basic numerical results for the quench spacetime and extremal surfaces therein. The key statements regarding the behaviour of entanglement entropy in a locally quenched CFT are then extracted in §4, where we describe the growth velocity v E and the return to equilibrium. We end with some open questions in §5. Some details of the numerical methods are collected in the Appendices. 4 There is a somewhat annoying fact that the tsunami velocity v * E (3) = 0.687 in three spacetime dimensions is mraginally lower than the speed of sound v s = 0.707, making it somewhat hard to convincingly point to precise origin of the effect. We nevertheless feel confident that the analogy with the global quench points to the tsunami velocity being the operative feature. 5 This constraint arises because our numerical solutions only determine the geometry to the exterior of the apparent horizon. For small regions A the extremal surfaces stay in this domain, but for larger regions, they do penetrate the apparent horizon -see [9,53].

Preliminaries: Holographic Local Quench
We are interested in the behaviour of entanglement entropy in a 2 + 1 dimensional field theory that has been driven out of equilibrium locally by an inhomogeneous relevant scalar operator. Holographically, this amounts to solving the gravitational dynamics of a 3 + 1 dimensional asymptotically AdS spacetime and its consequences for the area of extremal surfaces anchored on the boundary.

Metric Ansatz
In order to dynamically evolve a spacetime geometry following a local quench, it is convenient to choose our metric ansatz to be a generalization of the infalling Eddington-Finkelstein coordinates for black holes. We choose to work in an asymptotically AdS 4 spacetime, dual to a 2 + 1 dimensional CFT, where r denotes the radial bulk coordinate, with the boundary lying at r = ∞, and t is a null coordinate that coincides with time on the boundary. We have chosen our quench to be localized in the x-direction and translationally invariant in the y direction. Hence all the fields appearing above {A, χ, F x , Σ, B} depend only on the coordinates {r, t, x} with ∂ y being an isometry. This choice for the metric has many advantages: it provides us with coordinates that remain regular throughout the entire domain as the spacetime equilibrates, it leads to a characteristic formulation of our gravitational infall problem, and it comes with a residual radial diffeomorphism that is of great computational help [54]. Indeed, the metric (2.1) remains invariant under radial shifts, 6 On physical grounds, we anticipate that the black hole's horizon will grow locally as the effects of matter from the boundary are felt in the interior of the bulk. Hence a sensible gauge choice is to dynamically determine λ so that the coordinate location of the black hole's apparent horizon 7 remains fixed. This keeps the calculational domain simple.
Einstein's equations in the presence of a scalar field are given by For simplicity, we restrict our attention to m 2 2 AdS = −2 so that the asymptotic expansion of the scalar field near the boundary is analytic in powers of 1/r: We note that since t is a null coordinate, φ 1 (t, x) will have contributions coming from both the source and the response of the scalar field, as will be explained below.

Asymptotic Geometry
In a theory of gravity on asymptotically AdS spacetimes, asymptotic analysis alone is not sufficient to determine the bulk metric [55]. Indeed, the missing piece in the asymptotic analysis is the boundary stress tensor, determined by solving the full bulk equations: µν is the part of the metric undetermined by the equations of motion for d = 3. While our infalling coordinate chart (2.1) differs from the standard Fefferman-Graham chart typically used for asymptotic expansions, it is a straightforward exercise to carry out an asymptotic analysis. Demanding that the field equations are obeyed in the near-boundary r → ∞ domain we find (2.10) One may also show that the explicit map to the Fefferman-Graham coordinates {τ, ρ, ξ} takes the asymptotic form Additional care needs to be taken when dealing with scalar fields in a theory of gravity formulated in terms of null coordinates. Indeed, the falloff of scalar fields with m 2 2 AdS = −2 is known to behave in Fefferman-Graham coordinates as: as we approach ρ → ∞. By using the coordinate expansion above, we obtain (2.15) thus confirming our earlier claim that φ 1 = φ response + ∂ t φ 0 − λ φ 0 mixes the source and the expectation value of the scalar.

Boundary Stress Tensor
In order to solve Einstein's equations as efficiently as possible, we found it useful to use the boundary stress tensor and its conservation equations to find and propagate the undetermined fields a (3) and f (3) accurately in time (in our scheme, b (3) and c (3) need to be read off from the solutions directly). For asymptotically AdS 4 spacetimes, the boundary stress tensor in the presence of a scalar field of mass squared m 2 2 AdS = −2 can be expressed in the Brown-York form as where we have introduced some boundary data: γ µν is the induced metric on the boundary, K µν , K ≡ γ µν K µν its extrinsic curvatures, and γ R µν , γ R its intrinsic curvatures. Explicitly in terms of the asymptotic expansion coefficients we find that the energy-momentum tensor takes the form while the conservation equations in the presence of the scalar source φ 0 (x, t) read We take our initial state to be in thermal equilibrium, which translates to an initial condition on the bulk metric, which is then the planar static Schwarzschild-AdS 4 black hole spacetime with temperature The initial boundary stress tensor is then simply T µ ν = diag{1, 1 2 , 1 2 }. To model our local quench, we simply need to specify a source function φ 0 (t, x) and let the system evolve according to the Einstein equations, all while making sure that λ is gauge-chosen to fix the location of the apparent horizon.

Holographic Entanglement Entropy
Once we have obtained solutions for the local quench, we can study the subsequent dynamics of the entanglement entropy of a region A on the boundary using the covariant holographic entanglement entropy prescription [4]. The latter requires us to determine extremal surfaces anchored on the entangling surface on the boundary.
For simplicity, we exploit the translational invariance, and restrict our attention to a strip-region The extremal surfaces E A anchored on ∂A are straightforwardly determined by solving a set of ODEs. Using coordinates adapted to the ∂ y isometry, we parameterize the surface by coordinates y, τ . Consequentially, E A is then obtained by solving the geodesic equations in an auxiliary three dimensional spacetime with metricg M N dX Equivalently we solve the Euler-Lagrange equations obtained from the Lagrangian L = g yy g M NẊ MẊ N . While we have phrased the determination of E A as a boundary value problem, it is practical to switch to an initial value formulation. We parameterize the solutions by specifying the turning point, or tip, of the geodesic in the bulk, X M * (τ ) = {t * , r * , x = 0}, and evolve towards the boundary using an ODE solver (for instance the Matlab solver ode45 ) until both ∂A and a specified UV cutoff are reached.
-9 -To this end, we have chosen to transform our system of 3 second order ODEs into a system of 6 first order ODEs in the variables (2.23) With these new variables, 8 The boundary conditions at the turning point are The conditions on P t and P + are a consequence that, because of symmetry, we expecṫ t =ṙ = 0 at X * , whereas the condition for P x has been chosen to normalize the action by setting L = 1. The sign determines whether the geodesic will go towards the positive or negative x-axis.
To translate from the length of the geodesic to the actual entanglement entropy S A we pick an IR regulator L y along the translationally invariant direction and a UV cutoff . We choose to present the results for the regulated entanglement entropy by subtracting off the corresponding answer in the unperturbed theory. There are two natural regularizations we can use: Regulator 1: We subtract the entanglement in the 'instantaneous thermal state' obtained by taking the Schwarzschild-AdS 4 metric with a horizon located at r + (x, t) = M 1 3 + λ(t, x). This choice allows clean matching of the asymptotic coordinate chart. Regulator 2: We alternately can choose to subtract of the vacuum entanglement entropy for the same region, with a dynamical UV cut-off vac (x, t). This gives (2.25) The two regulators differ by a finite amount that is invariant temporally, allowing us to cross-check our numerical results. In what follows we will simply quote ∆S A normalized by L y .

The Quench Spacetime and Extremal Surfaces
We now turn to describing the results of solving Einstein's equations sources by the scalar field boundary condition. We then describe properties of the extremal surfaces of interest in these geometries.

Numerical Solutions
We use the characteristic formulation of Einstein's equations resulting from the null slicing of spacetime outlined in [54] to numerically find the geometry. Even though we start with a complicated set of PDEs, the characteristic formulation simplifies the equations of motion into two categories: the equations for the auxiliary fields that are local in time and reduce to a nested set of radial ODEs, and the equations for dynamical quantities that encode the evolution of the geometry.
To numerically integrate the Einstein and Klein-Gordon equations, we discretize the radial direction using a Chebyshev collocation grid. This choice of discretization for the extra dimension is particularly well suited to find smooth solutions to boundary value problems while ensuring their exponential convergence as the grid size is increased. We opted to choose a rational Chebyshev basis to deal with the non-compact spatial direction. The main advantage of working with a rational Chebyshev grid is that the boundary conditions at x = ±∞ are already implemented behaviourally; as long as the solution decays at least algebraically fast or asymptotes to a constant, we can avoid specifying the boundary conditions explicitly [56]. We use a grid of 41 points in both directions. To propagate in time, we use an explicit fifth-order Runge-Kutta-Fehlberg method with adaptive step size. We also avoid aliasing in both the radial and spatial directions by applying a low-pass filter at each time step that gets rid of the top third of the Fourier modes.
We chose the source function to be φ 0 (t, With it, we can ramp up the scalar field to reach its maximum value α at time t = t q ∆ before it vanishes again. The parameters {s, t q , ∆} are chosen to facilitate the numerics, whereas σ determines the width of the perturbation. In practice, we found s = 0.15, t q = 0.25 and ∆ = 8 to give us satisfying accuracy for the late-time behaviour of the scalar field while preserving a nicely localized shape for the pulse. So we therefore study the quench protocols parametrized by two parameters: an amplitude α and a width σ. Along with the initial temperature of the system which we take to be parametrized by M , we have three parameters at our disposal. The evolution of the spacetime following our quench is fairly simple. The injection of local excitation results in hydrodynamical evolution almost from the very beginning (cf., [5,6] for analogous statements with spatial homogeneity). Since our perturbation excites the sound mode of the system, we have the initial energy-momentum perturbation dispersing at the speed of sound. The presence of shear viscosity results in entropy production, manifested in the solution by the local growth of the horizon area element. In Fig. 1a we display the spatial and temporal profile of the function λ(x, t), related to the area element of the horizon. We see that the initial perturbation indeed results in entropy production, as expected. Curiously, the initial perturbation splits to -12 -two localized perturbations after some time; those follow the expected hydrodynamic evolution. Fig. 1b shows the equivalent evolution of the energy density for the same set of parameters. Finally, Fig. 2 shows that following the conclusion of the quench the total energy is conserved. These features verify the intuitive picture of hydrodynamical evolution following a local excitation of the system.
To quantify the entropy production, we can monitor the growth of the area of the apparent horizon as a function of time. In order to express the result in physical units, we need to convert from the natural time scale on the horizon to the time measured in the boundary. Recall that our solutions for the metric components are obtained on a slice of constant ingoing time coordinate t. We could, following [52], map the horizon data along ingoing null geodesics to the boundary. We will refrain from doing so explicitly and instead work directly in the chosen coordinates leaving implicit this translation. 9 We also overlay the plot for the total energy T 00 dx produced by quenching the system in red for direct comparison.
Using the induced metric h ab on a constant t slice we obtain the area element on the horizon which can be integrated directly. Since the naive answer is infinite, we regulate it by removing the contribution from the initial equilibrium state (i.e. subtract off the static Schwarzchild-AdS answer) to obtain: The numerical results are expressed in Fig. 3, where we also show the total energy for comparison. Notice the striking resemblance of the horizon's area evolution with that of the total energy injected into the system by the quenching scalar field. This seems to indicate that the growth of the horizon is dictated by processes governed by the speed of sound, such as energy and momentum transport. This is indeed the intuition we would have from the hydrodynamic regime of slow variations and it is a reassuring check of the set-up that this indeed is upheld.

Extremal Surfaces
Having the solution at hand we can compute the extremal surfaces as described in §2.4. In Fig. 4 we display the radial depth of the turning point for the extremal surfaces, as function of (boundary) time. Different points correspond to different extremal surfaces, which contribute to entanglement entropy of surfaces of varying lengths. We have plotted the radial depth both in the computational coordinate (in which the horizon is at fixed radial distance) and in coordinates in which the horizon grows.
Since our calculational domain ends at the apparent horizon, we cannot probe extremal surfaces that extend past into the trapped region. These are known to exist in various explicit simulations (cf., [53] for a comprehensive survey in Vaidya-AdS spacetimes). Pragmatically, this restricts our attention to small regions A. We will nevertheless see that despite this restriction we can still extract interesting physical features of S A using surfaces that lie outside the apparent horizon.
One of the interesting features to notice from Fig. 4 is that the geodesics never go beyond their initial depth in the bulk when we consider their position in the ungauged radial coordinate, i.e., where the radial depth is with u * = 1/r * being the radial position of the tip in the coordinate system where the apparent horizon is at a fixed coordinate locus.

Propagation of Entanglement Entropy
Armed with the numerical results for the spacetime geometry and the extremal surfaces therein, we are now in a position to extract some physical lessons for the evolution of -14 -entanglement entropy following a local quench. We restrict our attention to regions A centered around the source of the initial excitation which is taken to be w.l.o.g. at x = 0. We will examine the behaviour of ∆S A as a function of the width L of the strip and time t after the quench. We note that the region of parameter space that we can explore numerically is limited. The amplitude α of the scalar field cannot be too large, otherwise the timeevolution of the quench solution does not converge. Similarly, the evolution code becomes unstable if the spatial discretization falls below a critical grid size, which has for consequence that we cannot resolve quenches with width σ below a certain threshold. The width L of the entangling surface is in turn constrained by the initial values we can pick for M , which determines the position of the event horizon of the initial configuration: if M is taken to be large, then we cannot find extremal surfaces that go deep enough in the bulk to probe larger regions A, whereas if M is taken too small, then it becomes increasingly harder to quench the spacetime with a scalar source. We found that using quenches with width σ = 2, together with M ranging from 0.005 and 0.2 and α between 0.1 and 0.5, yielded interesting results that remained mostly the same, albeit delayed in time, as those with σ chosen larger.
Before proceeding we remind the reader that for regions A which are much wider than the width of the quench source profile, there is a time delay before the entanglement entropy starts to change. This is consistent with the causal properties one would The blue data points represent the radial depth u * in the fixed, gauged coordinate system, whereas the red data points represent the ungauged radial depth U * .
-15 -required of entanglement. Only when the quench can affect both the region and its complement (by being in the past of the entangling surface) would we expect a change in the entanglement for A. This is clearly borne out in our simulations and is used to benchmark that we are on the right track.

An Emergent Light-cone
We first note that the entanglement generated by the local quench is linearly dispersing, i.e., it traces an effective light-cone. This is quite reminiscent of the Lieb-Robinson bound [57] in non-relativistic theories, where correlations follow an effective information light-cone. The speed of entanglement propagation is then denoted by v E below.
The velocity v E we find is bounded from below. A-priori one might guess whether the lower bound is given by the speed of sound, which is the speed in which the initial pulse spreads, thereby further exciting the system and generating additional entanglement on larger scales. The true speed is however a bit lower, as we shall see, suggesting that the mechanism of entanglement propagation differs from that which drives physical transport of energy and other conserved charges in the system. 10 We therefore interpret the velocity v E as the speed in which the initial entanglement, generated locally by the quench, propagates in time. The entanglement velocity can be extracted from the emergent light-cone defined along the curve where ∆S A (t) reaches a maximum for every L in the L − t plane. We remark that unlike the results of [23], the height of this peak does not remain constant in our setup. Instead, we find that the maximum value of S A (t) increases as we increase L.
This behaviour of the entanglement entropy can be quantified rather explicitly. We find that dependence is strongest when the amplitude of the scalar field is varied. For small sizes L, the maximum of S A increases linearly with L. If we denote the slope of these curves by s, then we find the interesting relation 10 A-priori this statement statement appears reasonable, since the propagation of energy in the system is governed by the ability of the system to homogenize, which per se is not the same as becoming quantum entangled. There is thus far no clear mechanism for intuiting entanglement transport in quantum field theories, though the attempts of [47] suggest potentially interesting mechanisms for the same. In the first case, the linear behaviour is shown in Fig. 5. In the second instance (not pictured), while the linear nature breaks down when L is large, the slopes for small to intermediate regions still depend quadratically on the amplitudes. The dependence on temperature is less interesting. When the temperature M changes, the maximum of the entanglement entropy shifts slightly, as can be seen in Fig. 6.
For general values of parameters, the entanglement velocity v E changes with parameters, always bounded from below by the tsunami velocity (4.2), and above by the speed of light. We do however find two universal results which we now turn to.

Universal Behaviour at High Temperature
In the limit of an approximate global quench where the region A is contained within the local quench, i.e., L σ, and at high temperatures, we find a universal light-cone velocity v E = 1 (to very high accuracy), regardless of the amplitude of driving scalar field (including values well within the non-linear regime) 11 . This is depicted in Fig. 7. We note that for some values of parameters, this universal behaviour can be affected by edge effects of the local quench, and is seen for small enough surfaces only. As we decrease the black hole temperature, the velocity at the small surfaces becomes lower than 1. This confirms that v E = 1 is a high temperature effect only.

Wide Quench Profiles
An interesting feature of the emergent light-cone is the abrupt change of velocity as the width of the region A, L, is increased. When the size of the region A becomes of the same order as the width of the local quench, the curve traced by the peak of the entanglement entropy goes from one linear regime to another, as shown in Figs. 8a, 8b, and 8c.
Interestingly, for the first two data sets (for which α = 0.1, M = {0.005, 0.01, 0.02}, σ = 2), the light-cone velocities of v E = {0.678, 0.688, 0.706} are very close to the tsunami velocity of a Schwarzschild-AdS d+1 black hole found in [28], given by We note that temperature does not seem to have an effect on v E , which is consistent with the above formula. For these parameters, the evolution is described by linear -18 -response to good approximation, and in that regime the tsunami velocity seems to capture the spatial propagation of entanglement to very good accuracy.
This behaviour should be anticipated on physical grounds. When the region A is completely immersed in the quench source, we are back to the case where we may approximately think of the situation as a global quench problem. The fact that the source is not homogeneous in A c is irrelevant because all that matters is that the excitations produced by the quench are in the causal past of the entangling surface ∂A. With this in mind we immediately anticipate that the results for the Vaidya quench explore in [28,30] should apply and one see a linear growth with the tsunami velocity.
The story of the local quench however should be a lot richer than the homogeneous global quench. For one, we can encounter an interplay between the size of A and the width of the pulse. We also expect that the non-linearities of gravity will play a role as we try to increase the amplitude. Indeed we see that velocity v E increases as we increase the strength of the non-linearities in the bulk evolution -this is illustrated in Figs. 9a and 9b (where the scalar field amplitude was doubled from 0.1 to 0.2). This goes against the idea of the tsunami velocity as an upper bound on the speed propagation of the entanglement propagation, at least when that evolution is spatially resolved. Coupled with the earlier observation regarding the upper bound on v E ≤ 1, we find it natural to conjecture that The details of deviation from the two extreme limits appear to depend on various effects -20 -which we have not yet disentangled. While the upper bound follows form causality, it is unclear at present whether the tsunami velocity encountered (herein and before) is a fundamental bound on information processing in strongly coupled systems. It would be interesting to come up with a model which allows us to explore the different propagation velocities perhaps along the lines of [47].

Entanglement Decay
The process of return to equilibrium is characterized by universal behaviour and critical exponents. Therefore, an interesting quantity in our model is the decay of the entanglement entropy after it has reached a local maximum. To our knowledge this is the first time this decay has been calculated in either holographic theories or in higher dimensional conformal field theories. From our numerical data we find that the profile for the decay is best fitted by an exponential damping where the parameters a i depend on the specifics of the sources chosen to implement the quench protocol. In Fig. 10 we depict the behaviour for a particular simulation -21 -(parameters in the caption). Note also the time delay in the initial growth, which illustrates the causality feature discussed earlier.
It is interesting to contrast our result for the exponential return to equilibrium against a more slow return seen in some spin chain models. For instance, in [58] the authors study free electrons in a half-filled chain and determined the growth and decay of the entanglement entropy after a local quench. In that set-up they find a very slow return to the unperturbed value. In two dimensions the decay is characterized by S A (t) ∼ a 1 log(t)+a 2 t as t → ∞. The parameters a 1 , a 2 .are again obtained by fitting and depend on the specific details of the quench.
It is somewhat intriguing that the holographic computations relax much faster. This is reminiscent of features of scrambling in black hole physics, which we comment on in our discussion §5.

Conclusions and Future Directions
The main focus of the present paper was to describe the dynamics of the holographic entanglement entropy following a local quench. While this problem has been studied in the past using various known exact solutions to model the quench, we have carried out a full numerical simulation of Einstein's equations in the presence of a perturbing external source on the boundary of AdS. Given the explicit numerical solution to the quench geometry, we can study the dynamics of entanglement entropy by exploring the behaviour of extremal surfaces that are anchored on the boundary.
The upshot of our analysis was a clear signal that entanglement entropy disperses linearly, in a manner reminiscent of the Lieb-Robinson light-cone. The dispersion velocity appears to depend on the details of the quench, though we were able to bound the result between two interesting bounds that have been discussed in the literature earlier. On the one hand we found that for wide quench profile, the propagation speed saturated a putative lower bound, given by the entanglement tsunami velocity obtained by [28] in the context of global quenches (modeled using the Vaidya-AdS spacetime). On the other hand well localized quenches appear to propagate entanglement at the speed of light. It is rather curious that we have results very similar to the Vaidya-AdS quench, for the geometry we construct is not the same. This lends support to thesis of [28,30] that the holographic tsunami velocity ought to be a generic phenomenon.
The second aspect of holographic entanglement entropy which is interesting in our study is the rather rapid reversion of result to the equilibrium value. In various simulations we have tested, the reversion is exponentially fast, in contrast to the much slower logarithmic decay seen in spin models. This suggests again, as has been suspected in the past, that black holes are very efficient at information processing, cf., [59,60].
There are many other interesting areas for further investigation. It would be interesting to study other quench protocols and other theories, including massive models, primarily to extract a more detailed dependence of the entanglement velocity and the rate of equilibration. A particularly interesting direction is the study of (global and local) quenches past critical points, generalizing the results of [61] to higher dimensions. It would also be interesting to study other non-local measures besides the entanglement entropy, which are more sensitive to the spatial structure of entanglement in quantum field theory, and to the differences between strongly coupled holographic CFTs and CFTs of small central charge. In particular, the mutual information of disjoint intervals would be interesting to calculate in our setup for local quenches. Finally, one can make a direct connection to the study of entanglement entropy following a local quench in two-dimensional CFTs, for which we have analytic results to explain the behaviour at large central charge [50]. We hope to report on these results in the near future [62].

A Apparent Horizons
In §2, we used the residual radial reparametrization freedom of the metric (2.2) to fix the coordinate location of the black hole's apparent horizon. Here we provide some details on the process we used.
The notion of apparent horizon depends on the existence of trapped surfaces, which in turn depend on a chosen foliation. Given a spacelike surface S, a trapped surface on S corresponds to the region where both ingoing and outgoing future-directed null geodesic congruences orthogonal to S have non-positive expansions. The apparent horizon is then defined as the boundary of this trapped region, on which the geodesic congruences have vanishing expansions.
In our case, a spacelike surface can be parametrized by the two orthogonal vector fields spanning the x and y direction: We now construct future-directed null geodesic congruences orthogonal to both e x and e y . Ingoing geodesic congruences can be parametrized by the tangent null vector field k M = (0, −1, 0, 0), whereas outgoing geodesic congruences have The normalization is chosen such that g M N k M N N = −1. Since we are interested in the rate of change of the cross-sectional area of null geodesic congruences along their transverse directions, we need to define the transverse metric With this in hand, we can calculate the expansion θ ≡ h M N ∇ M k N . 12 This yields a condition on the dynamics of the field Σ: In addition, taking a time derivative of this relation yields a stationarity condition that ensures that the horizon condition holds for all times. One can show that the resulting constraint can be expressed as a second order spatial ODE that determines the value of A(r, t, x) at the apparent horizon.

B Numerical Details: Integration Strategy and Boundary Conditions
We use the characteristic formulation of Einstein's equations resulting from the null slicing of spacetime outlined in [54] (see also [63]) to numerically integrate our solution.
The clever idea behind this scheme is that both A(r, t, x) and time derivatives disappear completely from the equations of motion if we replace the latter with d + = ∂ t +A ∂ r , the directional derivative along outgoing null geodesics. The equations thereby obtained reduce to a set of nested radial ODEs that is much easier to tame than Einstein's equations in all their glory. For numerical purposes, we need to change variables to u = 1/r to make the domain compact, and redefine the fields appearing in the metric by subtracting the known divergent pieces as u → 0. We do so as follows: We remark that the presence of a scalar source forces us to shift the field Σ by an appropriate function of φ so that it satisfies our asymptotic analysis. Given φ, λ, and b, we can proceed to solve for c, f x ,Σ,Φ, andB -in that order. The prescription to find appropriate boundary conditions for these new fields is to expand the equations of motion near u = 0 and verify their agreement with the asymptotic -25 -analysis conducted in §2.2. We have, after an appropriate normalization of the ODEs: ∂ u c(u, t, x) + c(u, t, x) = g c (t, x) ⇒ c(u, t, x) = C c (t, x) u 3 + g c (t, x), fx (t, x) u + g fx (t, x), −u ∂ uΣ (u, t, x) +Σ(u, t, x) = gΣ(t, x) ⇒Σ(u, t, x) = CΣ(t, x) u + gΣ(t, x), −u ∂ uΦ (u, t, x) +Φ(u, t, x) = gΦ(t, x) ⇒Φ(u, t, x) = CΦ(t, x) u + gΦ(t, x), u 2 ∂ uB (u, t, x) +B(u, t, x) = gB(t, x) ⇒B(u, t, x) = CB(t, x) u 2 + gB(t, x) .
• For c, we find that agrees with the asymptotic expansion for Σ; spectral methods take care of making the non-analytic part vanish: C c (t, x) = 0.
• For f x , we find that g fx (t, x) = O(u 2 ). As a result, we need to specify ∂ u f x (u = 0, t, x) = f (3) (t, x) as a boundary condition, which in turn determines C  which agrees with the asymptotic expansion for d + Σ. However, we do things a bit differently to determineΣ on the entire domain [63]. Indeed, in order to ensure that the computational domain remains on a fixed rectangular grid, a condition (A.4) on d + Σ was derived in Appendix A. We use this horizon condition as a boundary condition forΣ(u = u h , t, x). As a consistency check, one can verify that is indeed satisfied at every time step.
• ForΦ, we find gΦ(t, x) = − forΦ to be in agreement with the asymptotic analysis for d + Φ.
• ForB, things are a bit trickier. The shift above is the only one that reconciles the expansion of d + B near the boundary with the value of gB(t, x) obtained from its redefinition; every other choice leads to a contradiction between the equation of motion and its expected behaviour, which requires With those solutions in hand, the next step is to calculate ∂ t λ. Since λ determines the position of the apparent horizon, it makes sense to solve for its dynamics using information about the horizon for increased accuracy. To proceed, we need two equations: the horizon condition (A.4), which determines a condition on d + Σ, and the stationarity constraint, which ensures that the horizon condition is satisfied at all times. Rather than using the field redefinition of d + Σ in the horizon condition, we use the definition of d + to express d + Σ as As a result of our knowledge of a hor (t, x) from the stationarity condition, we can calculate ∂ t λ, which in turn enables us to solve for a everywhere by using the last relation above. Now that we have solved for all the fields on a particular time slice, all that is left to do is to propagate λ, T 00 , T tx , φ and b forward in time (the last two from our knowledge of d + Φ and d + B), and reiterate the procedure.