Nonlocal probes of thermalization in holographic quenches with spectral methods

We describe the application of pseudo-spectral methods to problems of holographic thermal quenches of relevant couplings in strongly coupled gauge theories. We focus on quenches of a fermionic mass term in a strongly coupled N=4 supersymmetric Yang-Mills plasma, and the subsequent equilibration of the system. From the dual gravitational perspective, we study the gravitational collapse of a massive scalar field in asymptotically anti-de Sitter geometry with a prescribed boundary condition for its non-normalizable mode. Access to the full background geometry of the gravitational collapse allows for the study of nonlocal probes of the thermalization process. We discuss the evolution of the apparent and the event horizons, the two-point correlation functions of operators of large conformal dimensions, and the evolution of the entanglement entropy of the system. We compare the thermalization process from the viewpoint of local (the one-point) correlation functions and these nonlocal probes, finding that the thermalization time as measured by the probes is length dependent, and approaches the thermalization time of the one-point function for longer probes. We further discuss how the different energy scales of the problem contribute to its thermalization.


Introduction
Quantum quenches are processes where an isolated system is driven to a far-fromequilibrium state by rapidly varying some control parameters. It has been possible to produce and study such processes in laboratory experiments in recent years, in particular, with ultra-cold atomic gases [1,2]. This experimental progress has provided a great impetus to improve our theoretical description of quenched systems. Certainly theoretical progress made with investigations within a variety of different frameworks, e.g., two-dimensional conformal field theories [3], (nearly) free field theories [4,5] and integrable models [4,6], as well as some results applying to weakly interacting relativistic quantum field theories in higher dimensions [7]. It remains a challenge to find broadly applicable and efficient techniques, as well as extracting insights into general organizing principles for the behaviour of far-from-equilibrium systems.
A new theoretical tool allows for the investigation of quenches for (certain) strongly coupled field theories is gauge/gravity duality [8]. Assuming the robustness of this holographic duality in non-equilibrium situations, as studied in e.g., [9,10], it is possible to study the behaviour of the boundary field theory, either when it is perturbed, or far from equilibrium. Initial work by [11,12] has lead to a large body of work in the field of quantum quenches of field theories at strong coupling, including [13][14][15][16][17][18][19][20]. In the gravity dual of the quantum field theory, the quench usually has a simple geometric interpretation and is introduced, e.g., in the form of a gravitational shock wave collapsing into a black hole and a collapsing shell of matter described by the Vaidya metric, collapsing into a black hole [10, [21][22][23]. Applications of holographic quenches include quenches across phase transition points in the field theory [24], and to model hadron collisions in particle accelerators such as RHIC [25].
In an ongoing research program including [14], [15] and [16], we study the response of a strongly coupled N = 4 supersymmetric Yang-Mills thermal plasma, quenched by a relevant operator, using the holographic duality. Having previously studied such quenches, we now apply more powerful numerical techniques to find the full timedependent profiles of the perturbations of the metric and scalar field in the dual AdS spacetime. This allows us to utilize nonlocal probes such as two-point functions and entanglement entropy to better understand thermalization at various distance scales.
The quench that we study here for the Yang-Mills plasma in four spacetime dimensions, is that of switching on a fermionic operator in a smooth manner, by giving it a time-dependent mass. This is dual to a radially collapsing scalar field in an AdS black brane geometry in five dimensions. Since the mass is turned on in a homogeneous manner in the boundary theory, we are studying a global quench. Further, as implied in our description of the gravitational dual, we are studying a thermal quench where the theory begins in a thermal state, rather than in the vacuum -the latter simplifies the analysis, as we will describe below.
It is straightforward to find the solution for a static scalar field on an AdS background containing a planar black hole -see [26] for the solution in N = 2 * theory, and [27] for general operators satisfying the constraint 2 ≤ ∆ < 4. However, once the scalar field is given a time-dependent source, the nonlinear Einstein and Klein-Gordon equations become highly nontrivial to solve. Treating the scalar field as a perturbation backreacting on the metric only at second order, the problem becomes more tractable, since the scalar field and metric components decouple at leading order in the Klein-Gordon equation. The solution to the metric then becomes that of the static background, plus a time-dependent contribution which is second order in the amplitude of the scalar. Despite this simplification, in the asymptotic series for the scalar field, one (time-dependent) coefficient remains undetermined, and can only be solved by evolving the scalar field forward in time from a known initial configuration.
In [14] and [15], a finite difference method was employed, which is computationally quite costly. These studies were limited to first order in the amplitude of the scalar, meaning that the normalizable coefficient in the asymptotic expansion of the scalar field was calculated, as well as some terms in the metric's asymptotic expansion which could be directly calculated from this normalizable coefficient. However, the full second order profile of the metric could not be determined in this way, making the calculation of nonlocal probes in the geometry impossible.
Chebyshev spectral methods are powerful methods for solving systems of differential equations [28]. Representing the solution to the equations by a series of Chebyshev polynomials, we can approximate the full radial profile of the solutions to a high degree of accuracy. In the present paper we will apply these methods to the problem of solving for a massive scalar field in a five-dimensional AdS spacetime, as well as the time-dependent profiles of the metric perturbations.
The organization of the rest of our paper is as follows: we first introduce the physical setup of the scalar field on the AdS-black brane spacetime in sections 2 and 3, as well as the coordinate system used. We then go on in section 4 to show the calculation of the thermalization of the system, by studying different nonlocal quantities that can be calculated in this spacetime. First we examine the evolution of the apparent and event horizons. We then go on to calculate the two-point correlation functions.
Finally we calculate the entanglement entropy of a strip on the boundary. These calculations require knowing the full time-dependent geometry, i.e., the full profile of the second order metric components, and therefore rely on our numerical simulations of the evolution. We find in particular that the apparent and event horizons thermalize much sooner than the local one-point functions of the quenching operator, and that wider separations in the two-point correlator and entanglement entropy thermalize later than for narrower regions. Furthermore, two-point functions and entanglement entropies that are wide enough will thermalize later than the one-point function as well.
In section 5 we investigate the thermalization behaviour of the previous section in closer detail. The entanglement entropy and two-point functions are dual to minimal surfaces and geodesics extending into the geometry of the spacetime, respectively. Since the radial direction in the AdS geometry is related to the energy scales in the field theory, we can see the thermalization as happening due to the interaction of a range of energy scales, rather than a scalar quantity equilibrating over time. We end this section by discussing how the different scales of the problem contribute to the thermalization.
In appendix A we discuss a method for solving the perturbative problem using interpolating Chebyshev polynomials, which is far less costly computationally than the finite difference methods used in [14,15]. This method leads to equivalent results compared to the finite difference method used in the above-mentioned papers. We also discuss convergence properties of this method.
In the present paper we will focus on a bulk spacetime of dimension d + 1 = 5, with a scalar field of mass 1 m 2 = ∆ (∆ − d) = −3. We emphasize that we specialize to these cases by way of example, and the methods and algorithms described in this paper can easily be adapted for different d and ∆ (with the restriction d 2 ≤ ∆ < d). This will serve as a prelude to the new solution of the full non-perturbative backreaction of the scalar field on the AdS-black brane geometry. 1 We set the AdS radius to 1. 5 2 The physical setup The physical system we would like to study is that of a scalar field φ on an AdS-black brane spacetime. The evolution equations of the metric and scalar can be found by varying the five-dimensional Einstein-Hilbert action Of particular interest to us is the case m 2 = −3, since the scalar field is then dual to a fermionic mass operator with ∆ = 3 in a thermal N = 2 * gauge theory living in four flat spacetime dimensions [29][30][31]. We use the background ansatz of an infalling Eddington-Finkelstein metric and a scalar field which depend only on the radial and time directions in the spacetime, while being isotropic in the four transverse directions: In (2.2) r is the light-like radial coordinate of the spacetime, v is the time coordinate and y are the coordinates corresponding to the spatial directions on the conformal boundary. We would like to send in a scalar field φ (v, r) from the boundary of this spacetime at r = ∞. Varying the metric and scalar field in (2.1) leads to the equations of motion [14] 0 =Σ ∂ r (Σ) + 2Σ ∂ r Σ − 2Σ 2 + 1 12 whereḣ for any h.
Setting the scalar field to zero, there is no longer a source for dynamics in the spacetime. We are then left with a static, planar black hole metric which can be described by the line-element [14,15] ds 2 = −(r 2 − µ 4 r 2 ) dv 2 + r 2 (d y) 2 + 2drdv, (2.5) 6 where r = µ is the position of the event horizon. Of course, this parameter also sets the temperature of the corresponding plasma in the boundary theory, i.e., T = µ/π. Since the scalar field is initially zero when we turn on the quench in the asymptotic past, the above static spacetime is the initial equilibrium configuration of our system and µ sets the initial temperature in our thermal quenches. As the mass coupling of the fermionic operator is switched on, the changing boundary conditions excite the scalar field in the AdS-black brane background, collapsing into the black hole. The scalar field excitations evolve and backreact on the metric, and the bulk fields reach a different equilibrium configuration in the asymptotic future. In this final static configuration, the metric will be modified from its initial form in (2.5). In particular, the black hole would have grown due to the energy it absorbed from the infalling scalar field excitations. The scalar field will also have a nonzero static profile because of the new boundary conditions imposed at asymptotic infinity.
Beginning with an initial state allows us to simplify the analysis of the quenches.
In particular, the initial state is provides an energy scale, i.e., the initial temperature T i , and we will only study quenches where the final mass m f of the fermionic operator is small compare to that scale. That is, we only consider quenches where m f /T i ≪ 1, following [14,15]. In the dual gravitational description, this choice corresponds to treating the scalar as a perturbation on the background geometry. In other words, we assume that the black hole is very large so that it is possible to perform an expansion in the amplitude of the scalar field in equations (2.3), as in [14,15]. To leading order in its amplitude the scalar field equation becomes the equation of a scalar field on the static background metric (2.5). In references [26,27], the analytic profile of the static perturbative scalar field was found for both particular and general values of m, respectively. In the special case where m 2 = −3 that we will be considering, the scalar field was found to have the profile [14] φ where ℓ parameterizes the amplitude of the bulk scalar and also the value of the mass coupling in the boundary theory. Regardless of the dynamics during the evolution of the quench, the scalar will relax to profile of the above form in the final equilibrium configuration. In [14,15], it was still necessary to know the full evolution of the bulk scalar in order to calculate late-time quantities such as the change in the stress-energy tensor of the dual field theory and the change in temperature, which were determined 7 in terms of integrals of the normalizable mode over time.

Dimensionless coordinates
We introduce the dimensionless coordinates by scaling out a factor of the black hole horizon position µ (which has units of energy) After scaling out the appropriate factor of µ, the warp factors become Despite having dimensions of length, a more careful analysis shows that the scalar field should not be rescaled by a factor of µ −1 . The new radial coordinate ρ is particularly useful, since now the conformal boundary of the spacetime is located at ρ = 0. As noted above the dimensionful constant µ can be interpreted as the initial position of the black brane horizon. Therefore, the initial horizon is now situated at ρ = 1. With these redefinitions, the field equations (2.3) become the dimensionless equations These are the evolution equations for the scalar field and the two warp factors in the metric. Along with these, the Einstein equations provide two constraints, namely 8 Used in combination, (3.6) and (3.7) determine the response of the warp factor a up to an arbitrary integration constant [14].

m 2 = −3
Specializing equations (3.3) -(3.7) to a scalar with mass m 2 = −3, we find an asymptotic solution to the scalar and warp factors as ρ → 0 of [14] φ =p 0 ρ + p ′ 0 ρ 2 + ρ 3 p 2 + ln ρ where p 0 , p 2 and a 2 are functions of τ , a prime here denotes a derivative with respect to τ . The coefficient p 0 is the so-called 'non-normalizable mode' or the source coefficient [32]. Here we will choose this coefficient to have a time-dependent profile, implying that the scalar field is sourced at the conformal boundary of the AdS spacetime and excitations are sent into the bulk geometry in a time-dependent manner. The coefficient p 2 is the so-called 'normalizable mode', or the response coefficient. This is the coefficient which is to be determined given a source p 0 . While analytic solutions of p 2 are known when p 0 varies very slowly from time τ = −∞ to τ = +∞ [14,15], as well as for p 0 made time-dependent abruptly and over a very short period of time [16], no analytic solutions for the normalizable mode are currently known for a source with general time-dependence.
The solutions presented in (3.8) are for the full nonlinear equations. There is an additional constraint on a 2 coming from eqs. (3.6) and (3.7) [14]: The full warp factors can in principle therefore be determined completely given p 0 and The nonlinearities in the equation determining the scalar field make it challenging to extract the response coefficient p 2 . For this reason, in [14,15], the scalar field was treated as a perturbation on the spacetime, linearizing the Klein-Gordon equation (3.3). It then becomes a simple procedure to numerically determine the response. In the following, we also carry this amplitude expansion to second order in the metric coefficients in order to determine the leading-order backreaction of the scalar field on the background.
In the following subsection we will describe the asymptotic solution of the scalar field and warp factors in this perturbative regime. In the appendix we show how to solve the system using Chebyshev interpolation methods, which allow us to find the full profile of these metric perturbations, rather than single terms in the asymptotic expansion.

Leading-order backreaction
Since in our analysis, the scalar field backreacts only perturbatively on the spacetime, we are implicitly probing the limit of a very large black brane in the AdS spacetime.
This is the gravitational dual of switching on an operator in a thermal plasma at a very high temperature. In the boundary theory, the dual of the expansion in the amplitude of the scalar field is the expansion in m f T 1 ≪ 1, where m f indicates the mass of the fermionic operator and T i is the temperature of the initial thermal state, as described above. As discussed in detail in [14,16], the long-time evolution of the warp factors found in the perturbative regime is questionable in the very short quench limit. In this limit we find that the change in radius of the black brane becomes significant.
Nonetheless, we were able to find fully general results for certain questions, namely the scaling of the response (at early times) and the energy injected into the system, as described in [16]. This justifies the perturbative approach at all quenching rates in [14,15].

Rescaling the parameters
In this paper we give the scalar field source p 0 the time-dependent profile p 0 (τ ) = 1 2 1 + tanh τ α , (3.19) where α is the characteristic timescale on which the quench takes place. It will be useful to rescale our coordinates and fields such that we can compare different quenches on the same time and length scale, which make it easier to see how quantities behave in the fast quench limit.
As discussed in [14], the required rescaling is ρ → αρ, τ → ατ and y → α y. The fields rescale as 2 a → a/α 2 , s → s/α and φ → αφ. In these rescaled coordinates the horizon of the black hole will be located at ρ = 1/α, and the source of the scalar field would be p 0 (τ ) = 1 2 (1 + tanh (τ )) . (3.20) 2 More correct is to say that this is a leading order fast-quench rescaling. The rescaling to second order is a → 1 In other wordsâ is not rescaled

11
The expression for the metric remains unchanged:

Probes of thermalization
Knowing the profiles of the metric coefficients and the response coefficient in the asymptotic expansion of the scalar field, we would like to obtain a meaningful measure of the thermalization time of the field theory following the quench. Since the geometry fluctuates, but then returns to a static configuration after some time, we can conclude that the gauge theory plasma does (effectively) thermalize. An interesting question to ask then is whether the broken conformality of the theory introduces different scales for which thermalization occurs at different rates. Indeed, in [21][22][23] it was observed for Vaidya-type metrics that the theory thermalizes at the UV (short distance) range before thermalizing in the IR (large distance) range.
The Vaidya approach [21][22][23] considers a thin planar collapsing shell of null dust in AdS spacetime (an expanding shell in its original construction [33]), which produces a metric outside the shell equal to that of an AdS-black brane, and leaves the inside of the shell to be that of empty AdS spacetime. While this may seem an exotic form of bulk matter to consider, these constructions can be related to a collapsing thin shell of a massless scalar field in AdS [10]. In any event, the gravitational picture suggests that the dual field theory thermalizes instantaneously at the higher energy scales, while working its way down from the UV to the IR scales. This type of setup certainly describes exceptional quenches of the dual field theory, and it is not clear to what extent the lessons learned from these studies extend to general quenches. Hence in the present case, as in [14,15], we are considering rapid but smooth quenches, where in the gravitational dual the bulk scalar field evolves smoothly in space and time throughout the background geometry.
In [14,15], the thermalization time was approximated by observing when the response coefficient of the scalar field was within 5% of its final equilibrium value. 3 Of course, one limitation of this method is that it is essentially measuring the thermalization time using the one-point correlator of the quenching operator O ∆ . While in principle, the response of the one-point function depends on the whole range of energies from the IR to the UV, it cannot be used to distinguish between the different 3 In this paper, we will use a stricter 2% criterion -see below.
12 contributions from the different scales. It is therefore be interesting to employ nonlocal probes, as in [21,22], to study the thermalization process more carefully.
An important step in this direction was made by [19], in which the authors probed the thermalization of a periodically driven quench using holographic two-point functions and entanglement entropy. We will now extend their results for a non-periodic quench. An important difference between the current paper and [19] is that the source we use is not periodic, and can be tuned to be a step function in the case of an instantaneous quench. Periodic quenches are not truly realizable in the perturbative regime we are considering, since after a finite time the full nonlinear backreaction of the scalar field on the background must be considered.
In this section we first describe the analytic and numerical methods used to calculate the perturbation of the apparent and event horizons of the black brane. We then go on to discuss the calculation of two-point function and entanglement entropy in the field theory using holographic methods. We then show that our results agree with the results in the literature [21,22], namely that wider probes have longer thermalization times than narrower ones.

Evolution of the apparent and event horizons
As the scalar excitations are sent into the bulk geometry, and fall onto the black hole, the black hole will necessarily grow. While one would need a fully nonlinear evolution of the spacetime to see the full reaction of the geometry, it is still possible to probe the growth of the black hole horizon in the perturbative regime, as per [19].
When one speaks of the black hole horizon, it can mean either the apparent or the event horizon. In the case of a static black hole, the two horizons necessarily coincide.
In the dynamical case, they can evolve at different rates, with the condition that they coincide again once equilibrium is reached.
The apparent horizon is located at the radius where an outward pointing null geodesic stays at constant radius at that moment in time, i.e., it is the trapped surface of null geodesics. The event horizon is the surface outside which a light ray must be in order to escape to infinity. Intuitively a light ray may be able to move toward the outside of the black hole, but if it is inside the event horizon, then the apparent horizon, which is also growing outward, will eventually catch up with the light ray and cause it to fall into the black hole.
Locating the apparent and event horizons is useful, because it gives a nonlocal 13 measure for the thermalization of the quenched system. It further is a good consistency check of our numerics, since these properties of a spacetime is well understood. If our numerical methods correctly evolve the metric, we would expect the apparent horizon to always be located inside the event horizon. We would also always expect the area of the event horizon to grow monotonically as we pump energy into the black hole.
The apparent horizon is located at the radius where the expansion θ of a congruence of outward pointing null vectors vanishes (i.e., it stops expanding outwards). Working in the coordinates of equation (2.2), we we characterize such a congruence with the null The null vector k points toward the boundary of the spacetime outside of the initial stationary black hole, and points inward inside the initial horizon.
Following [34], the expansion of a congruence of affine parameterized null vectors n is given by (4.1) However, it turns out that k β ∇ β k α = 1 2 A ′ k α , i.e., k is not affine (the prime meaning the derivative with respect to r). To remedy this, we rescale k by exp{− 1 2 A ′ dλ}, where λ is the parameter along which the congruence k evolves. This ensures that the rescaled null vector satisfies the geodesic equation with λ as an affine parameter.
reference [34] then gives the expansion of k to be Substituting in for ∇ α k α , we see that θ = 0, when where the prime represents a derivative with respect to r, and the dot represents a derivative with respect to v.
In order to solve the equation, we change coordinates to the rescaled coordinates ρ and τ , in which the unperturbed event horizon is located at ρ = 1 α , α being the quenching rate. Equation (4.3) then gets modified to be α 2 ρ 2 a s ′ − 2ṡ = 0, (4.4) where the equation is now in terms of the new radial and time coordinates ρ and τ .
Expanding a and s in terms of the perturbation parameter ℓ, and using the ansatz that the time-dependent position of the apparent horizon is 1 α + ℓ 2 ρ a (τ ), we see that to zeroth order in ℓ, (4.4) is trivially satisfied. However, at order ℓ 2 , (4.4) gives an expression for ρ a , namely The entropy of a black hole at equilibrium is related to its horizon surface area by the Bekenstein-Hawking entropy formula [35] S = A hor 4 G , (4.6) G being Newton's constant. Of course, for the planar black hole under consideration, the area of the event horizon is infinite and hence we consider instead the area density of the horizon. That is, we can calculate the measure which would be integrated over the (spatial) gauge theory directions to evaluate the total area of the horizon. In a static configuration, i.e., at equilibrium, this area density can be related as above in (4.6) to an entropy density which is dual to the thermal entropy density of the corresponding plasma in the dual field theory. It was proposed, e.g., in [36], that this entropy density of the apparent horizon should have the same interpretation as the dual entropy density in the boundary theory even in dynamical situations.
As above, we use V to denote the area density or volume element of the black hole horizon. The full (dynamical) area density is then given by Therefore the perturbation of the area density of the apparent horizon is given by Note that at late times, after the system has equilibrated,ḃ = 0, and therefore at equilibrium. We expect that the area density of the black hole should increase due to the energy it absorbs, and therefore that the above perturbation should be positive.
However, this positivity is not immediately obvious, e.g., we do not have an analytic proof that δV a > 0. It is therefore a useful test of our numerics that this quantity comes out to be positive at equilibrium.
We would also like to calculate the location and area of the event horizon of the black hole. This is a more involved calculation, since it is a global property of the spacetime, and therefore cannot be read off from the fields at any one moment in time.
The equation satisfied by the event horizon can be obtained from the line-element.
The position of the event horizon is the outermost radius at a point in time from which a null ray cannot escape to infinity. Since the event horizon is an expanding null surface, an outward-pointing null ray lying on the event horizon will move outward with the event horizon and will stay at the same radius as the event horizon throughout the evolution, until it becomes stationary again when equilibrium is reached. In other words, the event horizon follows a null trajectory.
Working in the rescaled coordinate system, and working at fixed position in the transverse directions, we can therefore set the proper-time along a radial geodesic situated at the event horizon to zero: Substituting for a = 1 α 2 ρ 2 − α 2 ρ 2 + ℓ 2â , and ρ = 1 α + ℓ 2 ρ e , and dividing by dτ 2 , equation (4.11) simplifies to [19] dρ e dτ = 2αρ e − 1 2â . (4.12) Notice that at late times when dρe dτ = 0, the equation has the solution ρ e = 1 4αâ , therefore coinciding with the radius of the apparent horizon (4.5) at late times. The differential equation (4.12) has a general solution where ρ i is the radius of the event horizon at initial time τ i . Taking the limit for late times τ → ∞ in equation (4.13), we see that the prefactor exp (2ατ ) diverges. In order that the null ray (lying on the event horizon) does not shoot off to plus or minus infinity, and for the solution to make physical sense, we need This gives the position of the event horizon at time τ i as To calculate this quantity numerically, it is somewhat easier to integrate the interval −∞ < τ ≤ τ i than τ i ≤ τ < ∞, since one only needs to know the past evolution of the system (as well as the overall evolution). The position of the event horizon can then be expressed as To implement the evolution of the event horizon, we calculate the function (4.17) at each timestep. At the end of the numerical evolution of the spacetime, we can calculate ρ e as The value of ρ temp (∞) is determined by numerically taking the limit lim τ →∞ ρ temp (τ ).
The result is accurate, because we calculate the function up until late times, after the evolution has reached equilibrium.
Similar to the case above, the area density of the event horizon is given by We are now ready to calculate and compare the evolution of the area density of the apparent and event horizons. See figure 1 for the compared evolution of the apparent and event horizons for various quenching times α = {1, 1 2 , 1 4 , 1 8 }. As we argue in appendix A.4, α = 1 8 essentially corresponds to abrupt quenches. Notice that the area of the apparent horizon can decrease, but the event horizon necessarily increases monotonically with time. In fact, the apparent horizon evolution is already non-monotonic for quenches occurring at a thermal time-scale (see figure 2 for more details). Also note that the perturbation of the event horizon always has a larger area density than the apparent horizon, as expected. As α decreases, i.e., the quenches become faster, we know that more energy gets pumped into the geometry [16] and the final area density of the perturbed horizon also grows. Finally, both apparent and event horizons equilibrate to The plots are (from left to right, top to bottom) for α = 1, 1 2 , 1 4 and 1 8 , respectively. Note that the area of the apparent horizon can decrease, but the event horizon necessarily increases monotonically with time. the same area density (i.e., the same radius) towards the end of the evolution. This, in addition to the convergence tests of the code (see appendix A.3) gives us confidence that our numerical evolution captures the correct evolution of the radial profile of the metric perturbation's evolution.
Intuitively, the evolution of the perturbed horizons of the black hole should provide us with a measure of the time required for the system to return to thermal equilibrium after the quantum quench. However, to produce quantitative results, we need to provide a precise measure with which we can extract the thermalization time. Hence we define our thermalization measure for a general dynamical quantity f (τ ) as, 4  (4.20) for the apparent and event horizons (δV a in purple and δV e in blue) with the thermalization measure for the normalizable mode p 2 of the scalar field (orange), as well as its nonnormalizable mode p 0 in red. The plots are (from left to right, top to bottom) for α = 1, 1 2 , 1 4 and 1 8 , respectively. In all cases, the horizon thermalizes before the onepoint function, and this becomes more noticeable for smaller α.
which we will apply throughout the following, i.e., both here in examining the horizon behaviour and also in considering various nonlocal probes in the following sections.
From the above definition of the thermalization measure above, we see that f th (−∞) = −1, and f th (∞) = 0. Throughout the following, our criterion for saying that a quantity has thermalized will be that the corresponding measure comes within 2% of its final value, i.e., the thermalization time τ th will be defined with |f th (τ )| ≤ 0.02 for τ ≥ τ th . Figure 3 shows the thermalization measure (4.20) for δV e and δV a , i.e., the entropy densities on the event and apparent horizons, as a function of time. For comparison, we also plot the thermalization measure for the expectation value of the fermionic mass operator O 3 , i.e., for the coefficient p 2 in the bulk scalar field. The corresponding thermalization or equilibration times determined with our 2% criterion are given in table 1. In figure 3, we also shown the result of applying equation (4.20) to the source coefficient p 0 . In the figure, we see that p 2(th) makes excursions far beyond (−1, 0) while δV e and δV a remain within this range at all times. However, we note that although p 2 fluctuates much more than the horizon position, it reaches small values compared to its extrema before the horizons equilibrate. Nonetheless we see in table 1 that the thermalization time is slower for the expectation value than for the equilibration of the horizons. This relative difference becomes more pronounced as we make α smaller. Note that figure 3 is plotted in terms of the (dimensionless) rescaled time, as in eq. (3.20), and so the thermalization times in table 1 are measured in terms of the same rescaled time. The physical thermalization times 5 would carry an extra factor of α, i.e., τ physical = α τ rescaled . In the table, we see that the equilibration times of the horizon become approximately constant for small α, in terms of the rescaled time τ . Hence, in the physical time, the equilibration of the horizon perturbation therefore scales approximately as α. In contrast, as shown in the table, the equilibration time of p 2 becomes approximately constant when measured in the physical time, as previously noted in [15,16].

Analytic expression for the correlator
We now consider two-point correlators as probes of thermalization in the field theory. More specifically, we mainly consider perturbations to the equal-time two-point correlator due to the quench. This is because the mass coupling of the quenching operator in the field theory is small compared to the thermal scale, and therefore only perturbations of the correlator will be time-dependent and contain information about thermalization (as noted in footnote 4).
For ease of computation on the AdS side, we will consider the correlator of an operator with large conformal dimension (i.e., not the quenching operator). The correlator of such an operator can be calculated in the geometric optics limit by the length of a boundary-to-boundary spacelike geodesic [37,38].
Because it will turn out that the perturbations of the length of the geodesic remain finite, we needn't concern ourselves with the regularization of the static geodesic length.
As a reminder to the reader, the line element of the spacetime in dimensionless, rescaled coordinates is given by where we have defined y = µ x/α as the dimensionless boundary spatial directions.
To calculate the two-point correlator, we will calculate the length of a spacelike geodesic with endpoints at (τ = τ * , y 1 = −y m , y 2 = 0 = y 3 ) and (τ = τ * , y 1 = y m , y 2 = 22 0 = y 3 ) (i.e., with endpoints at equal times, and symmetric in the y 1 axis.) If we allow the geodesic to extend into the bulk, it will have both ρ and τ profiles that depend on y 1 . The length of the geodesic is given by where a prime denotes a derivative with respect to y 1 . The geodesic can then be viewed as the solution of the Euler-Lagrange equation obtained from (4.22) when treating L as an action.
Expanding the metric coefficients in the perturbative parameter in ℓ 2 as given in equations (3.12) and (3.13), and the time and radial profiles of the geodesic as and a perturbation of that length given by If we perform integration by parts, we can change the term involving ρ ′ 2 to a term involving ρ 2 plus a total derivative term. In the case of a geodesic, this total derivative term vanishes when integrated. We are then left with terms involving ρ 2 and τ ′ 2 . It turns out that since we are perturbing around an extremal trajectory, these two terms vanish by the equations of motion of ρ 0 and τ 0 , and we therefore needn't consider perturbations of the radial and time profiles of the geodesic in order to calculate the perturbations of its length [19]. Since the perturbations on the shape of the geodesic τ 2 and ρ 2 play no role in the calculation, we will for simplicity refer to τ 0 and ρ 0 as τ and ρ, respectively. Because L 2 depends on the unperturbed profile of the geodesic, we must first solve for ρ and τ . Since y 1 is an arbitrary transverse direction, we will simply refer to it as y.
As it turns out, it is useful to solve the problem by choosing ρ as our independent parameter, and τ and y as our dependent parameters 7 . We can find a closed form solution of τ (ρ). The independence of the integral in (4.24) on constant shifts in τ and the condition that the geodesic be smooth at y = 0 lead to the equation [19] ( Dividing equation (4.27) by ρ ′ , and using the chain rule, the equation becomes with a general solution of .
Because of the fact that the perturbations of the geodesic shape do not contribute to the twopoint correlator at order ℓ 2 , we only consider the static geodesic in our calculations. As such, we can parameterize the geodesic with either y or ρ, using the fixed endpoints ± ym α and ρm α , respectively. If τ 2 or ρ 2 do contribute, as it does in the case of the entanglement entropy, that integral must be evaluated in y-coordinates, since ρ m would change at order ℓ 2 , while we would have to make the choice of keeping y m fixed.
It turns out that the change of coordinates relates ρ = z, and where the Poincaré time will agree with the EF boundary time τ * . Replacing t with τ * and z with ρ, the above expression is identical to the geodesic contour in time given in equation (4.29). The τ (ρ) of the geodesic therefore corresponds to a constant time slice in Poincaré coordinates on the static background. Inverting equation (4.31), expressing t as a function of τ and z, we see that constant τ corresponds to an infalling null ray in the static Poincaré geometry, and constant τ in EF coordinates is actually the path of a lightray falling into the black hole from the spacetime boundary.
Since the integrand in (4.25) has no explicit dependence on y, we know that the "Hamiltonian" will be constant in y. The equation simplifies to [19] , the maximum value of ρ on the geodesic. Substituting in for τ ′ (y) in D(τ, ρ) from (4.27), (4.33) can be simplified to [19] Changing the integration variable from y to ρ in (4.24) and (4.25), and substituting in from equation (4.33), the new expressions are , The expression for τ ′ is given by (4.27), i.e., an expression depending on ρ and ρ ′ . The ρ ′ terms in the integral can be substituted by an expression depending on ρ and ρ m from (4.34). The expression for the integral therefore has no explicit y dependence, and other than knowing the solution for τ in terms of ρ, we need not know the profile of y in terms of ρ at all! The additional factor of 2 in (4.35) comes from the fact that the 25 integration limits correspond only to half the y-interval, [0, y m /α]. Making the above substitutions give us the perturbation of the two-point function as being (4.36)

Numerical calculation of the perturbed two-point function
There is still some difficulty with integral (4.36), namely that the integrand diverges near ρ = ρ m . Although it is a one-over-square-root divergence, and the integral itself will still be finite, this does pose a problem numerically. In order to avoid integrating a divergent quantity, we introduce a second change of variables: Therefore the interval in ρ, [0, ρ m ] corresponds to the reverse of the interval [0, 1] in q.
This transforms the integral into its final form This new integrand in terms of q contains no divergences, and can easily be integrated numerically.
With the numerical evolution of the scalar field-metric system, we found the profiles of the metric componentsâ and b, as described in appendix A. This was done by solving at each timestep τ i for the coefficients to the Chebyshev polynomials, and interpolating the values of the functions at the collocation points (see appendix A.2):

26
where L ρ is a numerical domain of the radial coordinate defined in (A.28), and the fixed coefficients C j,s are given by (A.67). N is the number of collocation points at which we choose to solve the functions, and also the number of Chebyshev polynomials we choose to use to model the functions at each timestep. By recording the coefficients F j gc and F j bc of the polynomials, we can calculateâ and b for any values of τ , by simply interpolating between the values of the coefficients for intermediate times. In equation (4.39) we can therefore replace τ i → τ .
In order to calculate the value of (4.38), we discretize the q-interval, and simply integrate the interpolating function calculated by Mathematica. Doing this for multiple values of τ * , we can see the full time evolution of the perturbed two-point function.
Doing this for multiple values of ρ m and the quenching parameter α, we can see how the system thermalizes at different length scales at different quenching rates.
In figures 4 and 5, we plotted the thermalization measure L 2(th) of the two-point function (as defined in (4.20)) for various values ρ m = ρ h × {0.1, 0.5, 0.9, 0.99, 0.999} (ρ h = 1/α being the horizon position in the current coordinates) of the depth that the geodesic extends into the geometry. It is straightforward to convert ρ m into the corresponding separation of the end-points in the two-point function, i.e., 2y m , by making the replacement ρ = zρm α in equation (4.34) and then integrating:  Of course, as we vary the rate of the quenches, α provides a natural scale with which to compare these separations. In particular, we examined α = {1, 1 2 = 0.5, 1 4 = 0.25, 1 8 = 0.125}. Alternatively, we can associate an energy with the two-point correlators using E 2pt = 1/y m , which is roughly the minimum energy scale to which these nonlocal probes are sensitive. For the different quenches, we might then compare E 2pt with the quenching rate 1/α. In each plot, the thermalization measure is plotted for a range of quenching times α. Since the two-point functions we calculate are for points on the boundary of the spacetime ρ = 0, we plot the thermalization measure against the boundary time τ * . We see in figure 4 that compared to the time-scale α set by the quench, the faster the quench is, the longer the two-point function takes to equilibrate. We notice in figure 5 that faster quenches still equilibrate faster in the unrescaled "physical" time ατ * . That is to say, as we increase the rapidity of the quench, α decreases faster than the thermalization time τ therm for a correlator with      (4.20), while on the right we show the same plot, but for the unrescaled equilibration time ατ (th) . The blue, purple, yellow and green curves correspond to ρ m = 0.9ρ h , 0.99ρ h , 0.995ρ h and 0.999ρ h respectively. Notice how the trends change sign from the left to the right plots.
fixed width 2 ym α . This behaviour is more accurately reflected in figure 7, where we plot the thermalization times (both rescaled and unrescaled) for different ρ m , as a function of α. We see that as α decreases ( 1 α increases) the slopes of the monotonic curves change sign.
In figure 6 we plot the thermalization measures of the perturbation of the twopoint function for α = 1 2 and ρ m = 0.1, as well as the thermalization measure of the non-normalizable mode-squared. They very closely match each other, as one might expect for small separations. For such small separations, the geodesic does not dip very far into the bulk geometry, and only "sees" the near-boundary metric, and its perturbations. Since in the near boundary limit, the metric perturbationsâ and b are proportional to p 2 0 , so will L 2 be, as can be seen from equation (4.35). In the dual field theory picture, one would say that for small separations the two-point function probes the UV.
In figure 8, we compare the thermalization times for two-point functions with different separations, but the same quenching parameter α. We notice that the larger the separation of the two points, the longer the two-point function takes to equilibrate.
Although not shown in the above plots, it is possible with wide enough separations for the two-point function to have longer equilibration times, than for O 3 , as we discuss

Analytic expression for the entanglement entropy
Another useful scale-dependent probe of thermalization is entanglement entropy (EE).
An elegant method was proposed by Ryu and Takayanagi [39,40] to calculate EE for holographic theories. In particular, The Ryu-Takayanagi prescription involves evaluating the Bekenstein-Hawking formula (4.6) on all bulk surfaces γ which are homologous to the entangling region on the boundary of the bulk spacetime. The holographic EE is then given by extremizing over all such bulk surfaces: Note that this prescription was originally proposed for static situations but has extended to consider dynamical bulk geometries in [41]. We will simplify our calcula- For simplicity, we consider a boundary entangling region Σ with a strip-geometry.
That is, the region is three dimensional, and is infinite in the directions y 2 and y 3 (regulated by K), but has a finite width in the y 1 -direction. The metric for the bulk spacetime can be expressed as A surface γ in the bulk spacetime connecting to the boundary of Σ, has a surface area given by In (4.43) the factor of 2 comes from us only integrating over half the interval of y 1 (renamed to y) in the second line, since γ is symmetric about y = 0. Following the Ryu-Takayanagi prescription [39,40] described above, the appropriate surface γ for calculating the EE is then the one that minimizes S Σ . Once again, the quench is treated as a perturbation on the spacetime, and as in the case of the correlator, this splits the entropy into the static part plus a perturbation: 8 In the rescaled coordinates, the entropy has a time-independent zeroth-order con- (where D was defined in (4.26)) while the time-dependent perturbation of the EE, 2) , is given by Once again we will drop the subscripts on the coordinates. The perturbations on the shape of the surface τ 2 and ρ 2 do enter at order ℓ 2 of the EE (although only ρ 2 actually contributes to the entropy), but we will be explicit when referring to them.
Notice that as in the case of the two-point function, the expressions for the zeroth and second order entanglement entropy have no explicit y-dependence. This is due to the simple choice of geometry of the entangling surface. Treating (4.45) as an action, we can find the conserved charge from time translation invariance, as well as the Hamiltonian like we did in the case of the two-point function.
By dividing by ρ ′ , the chain rule again leads to the solution for the time-profile of the minimal surface of The Hamiltonian to our action (4.45) will be constant because it has no y-dependence, and leads to the identity D(τ, ρ) α 6 ρ 6 = ρ 6 m . (4.49) Using expression (4.26) for D, and substituting in for τ ′ from (4.47), we have find the equation for the radial profile . This gives us the additional equation . (4.51) Using equations (4.47), (4.49), (4.50), and (4.51), we can express (4.45) and (4.46) as integrals over ρ only. This makes the calculation much simpler since we do not need to solve for the y-profile of the surface γ (e.g., we would need to solve for the y-profile in the case of a spherical entangling surface). The expression for the EE at zeroth and second order in perturbation theory then becomes and , (4.53) respectively. In (4.52) and (4.53), ǫ is the near-boundary cut-off, which we introduced since both the leading order and perturbative EE have UV divergences.
Both S Σ(0) and S Σ(2) have divergences close to the boundary of the spacetime [19].
We can identify these divergences by using the perturbation series ofâ and b near the boundary of the spacetime (i.e., in small ρ): We then substitute these series into the expressions for the integrands in (4.52) and (4.53), and expand the integrand close to ρ = 0 to find its divergent parts. Upon integrating, the divergence of the EE in terms of the cut-off is It is worth noting that the logarithmic term is universal, as pointed out in [42,43].
This divergence can be recast in the form where χ is a universal numerical constant, A Σ is the area of the entangling surface Σ on the boundary, and m f is the mass of the fermionic operator O 3 (our quenching operator) in the boundary CFT. Comparing our divergence to the desired form, we see that 2α 2 K 2 is the surface area (density) of Σ (the factor of 2 coming form the fact that there is a surface of area density K 2 on either side of the strip), ℓ p 0 is the non-normalizable mode of our scalar field, and by the holographic duality (and our conventions) is equal to m f . This means that the remaining numerical constant in the divergence 1 12 times some proportionality factor will be equal to the universal constant χ. In fact, this proportionality factor agrees with the above two references, and therefore the factor 1 12 is our universal coefficient and must agree for the logarithmic term for any entangling surface Σ. As a further test that we have the correct constant, we calculate the logarithmic term for a spherical entangling surface Σ. Showing the final result, we obtain a logarithmic divergence in ρ of 4πR 2 ℓ 2 p 2 0 12 log(ǫ), (4.58) R being the radius of the sphere on the boundary. We see that since 4πR 2 is the surface area of the spherical entangling surface, we indeed obtain the same universal coefficient of 1 12 . We should note that in our calculations of e.g., the thermalization time associated with the entanglement entropy, we will simply discard the divergent contributions in equation (4.56) and work only with the finite part of the entanglement entropysee equation (4.83) below. In fact, given the definition of the entanglement measure in equation (4.20), the area law divergence in equation (4.56) will drop out, since these agree for the EE at all times. However, the logarithmic divergence shown there will not cancel, since it is proportional to p 2 0 (τ ), and are therefore different at early and late times. Hence one might worry that the precise results will be sensitive to cutoff redefinitions. However, we do not expect this issue will effect the qualitative features determined in the following. This matter could be avoided altogether by using a renormalized version of the entanglement entropy, e.g., y m ∂ ym S Σ [44,45]. We hope to return to this approach in future work.

Contribution from surface perturbation
As mentioned, the perturbation of the static surface γ, namely ρ 2 and τ 2 does not contribute to the entropy at order ℓ 2 , with the exception of the boundary term of the EE. Expanding S Σ(0) in (4.45) in terms of ρ 2 and τ 2 as S Σ(0,0) depending only on the unperturbed profile ρ 0 and τ 0 , we obtain the integral Notice that the numerator of the factor multiplying τ ′ 2 is the left hand side of equation (4.47), and is therefore zero on-shell. We therefore only care about the terms involving ρ 2 and its derivative. Performing integration by parts on the second term, we obtain a term multiplying ρ 2 plus a total derivative term. The term involving ρ 2 combines with the first term in (4.60) to give the equation of motion for ρ, and is therefore zero on-shell. All that remains is the total derivative term, which we can integrate to evaluate at the limits of integration: Solving perturbatively for ρ(y) in (4.50) in small y, we see that factor of ρ 2 in the above equation vanishes when y = 0. We can solve (4.50) using the new coordinate In that case we can evaluate the coefficient of ρ 2 near x = 0 (y = y m /α). We see that the coefficient diverges as where δ is the cut-off in the x-direction. Since this factor contains a divergence, it is necessary to find the small-x expansion for ρ 2 in order to find potential divergent and finite contributions to the EE from the boundary term.
In order to solve for the perturbation ρ 2 of the surface γ, we derive its Euler-Lagrange equation from S Σ in equation (4.43) for both ρ 2 and τ 2 . The equations yield the equations of motion for ρ 0 and τ 0 at order ℓ 2 , and the coupled linear equations of motion for τ 2 and ρ 2 at order ℓ 4 . These equations involveâ, b, τ 2 and ρ 2 and their derivatives up to second order. These full equations also contain nonlinearities in τ 0 and ρ 0 and are too formidable to be explicitly included in this paper. Nonetheless, we can find the leading order expansions for ρ 2 and τ 2 . We do this by substituting in for the perturbation series ofâ and b in small ρ 0 and the asymptotic series for τ 0 and ρ 0 in small x.
We find the leading order degenerate solutions in terms of the boundary time τ * to be τ 2 (x) = m(τ * )x 3/4 + . . . , (4.66) where m(τ * ) + n(τ * ) = − √ 2 9 α 2 ρ 9/4 m p 2 0 (τ * ). (4.68) This means that to leading order ρ 2 has the right behaviour in x to cancel the divergence in (4.63) and have the boundary term make a finite contribution to the entanglement entropy. Although the equation above does not tell us the exact value of n, we can reasonably expect it to be of the form n(τ * ) ∝ α 2 ρ 9/4 m p 2 0 (τ * ), (4.69) with a purely numerical factor missing.
In order to solve for the numerical factor in (4.69), we solve for ρ 2 and τ 2 again in a static background. That is, we solve the system at large times, after it has fully equilibrated, and p 0 (τ * ) = 1. In that case, we can find a simpler equation of motion where κ is a constant. By similar arguments as for the unperturbed time-profile of γ in the previous subsection, we can set K = 0. The resulting equation is much simpler than we obtained in the time-dependent case, and contains only terms either independent of τ 2 , or terms that are linear in τ ′ 2 . It is therefore possible to solve for τ ′ 2 as In the above equation, F is linear in ρ 2 ,â, b and their derivatives, but is nonlinear in the unperturbed profile ρ 0 and τ 0 of the minimal surface.
Substituting in for τ ′ 2 (and τ ′′ 2 ) in the static Euler-Lagrange equation for ρ 2 , we have an equation where the only unknown function is ρ 2 . Solving again for ρ 2 perturbatively in x, we can find an exact solution for its leading coefficient. In the static case Substituting back for the leading-order solution of ρ 2 evaluated at the cut-off δ, we see that the boundary term has a finite contribution to the EE of δS = α 2 ℓ 2 K 2 5 36 p 2 0 (τ * ) .

(4.75)
We note that it is necessary to work with ρ 2 (y) rather than y 2 (ρ) in the treatment above. To leading order in ℓ the profile of γ remains static, and we can parameterize it with either ρ or y. The expressions (4.45) and (4.52), as well as (4.46) and (4.53) are therefore related by a change in coordinates. However, when dealing with the perturbations of the minimal surface, depending on whether one lets it fluctuate in the y-direction or ρ-direction, either y m or ρ m will be corrected for by the perturbations.
We must therefore make a choice of which of these parameters to keep fixed. Since y m is the field theory observable, i.e., it determines the width of the strip in the boundary, we choose the perturbations of ρ and τ to be functions of y, allowing for ρ m to be adjusted at order ℓ 2 , though not affecting the results calculated on the static surface γ.
We further note that this is a very unexpected result, since in the case of unperturbed holographic EE, the boundary term is typically ignored as it is assumed to vanish from the equations of motion -see for example [17,19,42,43], in which the boundary term is implicitly set to zero in the presence of a relevant perturbation. We point out that in the case of a spherical entangling surface on the AdS boundary, the boundary term above contributes both a log-divergence as well as a finite contribution. While the divergence in that case is readily obtainable by the same perturbative treatment as we did for the strip, the boundary term also has a finite contribution depending on a normalizable mode, which would require knowing the full profile of ρ 2 in y in order to be extracted. That is outside the scope of this paper.

Regularization of the entanglement entropy
In equations (4.52) and (4.53), the integrands have inverse square-root divergences near ρ = ρm α . This is not a problem mathematically, but since we have to numerically integrate these expressions, our results would be more accurate if the integrands didn't diverge at all.
In order to separate the EE into a finite and divergent part, we use the perturbation series of the metric perturbationŝ since only the leading-order terms in these series contribute to the divergence of the EE. We therefore have a finite part of the EE, , (4.78) as well as the divergent part . (4.79) In order to regularize S Σ(2)(div) , we add the counterterm There is some ambiguity in the counterterm (4.80), namely that it need only have the right asymptotic behaviour to cancel the divergence of the integrand in (4.79). That means we could e.g., use p 2 0 (τ (ρ)) instead of p 2 0 (τ * ) in (4.80), and it would still cancel the divergence of the EE, but yield a different finite result for the regularized EE. In order to circumvent this ambiguity, we need to add back the finite contribution that gets subtracted in the counterterm, that is, the non-divergent limit of the integral. Therefore, we need to add back the finite contribution to obtain a finite EE that is invariant under this particular regularization scheme.
In order for these equations to be accurately integrable numerically, we make the same change of coordinates (4.37) as we did for the two-point function, namely This change of equations yields the new full expression for the regularized EE at order δS being the boundary term defined in equation (4.60), and its final expression before change of coordinates shown in (4.75).

Numerical calculation of the entanglement entropy
Using equation (4.83), we can calculate the evolution of the entanglement entropy for different quenching rates α and for entangling surfaces with different widths y m (corresponding to different surface heights ρ m ). The procedure is very similar to the one described for the two-point function, and we give only a brief overview here.
We calculate the EE by discretizing the integrand in the first term of (4.83) in q, and then integrating the interpolating function instead. This shows a speedup in the  (4.20)) for different-sized entangling regions. The evolution is a function of the rescaled boundary time τ * . The plots are, from left to right, top to bottom, for ρ m = 0.1ρ h , 0.5ρ h , 0.9ρ h , 0.99ρ h and 0.999ρ h . In each plot the thermalization measure is shown for quenching parameters α = 1 (blue), α = 1 2 (purple), α = 1 4 (brown) and α = 1 8 (green). Note that the smaller α is, the longer thermalization takes, in this rescaled boundary time.
numerical calculation, without noticeable loss of precision. The other terms in (4.83) do not require such a discretization procedure, since they do not involve the numerical metric componentsâ 2 and b 2 .
The metric components are calculated from their numerically calculated Chebyshev coefficients, as described in the appendix.
The evolution of the perturbation of the EE is seen by calculating it in a range of boundary times τ * .
We plotted the thermalization measure of the regularized EE at different quenching parameters α for the different sizes of the entangling surface Σ in figures 9 and 10, 9 as we did for the two-point functions in section 4.2.2. We see a similar behaviour as we saw for the two-point functions, namely that the faster quenches have longer equilibration times as measured by the rescaled boundary time τ * than the slower 9 Recall that the thermalization measure is only calculated for the finite part of the entanglement entropy in (4.83).  . On the left we show the rescaled thermalization time τ (th) , while on the right we show the same plot, but for the unrescaled thermalization time ατ (th) . The blue, purple, yellow and green curves correspond to ρ m = 0.9ρ h , 0.99ρ h , 0.999ρ h and 0.9999ρ h respectively. Notice how the trends change sign from the left to the right plots. The plots are (from left to right, top to bottom) for α = 1, 1 2 , 1 4 and 1 8 , respectively. Each figure contains the plot for a minimal surface of height ρ m = 0.1ρ h , 0.5ρ h , 0.9ρ h , 0.99ρ h and 0.999ρ h , respectively. The plots for 0.9ρ h 0.99ρ h and 0.999ρ h are orange, bright blue, and red, respectively. We also plotted p 2 (th) in dashed lines for comparison.
We can see that the larger the entangling surface Σ (i.e., the depth ρ m ), the longer the thermalization time of the entanglement entropy is in each case. quenches for each surface size. We also see that, as in the case of two-point functions, the faster quenches have faster equilibration times when we measure the thermalization in unrescaled boundary time ατ * . We also plot these opposite trends in figure 11 as we did in the two-point function case.
We also plotted S Σ(2) (th) for each quenching parameter α separately in figure 12, but for the different sizes of the entangling surface (measured by the depth that the minimal surface γ extends into the bulk). We again obtain similar results as for the two-point functions, namely that the EE of the larger entangling regions equilibrates slower, but that the thermalization time for fixed ρ m decreases at a slower rate than α. Since it is possible to have longer thermalization times for the EE than for the one-point function for larger α, we believe that it may be possible to obtain larger thermalization times for arbitrarily small α if we let the entangling surface Σ be large enough.

