On entanglement spreading from holography

A global quench is an interesting setting where we can study thermalization of subsystems in a pure state. We investigate entanglement entropy (EE) growth in global quenches in holographic field theories and relate some of its aspects to quantities characterizing chaos. More specifically we obtain four key results: 1. We prove holographic bounds on the entanglement velocity $v_E$ and the butterfly effect speed $v_B$ that arises in the study of chaos. 2. We obtain the EE as a function of time for large spherical entangling surfaces analytically. We show that the EE is insensitive to the details of the initial state or quench protocol. 3. In a thermofield double state we determine analytically the two-sided mutual information between two large concentric spheres separated in time. 4. We derive a bound on the rate of growth of EE for arbitrary shapes, and develop an expansion for EE at early times. In a companion paper arXiv:1608.05101, we put these results in the broader context of EE growth in chaotic systems: we relate EE growth to the chaotic spreading of operators, derive bounds on EE at a given time, and compare the holographic results to spin chain numerics and toy models. In this paper, we perform holographic calculations that provide the basis of arguments presented in that paper.


I. INTRODUCTION AND SUMMARY OF RESULTS
A global quantum quench -unitary time evolution from a translation invariant short-range entangled initial state -is an interesting setting in which we can study thermalization in closed systems and probe their chaotic dynamics. Studying the real time dynamics of strongly interacting many-body systems or field theories is a formidable task. Holography has the potential to provide invaluable insight into the dynamics of these systems, by mapping the quantum problem into a higher dimensional classical geometric one [2][3][4][5][6][7]. In conformal field theories in two spacetime dimensions, global quenches have been thoroughly explored [8][9][10]. In two spacetime dimensions one can also study the problem in integrable lattice models [8,[11][12][13], and using numerics in nonintegrable chains [1,14,15]. One can also study free theories in higher dimensions [16,17]. In this paper, we will concentrate on higher dimensional chaotic systems.
Many aspects of the questions we ask in this paper have been already understood in prior work [18][19][20]. See also [21][22][23][24] for early work on quenches in holography. The importance of the problem however warrants further scrutiny, and the results of this paper should provide valuable data for efforts into understanding quenches in strongly interacting chaotic systems. In a companion paper [1], we use the results of this paper along with numerical results from spin chains, to propose that the entanglement entropy in a quench in a chaotic systems is close to saturating a combination of two constraints: one that follows from recent insight into quantum chaos [25][26][27][28] and the positivity of relative entropy [29], and another bounding the rate of growth of entropy. [1] can be read as putting the results of this paper into context, relating them to the picture of the chaotic growth of operators, and analyzing them from the point of view of toy models and bounds.
The outline of the paper and the summary of the results is as follows. In Sec. II we first review the gravity duals of global quenches. Second, we present a new computation of how close Ryu-Takayanagi (RT) surfaces [5,6] approach the horizon of a static black brane with metric with horizon at z h = 1. We work in the limit where the characteristic size R of an arbitrary shaped boundary entangling surface Σ (on which the RT surface is anchored) is large compared to inverse temperature β: where δz ≡ 1 − z, R insc is the radius of the largest ball inscribable in Σ, and µ is defined in (2.10).
It is explained in [1] how this result can be used to compute the butterfly effect speed v B that characterizes the growth of operators in a thermal state [25,26]. Third, we turn to the time evolution of entanglement entropy in a global quench. It was derived for holographic theories in [18][19][20] that at early times the entropy grows linearly 1 where β is the effective inverse temperature associated to the system, s th is the thermal entropy density, A Σ is the area of Σ, andŜ(t) ≡ S(t) − S vacuum is the entropy with the vacuum contribution subtracted. 2 We will sometimes refer toŜ as the extensive piece of the entropy, as for t ∼ R it scales with the volume. This equation is expected to hold in any chaotic system. 3 We review the holographic derivation of (1.4) for the (simplest) strip geometry in Sec. II using Hubeny-Rangamani-Takayanagi (HRT) surfaces [7]. Through this example, we introduce the entanglement velocity v E controlling the early time linear growth of entropy [18][19][20]. We then move on to our main results, which are organized in an order of increasing of geometric complexity: 1. In Sec. III we discuss several inequalities involving v E and v B valid in holographic theories obeying the Null Energy Condition (NEC). We prove that in holographic theories This inequality is proven in [1] for any unitary quantum system, so (1.5) is a consistency check for holographic theories in the dynamical regime. 4 The two velocities are controlled by different parts of the geometry, v E by the region behind the horizon and v B by the near horizon geometry, hence it is somewhat surprising that one can prove a relation between them.
1 See [30] for a recent study of small regions. 2 The vacuum contribution is a divergent area law. 3 It is also obeyed in the quasiparticle model of entropy growth [31] and in free scalar theory [17]. 4 That vE should be smaller than vB was also discussed recently in [32].
We also prove several inequalities valid in holographic theories: , (1.6) where the (S) superscript refers to the Schwarzschild black brane value given in (2.17) and (2.30), which only depends on the number of dimensions. The first of these inequalities was conjectured in [19,20] based on the evaluation of v E in many examples. The Schwarzschild black brane is the holographic dual of a conformal field theory at finite temperature (or an eigenstate with energy density). In a given dimension we can consider quenches that create not only energy density, but also charge density, or we can investigate non-conformal systems. These setups lead to final state black branes different from Schwarzschild. The inequalities (1.6) imply that these complications slow down the spread of entanglement and operators, and also how close entanglement growth in the linear regime can come to saturating (1.5).
2. In Sec. IV we discuss the entropy growth for spherical regions of radius R. In the limit (1.2), we show that the extensive part of the entropy is the same irrespective of the details of the quench. While this finding conforms with the field theory intuition that all short-range entangled states should have the same dynamics in a strongly interacting chaotic system after the local thermalization time t loc ∼ β, the disparity in geometry between the end of the world brane quench model [18] and the infalling shell model [19][20][21][22][23][24] makes this universality a nontrivial new result.
We first use numerics, and motivated by the numerical results we take a scaling limit of the extremal surface equations which leads to major simplification. We are able to solve these equations analytically. For the black brane with metric (1.1) the entropy as function of time is given by: , (1.7) where 1 ≤ z f ≤ z HM , and z HM is the locus behind the horizon, where −f (z)/z 2(d−1) is maximal.
The equations are a bit complicated, but are completely explicit. On Fig. 1 we show the evaluation of the integrals for the d = 4 Schwarzschild black brane. determined by (1.7) may be non-single valued, leading to the HRT surface changing discontinuously and t S > R/v B . We will refer to this situation as discontinuous saturation; although the entropy as a function of time remains continuous, its derivative is discontinuous.
We analyze necessary conditions for discontinuous saturation in Sec. IV D by analyzing (1.7).
We prove that in d = 3 the entropy always saturates discontinuously, but in higher dimensions Schwarzschild black branes (e.g. in d = 4 shown in Fig. 1) give continuous saturation and In [1] we give an alternative proof of this equation for any black brane that gives continuous saturation based on entanglement wedge subregion duality. The proof applies to any theory of gravity, and in [1] we verify the relation (1.8) explicitly in four derivative gravity.

3.
In [19,20] an "entanglement tsunami" picture for entropy growth was suggested: there would be a sharp wave front propagating inward from the entangling surface Σ. Behind this wave, degrees of freedom become fully entangled with the outside, while the degrees of freedom yet to be reached by the wave remain unentangled. At early times the wavefront moves with v E giving the relation (1.4).
In Sec. V we probe this proposal in the thermofield double state [33], where degrees of freedom in the right copy (R) are strongly entangled with their partner in the left copy (L). If we time evolve only R, this high degree of local entanglement spreads in R. The technology developed in Sec. IV can be used to compute the two-sided mutual information [18,34] between two (concentric) spherical regions A L (0) and B R (t), where the sphere B in R has been time evolved for time t with respect to A L (0): We regard L as auxiliary system that keeps track of the movement of information in R, and ] as the "mutual information in time" [35] between A R (0) and B R (t).
Let us fix R A and some time t. The "entanglement tsunami" picture would predict that information in R A only gets scrambled in a region of size R A + v E t at least for t R A . However, in Sec. V we find that as a function of R B the mutual information only saturates at which is valid for any time t. More discussion of the two-sided mutual information and the interpretation of our results can be found in [1], where we advocate for this quantity as a useful probe of information spreading in any quantum system. In this paper we concentrate on the computational aspects of obtaining the two-sided mutual information for a holographic theory.
There is another result that may help in ameliorating the tsunami picture. In [1] using ideas from [29], we derive an upper bound on the entropy in quencĥ where the condition on time comes from having taken the scaling limit. It was suggested in [19,20] that this equality may hold in holographic systems. (1.13) was proved in the quasiparticle model in [31]. It plays an important role in the discussions of [1]. There are rigorous bounds of the form (1.13) but with a coefficient that depends on the Hilbert space dimension and operator norms [36,37]. It would be interesting to establish (1.13) using field theory techniques. 6 The same methods that allowed us to establish (1.13) give us an intriguing equation for the rate of growth of entropy. It is easiest to state our result in the boundary state quench model [8], whose holographic dual is the eternal black hole cut in half by an end of the world brane [18].
The extremal surface determining the entropy at a given time, is a tube connecting the entangling surface Σ to its image on the brane Σ im . The rate of growth of entropy in the scaling limit is where the area is measured in the field theory coordinates x defined in (1.1). We apply this technology to develop an early time expansion for the entropy, and obtain the first subleading correction to (1.4): where a [Σ] is given in (6.25). 7 (1.14) together with the bounds discussed in [1] is very suggestive of a tensor network interpretation, but we have not been able to make this connection precise.
The results in this paper provide us with a wealth of information about entropy growth in a quantum quench in holographic systems. The key technical result of our paper is the determination HRT surfaces in the limit (1.2). We find hints that various geometric objects, the largest inscribable ball (appearing in (1.3) and (1.12)) and the image on the end of the world brane (appearing in (1.14)) play an important role in the dynamics of entropy. We provide an improved understanding of when the saturation is continuous and when it is not. We study the "entanglement tsunami" in detail, and find that it propagates with the butterfly effect speed v B . In [1] we interpret the results of this paper. As a model for the time dependence, we suggest that the entropy is as large as possible given the two constraints (1.11) and (1.13). We also construct a microscopic toy model based on the chaotic growth of operators that saturates these bounds. It would be very interesting to compute the time evolution of entropy for shapes different from the sphere and the strip, and to see if the bound on saturation time (1.12) is saturated for shapes other than spheres (at least for certain black holes). Another interesting direction is to refine the bounds and the microscopic models presented in [1], to bring them closer to reproducing the entropy curve obtained from holographic and other chaotic systems.