Scaling of the thermalized correlator and entropy
After reaching equilibrium, we expect our Yang-Mills plasma to satisfy equilibrium thermodynamics. At the level of the thermal entropy, we know that [15] 4.84) meaning that up to constant prefactors, the above relation gives the equilibrium behaviour for the system. Here T f is the final temperature of the system, and λ is the field theory coupling of the quenching operator. It should also be noted that λ/T f ∝ ℓ relates the small parameter in the AdS picture to the coupling λ in the field theory picture.
For wide entangling regions, as the ones we considered in the previous subsection,the minimal surface γ will become wide, with the largest contribution coming from the part deep in the bulk. As the surface becomes wide in the transverse y-direction, more of it will lie close to, and parallel with the horizon of the AdS black brane. Most of its area will come from a surface that almost coincides with a part of the horizon of roughly the same width. In the dual picture, since the entanglement entropy is proportional to the area of γ, an entangling region will have the largest contribution of its EE be proportional to the thermal entropy. In the limit of infinitely wide surfaces, the EE is simply equal to the thermal entropy [46]. On the right, we plot the perturbation of the EE for different quenching rates, namely α = 1, 1 2 , 1 4 and 1 8 corresponding to the coloured plots blue, purple, yellow and green, respectively. We show the perturbations of the entropy for different values of y m , as well as the best-fit straight lines through the data. The data in each case is clearly well approximated by straight lines.
Equation (4.84) now implies that both the zeroth order and and second order EE (expressed as O(ℓ 2 )) should be proportional to T 3 f . In this paper we have thus far kept the dependence on the temperature hidden, by setting the black hole radius µ (in unrescaled r-coordinates) to 1. It happens that the temperature is proportional to µ, so we should reintroduce µ, as well as the AdS radius L, to see what behaviour to expect from our EE.
It is easy to see that by reintroducing L into our equations, that S Σ ∝ L 3 , which already has the correct units for the area of γ. Therefore we should introduce a factor of µ for each other factor with units of length. Since the profile of ρ in y becomes proportional to y m for wide surfaces (γ becomes a slab shape), we expect S Σ(0) to scale as L 2 y m (as we verified). Therefore we need a factor of µ 3 to give the entropy the correct units.
The second order EE, S Σ(2) also has a factor of L 2 . We should expect it to also scale linearly with y m in order to balance the factor of µ 3 , as required by the arguments above. In figure 13 we show that both S Σ(0) , and S Σ(2) scale linearly with y m , as we would expect it to. Moreover, S Σ(0) has the correct slope of 2 which we would expect because of y 2 being exactly half of the width, and therefore being proportional to 1 2 of the area of the minimal surface.
We therefore see that the EE scales as µ 3 ∝ T 3 f , as predicted from equation (4.84). Note that the additional scaling of T −2 f in the EE at perturbative order is contained in the perturbation parameter ℓ 2 .
By similar horizon arguments we can predict that the two-point correlator should also scale linearly with y m for wide separations. In figure 14 we see that both L 0 and L 2 scale linearly with y m for wide separations, and moreover that L 0 scales with the correct slope of 2.