A. Holographic quench models
In a field theory model of a global quench we want to create an initial state that is short-range entangled, translation invariant and has finite energy density. One way to create such a state is to dump in energy in an uncorrelated manner by smearing a local operator over the whole system and acting with it on the vacuum. This is the setup that Liu and Suh considered holographically [19,20].
Another way to model a quench that is more convenient for CFT computations is to consider a 7 It was suggested in [32] based on tensor network intuition that the early time expansion would involve these powers of t.
conformal boundary state in CFT as the initial state [8]. 8 The holographic dual of this setup was considered by Hartman and Maldacena [18]. The two setups are shown in Fig. 2.
Left: Vaidya quench model. The infalling null shell is drawn by an orange arrow. Below the shell the geometry is pure AdS, above the shell it is a static black brane. We get time evolution of the entropy in the boundary theory because the HRT surface lives in both parts of the geometry and passes through the null shell. Right: The boundary state model of a quench is dual to an end of the world brane cutting the eternal black hole in half. The lower portion of the figure illustrates how the Euclidean path integral prepares a short-range entangled initial state from the boundary state [18]. The time evolution of the entropy comes from the HRT surface entering the black brane horizon and ending on the brane.
We analyze both setups below for large regions and times R, t β, and find that the entire time evolution of entanglement entropy is universal: it does not depend on which setup we consider.
It has been already been demonstrated in [18][19][20] that in the linear regime (1.4) the Vaidya and the end of the world brane models give the same result. Because these setups differ significantly on what conditions the HRT surface has to obey in the bulk, it comes as a surprise that even the saturation behavior agrees between the two models. We discuss the details below.