Thermalization
We have so far discussed the different probes of the thermalization of the system. In this section we explore the mechanisms behind the thermalization behaviour seen in the two-point correlator and entanglement entropy.
We first discuss the thermalization times for the different probes introduced, before going on to examine how the different scales of the problem contribute to the observed thermalization. The correlator and entropy are integrals over the radius of the AdS Note that the thermalization time for the one-point function is nearly the same for α = 1 4 and α = 1 8 .
spacetime, and different parts of the profile make different contributions. We compare these contributions with the thermalization times of the integrands at fixed radii. We then go on to see how the profile of the scalar field and different components of its stress-energy tensor equilibrate. We end this section by bringing all these observations together, and speculate about the cause of thermalization at the different scales.

Thermalization times of the entanglement entropy and two-point correlator
We can ask how long the two-point function and the entanglement entropy take to thermalize for different separations of the points, or widths of the strip, respectively. Here we show the plots of the thermalization times of the EE and correlator as a function of the width of the surface and separation of the points, respectively. The thermalization time is determined by applying equation (4.20) to the EE and correlators, and choosing a thermalization threshold of 2% of its final equilibrium value.
We plotted the thermalization times of both the correlator and EE for various values of α in figure 15. As one can see for narrow surfaces, the increase in the thermalization time is not monotonic. This occurs due to the fluctuations that occur in the quasinormal modes, which are large compared to the size of the EE and correlator for small widths. For wider surfaces we see a linear growth of thermalization time with the width of the surface. Although these thermalization times observed here are smaller than that for the normalizable mode i.e., the one-point function which are also shown in figure 15 (at least for faster quenches), its monotonic nature, and its linear nature, indicates that for wide enough separations and widths, the two-point correlator and EE should have longer thermalization times than the normalizable mode.