B. Geometry setup
The static geometry of the most general translation and rotation invariant asymptotically AdS black brane is where the boundary is at z = 0, x is the field theory spatial coordinate, and writing g tt in this form is for convenience. To get an asymptotic AdS spacetime we need to impose the boundary conditions and we use the freedom of scaling the coordinates to set the horizon radius z h = 1, leading to the other boundary condition: A useful example to bear in mind is the one parameter family of Reissner-Nordstrom (RN) black branes with emblackening factor where the parameter 0 ≤ q ≤ 1 characterizes the proximity of the black brane to extremality, q = (Q/Q ext ) 2 . We get he Schwarzschild geometry by setting q = 0.
It will be beneficial to work in infalling Eddington-Finkelstein coordinates: .
We can extend the t coordinate inside the black brane using this equation. If we choose an integration contour that goes above the pole at z = 1, behind the horizon we have We can obtain a Vaidya geometry by gluing an empty AdS space to the above described black brane along the v = 0 null plane. The geometry can be written as In the end of the world brane setup the brane is located at the plane of time reflection symmetry, t I = 0, hence in the coordinates (2.5) its position is

C. Static surfaces
We first want to ask how close RT surfaces anchored on large boundary regions approach the horizon of the static black brane (2.1); the answer will have applications throughout the paper. We take the surface to be parametrized by boundary polar coordinates, z(ρ, Ω). The area functional in the geometry (2.1) is: where ∂ Ω is the gradient on the unit S d−2 . As explained in more detail (for symmetric surfaces) in [1,39] the minimal surfaces that we are interested in will for mostly lie flat skimming the horizon and then shoot out exponentially to the boundary. Then we can focus on the near horizon region to understand the important physics. It is convenient to define a somewhat peculiar new variable where ∇ 2 x is the Laplacian in boundary coordinates, and µ 2 > 0 because f (1) < 0. This equation has solutions blowing up as s ∼ exp (µ x) corresponding to the fast shooting out to the boundary explained above.
One can get a good approximation to the minimal surface by sourcing (2.10) uniformly around the entangling surface Σ. Using the appropriate Green's function of (2.10) this approximation amounts to where K ν (x) is the modified Bessel function of the second kind, and we will fix the constant C below. Note that (2.11) defines two disconnected surfaces (for a connected Σ), one that we are looking for and one that runs off to infinity that we should discard. The solution has the property that it blows up near Σ, which corresponds to the RT surface leaving the near horizon region and arriving to the boundary. It decays exponentially away from the entangling surface, which corresponds to the surface lying flat close to the horizon. Where s( x) gets big, we should not trust (2.11). One can actually solve for the full surface analytically in a double expansion [39]: it should involve adjusting the source in (2.11) so that it is smeared an O(β) amount around Σ of characteristic size R. The precision of (2.11) is however enough for our purposes.
The location where s( x) is the smallest corresponds to the point x which is farthest from all points on Σ. For a convex Σ this is the center of the largest inscribable ball, and this is where we put the origins of the polar coordinate system. To satisfy the boundary condition (2.9) we require that C ∝ e µ R insc , (2.12) where we used the large x asymptotics of K ν (x), and we have not written out the prefactor involving powers of R insc . In order for s( x) 2 1 and for the minimal surface to reach the boundary around Σ we need to take This is the closest approach of the minimal surface anchored on a convex Σ to the horizon.
Two applications of this result are explored in [1], which we briefly describe here. The reader should refer to that paper for a more complete treatment. First, acting with a smeared operator on top of the thermal state is dual to creating a particle near the boundary of AdS, which then falls into the black brane. Let us define δz ≡ 1 − z. Its trajectory at late times can be approximated by where we have to remember that f (1) < 0. In the field theory, the support of the information about the insertion of the operator gets scrambled inside a ball of radius v B t, and we need to access the density matrix of this entire ball to be able to reconstruct the operator insertion. The dual perspective is that we need to find an entanglement wedge that contains the infalling particle. We found above that two entanglement wedges whose boundary entangling surfaces share the same largest inscribed ball reach equally deep into the bulk, hence the minimal choice of entanglement wedge corresponds to a spherical Σ concentric with the operator insertion. In order for the entanglement wedge to contain the infalling particle we need .
This result agrees with prior determinations of v B from entirely different computations as discussed in [1]. For an RN black brane (2.4) Second, in the Vaidya setup if the saturation of entropy is continuous, saturation happens when the HRT surface climbs out from behind the shell and its tip is barely touching it [19,20]. If the saturation is discontinuous, then the time when this happens provides a lower bound on saturation time. The shell is following a null line v = 0, which in Schwarzschild coordinate is at where we went to late times or close to the horizon. (Note that before the shell the spacetime is pure AdS, but we do not need to use this fact.) This formula is identical to (2.14). The (lower bound on) saturation time is obtained by plugging δz = into (2.18) giving circumstances we get continuous saturation as explored in Sec. IV D, while for strips the saturation is always discontinuous. One of the most important questions left open in this paper is whether shapes other than the sphere give continuous saturation.

D. Review of the strip
In the end of the world brane setup the analysis of the strip is particularly simple, the results below are a review of [18] in our set of coordinates. Before saturation the two sides of the strip are disconnected, and we can concentrate on one of them. Because of symmetry, the HRT surface does not move in the boundary spatial directions, and it is just determined by the function z(v). 9 The action is where z c is the point where the surface ends on the brane. Because the action does not depend on v we have a conserved quantity: which we can solve for z . We can also determine E as a function of z c using z brane = 0. Plugging back into the area functional and changing integration variables to z we get: where we used (2.7). One can evaluate the above integrals numerically, and obtain the entropy curve (t(z c ), A(z c )). Instead, we proceed analytically.
Let us define z HM (after the authors of [18]) as the maximum of −f (z)/ h(z) z 2(d−1) . We will use frequently in the rest of the paper that at this point An important observation in analyzing the integrals is that as z c → z HM both t and A diverge [18].
This is exactly the regime that we are interested in, if we want to consider R, t β in the field theory.
We do not even need to obtain (2.22) to reach this conclusion. It is easy to see that if (2.21) gives z (v) = O(1) than the HRT surface will reach the boundary in time t = O(β). In order to go to long times we need to suppress z (v). Solving (2.21) for small v we obtain: where a = O(1) depends on the details of f, h and we do not write down its explicit expression. We see that z c → z HM indeed suppresses z (v), and should correspond to late times in the boundary theory. We can actually estimate how small δz c ≡ z HM − z c has to be. Let us define 25) and assume that δz c , δz(v) 1. (Note that this δz is not the same as the one used in Sec. II C.) Expanding (2.21) in these small quantities and solving the resulting differential equation we get that The solution looses its validity once The exponential growth of solutions, and that we have to choose the control parameter δz c to be exponentially small is analogous to the exponentially close approach of RT surfaces to the static black brane horizon discussed in Sec. II C, and it will appear throughout the paper.
In summary, we have to choose δz c exponentially small in R to make the HRT surface skim z HM for a time of O(R). The HRT surface then steeply shoots out to the boundary. It is not hard to convince ourselves that this shooting out part gives subleading in 1/R contribution to both t and the subtracted extremal surface areaÂ We thus conclude that to obtain the entropy as a function of time, we do not need to know the details of how the extremal surface reaches the boundary. We can approximate it with a surface lying on z HM for v ≈ t, and neglect the (subleading) area law contribution toÂ(t), and an O(1) correction to t coming from the shooting out part of the HRT surface. In the rest of the paper we will use a similar philosophy for more complicated geometries.
The entropy is given by the volume element on the surface z HM [18][19][20]: where we used that s th = 1/(4G N ) with our choice of z h = 1. The value of v E for a Schwarzschild black brane is This behavior is valid for any entangling surface at early times [18][19][20], but for the strip it lasts until saturation. We plot the entropy curve in Fig. 3.
We note that for the d = 2 BTZ black hole there does not exists a finite value of z HM , because −f (z)/z 2 is monotonically increasing. Hence, the computations above are not valid. Nevertheless, plugging in d = 2 into our formulas gives the correct result for d = 2 [8,23,24].
We have only analyzed the end of the world brane quench model so far. It was shown in [19,20] that the same results apply in the Vaidya setup. They can be derived from integrals very similar to (2.22) explicitly, but can also be understood just from a careful analysis of the HRT surfaces.
First, the HRT surface passes through the null shell, and caps off in the AdS region of the geometry.
This part of the surface is a chunk of an RT surface that we would get in the vacuum, hence only Entropy growth for a strip of width 2R from a d = 4 Schwarzschild black hole in the limit of large region sizes. The above rescaled curve is exactly linear with slope v E , which gives a saturation time contributes to the area law pieces of the entropy that we discard. In the black hole region of spacetime the equations of motion are the same as in the end of the world brane case, thus the only way that the surface can stay inside the horizon for a long time, is to skim z HM , as we argued around (2.24). Indeed, we can see from the plots of [20] that after crossing the shell the HRT surface reaches the vicinity of z HM exponentially fast, and stays there for a long time, before departing to finally reach the boundary. We can argue in the same way as for the end of the world brane case that only the part of the surface skimming the z HM surface is relevant in the limit R, t β. We conclude that the extensive part of the entropy (in the large strip size limit) is the same in the two holographic setups considered above. We will refer to the independence of the quench process as universality, and discuss it in more detail at the end of Sec. IV C.