Equilibration of the correlator and entropy profiles
We would like to know how the two-point correlator and the entanglement entropy thermalize. The thermalization time of the previous subsection is informative, insofar as it tells us that wider surfaces have longer equilibration times than narrow surfaces or separations, as well as the limiting behaviour for wide surfaces. The regions with wider separations have minimal surfaces or geodesics that probe deeper into the bulk geometry. This provides us with a clue as to what may be causing the observed difference in thermalization time, namely that the part of the surface deeper in the geometry equilibrates later than parts near the boundary.
In this subsection we will show how the thermalization of the EE and two-point correlators depend on different parts of the dual minimal surfaces or geodesics at different depths in the AdS-geometry. First, it turns out that the parts of the integrands of the correlator or EE integrals corresponding to larger ρ in the regularized (i.e., finite) version of integral (4.53) make larger contributions to the full integral. This makes sense for wide surfaces, since most of the area is near ρ m , close to the black brane horizon. Secondly we show that it is at larger ρ that the integrand thermalizes last.
In figure 16, we show the fractional contribution to the regularized entanglement entropy (4.53) That means that for ρ m = 0.9999ρ h , approximately the last 15% of the integration interval contributes 80% of the total regularized entropy.
Next, in figure 17, we plot the excitation (i.e., when the equilibration measure (4.20) of the integrand at a particular radius is more than 2% of its final value away from its initial equilibrium value) and equilibration boundary time τ * of the EE integrand (for various α and for a wide surface) as a function of its radial position, for the minimal surface where ρ m = 0.999ρ h . Note that the we say the profile "equilibrates", rather than thermalizes, since the integrand of the correlator or EE at a particular radius is not a physical quantity in the boundary theory that can thermalize. Rather, it comes to rest in some equilibrium, after which it is equilibrated.
We can conclude from these plots that the parts of the surface that lie deeper into the geometry are also generally the ones that thermalize the latest (note that is fig. 17, the equilibration curve for α = 1 is not significantly later deep in the bulk than near the boundary, while the effect becomes more pronounced for the smaller values of α).
It is precisely this part of the surface that contributes the most to the EE. Although not shown, we see a similar behaviour in the case of the correlator.
In the next subsection, we show why it may be that these deeper parts of the geodesics and minimal surfaces thermalize later than the near-boundary part.