III. HOLOGRAPHIC BOUNDS ON THE SPEED OF INFORMATION SPREAD
The Null Energy Condition (NEC) imposes two conditions on the functions f (z) and h(z) appearing in (2.1): where d is the number of boundary spacetime dimensions.
The idea that we are going to use is to turn these inequalities into differential equations: where a(z) and b(z) are positive functions that we take as given. For example, for the RN black The differential equations (3.2) can be solved by straightforward integration, and one can fix the three integration constants that arise using the three boundary conditions (2.2) and (2.3). Naively the solution for f (z) involves a double integral over b(z), but we can get rid of one integral using integration by parts. We can bring the solution to the nice form: where we introduced the function It is easy to check that these indeed solve the differential equations (3.2), and they manifestly satisfy the boundary conditions (2.2) and (2.3).
Now we are in good shape to prove bounds on the speed of information spread form holography.
In Sec. II we have introduced the quantities v E and v B and determined them in terms of the where f (S) (z) = 1 − z d is the Schwarzschild emblackening factor. From the definition (3.4) it is clear that H(z) is a nonnegative monotonically increasing function of z, and because b(z) is also nonnegative, the last two terms are readily seen to be positive behind the horizon, z > 1. We now focus on the first term, and show that it is also nonnegative. Again, from the definition (3.4) and using that h(z) is a positive monotonically increasing function of z we see that where in the only nontrivial step we increased the numerator by using 1/ h(z ) ≤ 1/ h(1) valid for z ≥ 1, and decreased the denominator by using Using the fact that if a function is greater than another one than this is true for their maxima as well, we conclude that with the Schwarzschild value given in (2.30).
Bounding v B is even more straightforward. Using (3.3) we get the more explicit formula for v B : , which is just the Schwarzschild value (2.17). In the only nontrivial step above, we increased the numerator and decreased the denominator again. We conclude that In [1] we argue that in any quantum system v E ≤ v B . We can test that field theory argument To prove this, we first take the explicit formulas from (3.3) and (3.8) . (3.11) We decrease the numerator by dropping its second term, and we get major cancellations between 21 the numerator and denominator afterwards: where in the second line we used the definition of H(z) (3.4) and that h(z) is a monotonically increasing function. We now plug in the maximum value of the function ( the value of the function at z = z HM is necessarily smaller than this value: where finally we noticed that the bound that we obtained is equal to the ratio of speeds for the Schwarzschild black brane (2.17) and (2.30), which is smaller than 1. Thus, we conclude that this ratio is maximal in an uncharged quench. We also take this result as a strong confirmation of the field theory argument for v E ≤ v B given in [1], which perhaps is the most important result that the NEC gives us.
One may wonder whether there is a lower bound in v E /v B , which would be interesting from the perspective of the bounds analyzed in [1]. There is no such bound however, as for a near extremal RN black brane we have hence the ratio can be arbitrary small.

IV. SPHERICAL ENTANGLING SURFACES
A. Equations of motion for the HRT surface We will look for the extremal surface parametrized as (v(ρ), z(ρ)) in the spherical case. In the Vaidya setup in the AdS region the surface is described by a hemisphere that reaches the shell at while in the end of the world picture the extremal surface has the topology of a cylinder and ends on the brane on a sphere of radius ρ c . 10 The equation of motion for the HRT surface in black brane background can be obtained by extremizing the area functional from which v (ρ) can be expressed. The equation of motion of z(ρ) is Plugging (4.3) into this equation, we get a single second order ODE for z BH (ρ), which can be solved numerically: where the last equation is obtained by integrating the field equations across the null shell [19,20]. 10 We are using the same symbol ρc for the two radii because they play similar roles, but the two setups are different.
For the end of the world brane setup the near brane behavior is dictated by the condition that the HRT surface has to end on the brane perpendicularly. (4.7) Note that v BH (ρ) diverges as ρ → ρ c . From these boundary conditions we can read off the value of E:  On the left we plotted the HRT surfaces (v(ρ), z(ρ)) in coordinates that correspond to the Penrose diagram shown on the right. On the right only the part of the Penrose diagram is shown that is covered by the coordinates (z, v). Top: HRT surfaces in the Vaidya setup, with the null shell shown with green and the horizon with blue. Darker color correspond to earlier times. At very early times, much of the surface is in the pure AdS region, and is an almost perfect hemisphere in the coordinates (2.1), deformed by the conformal mapping that gives the Penrose diagram on the right. As time evolves the surface goes behind the horizon, most of it lies on z HM , and near saturation it climbs out from behind the horizon and the entropy saturates, when the surface is only barely touching the shell. Bottom: HRT surfaces in the end of the world brane setup. At early times the HRT surface is a tube connecting the boundary theory entangling surface Σ to the brane, and the image on the brane is Σ im ≈ Σ. The linear regime of entropy growth takes place, when Σ im migrates up to z HM . We can clearly see that as we go to later times (lighter color) Σ im is shrinking, but staying at z HM . Finally Σ im migrates towards the bifurcation surface, but this is an O(β/R) effect. For similar figures, see [40]. . Finally, the red curve corresponds to the analytic prediction for R β given in (4.18). The red dashed lines show how to read off v E and v B from the plot. Because in the end of the world brane setup the initial state is different from the CFT vacuum, subtracting the vacuum contribution gives negativeŜ(t) for early times. Because this contribution is area law, for larger R this effect goes away. The data points for R = 4 come from the HRT surfaces plotted on Fig. 5. The data points for R = 15, 25 come from the HRT surfaces plotted on Fig. 4, and are somewhat hard to distinguish due to very accurate overlap. The main message of this plot is that as we increase R the data points collapse onto the red curve (4.18).

C. Scaling limit and explanation of universality
Based on our experience from the numerical solutions, we make an attempt to capture the limit of large R by making the scaling Ansatz: In both setups in this limit we expect the following scaling behavior: hence to leading order (4.5) becomes an algebraic equation: where for transparency we have introduced the scaled version of E ≡ ρ d−2 c E and ρ ≡ ρ c r.
Let us take the end of world brane setup first. Plugging in (4.8) into (4.11) we get This equation in fact determines what Z(1) = z c is, we simply have to set r = 1 in (4.12) and use (4.8) to get: which is exactly the equation (2.23) that z HM satisfies. 12 We see that contrary to the strip case, the extremal surface does not like to stay on z HM even in the large R limit, instead it moves away from it according to (4.12), see also Fig. 7. For r → ∞ we get Z = 1, so the scaled extremal surface does not get out from behind the horizon. The scaling surface is shown on the Penrose diagram on Fig. 8, which summarizes the information on Fig. 5 and Fig. 7. The scaling solution we find can be viewed as the precise realization of the "critical surface" envisioned in [19,20].
We note that (4.12) can be corrected order by order in 1/ρ c . If we wrote Z = Z 0 + 1 where the leading term is z HM for the d = 3 Schwarzschild black brane. The formula (4.14) agrees with what we find in numerics as demonstrated on Fig. 4. The most important observation is that at all orders in 1/ρ c the corrections decay as we go to the horizon, and the scaled surface stays inside the horizon. Then it must be that from this perspective nonperturbative exp(−ρ c ) effects cause the surface to exit from behind the horizon. This is reminiscent of what we found in the strip case in Sec. II D. To avoid clutter, we specialize to the Schwarzschild black brane to get: For r → ∞ the latter equation simplifies to: so we get an exponentially growing solution that has to be tamed by the exp (−#ρ c ) coefficient.
We have not analyzed (4.15) in detail, but it is worth noting that the coefficient function ofz where Z = Z(r) is the solution of the algebraic equation (4.12). Becausez(ρ) shoots out exponentially fast to the boundary, we can drop its contribution to both the extensive part of the entropy where we introduced r * ≡ R/ρ c ,Â was defined in (2.28), and we introduced a (d − 1) factor for convenience. Below in Sec. IV D we simplify the expression (4.18) somewhat, but for the purposes of this section this formula suffices. By changing r * from 1 to ∞ we obtain the parametric curvê

S(t).
There is one additional subtlety: it can happen that for a given time τ there are multiple extremal surfaces, and the holographic entanglement entropy formula requires us to choose the one with the minimal area. In the end of the world brane setup the static surface is available at all times, while for the Vaidya setup it becomes available for τ ≥ 1/v B . It can be shown (see Sec. IV D for details) that and a natural minimal assumption is that for τ < 1/v B we take the extremal surfaces that probe behind the horizon and for τ ≥ 1/v B we take the static surface. However, it can happen that the parametric curve (τ (r * ), A(r * )) itself does not give a single valued A(τ ). Then for a given τ , we have to choose the smallest A. On Fig. 6 we show the resulting curve (τ (r * ), A(r * )) against the results of numerical calculations in d = 3 for the Schwarzschild black brane and find perfect agreement. 13 In Fig. 9 we plot the resulting curves for RN black branes in d = 3, 4, and we find that the multivalued behavior discussed above does happen in some cases. In Sec. IV D we analyze analytically when such behavior takes place. We note that from the rate of growth bound (1.13) (to be proven in Sec. VI) and the result (4.19) it follows that black branes with v E v B < 1 d−1 will necessarily give a non-single valued entropy curve, as the saturation time in this case has to be t S > R v B . See [1] for related discussion. 13 There is an incredibly subtle effect for the d = 3 for the Schwarzschild black brane. The saturation is not continuous, but this is not visible on Fig. 6, as it happens for (tS − t)/tS ∼ 10 −5 . We zoom in on this neighborhood on Fig. 10. and t, we do not need to know its details. To get an identical equation to (4.12) we need the two energies, E Vaidya and E brane to agree. From this requirement and (4.14), we can derive that in the Schwarzschild case We confirm this prediction with numerical computations on Fig. 4.
We conclude that the entropy growth in the global quench does not depend on the details of the quench process, the leading order entropy is the same both in the end of the world brane and the Vaidya setup. The analytic argument is backed up by numerical data summarized on Fig. 6. There are many other quenches that we may consider. We can form a black brane from collapsing massive instead of null matter. We can also smear out the process of a quench in time by a small amount.
Because black branes form quickly, which is dual in the field theory to fast local thermalization in time t loc ∼ β, the details of the quench should not effect the leading part of the entropy. Indeed, in any quench setup we would find that for times t β the extensive part of the entropy is determined by the scaling solution analyzed above. The details of the quench only determine a small part of the surface before it reaches the scaling surface. We demonstrated this phenomenon in the Vaidya case, see Figs. 4 and 7. It is a quite satisfying finding that entropy growth for short-range entangled initial states is universal, it does not depend on the state chosen or more generally on the quench process. We found the same universality for strips, and it continues to hold for entangling surfaces of arbitrary shapes as will become clear from the discussion in Sec. VI.

D. Continuous and discontinuous saturation
Let us examine (4.18) for large r * . For this analysis it is more convenient to regard r as the dependent variable. From (4.12) we get . (4.21) We rewrite the integrals as where we defined z f such that r(z f ) = r * , and from now on we will parametrize the A(τ ) curve as (τ (z f ), A(z f )). Remarkably, we can reduce the number of integrals we need to perform, as A(z f ) can be rewritten as follows: 14 From this equation an interesting relation follows. To derive it we note that from (4.17) it follows Combing this relation with (4.23) and the definition of τ, A in terms of t, A (4.18), we obtain the derivative relation From the experience with RN black branes, discontinuous saturation is correlated with the near saturation behavior of τ (z f ) (and hence A(z f ) using (4.23)), thus we investigate the behavior near saturation, z f → 1. We emphasize that it is logically possible, to have more complicated behavior than we analyze below, and it would be interesting to understand whether the criterion that we give for discontinuous saturation is necessary or only sufficient. In the limit z f → 1 both the integrals and r(z f ) diverge, but τ, A have finite limits (4.19).
We first examine r(z f ): where a i are determined by the emblackening factor f (z) and d, and we will avoid writing the explicit expressions down to avoid clutter. We now focus on the integrand I τ near z = 1, where it diverges, and write down the potentially divergent terms:  Subtracting these terms from the integrand and adding them back, we obtain the integral in a . . . , (4.28) where in the second line we evaluated the remaining integral We put all this back together to obtain an expansion for τ (z f ). In d = 3 we get: where c 3 involves a 5 defined in (4.28), and hence requires numerical integration. The leading behavior is however determined by the coefficient of the δz f log δz f term: combining the Null Energy Condition (3.2) evaluated at z = 1 with f (1) < 0 implies that Near saturation τ (z f ) grows with δz f , thus saturation is discontinuous for any black brane. In the caption of Fig. 10 we give c 3 for the black branes analyzed there, and find perfect agreement with the numerical evaluation of (4.22).
In higher dimensions we get: where c d has a complicated expression that involves a 5 defined in (4.28), and hence requires numerical integration. Because the nonanalytic term with coefficient c d gives the leading correction, we are unable to determine analytically whether the saturation is continuous or not. If it is continuous, we near saturation we get a power law behavior This behavior is valid for 1 v B − τ 1. This is the same power law that the upper bound (1.11) or the combined bounds discussed in [1] give, but the overall coefficients do not agree. In [19,20] instead of this regime, they zoomed in on times 1 v B − τ β/R, and found a different power law in the Vaidya setup. For such time differences the details of the quench matter, as can also be seen from their analysis. Here instead we are focusing on universal features.

E. Early time growth and some worked out examples
The early time growth of entropy is significantly easier to work out that the late time behavior: we have to expand r(z) and the integrand I τ (z) in (4.22) for z = z HM − δz, integrate, and combine to get an expression for τ (z f ). We plug into (4.23), obtain A(z f ) as a power series, finally invert the relation τ (z f ) perturbatively, which gives It was suggested in [32] based on tensor network intuition that the early time expansion would take this form. This expansion could provide clues to how to improve the field theory bounds discussed in [1], which give exact linear growth at early times.
Next we investigate two limits that give some simplification. In Fig. 9 we observe that as q → 1 the A(τ ) becomes more and more linear. Motivated by this observation, we investigate this limit in some detail in an expansion in δq ≡ 1 − q. We obtain the simplest formulas in d = 3, and hence we will restrict to this case. Because in the extremal limit z HM → 1, it is easier to work in terms of a new bulk radial coordinate, z ≡ 1 + (z HM − 1)(1 − ξ). In the limit δq → 0 we find that the integrand I τ (ξ) is dominated by the small ξ region: showing the cuspy behavior seen on Fig. 9, but magnified. Top: d = 3 Schwarzschild black brane, which was the primary example discussed in this section. Note that the extremely small values on the axes is the reason why we were not able to see discontinuous saturation on Figs. 6 and 9. The curve plotted on the first plot is v B τ − 1 = − 1 2 δz f log δz f − 4.2598δz f . The zero of this curve is hence approximately at δz f = e −8.52 explaining the small scale of the graph. The rest of the curves are easy to obtain from this formula. There is a perfect agreement between the expansion and the numerical evaluation of the integrals (4.22). Middle: d = 3 RN black brane with q = 3/4, for which v B τ − 1 = − 13 2 δz f log δz f + 0.03822δz f + . . . . We see that the scale of the curves is a lot bigger than what we got for the Schwarzschild case, and the series expansion only matches the numerical data for small values of δz f . Bottom: d = 4 RN black branes with q = 0, 1/2, 2/3, 3/4 plotted with blue, orange, green, and red respectively. The values of c 4 defined in (4.31) are c 4 = −1.9955, −0.5164, 1.0005, 2.5126 respectively, and the sign of c 4 determines whether saturation is continuous or discontinuous: for q = 0, 1/2 we get continuous saturation, while the other two cases lead to discontinuity. The solid curves include the subleading term from (4.31) as well.
Using (4.23) we finally arrive at the parametric form of the entropy curve: This curve looks extremely straight, as seen in Fig. 11, and saturates at the value ζ f = 0.2337.
It has been recently observed that general relativity simplifies in the large d limit [41,42].
Motivated by this observation, we investigate the large d limit of the entropy curve for the Schwarzschild black brane. In this limit, the black brane becomes a membrane and z HM → 1, as in the q → 1 limit discussed above. Hence, it is convenient to change the bulk radial coordinate to z ≡ 1 + (z HM − 1)(1 − ξ), and expand in 1/d. A straightforward computation gives:

V. TWO-SIDED MUTUAL INFORMATION
As explained in the Introduction, we would like to determine the two-sided mutual information to probe the tsunami picture of [19,20]. We are looking for two-sided surfaces with R L , t L = 0 and R R , t R . We will calculate the mutual information between these two regions. As in the one-sided setup, the area and time of the extremal surface is dominated by the (doubled) scaling surface.
Here, because we do not have the constraint of Z 2 symmetry, we are allowed to boost this (doubled) scaling surface, which amounts to a shift in its v coordinate. In this section we will only discuss in detail the case, when the black brane under consideration gives continuous saturation in the one-sided setup. We will comment briefly on the discontinuous case, which is more complicated to analyze. The scaling surface has a size on z HM that we will keep denoting by ρ c , and we will use the rescaled field theory radial coordinate r defined below (4.11). The left and right boundaries are connected to the scaling surface at some r * ,L = R L /ρ c and r * ,R = R R /ρ c respectively. Depending on where these connections are compared to where the scaling surface touches z HM , we have three distinct possibilities. We explain these possibilities below, and illustrate them on Fig. 12.
(a) The midpoint is to the left of both shooting out points, see Fig. 12a. Using the quantities introduced in (4.18) the equations governing this situation are: where v shift quantifies the boost we made. We found the easiest to think about the left side by We introduce the natural rescaled variables where 0 ≤ I ≤ 2 is the two-sided mutual information in the focus of our interest. Using the relation 3) (5.1) can be rewritten in terms of the rescaled variables as: where we suppressed the L subscript on r * . Henceforth, we will always use (5.3) to eliminate r * ,R .
In the following we want to fix T and determine the function I(R). To do this we have to determine the relation r * (R) for fixed T . This is easily done numerically. There is one subtlety, the domain D of the function r * (R) is bounded. The function τ (r * ) is only defined for r * ∈ [1, ∞) 15 and takes values in [0, 1/v B ). First, let us plug in r * = 1 into the first equation of (5.4), and define R 1,T through it: Second, taking r * → ∞ in the first equation of (5.4) gives From (5.5), (5.6) and the blue curves on Fig. 13, we conclude that the domain of the function At the endpoints of this domain r * = 1 and r * = ∞ respectively. What we said above remains true even in the case of discontinuous saturation.
(b) The next case we consider is when the midpoint of the scaling surface is between the left 15 We note that r * ≥ 1 from its definition r * = R/ρc.
Again, the domain of r * (R) is a bit subtle. Let us define R 2,T as the root of the equation The equation has a solution only if v B T < 1. By a similar logic as above, and analyzing the orange curves on Fig. 13, we conclude that the domain of the function r * (R) is .
The last equation is not true in the case of discontinuous saturation, and it has important consequences.
(c) Finally, the midpoint of the scaling surface can be to the right of both shooting out point, see Fig. 12c. This immediately implies that R < 1 the same way as the setup in (a) implied that R > 1. Again there are only some sign changes: (5.11) The domain of r * (R) can be read off from the green curve on Fig. 13, and we get There are two notable points: the values of R where the mutual information starts to be nonzero, and where it saturates Multiplying the second equation by R L and recalling the definitions of R, T given in (5.2), we get This result was quoted in the Introduction (1.10). In the discontinuous case the mutual information becomes nonzero at an earlier time, but saturates, where (5.15) says so, see Fig. 15. 16 This is quite interesting, because in the discontinuous case we do not know how to extract the value of v B from the entropy growth in a quench. 17 The two-sided mutual information can reveal the value of v B .
Similar comments apply to spin chains [1].
We now comment on how these results refine the "entanglement tsunami" picture proposed in [19,20], see also [1] for further comments. It was suggested that at early times only degrees of freedom behind the tsunami wavefront spreading with the entanglement velocity v E would become entangled. This picture would imply that for fixed R L , t the radius R R at which the mutual information saturates is R L + v E t, which is smaller than the result (5.15). 18   . Two-sided mutual information for the d = 3 RN black brane with q = 3/4, which does not give continuous saturation in the one-sided setup. On Fig. 9 it corresponds to the red curve with the largest cusp. Besides (5.13), we have plotted I for T = 4, and marked the points (5.14) with dashed red lines. While I saturates at the point R (saturates) = v B T + 1, it starts is nonzero earlier than what we get for black branes with continuous saturation. Also, while the I curve starts with zero slope on Fig. 14, here the slope is nonzero.
These give the coordinates for points on the R side. To get the points on the L side, we add iπ to the time coordinate. This switches the sign of T 1 and X 1 . The length of a geodesic connecting two points is given by We are interested in the mutual information between intervals on the L and R side. We want to connect the endpoints of intervals separated in time and centered at the origin on the two sides, thus we are looking for a geodesic connecting the point (t L = 0, R L ) on the L side, with (t R , R R ) on the R side. Note that the length of the intervals is 2R L,R . In the computation below we assume that R R > R L , the R L < R R is easily obtained from the result.
Subtracting a universal log(z cutoff ) term from the length, we get To compute the mutual information, we also need to know the length connecting the ends of the intervals to themselves. This is e.g. d LL = log [2 cosh(2R L ) − 2]. The mutual information is then

(5.20)
This formula is valid as long as the answer is positive. If it is negative, a disconnected surface dominates and the correct answer is zero. The rescaled mutual information (5.2) is plotted in Fig. 16.

A. Scaling limit for arbitrary shapes
In this section we will restrict to d = 3 in order to alleviate the notation; all our results generalize straightforwardly to higher dimensions. It will be advantageous to consider the infalling time v and the angular coordinate θ as independent variables. The extremal surface is then given by two functions ρ(v, θ), z(v, θ), and we can write down a complicated action in terms of these variables. Instead of doing this, using our experience with spherical surfaces, we do a rescaling: where R is an arbitrary scale associated to the entangling surface Σ. We will consider the end of the world brane setup only, but the arguments at the end of Sec. IV C imply that we get the same extensive part of the entropy for any quench for arbitrary shapes. The scaling solution is determined by an image on the brane r im (θ) and the boundary shape r bdy (θ). The boundary conditions on r(V, θ) hence are It may be instructive to compare to our treatment to the discussion in Sec. IV C of the spherical case. The approach used in this section could also be used to derive the results for the spherical case in a streamlined fashion. In Sec. IV C we scaled with ρ c instead of R, which meant that we fixed r im (θ) = 1, obtained one scaling solution that we connected to the boundary at some r * , finally we had to scale the answer with the ratio r * = R/ρ c . Here we find it more convenient to scale with R, which has the consequence of r bdy (θ) staying fixed, but the image on the brane r im (θ) changing as a function of time.
In terms of the scaled variables (6.1) the action is: Note that the large R limit brings major simplification: no derivatives of Z appear in the action.
This implies that the equation of motion for Z is algebraic. We write it in a suggestive form: The other equation of motion is more complicated: Using (6.4) we see that some terms in this equation simplify: and we arrive at the following final form of the equation: (6.7) Note that we have eliminated ∂ V r. One can reduce the number of equations from two to one, by expressing Z as a function of r and its derivatives from (6.4), but the resulting equation is very complicated even for the Schwarzschild black brane. Before making some comments about the behavior of this complicated equation, we analyze some simple cases.

B. Simple cases
Let us first take the example of a spherical surface. (6.4) and (6.7) simplify drastically to (6.8) Let us see how these agree with the result of our analysis in Sec. IV. There we were using r as the independent coordinate, so we have to account for the change of independent coordinate. The first equation above is the same as the first equation of (4.17), if we use that dr dV = 1/ dV dr . The second equation above is the derivative of the first, if we use the relation between r and Z (4.21).
Second, we can find a special solution to the equations. Setting Z = 1 solves (6.7), and (6.9) becomes ∂ V r = v B 1 + (∂ θ r) 2 r 2 , (6.9) a true tsunami equation, with v ts = v B (2.16). From our experience with the spherical case in Sec. IV we know that this HRT surface is only relevant, if saturation is continuous. We do not know, if saturation can be continuous for any shape other than the sphere. If it can be, then the saturation time is immediately seen to be t S = R insc /v B , the time it takes the tsunami (propagating from Σ inwards) to cover the interior of the region. This is the same time it takes for the static surface to become available in the Vaidya setup, as discussed around (2.19). If the saturation is discontinuous, this solution still provides a lower bound on the saturation time: (6.10) The above argument is a second holographic proof (the first being (2.19)) of the same inequality proven in field theory in [1].
C. Comments on the general case which agrees with the result in Sec. II D that for the strip the scaling surface stays on z HM forever.
Let us fix r im (θ) for the following argument. From (6.12) it follows that for angles θ for which r im (θ) is more curved will depart form z HM faster, while flatter parts hang around longer. At least initially, the surface becomes more wiggly. Eventually, we expect that every bit of the surface reaches the vicinity of the horizon, where we can perform another expansion.
For large V we have only found one admissible behavior, which implies that the surface approaches the horizon. 19 The form of the near horizon expansion is: where a(θ), b(θ) are arbitrary functions, while r i (θ), z i (θ) are functions of a(θ), b(θ) and their derivatives that we have determined, but do not write down here. As usual in an asymptotic analysis, the freedom of choosing a(θ), b(θ) can be used to obtain a regular solution near the brane at z HM with an arbitrary r im (θ). Note that for fixed r im (θ) the surface becomes spherical for large V according to (6.13). It would be a very interesting to understand the details of how the HRT surface interpolates between the two regimes (6.11) and (6.13). A straightforward numerical solution of (6.4) and (6.5) is prevented by the formation of singularities for intermediate times. We understand the reason for the formation of cusps from two perspectives: (6.11) leads to a wigglier wave front as V grows, while evolving backwards in V in the near the horizon where (6.13) holds, should also lead to singularities according to what was explained around (6.4) and in Fig. 17.
For the above discussion we have fixed r im (θ), but we are actually interested in keeping r bdy (θ) fixed, which requires adjusting r im (θ) as time evolves. For early times, the expansion (6.11) relates the two, and this is explored further in Sec. VI E. We speculate that for intermediate times we have to start with an almost spherical r im (θ) to reach a more deformed r bdy (θ). For late times, presumably we start with some (possibly cuspy) r im (θ), which under evolution in V becomes more wiggly and can develop cusps, finally it gets smoothed out to end as r bdy (θ) at V = V f . It would be important to understand, if this is indeed what happens. While our understanding is incomplete, we have enough control to prove an important result below. 19 The scaling solution for the strip is an exception to the story sketched here, as it lies on zHM for all V .
−H(V f ) ≤ v E dθ r 2 + (∂ θ r) 2 = v E area[r bdy (θ)] . (6.19) Undoing the rescalings (6.1), we get a bound on the rate of growth where A Σ is the area of the entangling surface. Recall that s th = 1/4G N is what converts area to entropy. Note that at late times, when Z(V f ) → 1, this is a huge overestimate, as the prefactor in (6.18), −f (Z)/Z 4 → 0.
There is an interesting result that we can derive along the same lines. Unfortunately, we have not found a clear field theory interpretation of it. 21 Note that in (6.6) we have obtained a simple expression for Q valid on-shell that we have not used above. Using that expression we get: = v E area[r im (θ)] , (6.22) where we used that Z(V i = 0, θ) = z HM according to (6.11), and that (2.23) relates f (z HM ) to f (z HM ). Note that the area of the image on the brane is measured in the field theory coordinates x defined in (2.1). This is also a somewhat different proof of the rate of growth bound (6.20), as area[r im (θ)] < area[r bdy (θ)]. The higher dimensional generalization is straightforward, and gives the area of the image on the brane on the right-hand side of (6.22).
In the spherical case, by manipulating integrals, we found a mysterious relation (4.25). Now we 21 Conversely, we do have not found a holographic proof of the inequality (1.11) proven in [1]. This situation may be analogous to the case of the monotonicity of renormalized entanglement entropy [43], where there exists a field theory proof [44], but no holographic argument.