Equilibration profile of the scalar field and its stress-energy
The scalar field encodes both the source and response of the field theory to the quench.
For this reason, we will consider the scalar field and its stress-energy as an indicator of how the energy of the quench enters the interior of the AdS bulk. It should be remembered that 1 ρ is proportional to the energy scale of the field theory. Therefore the propagation of energy into the bulk is in a sense dual to the energy of the quench being distributed through the different energy scales of the field theory -from the UV down to the thermal scale.
We show the contour plot with excitation and equilibration curves of the scalar φ ρ in figure 18. We show φ ρ rather than φ, because φ ρ is the natural quantity that was calculated in our numerical simulations. Also note that in this figure, as well as in 50 figures 19 and 20 we show a contour plot of the fields' profiles, where its values are the contours shown in the plot, while the solid coloured regions between the contours have intermediate values. We remind the reader that lines of constant τ in these plots are null rays infalling into the black brane, rather than constant time slices, as also explained in section 4.2.1 after equation (4.31).
We will also plot two of the components of the stress-energy of the scalar field, T φ 00 and T φ ρ ρ , because it is the "matter" stress-energy which sources the backreaction of the metric in the Einstein equations. The stress-energy is given by In the first line above, S φ is the part of the bulk action (2.1) containing only scalar field terms, i.e., the matter action. We show the contour plots for the these two components of the stress-energy in figures 19 and 20, respectively.
We see several discontinuities in the equilibration curves of these quantities. This however is not showing some novel physics, but is rather a remnant from the strict 2% cut-off, as seen in figure 21. That is to say that points on either side of such a discontinuity does not have very different behaviour in time, but rather one would have a slightly higher amplitude, which allows it to cross the 2% threshold at a much later time than one with a slightly smaller amplitude, giving the discrete jump in equilibration time. What is interesting in each of the plots 18 -20, is that the equilibration time deep into the bulk is much later than near the boundary. We have also included the τ profile of a minimal surface γ, corresponding to a wide entangling region at the thermalization time of the corresponding entanglement entropy. As can be seen in the three figures, the profile is mostly outside of the spacetime regions where the scalar field and the stress tensor fluctuate most. We can therefore think of these contour profiles as indicating the level of disturbance the EE (and correlator) experience at a certain boundary time τ * , from how much of the τ -profile extends into these regions.
The minimal surface or geodesic can be seen as being dragged through this contour plot of φ ρ and T (φ) , exciting the entropy and two-point correlator, until most of this profile has passed through the disturbed region and is deemed thermalized. Because the time-profile of the minimal surfaces/geodesic can stretch infinitely far into the past as ρ m → ρ h , we can expect that the thermalization time of the entanglement entropy or two-point correlator could be made arbitrarily long.  Figure 19: (Colour online) A contour plot of T φ 00 with α = 1 8 as a function of ρ and τ * . We add in the excitation and equilibration curves for the tensor component in green, as we did for the scalar field in figure 18. We also show the same contour for the thermalized entanglement entropy.  (4.20)) for the scalar field for particular radial values on either side of the discontinuity seen at ρ ≈ 3.8 in fig 18. The red curve is for ρ slightly smaller than 3.8, and the blue curve for ρ slightly larger. The blue one crosses the dashed line representing the equilibration threshold, and will therefore have a much later equilibration time than the red curve, although the behaviour of the function is very similar at the two radii.

Heuristics of thermalization
As we have seen, the nonlocal probes that thermalized most like the one-point function were those that had relatively small separations, and reflect the physics closer to the AdS boundary. Those that thermalized most slowly were those with larger separations.
The time scales here approached (and potentially exceed for wide-enough surfaces) the thermalization time of the one-point function.
For the local quantities in sections 5.2 and 5.3, we saw a wide range of equilibration times. For the quantities close to the boundary, we saw that they equilibrated on a time scale similar to the non-normalizable mode. This makes sense, since near the boundary, the dominant term in the respective asymptotic series is in fact the one containing the source term p 0 , or p 2 0 . However, deeper into the geometry the higher-order terms in the expansion containing the response coefficient p 2 will have an increasing contribution, so that we can expect longer thermalization times. This is exactly what was observed.
The only scales that we introduced in this system are the quenching parameter α (corresponding to the non-normalizable mode p 0 ), the emergent response p 2 , and the temperature of the system. The other scale that comes into play for the geodesic or entanglement entropy is the width of the probes. We should expect wider separations in the entangling surface and the points of the correlator to introduce an extra scale, since their boundaries are causally disconnected in the boundary spacetime. In figure   15 we see that the thermalization time of the two-point correlator and entanglement entropy grow (at least, roughly) linearly as a function of the separation in the function.
The thermal wavelength of the dual field theory is λ T ≡ 1 T = π, given our conventions [15]. One might expect that if the thermalization process occurs quasilocally in the field theory, correlation functions or entanglement entropies on scales larger than the thermal wavelength (i.e., involving points or boundaries separated by more than λ T ) should thermalize with approximately the same time. Nevertheless, we see monotonic, linear growth in these probes' thermalization times for separations 2 y m > λ T . In fact, in figure 15, the widest separation for our two-point function is approximately twice the thermal wavelength. This behaviour suggests that arbitrarily large regions will see arbitrarily long thermalization times set by the time for these points to come into causal contact. This is the behaviour observed by Calabrese and Cardy [4] in considering the entanglement entropy of an interval in a two-dimensional free field theory, and which lends itself to a simple quasiparticle interpretation. They saw a linear increase in the entropy after the quench, until such a time as the two ends of the region would come into causal contact (as though the information was carried by quasiparticles), after which the EE would quickly thermalize. As in our case, they did not see an upper bound on the thermalization time. For holographic calculations of two-point functions and EE in a two and higher dimensional boundary theories, similar behaviour was found 10 in studies using a Vaidya metric in the bulk [21,22] and the precise evolution by which the entanglement entropy is saturated was extensively studied in [23]. Hence our results are in agreement with these other holographic studies and hence it seems that the trend of longer thermalization times at larger length scales holds true in both strongly and weakly coupled field theories.
This discussion pertains only directly to perturbative quenches. However, we may expect to see similar behaviour for fully nonlinear quenches, since the response p 2 would be appropriately modified in the nonlinear regime.
10 In some Vaidya studies [21], the thermalization time at a particular length scale is estimated as the time when the entire probe is completely contained inside the collapsing shell. With this geometric definition, it is clear that the thermalization time can become arbitrarily long since in terms of the boundary time, the shell takes infinitely long to cross the location of the black hole horizon.

Conclusion
The standard toolkit of numerical relativity [34,47] faces challenges when confronted with typical problems in asymptotically anti-de Sitter spacetimes, motivated by the gauge theory/string theory correspondence [8]. The main challenge is that gravitational simulations in asymptotically Minkowski spacetimes mostly have a compact physical dependence domain; on the contrary, in AdS, control over the whole space-time, and especially near the boundary is crucial. The latter is emphasized in problems related to holographic quenches, where the temporal history of a quantum gauge theory coupling is encoded as a non-normalizable component of the gravitationally dual bulk scalar field near the boundary. In this paper we described the application of pseudo-spectral methods based on Chebyshev polynomial expansion to problems of holographic quantum quenches. We paid special attention to convergence and accuracy issues of the proposed spectral framework.
Our main physical application was the extension of the earlier work on holographic quenches in strongly coupled N = 4 supersymmetric gauge theory plasma induced by a time-varying coupling of certain dimension ∆ = 3 operator [14]. Here, having access to the full bulk metric (albeit only to the leading order in the gravitational scalar backreaction), enabled us to compare local and nonlocal probes of the ensuing thermalization process. Specifically, we compared the relaxation of event and apparent horizons, the equal time two-point correlation function of operators of large conformal dimension, the entanglement entropy of strip-shaped regions with the relaxation of the one-point correlation function of the quenching operator. The nonlocal probes of thermalization were discussed earlier in the literature [21,22]. In fact, our discussion of non-local probes parallel that of [19]. The important difference is that the authors of the latter work considered periodically-driven holographic quenches, in which the boundary gauge theory never reaches an equilibrium, thus making the comparison of various thermalization criteria impossible.
As a criterion for thermalization of a probe f we considered the quantity (4.20) . ] ∝ α 4−2∆ (∝ α −2 for ∆ = 3) in the limit of fast quenches, i.e., as α ≪ 1 [16]. Probably our most dramatic finding is that for all nonlocal probes discussed (and for wide range of probing energy scales, if applicable), f non−local th remains finite in the limit of fast quenches. As a result the thermalization of non-local probes, within criteria (6.1), appears to be faster than that of O ∆ . This effect becomes more pronounced as the quenching rate 1 α increases. While we focused in this paper on d = 4 boundary space-time dimensions and for ∆ = 3 of the quenching operator, we believe that this observation would extend to general d and ∆ insofar as d 2 < ∆ < d. We now finish with open problems. First, it is important to lift the restriction of 'leading order backreaction' -for this, one needs to extend the proposed numerical framework to the full nonlinear evolution. We believe that this does not pose conceptual or technical difficulties: all the numerical steps can be easily generalized, even our spectral uplift procedure from g c (τ, ρ) to a c (τ, ρ) (see appendix A.2). Likewise, using realistic dual scalar potentials, as in [29], does not pose a problem either. Another benefit of the spectral approach is that it relatively easy allows for a generalization to spatially non-homogeneous and non-isotropic quenches. We hope to report on the latter problem in a future work.
The φ log , a log and b log terms remove (subtract) logarithms 11 close to the boundary (ρ → 0) in the asymptotic expansion ofφ and the warp factors a and b, while staying bounded close to the horizon (ρ → 1 α ). Of course, there is a choice in selecting φ log , a log and b log . For φ log (τ, ρ) = φ log (p 0 (τ ), ρ) we choose where the coefficients F i (p 0 (τ )) are (uniquely) adjusted in such a way that the resultinĝ φ c are free from ln ρ up to terms O(ρ 9 log ρ). Explicitly, the first few coefficients F i (p 0 (τ )) are Note that the subtraction φ log remains bounded all the way to the horizon for fast quenches α ≤ 1. Similarly, we take Ideally, we would like to subtract as many log-terms near the boundary as possible; this would make spectral expansion of the functions more precise. It is possible to expand (A.2) to arbitrary order: for any i, F i depends on a source p 0 and its higher time derivatives, and thus is known analytically for our quenches, where The asymptotic expansions forâ and b contain log-terms with prefactors that, in addition to the functional source dependence, depend on response function p 2 (τ ) and its derivatives. Specifically, both B 1,i and A 1,i for i ≥ 6 depend on the derivatives of the response p 2 (τ ) up to order (i − 5). We can extract reliably p 2 (τ ) from the evolution of the φ c : however, we find that the errors in extracting derivatives of p 2 (τ ) does not justify truncating the (A.4) beyond the terms employed. Explicit expressions of the first few logarithm prefactors in (A.4) are given by We further define the functions that occur naturally in equations (3.14) -(3.16), namely π, β and g c : As in (A.2), we subtract the logarithmic terms of π(τ, ρ) near the boundary π(τ, ρ) = π c (τ, ρ) + π log (τ, ρ), We now present the equations which we separate into the evolution (containing time derivatives of the functions) and the constraint (without time derivatives of the functions) ones, evolution equations: constraint equations: One additional constraint equations is obtained combining (3.15) and (3.16). First, using the second equation in (A.8) we rewrite the latter equation as (A.17) Algebraically solving for β(τ, < r) from the first equation in (A.16), we can represent the remaining equation in (A.16) as (A.20) Note that given g c (τ, ρ), and using the definition ofâ in (A.1), we can always reconstruct a c (τ, ρ) as So far, we have not used (3.18) -this equation contains two τ derivatives, and so appears to be an evolution equation. It turns out that this is a momentum constraint, and should be imposed at a single spatial point, say ρ = 0; the equations (3.14)- (3.17) guarantee that (3.18) would then be true at any other point. The latter constraint determines a 2,2 (τ ) (see (3.12)) in terms of the source p 0 (τ ) and the response p 2 (τ ): In practice, we find it convenient to introducê a 2 ≡ a 2,2 + 5 36 : As explained in [14], there is no need to impose the boundary condition at the horizon provided we extend the radial integration past its location:

A.2 Numerical implementation
We use pseudo-spectral methods In practice we choose τ initial = −7.5, corresponding to the source value p 0 (τ initial ) ≈ 3 × 10 −7 ≪ 1 (see (A.5)); and τ f inal = 12.5 (or later for smaller α). In a nutshell, any function f (τ, ρ) we represent as truncated sum over Chebyshev polynomials T j (x), All constraints equations are then reduced to linear-algebraic equations evaluated at N-collocation points [28]. We use fourth-order Runge-Kutta method (RK4) to evolve functions in time.

A.4 Limit of abrupt quenches
Quantum quenches have two scaling regimes: the adiabatic one (α ≫ 1 ), and the regime of the abrupt quenches (α ≪ 1). The former one represents an expected slow, hydrodynamic response of the system to external forcing [14]. It was observed in [14] (and further studied in [16] ) that within a holographic framework a QFT exhibits a scaling response in the limit of abrupt quenches as well. The same scaling was observed outside of holography in a CFT deformed by a relevant operator [49]. Our code is ideally suited to study fast quenches since we use the α-rescaled scalar and metric variables [16]. Figure 25 illustrates the response function p 2 for fast quenches, {α = 1, 1 2 , 1 4 , 1 8 }. Note that the response becomes almost indistinguishable between α = 1 4 and α = 1 8 quenches. We take α = 1 8 to correspond to abrupt quench.