{'fraction_non_alphanumeric': 0.03651770098305289, 'fraction_numerical': 0.011097038779325734, 'mean_word_length': 4.436745177332095, 'pattern_counts': {'":': 0, '<': 19, '<?xml version=': 0, '>': 15, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Holographic Thermalization

Using the AdS/CFT correspondence, we probe the scale-dependence of thermalization in strongly coupled field theories following a quench, via calculations of two-point functions, Wilson loops and entanglement entropy in d=2,3,4. In the saddlepoint approximation these probes are computed in AdS space in terms of invariant geometric objects - geodesics, minimal surfaces and minimal volumes. Our calculations for two-dimensional field theories are analytical. In our strongly coupled setting, all probes in all dimensions share certain universal features in their thermalization: (1) a slight delay in the onset of thermalization, (2) an apparent non-analyticity at the endpoint of thermalization, (3) top-down thermalization where the UV thermalizes first. For homogeneous initial conditions the entanglement entropy thermalizes slowest, and sets a timescale for equilibration that saturates a causality bound over the range of scales studied. The growth rate of entanglement entropy density is nearly volume-independent for small volumes, but slows for larger volumes.


IV. Discussion
The observed nearly inviscid hydrodynamic expansion of the hot QCD matter produced in nuclear collisions at the Relativistic Heavy Ion Collider (RHIC) indicates that matter produced in these nuclear reactions is strongly coupled [1][2][3]. The inability of perturbation theory to account for this phenomenon has motivated studies of nonequilibrium dynamics in analytically tractable, strongly coupled gauge theories. The prototype of such theories is the maximally supersymmetric Yang-Mills theory in the large-N c , large 't Hooft coupling limit, which is holographically dual to the classical limit of a superstring theory on the Anti-de Sitter space AdS 5 ×S 5 background [4]. Hydrodynamical evolution with a minimal shear viscosity emerges naturally as late-time behavior in the longitudinal expansion of a thermal gauge plasma in this model [5,6]. In the dual supergravity theory, the thermal state of the gauge theory is represented by a black brane in the asymptotic AdS 5 space, and the near-equilibrium dynamics of the gauge theory giving rise to hydrodynamic behavior is mapped onto the dynamics of perturbations of the AdS 5 -black brane metric. A nice recent review on the holographic study of hot QCD matter, focused on equilibrium and near-equilibrium aspects, is [7].
While the near-equilibrium dynamics of the strongly coupled super-Yang-Mills (SYM) theory is well studied, the process of thermalization itself is still poorly understood. The RHIC data demand that the time scale for equilibration of matter is considerably shorter than expected in the framework of perturbative approaches to thermalization [8,9]. The dearth of other nonperturbative tools for the description of the short pre-hydrodynamic stage of the nuclear collision motivates the use of the AdS/CFT correspondence to study thermalization of strongly coupled plasmas. Phenomenologically, the central rapidity region of a high energy collision turns out to be boost invariant to a good approximation, so it is natural to consider boost invariant configurations [6,[10][11][12]. One interesting question is how thermal equilibrium is reached in such systems. This is the question we will address in the present paper, albeit in the case of translationally invariant plasmas.
We investigate the thermalization process not only in the strongly coupled (3+1) = 4−dimensional super-Yang-Mills theory, but also in the analogous lower (2 and 3) dimensional field theories. This has several motivations. The first reason is that the 2−dimensional version, which is dual to classical (super-)gravity on AdS 3 space, admits analytical solutions and thus allows us to explore thermalization over a wide range of parameters. A second motivation is that some version of the 2−dimensional dual theory, (1 + 1) dimensional super-Yang-Mills theory, can be formulated on a lattice and solved nonperturbatively [13,14]. Finally, the comparison of theories in different dimensions makes it easier to distinguish generic aspects of thermalization of strongly coupled quantum field theories from aspects that are special features of 4−dimensional field theories.
B. Holographic models of thermalization

Thermalization scenarios
What is appropriately called the thermalization time may not only depend on the probe of the state of the field but also on the initial field configuration that evolves toward thermal equilibrium. Here one can distinguish between two broad classes of scenarios. The first class considers small perturbations around thermal equilibrium and asks how equilibrium is reached. For a gauge theory with a holographic dual, this approach amounts to studying the decay of small perturbations of an AdS black brane geometry. This is conveniently studied in terms of the quasi-normal modes of fields propagating in the black brane background [15]. The imaginary parts of the eigenvalues of these modes describe the thermal relaxation rates of various excitations, such as anisotropic perturbations of the stress-energy tensor. In the context of relativistic heavy-ion collisions, this approach is relevant for the study of viscous corrections to hydrodynamics and other transport processes in the presence of a thermal gauge theory plasma (see e.g. [16] for a review).
The second class of scenarios considers the thermalization of an initial field configuration, which is generally far from equilibrium. Through the AdS/CFT correspondence, the approach to thermal equilibrium in the boundary gauge theory is related to the process of black hole formation in the bulk. In early work, before the motivation from RHIC experiments, the goal was to understand how black hole formation from a gravitational collapse of a shell of matter would be encoded in the field theory in the boundary [17][18][19][20]. Also in more recent work, the initial conditions have usually been defined directly in the dual gravity theory -in most, if not all, of these cases the precise form of the initial condition is not known in terms of gauge theory excitations. Examples of such configurations are colliding gravitational shock waves [21], sheets of fundamental strings that may provide an AdS model of densely packed flux tubes [22], and sudden perturbations of the metric near the boundary that propagate into the bulk [11,23]. Motivated by the properties of the initial state of a relativistic heavy-ion collision, which contains two highly energetic nuclei, these scenarios have in common that they inject energy into the AdS geometry at high momentum or short distance scales.
Some of these recent studies of gravitational collapse in AdS 5 started from translation invariant, but locally anisotropic perturbations of the metric near the boundary and followed their propagation into the bulk, ultimately resulting in the formation of an event horizon [11]. The initial anisotropy of the metric is dissipated by the black brane, asymptotically resulting in viscous, boost-invariant scaling hydrodynamics on the boundary [12,24]. Gravitational collapse in AdS d+1 induced by a scalar field perturbation was studied in [23], most explicitly for the case of AdS 4 , where the bulk equations of motion for gravity coupled to a massless scalar field were solved perturbatively for a small amplitude scalar perturbation that propagates in from the AdS boundary. The perturbation was taken to be translationally invariant along the spatial directions of the boundary and to vanish outside of a short time interval. An interesting technical result is that for this translationally invariant collapse in AdS 4 , the metric outside the infalling shell of matter coincides with a black brane metric to first non-trivial order in the perturbation [23]. This has the consequence that expectation values of local observables thermalize essentially instantaneously in the thermally quenched field theory dual to the infalling shell background.
Similarly to the approach of [23], in [25] the equilibrating field configuration was modeled by a homogeneous, infalling thin mass shell. The motivation for choosing such a configuration was that it naturally arises in the AdS dual description of the asymptotic limit of a transversely extended ensemble of flux tubes [66]. The latter represent the energy density deposited by the colliding leading gauge charges (i.e. the valence quarks of the colliding nuclei), and the AdS description involves a sufficiently simple geometry to admit analytic or semi-analytic solutions. Our model does not allow us to study the approach to local isotropy and hydrodynamical behavior, but it enables us to investigate the approach of the field to a thermal configuration in momentum space and to measure the growth of its entanglement entropy, as discussed in Sections II C and III C. In this sense, our present study is complementary to those mentioned above which focus on the approach to the hydrodynamical limit. Note that the work [25] considered a quasi-static approximation with metrics that correspond to stationary "snapshots" of the dynamic geometry containing the infalling mass shell. Here we study both the dynamic and the quasi-stationary geometries and show that there are some significant differences in the pictures that they yield of thermalization. Various aspects of the dynamic case of a thin infalling shell have been investigated in [26][27][28]. The effect of a dynamical boundary condition has also been studied in the background of a moving mirror [29]. New techniques for the evaluation of various Green functions in non-equilibrium backgrounds have been developed in [30], in addition to the complex contour approach by [31]. The counterpart to Chandrasekhar limit for a gravitational collapse of a degenerate star to a black hole in AdS/CFT has been investigated in [32]. Finally, thermalization after a quantum quench in gauge theory has also been modeled by a sudden change in couplings or other background fields than the bulk metric [33,34]. In probe brane approach such a change then can manifest as time-dependence in the induced metric on the brane.
Concretely, we focus on the (d + 1)-dimensional infalling shell geometry described in Poincaré coordinates by the Vaidya metric where v labels ingoing null trajectories and we have set the AdS radius equal to 1. The boundary space-time is located at z = 0 and x = (x 1 , . . . , x d−1 ) correspond to the spatial coordinates on the boundary. For constant m(v) = M , the coordinate transformation brings the metric (1) in the form On the boundary, the coordinates v and t coincide. We take the mass function of the infalling shell to be where v 0 parametrizes the thickness of the shell falling along v = 0. We will often be interested in the zero thickness limit v 0 → 0. The Vaidya metric describes a shell composed of tensionless null dust, which represents an analytical simplification of the tensionful shell models studied in [25] and of the model of [23]. The geometry outside the infalling shell (which corresponds to v > 0 in the zero thickness limit) is identical to that of an AdS black brane geometry with Hawking temperature T = dM 1/d /4π, while the geometry inside the shell (v < 0 when v 0 → 0) is the same as that of pure AdS geometry. The causal structure of the Vaidya space-time is shown in Fig. 1.
FIG. 1: The causal structure of the Vaidya spacetime (in the v0 → 0 limit) shown in the Poincaré patch of AdS space. In this presentation, the asymptotic boundary (vertical line on the right hand side) is planar, and the null lines on the left hand side of the diagram represent the Poincaré horizon.

Probes of thermalization
In order to probe the dynamics of the thermalization process, one can study a variety of observables of the boundary gauge theory. Expectation values of local gauge-invariant operators, such as the energy-momentum tensor and its derivatives provide valuable information about the applicability of viscous hydrodynamics, but they cannot be used to explore deviations from thermal equilibrium in detail. Nonlocal observables such as pair correlation functions, Wilson loop expectation values, and entanglement entropy provide much more information about progress towards thermalization. (They are also relevant to the physics probed in relativistic heavy ion collisions, e. g. through the jet quenching parameterq [35,36] and the color screening length.) To illustrate, consider spatially homogeneous states of a weakly interacting massless scalar field. The energymomentum tensor can be expressed as where n(k) = a † k a k denotes the occupation number of a momentum mode of the scalar field. (There is no spatial dependence of the stress tensor because the field configuration is homogeneous here.) It is obvious from (5) that this observable contains only limited information about the particle distribution. In particular, T µν cannot inform us whether thermal equilibrium has been reached; it only tells us whether the pressure is locally isotropic. This is not even a sufficient condition for the matching of the field theory to hydrodynamics, because the equation of state relies on the assumption of thermal occupation numbers for the various modes. It is also insufficient to answer many other questions, e.g. about the expected spectrum of radiation emitted by collisions among the scalar particles. The equal-time two-point function, on the other hand, is given by which allows us to extract detailed information about the particle distribution and thus probe its closeness to a thermal distribution in detail. Note that even the two-point function is insufficient to probe for the phase relationship between the occupation amplitudes of different field modes. This information can, however, be obtained by an analysis of the four-point function of the field.
In an interacting quantum field theory, the two-point function is determined not only by the mode occupation numbers n(k) but also by the spectral density function σ(k 0 , k). For the free field, σ(k) = 2πδ(k 2 − m 2 ), but more generally the spectral function is determined by the exact self-energy Π(k 0 , k) as follows: Since the self-energy is a function of the temperature, so is the spectral density function. For a weakly coupled field, the medium dependence is reflected, e.g. , in medium modifications of particle masses and widths. In the strongly coupled gauge theory, the vacuum spectral function has the generic form [37] σ where the exponent ν depends on the dimension of the considered field. This becomes a sigmoidal function at finite temperature. The mode occupation number n(k) and the spectral function σ(k 0 , k) combine multiplicatively to yield the true, physical excitation probability of a momentum mode k. The time evolution of the equal-time two-point function probes the approach of this excitation probability to thermal equilibrium. In the AdS/CFT correspondence there is a geometric intuition for why expectation values of local gauge-invariant operators are insensitive to details of the progress towards thermalization -they are only sensitive to phenomena near the AdS boundary. Thus they do not probe the details of phenomena occurring near the thermal scale. By contrast, nonlocal objects, such as two-point correlators of gauge-invariant operators and expectation values of Wilson loops, are dual to AdS quantities that probe deeper into the space-time and further away from the boundary. For example, at strong 't Hooft coupling, it is possible to approximate the path integral for the connected two-point correlator as a sum over all possible AdS geodesics connecting the two endpoints, which are placed on the AdS boundary [38]. (For certain subtleties regarding this statement, see [39][40][41].) The geodesics probe the interior of AdS space, which is dual to the statement that they are sensitive to a wide range of energy scales in the boundary field theory. As we will see, Wilson loop expectation values [42] and entanglement entropy [26,43,44] are also related to minimal lengths, surfaces and volumes of various kinds in AdS. These also extend into the bulk of AdS and hence probe a range of energy scales.

C. Overview of this work
We consider two-, three-, and four-dimensional field theories dual to gravity in asymptotically AdS 3 , AdS 4 , and AdS 5 spacetimes, respectively. The motivation for studying all these cases (rather than only the case of four-dimensional field theories potentially relevant for heavy-ion collisions) is that we are interested in generic conclusions on strongly coupled field theories with gravity duals. Lower dimensional theories are technically much simpler: analytic computations are possible in AdS 3 . Furthermore, our basic methods are similar in all dimensions -we calculate field theory observable that can be related to invariant geometric objects -minimal lengths, surfaces and volumes -in AdS spacetimes of different dimension. Our results show that many conclusions are not sensitive to the dimension.
Sec. II discusses three measures of thermalization -two-point functions, Wilson loop expectation values, and entanglement entropy. We compute these observables at thermal equilibrium, and, in later sections, use deviations from these results as probes of incomplete thermalization. Sec. III relates thermal quenches in two-, three-and fourdimensional theories, to collapsing shells of matter in AdS 3 , AdS 4 and AdS 5 that form black branes. The black brane formation represents the process of equilibration and thermalization of energy injected in the dual field theory. In these dynamical backgrounds we compute nonlocal observables as a function of time, and track the rate at which thermalization progresses for each observable at different spatial scales in the field theory. Sec. IV concludes the paper by summarizing our findings. An Appendix describes how the results would differ in a quasistatic approximation. A brief account of this work has appeared in [45].
One broad finding is that these models exhibit a "top-down" thermalization mechanism, see also [25]. We would like to stress that this is very much in contrast with the standard "bottom-up" paradigm [8] based on perturbative gauge theory. In the "bottom-up" scenario, hard quanta of the gauge field do not equilibrate directly by randomizing their momenta in two-body collisions, but they do so by radiating softer quanta, which gradually fill up the thermal phase space and equilibrate by collisions among themselves. In other words, infrared modes are the first to thermalize, and the thermalization proceeds gradually ("from the bottom up") to more energetic modes of the gauge field. The bottom-up scenario is closely linked to the infrared divergence of the splitting functions of gauge bosons and fermions in the perturbative gauge theory. This contrasts with the "democratic" splitting properties of excitations in the strongly coupled SYM theory, which favor an approximately equal sharing of energy and momentum [46,47]. We should therefore expect that the thermalization process proceeds quite differently in the strongly coupled gauge theory. The phase space for radiation and, at strong coupling, also the spectral weight (8) is largest for the most energetic quanta. They then divide their energy and momentum approximately equally among their siblings, leading to a rapid cascading down to a thermal distribution. This consideration suggests that the thermalization process more closely resembles a "top-down" scenario, in which energetic gauge field modes equilibrate first and soft modes last. Further broad conclusions appear in Sec. IV.
It is appropriate here to comment on the differences and overlap between the current work and the recent related work in the literature. In [26,27], geodesics (or equivalently, entanglement entropy) in AdS 3 Vaidya were studied numerically for a thin but finite shell. In contrast, we also study geodesics in AdS 3 Vaidya analytically for an infinitesimally thin shell. In [28], Wilson loops (or equivalently, entanglement entropy) were studied in AdS 4 Vaidya, for the case of circular and rectangular loops. We do not only confirm their results but also study geodesics in AdS 4 Vaidya, which have not been closely investigated in the context of holographic thermalization as far as we know. Moreover, we study AdS 5 Vaidya geometries for the first time by computing geodesics, circular and rectangular Wilson lines, and entanglement entropy for spherical volumes on the boundary.

II. PROBES OF THERMALITY
We are interested in processes where energy is injected into a strongly coupled field theory and equilibrates over time. Thus we need probes that take instantaneous snapshots of the state of the system which can be compared to the thermally equilibrated state. As discussed above, expectation values of local gauge-invariant operators are insufficient. Thus, in this section we describe three nonlocal probes -two point functions, Wilson loop expectation values, and entanglement entropy -which can be elegantly computed using geometric techniques in the AdS/CFT correspondence. We will show how these quantities computed in strongly coupled theory at thermal equilibrium differ from the corresponding quantities in the vacuum. In Sec. III we will use the same quantities as probes of thermalization following a quench, by comparing the instantaneous values of the probes with the values at thermal equilibrium.

A. Equal-time two point functions
We begin our exploration with the equal-time two-point functions, which can be analytically calculated in twodimensional conformal field theories at thermal equilibrium. We first remind the reader how the thermal occupation probabilities of momentum modes enter into the two-point function and how they are related to the spatial form of the equal-time Wightman function. We then explain how we generalize this approach to higher dimensional theories, where we can easily calculate the Wightman function in the saddle-point approximation.

Two-dimensional free scalar field in thermal equilibrium
We would like to understand how we can extract information about the thermalization process from the time dependence of the Wightman function where O is a local operator of dimension ∆. As a warmup exercise we consider the free massless scalar field in two dimensions with mode expansion We defined the operators a + k as dimensionless, which implies the commutation relation [a k + , a † k + ] = k + δ(k + − k + ), etc. Employing the free spectral function σ 0 (k) = 2πδ(k 2 ) = 2πδ(k + k − ), this can be written as Consider the dimension one operators ∂ + φ or ∂ − φ. The vacuum two point function of ∂ ± φ is time and space translation invariant, so we consider For nonzero temperature, the spectral function of the free field remains unchanged, and we find where n(k) = (e βk − 1) −1 is the standard (Bose) thermal occupation number. This can be rewritten as This is exactly the result expected from conformal invariance. The finite temperature computation is a computation on a Euclidean cylinder, which can be mapped to the complex plane using the exponential map, and we already know the answer on the plane. Undoing the coordinate transformation leads to the above two-point function.

Strongly interacting scalar field theory in equilibrium from AdS/CFT
For an interacting field theory the spectral function changes with temperature. As a preparatory example of an interacting two-dimensional thermal field theory defined by a gauge-gravity duality we consider the AdS 3 case, i.e. the BTZ black hole, and make use of the results presented in [48]. Writing ω = (k + + k − )/ √ 2 and k = (k + − k − )/ √ 2, the thermal Wightman function for a scalar field in such a theory is [48]: where G R is the retarded Green's function and −2 Im G R yields the spectral function. The precise form of the thermal retarded Green function depends on the dimension ∆ of the field in the operator. For ∆ = 1, we have [48] G R (ω, k) = 1 2π Using Im ψ( 1 2 + iy) = π 2 tanh πy, gives the spectral density of states as What can we learn from the equal-time Wightman function? One strategy is to use the vacuum subtracted Wightman function which at thermal equilibrium depends only on ∆t = t − t Using the Bose occupation number density n B (ω) = (e ω/T − 1) −1 , the equal-time Wightman function (∆t = 0) is then Inserting the explicit expressions, we obtain after some manipulations Introducing the notation cosh a = coth(k/2T ), and using the integral identity the final result reads Having derived the functional form of the equal-time Wightman function for a thermally equilibrated state in AdS 3 , we ask what we could learn from the same function at non-equilibrium background. In that case time translation invariance is lost. As in the above, we could define a vacuum subtracted Wightman function G > sub (t, k) now depending on the equal-time t = t , and expect that thermalization implies G > sub (t, k) → G > sub (k; T ). This gives us a way to estimate equilibration times for different momentum scales k, by analyzing how this limiting behavior is approached. The result (22), that the thermal contribution to the equal-time Wightman function of the two-dimensional field theory is proportional to the Bose distribution, justifies a more complete investigation of the equilibration properties of the Wightman function in two and more dimensions. For this purpose, we now turn to the methods made available by the holographic duality between strongly coupled supersymmetric gauge theories and gravity theories. We will also analyze thermalization times at different scales, in the spirit of the above motivation, but the details will differ somewhat.

The geodesic approximation
While the approach described above would give detailed information regarding scale-dependent thermalization, there is a technical challenge for strongly coupled field theories, even if they have a dual description in terms of gravity in an AdS space. The easiest case is two-dimensional field theory with a three-dimensional gravity dual. In this case, the thermal Wightman function in momentum space is known in closed form for operators of general conformal dimension ∆ (∆ = 1 was analyzed above) [48]. However, the integral with respect to ω that is needed to compute the equal-time Green function is sensitive to the ultraviolet completion of the theory and thus requires careful regulation. Thus, it would be more convenient to directly compute the equal-time thermal Wightman function from the AdS 3 theory. The next simplest case involves three-dimensional field theories with an AdS 4 dual. In this case there is not even a closed form expression for the thermal Wightman function. Thus we cannot pursue the approach discussed above analytically. Furthermore, we are really interested in computing the equal-time Wightman function in a non-equilibrium setting in which matter collapsing in AdS to form a black hole models the thermalization of energy injected into the strongly interacting field theory. In this situation, there are no known analytic solutions for the Green functions in any dimension, and we have to resort to approximation schemes and/or numerical analysis.
An insightful approach is to probe the non-equilibrium state of the strongly coupled field theory with a very heavy operator whose Wightman function can be approximated in terms of AdS geodesics as described below. According to the AdS/CFT dictionary, a massive scalar field ϕ(z, t, x) with mass m in (d + 1) dimensions is dual to an operator [4]. In coordinates where the boundary of AdS is at z = 0, the scalar field ϕ behaves like where the ellipsis corresponds to subleading terms in the z-expansion. In CFT language we interpret the nonnormalizable mode φ (0) (t, x) as the source to which O(t, x) couples and the normalizable mode φ (1) (t, x) as the vacuum expectation value of the renormalized operator O ren (t, x) [49], The bare operator O(t, x) and the renormalized operator O ren (t, x) are related as [67] where z 0 represents an IR bulk cut-off. The two-point function O(t , x)O(t, x ) is often computed using the on-shell (supergravity) action for ϕ, when the action is explicitly given and solutions to the wave equation can be readily computed. However, for our purposes it is more convenient to follow [38] and to observe that the equal-time Green function can be computed via a path integral as The first expression sums over all paths that begin and end at the boundary points (t, x) and (t, x ), with L(P) being the proper length of the path (L(P) is imaginary for space-like trajectories). The second expression is a saddlepoint approximation to the path sum as a sum over geodesics [38] (here L is the real length of the geodesic between the boundary points). The latter approximation is effective when the probe operator is heavy so that ∆ 1. It is easy, to check, for example, that in the zero temperature theory modeled by empty AdS space, this formula gives the correct conformally invariant equal-time two point function of the renormalized operator [38] Some care is necessary in evaluating Lorentzian correlators using the geodesic approximation, because careful consideration of steepest descent contours of integration is generally necessary [39][40][41].
The geodesic length L diverges due to contributions near the AdS boundary. Therefore, we define a renormalized length δL ≡ L + 2 ln(z 0 /2), in terms of the cut-off z 0 , by removing the divergent part of the geodesic length in pure AdS (see Eq. (51) below). The renormalized equal-time two-point function is The geodesic approach to computing equal-time correlators gives a clear intuition for how thermalization proceeds in a strongly coupled theory with an AdS dual (Fig. 2). Consider injecting energy homogeneously into the field theory -we model this in the AdS theory by dropping in a shell of matter from infinity. As the shell progresses inward into AdS space, the spacetime outside the shell will be well described by the AdS-black brane metric (except within the shell itself), while the metric inside the shell will be well described by the empty AdS metric. Now consider the two-point function calculation using geodesics. If the two points are close together on the boundary, the associated geodesic will not penetrate very much into the bulk space, and thus will not cross the matter shell. Thus, at these separations the correlation function should look thermal, because a thermal field theory is modeled by a black brane background in AdS. At larger boundary separations, the associated geodesic will penetrate the shell and be "refracted" by it (Fig. 2) leading to deviations from thermality. As time passes, the shell will penetrate deeper into the bulk and thus an ever larger range of spatial scales in the field theory will have associated geodesics that do not penetrate the shell and hence have thermal correlators. Thus we can come to a qualitative conclusion that thermalization proceeds top-down in this setup -ultraviolet, i.e. short distance, degrees of freedom thermalize first. In subsequent sections we will calculate the rates at which different spatial scales thermalize.

Two-point function for two-dimensional field theories in equilibrium: analytic computation
In order to gain intuition for the two-point function in the geodesic approximation, we return to the case of twodimensional boundary field theories. In order to obtain the thermal boundary-to-boundary Wightman function, we need to study space-like geodesics in the thermal black brane geometry in three dimensions. The metric is obtained by setting m(v) = M in (1), and is given by This is a black brane geometry with a horizon at r H ≡ √ M [68]. The Vaidya geometry in three dimensions with a constant m(v) = M can be put in this form. Parametrizing the geodesic by the geodesic length λ, the geodesic equations in the geometry (29) are where E, J are constants and˙≡ d/dλ. Plugging (30) and (31) into (32) and solving forṙ, Integrating this, we obtain or, equivalently, Here, λ 0 is an integration constant which we will henceforth set to zero, λ 0 = 0, by absorbing it into the definition of the parameter λ. Also, we defined the combinations The geodesic reaches the boundary r = ∞ (or z = 0) as λ → +∞. As we can see from (34) and (35), the signs of A ± , B ± determine the behavior of the geodesic. The relation between the values of the parameters E, J and the signs of A ± , B ± is shown in Fig. 3. Note that not all possible combinations of signs occur. We discuss different cases (i)-(iv) shown in Fig. 4 in turn below.
(i) If A ± , B ± > 0, the geodesic is entirely outside the horizon. r(λ) takes its minimum value r = r H [1 + (A 1/2 One sample profile of this case is shown in Fig. 4(i). (ii) If A + A − < 0 and B + B − < 0, the geodesic crosses the horizon r = r H and reaches the singularity r = 0. Since in this caseṙ (recall that we have set λ 0 = 0), r(λ) is a monotonically increasing function of λ and the geodesic crosses the horizon only once. One sample plot of this behavior is presented in Fig. 4(ii).
(iii) If A ± < 0 and B ± > 0, the geodesic crosses the horizon twice at but does not reach the singularity. One sample plot of this behavior is presented in Fig. 4(iii).
(iv) If A ± , B ± < 0, the geodesic crosses the horizon twice at (38) and furthermore hits the singularity twice at r is pure imaginary between λ O+ and λ O− . One sample plot of this behavior is presented in Fig. 4(iv). If one is interested in geodesics that start and end on the boundary of the BTZ black brane geometry (29), the only relevant cases are (i) and (iii). However, if one is interested in geodesics in the Vaidya geometry (1) for which the BTZ geometry (29) is only a part of the entire spacetime, then one should also consider cases (ii) and (iv) since they can appear as a part of geodesics in the full spacetime. r Substituting (34) into (30) and (31) and integrating the two equations, we find the expressions for t and x: where t 0 , x 0 are constants of integration corresponding the values of t, x at λ = +∞. Comparing (40) with (35), it can be seen that the quantity inside the absolute value of (40) vanishes or diverges at the horizon r = r H and thus t(λ) diverges there. This is due to the well-known fact that the Schwarzschild time t is not well-defined at the horizon. Note that t(λ) → ∓∞ as λ → λ H± and therefore λ H+ (λ H− ) corresponds to the past (future) horizon.
On the other hand, the Eddington-Finkelstein coordinate v defined by (see (29)) is well-defined across the future horizon λ = λ H− , because the log divergence in t(λ) gets canceled by the second term of (42) and v(λ) is finite across the future horizon. On the other hand, the divergence at the past horizon λ = λ H+ is not canceled and v(λ) diverges across the past horizon. The argument of the logarithm in (41), on the other hand, is always positive for r > 0; see (34). The above parametric expression of r, t, x, v in terms of λ is useful for understanding the behavior of the geodesic, but for computational purposes it is useful to eliminate λ and write t, x, v as functions of r. Inverting the relation (34), we obtain the following two branches for λ(r): Note that we have been setting λ 0 = 0. Substituting (43) into (40)-(42) and setting x 0 = 0, we get We will refer to the branch given by (t(r) + , x(r) + , v(r) + ) as branch 1 and the one given by (t(r) − , x(r) − , v(r) − ) as branch 2. Only by combining both branches can we recover the full geodesic described by (34), (40) and (41). From (40) and (41), we can compute the coordinate difference between the two boundary points on which the geodesics end: We see that geodesics connecting endpoints at the same time (∆t = 0) correspond to taking E = 0. If one is interested in such geodesics in the pure BTZ geometry, then we have only to consider the E = 0 case. However, in section III, we will be interested in geodesics in the Vaidya geometry, which are composed of two pieces stretching respectively in AdS and BTZ, and glued together at the shell location. In this case, one should also consider E = 0 geodesics since they can appear as part of geodesics between equal-time endpoints in the full spacetime. Now let us focus on the equal-time geodesics with E = 0 and compute the relation between the spatial boundary separation and the geodesic length L thermal in the black brane background. By setting E = 0 in (48), is computed as On the other hand, the geodesic length L thermal is computed from (43) as: where r 0 = 1/z 0 is the bulk IR cutoff and we dropped the O(1/r 0 ) quantity that vanishes as r 0 → ∞. In the last equality, we used (49).
As discussed in Section II A 3, the geodesic length (50) is divergent in the r 0 → ∞ limit and should be regularized by subtracting the corresponding divergent part of the quantity in pure AdS geometry. Because the pure AdS geometry can be obtained by setting M → 0 or r H → 0 in the BTZ geometry (29), we can simply send r H → 0 in various quantities to obtain the corresponding quantities for pure AdS. By sending r H → 0 in (50), we obtain the geodesic length in pure AdS: Subtracting the divergent part of this from (50), we obtain which allows us to compute the thermal renormalized two-point function defined in (28). Observe that the twopoint function result computed in the geodesic approximation coincides with the weak coupling result (14) for ∆ = 1. Actually this is a consequence of conformal invariance, because the boundary CFT lives in noncompact space (namely, the x direction is not compactified in our setting). We can also obtain the expression for the equal-time geodesic in pure AdS geometry by sending r H → 0 in (45). In doing this, we should send J → ∞ at the same time because (49) becomes J = 2/(r H ). In this limit, (43) and (45) give where we dropped constants. The second equation (54) says that the equal-time geodesic is a semi-circle on the z-x plane, Note that the pure AdS results (51), (53), (54) and (55) are correct not only for d = 2 but for general dimensions, because it is only the z-x plane that is relevant for such geodesics. In any dimension we can consider a Vaidya type background (1) with constant mass function m(v) = M . This is a (d + 1)-dimensional black brane geometry, and in the special case of three dimensions it can be put in the form (29). We consider space-like geodesics connecting the boundary points (t, x 1 ) = (t 0 , − /2) and (t , x 1 ) = (t 0 , /2) of such black brane geometries (for d = 3, we also have x 2 = x 2 and for d = 4, (x 2 , x 3 ) = (x 2 , x 3 )). If we parametrize the geodesic in terms of the coordinate x 1 , which we will simply denote x in the rest of the Section, the solution will be given by the profiles v = v(x) and z = z(x). Inserting the IR bulk cut-off z 0 , the boundary conditions read and are symmetric with respect to the z-and v-axis. The geodesic length reads with ≡ d/dx. The action is not explicitly dependent on x, which implies the existence of a conserved quantity. The conservation equation reads where z * is the value of z(x) at the midpoint. Notice also that v is a cyclic coordinate, implying a conserved momentum. From Eq.
Inserting the relation between v and z of Eq. (59) into the conservation equation (58) yields a first order differential equation in terms of z(x) only. We take the derivative of this equation with respect to x and solve the resulting second order differential equation starting form the midpoint, where z(0) = z * and z (0) = 0. We can then read off the boundary separation corresponding to a specific midpoint value z * via z( /2) = z 0 . By substituting the conservation equation into the expression for the geodesic length, we can simplify the expression to The renormalized length δL thermal is from which one obtains the renormalized equal-time two-point function through (28). We plot the renormalized length δL thermal as a function of the spatial scale in Fig. 5.

B. Space-like Wilson loops
Above we examined the Wightman functions for field theories at finite temperature in order to get a sense of what thermal equilibrium looks like according to this measure. A second non-local probe can be used to discuss the thermality of field theories is the (expectation value of the) Wilson loop. The Wilson loop is a gauge-invariant observable, constructed as the path-ordered contour integral over a closed loop C of the gauge field Wilson loops contain useful information about the non-perturbative behavior of (non-abelian) gauge theories, but are in general hard to calculate. In the AdS/CFT correspondence, the expectation value for the Wilson loop is related to the string partition function with a string worldsheet Σ extending in the bulk and ending on the loop C at the boundary where we integrate over all inequivalent string surfaces Σ such that ∂Σ = C at the AdS boundary and where A(Σ) corresponds to the string action. In the strongly coupled limit, we can consider a saddle point approximation for the string partition function and reduce the calculation of the Wilson loop expectation value to determining the minimal area surface of the (classical) string worldsheet whose endpoints trace out the desired Wilson loop C on the AdS boundary A(Σ 0 ) represents the area of the minimal area surface Σ 0 with boundary C. The surface Σ 0 is a solution to the equations of motion arising from the bosonic action of the string. Thus, the expectation value of a Wilson loop in the AdS/CFT correspondence is given by the area of a minimal surface in AdS space which is bounded by the loop C [42].
In this section, we will compute space-like Wilson loop expectation values in the vacuum and in the thermally equilibrated theory (the black brane background) for d ≥ 3-dimensional field theories. In Sec. III, we will compute Wilson loops in a strongly coupled theory following a quench and compare them to the thermal results to assess progress towards thermalization. We will study two loop shapes -circles and strips as sketched in figure 6 -to test how and whether the shape of the loop affects the rate of thermalization.

Circular Wilson loops
We first examine the circular space-like Wilson loop in a strongly coupled d ≥ 3-dimensional field theory at finite temperature. At zero temperature the dual geometry is pure AdS d+1 , and the associated minimal surface in AdS is simply a hemisphere [51]. This surprisingly simple solution to a complicated second-order non-linear equation of motion was obtained by taking an infinite straight Wilson line on the boundary, and applying a special conformal transformation to map this straight line to a circle. Since the conformal group of the boundary corresponds to reparameterizations of AdS we can apply the appropriate reparameterization to the AdS space, and it maps the original worldsheet -which was flat and extended straight into the AdS -to the hemisphere. This trick only works in pure AdS space where conformal invariance is preserved, so we are obliged to use numerical methods of solving the equations of motion when we want to calculate the Wilson loop for thermal or out-of-equilibrium states.
For now, we want to study the minimal hemisphere surface in a thermal black brane background, with metric given by (1) with m(v) = M , where M is the constant tension of the black brane. On the boundary we choose a two-dimensional plane (x 1 , x 2 ) in which the circular Wilson loop is set. To parameterize the Wilson loop we introduce polar coordinates (ρ, ϕ). The minimal space-like surface with the circular Wilson loop as basis also has an azimuthal symmetry. The tip of the surface occurs at (v, z, x) = (v * , z * , 0). The cross section at fixed z and v is a circle, and thus the surface is conveniently parameterized in terms of the radii ρ of these circles. The Nambu-Goto action for the string with circular symmetry in the (d + 1)-dimensional AdS black brane background is where we integrated out the ϕ-factor as both z and v are independent of ϕ due to the circular symmetry. This action has an explicit ρ-dependence, implying that the second order differential equation in z(ρ) cannot be integrated to a first order differential equation. The coordinate v is still a cyclic coordinate, which implies that the associated momentum is conserved. To solve for an equal-time Wilson loop, we set the conserved energy to zero, which implies a relation between v and z as given in eq. (59). Inserting this relation in the action, the resulting equation of motion reduces to a second order differential equation in z(ρ) One can easily check that the hemisphere is an analytic solution in pure AdS (obtained by setting M = 0 in the above). It is given by In case of a black brane background (M = 0) we have to resort to numerical means. We impose boundary conditions at the midpoint assuming reflection symmetry along the z-axis and read off the boundary radius z(R) = z 0 . Note that because the second term in the first bracket of (68) causes numerical issues at ρ = 0 (the numerator and denominator both go to zero, but the ratio should take a fixed value), we cannot construct the solutions starting from the midpoint. Therefore we choose to solve (68) in the neighborhood of the midpoint by expanding around ρ = 0 to quadratic order (odd powers vanish by the symmetry) The boundary conditions for the numerical solution are obtained by matching at a point ρ p close to the midpoint; z(ρ p ) = z p (ρ p ) and z (ρ p ) = z p (ρ p ). The logarithm of the expectation value of Wilson loop operators is approximated by the area (A thermal = α A NG ) of the minimal surface. The largest contributions to the area are coming from the near-boundary region, because of the diverging AdS volume there. We subtract the divergent piece, which is proportional to the UV-cutoff 1/z 0 . This procedure is interpreted in the dual CFT as removing the UV divergence of the Coulomb self-energy of a point charge. Subtraction of R/z 0 gives the regularized area The regularized area is depicted for AdS 4 and AdS 5 in Fig. 7.

Infinite Rectangular Strips
A less symmetric Wilson loop is the rectangular strip parametrized by the boundary coordinates (x 1 , x 2 ). Assume that the strip is translationally invariant along the x 2 -axis, such that the profile of the associated minimal surface in AdS d+1 is described by z(x 1 ) and v(x 1 ). In the following, we again denote x ≡ x 1 . We impose the boundary conditions As for the circular Wilson loop, the equation of motion for the minimal surface is obtained by minimizing the Nambu-Goto action. In this case, because of the symmetries, the analysis closely resembles that of geodesics with some small differences in the action and the equations of motion. The area for a segment stretching between x 2 ∈ (0, R) in Eddington-Finkelstein coordinates is given by The action does not exhibit an explicit x-or v-dependence, hence we can use the existence of two conserved quantities to simplify: these can be taken to be the midpoint of z(x), z * , and the momentum conjugate to v(x). The z * conservation equation reads while the conservation of the momentum conjugate to v(x) implies the same relation as in Eq. (59). Inserting this relation in the conservation equation gives a first-order differential equation in z(x). We then take the derivative with respect to x to obtain a second-order differential equation in z(x) We impose the boundary conditions z (0) = 0, z(0) = z * at the midpoint, and read off the boundary separation via z( /2) = z 0 . Inserting the conservation equation in the Nambu-Goto action we get the following formula for the area (A thermal = α A NG ) of the minimal surface .
The regularized area reads and is plotted in Fig. 8 for a four-and five-dimensional background. We have chosen to plot the results after dividing out R/π as we are concerned purely with the dependence of the area on the boundary separation .

C. Entanglement entropy
It is still an open question whether a suitable notion of "local entropy" density exists which is valid out of equilibrium and which satisfies some basics physical properties; in particular it must be a non-decreasing function of time. To define a local entropy density using gauge/gravity duality, one typically employs horizons since their area increases in time, and it has been suggested in [11,52] that apparent horizons provide more compelling notions of local entropy density then global horizons do. In our case, the entropy associated to the apparent horizon instantaneously thermalizes, in the limit of a zero thickness shell, and does not provide a useful probe of thermalization. Skepticism against the identification of entropy with apparent horizon area in full generality has also been raised in a number of other papers, including [27,52,53].
Another way of assessing thermalization at different spatial scales in an out-of-equilibrium theory is by measuring the entanglement entropy associated with volumes of different sizes and shapes. To review, consider dividing a quantum mechanical system in a state |Φ into two spatially disjoint parts A and B. The density matrix of the system restricted to A and its entropy are computed by taking the traces In a (1 + 1)-dimensional conformal field theory, where A is an interval of length , the entanglement entropy can be calculated analytically, yielding the universal result [54,55]: Here a is the UV cutoff of the field theory, c is the central charge, and β = T −1 denotes the inverse temperature. The entanglement entropy S(A) describes the amount of information loss associated with the restriction of an observer to the volume A. In the vacuum state for d > 2 dimensions, S 0 (A) is proportional to the surface area of A [56,57]; in d = 2 dimensions, as Eq. (80) shows, S 0 depends logarithmically on the length of the interval . At non-zero temperature, S(A) receives an additional contribution, which can be interpreted as thermal entropy [43,44]. In the limit T → ∞, the thermal contribution is proportional to the volume of the region A, just like the statistically defined thermal equilibrium entropy. Computing the time dependent entanglement entropy as a function of spatial scale and studying its approach to S T ( ) thus will provide a probe of scale-dependent thermalization.
There is a precise proposal for computing entanglement entropy in strongly coupled field theories with AdS duals [26,43,44], where we think of the field theory as living on the boundary of an AdS space. Specifically, consider the boundary ∂A of a connected region in the field theory whose entanglement entropy we wish to compute. For a two-dimensional field theory, ∂A is a pair of points, for a three-dimensional field theory ∂A is a closed curve, and for a four-dimensional theory a surface. Now construct the minimal surface σ A in AdS space that meets ∂A on the AdS boundary. For a two-dimensional field theory, σ A is a geodesic in AdS 3 that approaches the boundary points ∂A, while for a three-dimensional field theory, σ A is a minimal surface in AdS 4 with boundary ∂A. For a four-dimensional field theory, σ A is a minimal volume in AdS 5 which ends on the surface ∂A on the AdS boundary. The entanglement entropy of the region A of the field theory is then given by where G N denotes Newton's gravitational constant (and "Area" stands for the length of geodesics, the area of 2surfaces, and the volume of 3-surfaces). The authors of [26,43,44] showed that this formula precisely reproduces the universal entropy formula of two-dimensional conformal field theories in thermal equilibrium [54,55]. The restriction of the entanglement entropy to a finite spatial volume can be understood as a kind of coarse graining, by the discretization of the available momentum space modes, of the information contained in the quantum state [58]. Thus, to measure entanglement entropy in two dimensional field theories, we need to calculate space-like geodesics in the dual asymptotically AdS 3 background. This means there is an intimate relation between two-dimensional entanglement entropy and the equal-time Wightman functions computed in the geodesic approximation. Indeed, the two quantities can be interpreted in terms of each other, as was also observed in [27]. For three dimensional field theories, we need to compute space-like, two-dimensional minimal surfaces in a dynamical, asymptotically AdS 4 background. These coincide precisely with classical space-like worldsheets of fundamental strings with a fixed boundary, as previously studied in [28]. As we discussed above, the exponential of the action of such objects (essentially the area) is known to give the expectation value of Wilson loops in the field theory [42]. Thus, in these three-dimensional field theories, there is an intimate relation between entanglement entropy and Wilson loop expectation values, as we explore in this article.
For four dimensional field theories, the proposed measure of entanglement entropy for spherical regions involves minimal volume hyper-surfaces whose cross-sections with respect to the AdS radial coordinate are spheres. In particular, in AdS 5 we can look for a 3-dimensional volume with 2-sphere cross-section which terminates on the AdS boundary. This allows us to examine the thermalization of the entanglement entropy in the 4-dimensional field theory. The methods used to compute these volumes are very similar to those for the circular Wilson loop described in Section II B 1, with only some powers changing with the dimension. Briefly, the action for the general spherically symmetric p-volume in AdS d (p ≤ d − 1) black brane background is where A S p−1 is surface area of a unit S p−1 . Things continue to proceed as for the Wilson loops: the solution in the pure AdS (M = 0) background is now a hemi-hypersphere and still obeys (69). What changes is the volume evaluated on the solution, which we evaluate with a cut-off near the boundary and subtract to renormalize as before. The result in the vacuum (M = 0) is The regularized volume is calculated analogously to previous cases by subtracting the divergent part of the volume of a pure AdS 3-surface with terminates on the same sphere on the boundary, We plot the regularised volume δV thermal as a function of the radius R of the sphere on the boundary in Fig. 9.

III. SCALE DEPENDENT THERMALIZATION FROM ADS/CFT
In the previous section we discussed how the equal-time Wightman function (related to geodesics in AdS space) and Wilson loops (related to minimal surfaces in AdS space), and entanglement entropy (related to minimal lengths, surfaces and volumes in different dimensions) can probe thermal equilibrium. In this section we will use the same probes to study the dynamics of equilibration in strongly coupled field theories. The basic setup is to drop a shell of matter with vanishing rest mass ("null dust") into AdS space to form a black hole. This is dual to the homogeneous injection of energy into the dual field theory and its subsequent thermalization. We want to probe the rates at which thermalization occurs on different spatial scales. To this end we consider dynamical Vaidya-type backgrounds, again given by equation (1), but this time using the continuous mass function of equation (4). This dynamical background continuously interpolates between pure AdS and AdS with a Schwarzschild black brane. Recall that the parameter v 0 in the mass-function (4) dictates how "thick" the shell is. v 0 → 0 is the step function limit corresponding to an infinitely thin shell. In this limit, the infalling shell represents a shock wave. A. Equilibration of the two-point function

Two-dimensional field theories: analytic treatment
We first consider the d = 2 case where analytic computations are possible. In this case, the Vaidya metric (1) becomes In order to have analytic control over the solution, let us consider the thin shell limit v 0 → 0 of the mass profile (4), which gives where θ(v) is the step function.
Outside the shock wave, v > 0, the metric is the standard planar black brane metric (29) which we studied in detail in the previous section. Inside the shock wave, v < 0, the metric is the Poincaré AdS 3 : We would like to study geodesics in the AdS 3 Vaidya geometry which is (87) and (88) glued together across the infalling shell at v = 0. In particular, we focus on the equal-time geodesic which starts and ends at the same time t = t 0 on the boundary r = ∞. There are two possible cases: (i) the geodesic does not reach the shell and is entirely outside of it, and (ii) the geodesic crosses the shell. Because case (ii) is more involved than case (i), let us first discuss case (i) briefly and then turn to the discussion of case (ii).
In case (i), the geodesic is given by the equal-time geodesic in the pure black brane geometry, namely by (44)-(46) with E = 0. Therefore, the relation between the (renormalized) geodesic length δL thermal and the boundary separation is given by (52). Now let us turn to case (ii). In this case, the part of the geodesic that is inside the shell is given by the equal-time geodesic in the pure AdS geometry, namely by (54), which we write in the following form: where we set r * ≡ 2/ . On the other hand, the part of the geodesic that is outside the shell is given by the geodesic in the pure black brane geometry, namely by (44)- (46). In the present case, we should not set the parameter E to zero, because the geodesic gets refracted at the shell and it does not have to be in a constant t slice.
The two parts of the geodesics, the inside part and the outside part, should be connected so that the total geodesic length is minimized. Just like Snell's law, this can be stated as a refraction law for the angles entering and exiting the shell, as follows. Let us write the metric (87), (88) as and focus on the region very close to the shell at v = 0. Consider a point P in (P out ) just inside (outside) the shell, and let the coordinate difference between P in and P out be ∆X = (∆v, ∆r, ∆x). Take another point M on the shell v = 0, and let the coordinate difference between P in (P out ) and M be ∆X in = (∆v in , ∆r in , ∆x in ) (∆X out = (∆v out , ∆r out , ∆x out )), so that ∆X in + ∆X out = ∆X. Then the distance from P in to P out via M is We want to find the point M that minimizes ∆s. Extremizing this with respect to ∆r in , ∆x in , we find and therefore From these we obtain the following refraction law as required: By plugging the inside solution, (89), and the outside solution, (44)- (46), into the refraction condition (95), we obtain the relation between the parameters of the geodesic: where r sw is the value of r at which the geodesic intersects the shell, and how to determine it will be explained below. The first sign combination applies to branch 1 for r sw ≤ r H / √ 2 and branch 2 for r sw ≥ r H / √ 2, while the second sign applies to branch 1 for r sw ≥ r H / √ 2 and branch 2 for r sw ≤ r H / √ 2. Recall that branches 1 and 2 were defined in section II A 4, below eq. (46).
It is possible that the geodesic crosses the horizon (r = r H ) before reaching the shell (v = 0). When this happens, let us require that the geodesic goes to the boundary r = ∞ at λ = ∞ and, as we decrease λ, it enters the horizon with v staying finite. In the previous section, we saw that the geodesic crosses the horizon at λ = λ H± and, below (42), we observed that v is finite across λ = λ H− but not across λ = λ H+ . Therefore, we need that λ H+ < λ H− which means that E < 0. So, we choose the second sign combination in (96): This leaves branch 1 for r sw ≥ r H / √ 2 and branch 2 for r sw ≤ r H / √ 2. The equation (97) depends on r sw which is where the geodesic hits the shell. Because the shell is at v = 0, the value of r sw is determined by where in the second equality we used (46) and the ± signs correspond to branch 1 and 2, respectively. If we rewrite E, J in favor of r sw , r * using the relation (97), this becomes Note that the two possible signs in (98) lead to the same condition (99). To solve this, let us make the convenient definitions The range of the parameters is θ ∈ [0, π 2 ], ρ ≥ 0, and a ≥ 1. In terms of these, the relation (99) becomes with c ≡ cos θ. This can be solved for ρ, giving two branches of solutions: We can see that only ρ + is allowed, as follows. First, we can show that, for fixed c, the derivative of ρ(a, c) ± with respect to a is always negative (positive) in the range a ∈ (1, ∞), c ∈ (0, 1] [69], so that ρ(a, c) ± monotonically decreases (increases) as we increase a. Furthermore Therefore, for any finite values of t 0 , For r sw ≥ r H / √ 2 (i.e. ρ > 1/ √ 2), this means that we should take the ρ + branch. The case with r sw < r H / √ 2 (i.e. ρ < 1/ √ 2) needs some more work. We only need to consider geodesics that do not reach r = 0, since otherwise we cannot glue branch 1 to branch 2 at λ min before the geodesic reaches the singularity. Thus, when r sw < r H / √ 2, we need A ± < 0 and B ± > 0, that is J < 1 + E, for E < 0 and 0 < J < 1. Using (97), the last condition translates into the following condition for ρ: This means that we should take ρ + in (102) for all c ∈ [0, 1]. In Figures 10 and 11, geodesics in the (v, r) plane are plotted for different values of t 0 and . The geodesics start from the boundary at (r, v) = (∞, t 0 ), plunge in the bulk in the planar BTZ geometry, eventually refract at the shell at v = 0 and propagate in pure AdS. All profiles are symmetric under x → −x. The curve in black dashed is the apparent horizon.
In Figures 12 and 13, we display the boundary time and spatial separation dependence of the geodesics in the Carter-Penrose diagram of Vaidya spacetime. The transformations to the (U, V ) coordinates are given by for v < 0, for v > 0 and r > r H , and for v > 0 and r < r H . In these coordinates, the "singularity" is located at U + V = 1 and is depicted by a wavy line in Fig. 12 and 13, the shell is at V = 0 and the horizon is at U = 1/2 (dashed line). With all the relations above in hand, we can compute the geodesic length and spatial boundary separation in terms of the parameters of the geodesic. The geodesic length is the sum of the geodesic length for the inside part, which can be computed from (53), and the one for the outside part, which can be computed from (43). The result is Here, just as in (52), δL has been renormalized by subtracting an IR divergent quantity 2 ln(2r 0 ). The spatial boundary separation can also be computed by summing the contribution from the inside part (89) and the outside part (45), the result being Although δL and are written in terms of r sw and r * in the above expressions, we can write them in terms of t 0 and s ≡ √ 1 − c 2 using the definitions (100) and plugging in the solution ρ = ρ(a, c) + of Eq. (102) Explicitly, where s( , t 0 ) ∈ [0, 1] is parametrically defined by 2. d ≥ 2-dimensional field theories: numerical analysis As in the AdS 3 case, we consider a four-or five-dimensional dynamical metric which interpolates between pure planar AdS at early times and a Schwarzschild black brane at late times. The transition is induced by an infalling shell of null dust. The metric in (d + 1) dimensions was given in Eq. (1). We again consider the mass function (4), where the profile parameter is set to v 0 = 0.01 in the remainder. We consider geodesics with a boundary separation along x 1 , denoted x in the following, while the (d − 2) coordinates (x 2 , . . . , x d−2 ) of both endpoints are the same. This implies that z and v only depend on x, i.e. z = z(x) and v = v(x). We assume that x runs between − /2 and + /2. Then the length of the geodesic is given by: where ≡ d/dx. We notice that the integrand has no explicit x-dependence, implying the existence of a conserved quantity, similar to the black brane case. The conservation equation reads: The two equations of motion following from (114) are whereṁ(v) = dm(v)/dv. It is possible to show that (115) and (116) imply (117) after taking the derivative of (115) with respect to x. Therefore we restrict ourselves to solving (115) and (116). We construct solutions z(x) and v(x) with the x → −x symmetry which is already present in the equations of motion, and which satisfy the boundary conditions After constructing a geodesic for a specific choice of values (z * , v * ), information about the boundary separation and boundary time at which the geodesic is inserted follows from The on-shell length is obtained from Eq. (114), upon use of the conservation equation (115), This length should be regularized by subtracting the cut-off dependent part 2 ln(2/z 0 ), yielding δL( , t 0 ).

Thermalization
Having collected the geodesic lengths in the infalling shell backgrounds in AdS 3,4,5 , we can use these to describe the process of thermalization following a quench. To this end, we measure the approach to thermal equilibrium by comparing δL at any given time with the late time result δL thermal . In any dimension, this compares the logarithm of the two-point correlator at different spatial scales with the logarithm of the thermal correlator. However, we find it more revealing to considerL ≡ δL/ , where we divide by the spatial separation on the boundary (we will discuss the significance of this later in Section III C where we deal with entanglement and Kolmogorov-Sinaï entropies). In figure  14 we plot this measure for two-, three-and four-dimensional field theories. In all three cases we observe a delay in the onset of thermalization. The reason for this "delay" is simply that the effect of the medium only becomes fully apparent at distances of the order of the thermal screening length D ∼ (πT ) −1 . Although a very small volume of linear dimension D would appear fully thermalized after a time t 0 = /2, its contribution to the entropy of a large volume would be disproportionately small, because it does not support the momentum modes that constitute the thermal medium at large. As a consequence, the rapid linear increase of the logarithm of the two-point function seen in Fig. 14 thus only sets in after some delay. Obviously, this effect is more pronounced in higher dimensions. From the curves in Fig. 14 we can extract different thermalization times for any spatial scale: (1) The critical time τ crit at which the tip of the geodesic grazes the middle of the shell at v = 0. This can be computed by asking when a geodesic with a given boundary separation in the black brane geometry outside the infalling shell just grazes the latter: τ crit ( ) = z * z0 dz 1−M z d , where z * is determined by the boundary separation . (2) The half-thermalization time τ 1/2 : time that measures the duration for the curves to reach half of their equilibrium value.
(3) The time τ max : time at which thermalization proceeds most rapidly, thus for which the curves in Fig. 14 are steepest.
For two, three and four dimensional field theories all these times are plotted in Fig. 15. In d = 2 we can analytically derive the linear relation τ crit ≡ /2, as also observed in [27]. The linearity of τ crit ( ) in two dimensions is expected from general arguments in conformal field theory [55], and the coefficient is as small as possible under the constraints of causality following a quantum quench. The thermalization time scales τ 1/2 and τ max for three-and four-dimensional field theories (Fig. 15, middle and right) are sublinear in the spatial scale. In the range we study, the complete thermalization time τ crit deviates slightly from linearity, and is somewhat shorter than /2. These observations pose the question whether a rigorous causality bound for thermalization processes exists or not. We will come back to this question below, after examining other probes than two-point functions. One of the issues here is that the initial conditions are homogeneous -thus different domains will be independently coming to equilibrium at the same temperature possibly leading to apparent violations of causality.
General arguments for a sharp quench in a 2-dimensional field theory [27,55] predict a nonanalytic feature where thermalization at a spatial scale is completed abruptly at τ crit ( ). This is evident in Fig. 14 (left) as a sudden change in the slope at τ crit , smoothed out only by the small non-zero thickness of the shell, or equivalently, by the intrinsic duration of the quench. We find a similar (higher-order) non-analyticity for d = 3, 4 (Fig. 14, middle and right) and expect this to be a general feature for quantum quenches in all dimensions.   15 shows that complete thermalization of the equal-time correlator is first observed at short length scales, or large momentum scales (see also [25]), in contrast with the behavior of weakly coupled gauge theories. In the "bottomup" scenario [8] applicable to that case, hard quanta of the gauge field do not equilibrate directly by randomizing their momenta, but gradually degrade their energy by radiating soft quanta, which fill up the thermal phase space and equilibrate by collisions among themselves. This bottom-up scenario is linked to the infrared divergence of the splitting functions of gauge bosons and fermions in perturbative gauge theory. It contrasts with the "democratic" splitting properties of excitations in strongly coupled SYM theory that favor an approximately equal sharing of energy and momentum [46]. One might, therefore, have expected that thermalization proceeds fundamentally differently ("topdown") in strongly coupled gauge theories, and this is evident within the AdS/CFT paradigm with the natural initial conditions used here.

B. Equilibration of Wilson loops
Another nonlocal probe of thermalization is the Wilson loop expectation value. As we discussed in Sec. II B, this quantity is related in the AdS/CFT correspondence to the area of a minimal surface in AdS space that bounds the desired loop on the AdS boundary. In this section we compute such minimal areas in the infalling shell background.

Circular Wilson loops
We first consider minimal surfaces corresponding to circular Wilson loops in the Vaidya background. The corresponding minimal surfaces were analyzed for the vacuum (pure AdS) and finite temperature equilibrium (AdS black brane) situations in Sect. II B 1. This set-up was also analyzed in [28] in the context of studying entanglement entropy. The parameterization of the surface and boundary conditions are as given in Section II B 1, the only difference being that the mass function m(v) is now given by (4). The area functional is a slight modification of (66) (A = α A NG ) given by where ≡ d/dρ. The explicit ρ dependence again means there is no conservation equation. The equations of motion become quite involved and we omit them here. We solve the equations of motion numerically, and although we want our output as a function of the boundary radius R, for practical purposes we have to input boundary conditions at the tip of the surface at ρ = 0 where z = z * and v = v * . For each value of z * we consider, we find numerical solutions for the functions z(ρ) and v(ρ) for various values of v * , each of which reaches the boundary (which is cut-off at z 0 ) at ρ 0 (v * ). We then make an interpolation of ρ 0 (v * ), which allows us to find the exact value of v * such that ρ 0 (v * ) = R. For each value of z * , we now know the value of v * giving the surface with boundary radius R, and we calculate the area (using the functional above) and the boundary time t 0 = v(R). These are the quantities we finally plot. There are further numerical challenges with a divergence as ρ → 0 and with integrating our solutions near the AdS boundary, where contributions to the area are weighted much higher. Thus we used an expansion around ρ = 0, and different techniques to deal with a separate integral near the boundary, as necessary.
As in (72) we regulate the boundary by subtracting the cut-off dependent piece of the AdS area, We define δÃ ≡ δA/(πR 2 ) dividing by the area of the region on the boundary bounded by the loop. We then plot δÃ − δÃ thermal as a function of the boundary time, subtracting the black brane value to show the approach to thermalization. This is shown in Fig. 16 for Repeating the analysis performed for Wightman functions, we calculate the three thermalization times defined earlier (see Sec. III A 3) as a function of the loop diameter (Fig. 17). We use the diameter rather than the radius here as it is the analogue of the separation that we plotted for the geodesics. The complete thermalization time τ crit (D) is close to being a straight line of slope 1/2 for three dimensional theories over the range of scales that we study (also see [28]) (it would be unit slope as a function of the radius). But in four dimensions τ crit (D) deviates somewhat from linearity and is shorter than D/2. Overall, our thermalization times for Wilson loop averages are remarkably similar to those for two-point correlators . This suggests that higher-order correlators thermalize similarly to the basic Green function in these strongly coupled theories, so that any of these non-local probes fairly assesses progress towards thermalization.

Infinite Rectangular Strips
As a second example we consider an infinite strip, similarly to [28]. A key difference between the strip and the Wilson circle is that the size of the latter is set by a single scale, the radius. The strip could be regarded as a limit of a ellipsoidal loop with a highly elongated semi-major axis. Thus there are in effect two scales involved in the Wilson strip -the width and the (infinite) length. As before, at early times the background is AdS 4 , but evolves to a black brane at late times due to an infalling shell of null dust. The metric is at all times given by (1) with d = 3 or 4. Again, the mass function m(v) determining the evolution of the metric as a function of the bulk light-cone time v is modeled by the expression (4). We set the parameter v 0 = 0.01. At the AdS 4 boundary we consider a rectangular strip parametrized by the coordinates x 1 and x 2 , such that where R will be taken to infinity [70]. The coordinate x 1 is denoted x in the rest of the section. The area of the string surface with this rectangular base is given by when we consider the following embedding profile for the string: We notice that z and v only depend on x and that x itself does not appear in the action, similarly to the geodesic length in a dynamical Vaidya background. There exists a conserved quantity corresponding to the maximum value of z denoted as z * . The conservation equation then reads: From the action (124) we obtain the following two equations of motion: One can show that taking the derivative of (126) with respect to x and combinig with (127) leads to (128). This it is sufficient to solve (126) and (127). In order to do this numerically we use the symmetry of the string surface, i.e. we construct one half of the solution z(x) and v(x) and then use the reflection symmetry to construct the other half. Our boundary conditions are z(0) = z * , v(0) = v * , v (0) = 0 = z (0). After constructing a string surface satisfying these boundary conditions, we read off the boundary separation and the boundary time using z( /2) = z 0 , v( /2) = t 0 .
Using the symmetries of the embedded string surface and the conservation equation (126), we obtain the on-shell area as: We regulate the area by subtracting the cut-off dependent part We then compare the area in the shell background with the area in a thermalized background at different times for fixed scales in Fig. 18, where we defined δÃ ≡ δA/(R /π). In four dimensions (d = 3) it turns out that there is a range of times t 0 and Wilson strip widths for which there are three different minimal surfaces in AdS 4 . Thus, in this case, we observe a "swallow tail" in the Wilson loop thermalization curve for large enough boundary separations, as was also noted in [28]. However, in five dimensions (d = 4) the thermalization curves do not exhibit the swallow tail. In Fig. 19 we plot the three different thermalization times τ crit , τ 1/2 and τ max for d = 3 and d = 4 dimensional field theories. In three and four dimensions, the time for complete thermalization, τ crit (dashed green curve), exceeds the linear relation /2. The half-time τ 1/2 (orange curve) is sub-linear. The time τ max (dot-dashed blue curve) exhibits a different behavior in three and four dimensions. The time for complete thermalization for a strip of width is slower than the thermalization time of a circular loop of radius R. This may be because the long direction of the strip provides a second, larger scale.

C. Equilibration of the entanglement entropy
Finally we consider how the entanglement entropy equilibrates following a quench. Because it measures all contributions to the information loss caused by the restriction of the field theory to a finite volume, this quantity provides for a more comprehensive measure of equilibration than either the two-point function or the Wilson loop expectation value. The entanglement entropy can be considered as a special case of the standard coarse grained entropy of equilibrating quantum systems [58]. As discussed in Sec. II C, the entanglement entropy in a two-dimensional field theory is related to geodesic lengths in AdS 3 , while in three-dimensional field theory entanglement is related to minimal surface areas in AdS 4 . Both of these quantities have been computed above in the infalling shell background in the process of studying two-point functions and Wilson loops. It is interesting that the Wightman function in 2 dimensions and the Wilson loop in 3 dimensions are so closely related to entanglement.
For 4-dimensional theories, to study the entanglement entropy in spherical regions, we need to compute the volume of minimal three-surfaces with S 2 cross-section in the asymptotically AdS 5 infalling shell geometry. As with the pure AdS and black brane cases we discussed in Section II C, the method is a straightforward generalization of the circular Wilson loop case, so we do not repeat the details of Section III B 1. We replace the black brane tension M in equation (82) with the dynamical mass function (Eq. (4)). We can regulate the resulting volumes by subtracting the divergent part of the 3-volume of the 3-surface of the same boundary radius in pure AdS, which is given analytically in (83). We call this regulated volume δV. However, as in previous cases, to illustrate the thermalization process we find it more instructive to subtract the renormalized 3-volume of the solution of the same radius in the black brane background, δV thermal . We ultimately plot δṼ − δṼ thermal against the boundary time in the left panel of Fig. 20.Ṽ is defined as V divided by the volume of the region that it bounds on the boundary. As in previous cases we can also calculate the three thermalization times we defined earlier for the 3-volume, at different values of the diameter D = 2R. These are shown in the right panel of Fig. (20). Our thermalization times for Wilson loop averages and entanglement entropy seem remarkably similar to those for two-point correlators. Slightly "faster-than-causal" thermalization, possibly due to the homogeneity of the initial configuration, seems to occur for the probes that do not correspond to entanglement entropy in each dimension. For the latter, the thermalization time is linear in the spatial scale and saturates the causality bound. As the actual thermalization rate of a system is set by the slowest observable, our results suggest that in strongly coupled theories with a gravity dual, thermalization occurs "as fast as possible" at each scale, subject to the constraint of causality.
The average growth rate of the coarse grained entropy in nonlinear dynamical systems is measured by the Kolmogorov-Sinaï (KS) entropy rate h KS [59], which is given by the sum of all positive Lyapunov exponents. For a classical SU(2) lattice gauge theory in 4 dimensions, h KS has been shown to be proportional to the volume [60]. For a system starting far from equilibrium, the KS entropy rate generally describes the rate of growth of the coarse grained entropy during a period of linear growth after an initial dephasing period and before the close approach to equilibrium [61]. Here we observe similar linear growth of entanglement entropy density in d = 2, 3, 4 (leftmost panels Figs. 14, 16,20). For small boundary volumes, the growth rate of entropy density is nearly independent of the boundary volume (almost parallel slopes in the leftmost panels of Figs. 14, 16, 20 and nearly constant maximal growth rate in Fig. 21  left). Equivalently, the growth rate of the entropy is proportional to the volume -suggesting that entropy growth is a local phenomenon. However, in d = 2 where our analytic results enable study of large boundary volumes , we find that the growth rate of the entanglement entropy density changes for large , falling asymptotically as 1/ (Fig. 21 middle panel). Equivalently, the entropy has a growth rate that approaches a constant limiting value for large ( Fig. 21 right panel), and thus cannot arise from a local phenomenon. This behavior suggests that entanglement entropy and coarse grained entropy have different dynamical properties.
In summary, we have investigated the scale dependence of thermalization following a quench in 2-, 3-, and 4dimensional strongly coupled field theories with gravity duals. We found that the entanglement entropy sets a time scale for equilibration that saturates a causality bound. Our results raise interesting questions about the relationship between the entanglement entropy growth rate and the KS entropy growth rate defined by coarse graining of the phase space distribution.

IV. DISCUSSION
In this work, we explored the approach to thermal equilibrium of simple nonlocal observables in strongly coupled conformal field theories with a gravity dual. In particular, we investigated the behavior during thermalization of equal-time two-point functions, Wilson loops and entanglement entropy in boundary field theories holographically represented by asymptotically AdS geometries. The thermalization process was modeled in the bulk by the collapse of a spatially homogeneous thin shell of null dust that eventually forms a black hole. The thickness v 0 of the shell is related to the duration of the process of energy deposition in the boundary field theory. In the limit v 0 z 0 , where z 0 parametrizes the effective UV cut-off of the field theory, the energy deposition can be considered as instantaneous.
Our observables were calculated in a semi-classical approximation appropriate for high-dimension operators, in which they correspond to geodesics or minimal surfaces in the bulk. For the AdS 3 case, where the thermal limit is represented by the planar BTZ black hole, we were able to solve the problem analytically. Not only did this provide a check on our numerical approach, which is required for higher dimensional cases, it also allowed us to obtain results in regimes beyond reach of our numerical analysis.
In all cases we found that the thermal limit is reached after a finite time τ crit , which is a function of the geometric size of the probe in the boundary field theory, e.g. the separation of the two points of the equal-time Wightman function or the radius R of the circular Wilson loop. For those cases where the logarithm of the correlation function is proportional to the entanglement entropy of the enclosed area on the boundary, i.e. for geodesics in AdS 3 , minimal surfaces bounded by a circle in AdS 4 , and minimal co-dimension two surfaces bounded by a spherical shell in AdS 5 , we found that τ crit = /2 (AdS 3 ) and τ crit = D/2 (AdS 4,5 ), in the regimes we were able to study. Our result confirms a general rule for two-dimensional conformal field theories implied by causality (as also observed in [27]) and generalizes it to higher dimensions (as also done in [28] for AdS 4 ).
The conclusion we can draw from these results is that the decoherence and equilibration of instantaneously deposited energy generally propagates at the speed of light in the conformal field theory. This makes sense, because the process can be qualitatively understood as a cascade of gauge boson splittings, which transfers energy from all momentum scales into the infrared until all energy is distributed thermally. The speed by which this process evolves is determined by the speed of propagation of the gauge quanta, i.e. the speed of light. Alternatively, if the initial energy deposition were spatially inhomogeneous, one would expect the dissipation of the local inhomogeneities in the initial energy density to be constrained by the speed of sound, which is equal to c/ √ 3 in a conformal field theory. This argument invites the conjecture that local equilibration of energy, in the sense of hydrodynamics, is controlled by the speed of light in conformal field theories, whereas global equilibration requiring hydrodynamic transport of energy is governed by the speed of sound. On the other hand, we found that the naïve causality argument does not hold for the two-point function in 3d and 4d field theories and for Wilson loops in 4d field theories, which do not have an interpretation as entanglement entropy. In this case we obtained τ crit < /2. This is not as surprising as it might sound: we already know of other observables that thermalize faster than the entanglement entropy, for instance one-point functions.
We also found that the transition to full thermal equilibrium is abrupt and non-analytic, as expected from the causality argument. One might have objected that the non-analyticity of the transition to full thermal equilibrium is probably due to the geodesic approximation and will be smoothed out once we consider the exact Wightman function. Corrections to the geodesic approximation can be worked out by considering the exact first-quantized path integral representation of the propagator, whose saddle points are given by geodesics. Possible corrections are then due to fluctuations around the saddle point. However, if the geodesic saddle-point lies entirely within the AdS Schwarzschild geometry and does not graze the shell, such higher order corrections will not be sensitive to the presence of the shell, since the shell is invisible in perturbation theory around the geodesic. This remains valid when we couple the heavy bulk scalar field to other bulk degrees of freedom. At best, there could possibly be non-perturbative corrections due to other complex saddle points which are sensitive to the existence of the shell, but is unclear how to detect those in time-dependent situations. Therefore, for sufficiently small separations the two-point function appears to be exactly given by the thermal answer, which will then abruptly change to a non-thermal answer once the geodesic starts crossing the shell. As a result, we expect the non-analyticity to be valid beyond the geodesic approximation. The non-analyticity is probably a consequence of our treatment of the endpoints of the Wightman function, which are both located at the same bulk variable z 0 . The origin of a similar problem has been studied in the theory of boundary critical phenomena, where one typically introduces an extrapolation length τ 0 to tame the sharp UV cut-off of the conformal field theory [62,63]. This corresponds to the introduction of a smooth UV cut-off of the form exp(−Eτ 0 ), where E is the energy. The sharp cusp in the asymptotic result is then rounded over a region |t − /2| ∼ τ 0 . In our holographic setting, a similar approach would smear the two endpoints of the correlation function independently over some range of values z 0 . We have not studied the consequences of such a prescription.
One may wonder whether our results are in tension with those of [64], where instantaneous thermalization was found for two-point functions related to Brownian motion of a "quark" represented by a string stretching from the boundary to the horizon of an AdS 3 -Vaidya spacetime. While the two-point functions studied in [64] did not involve high-dimension operators, so that there is no reason to expect a geodesic approximation to be valid, one can get some intuition by computing geodesics in the induced geometry on the string worldsheet. We have done so, and find that geodesics connecting two points on the string outside the shell never cross the shell. So within a geodesic approximation, this would explain why the observable considered in [64] thermalizes instantaneously, in contrast to the observables studied in the present paper.
The fact that correlation functions of a small geometric size thermalize earlier than those of a large size implies that high momentum modes in the boundary field theory approach thermal equilibrium faster than long wavelength modes. In other words, the boundary field theory thermalizes first in its ultraviolet domain and later in the infrared or, as one might say, thermalization proceeds from the top down. In part, this is a consequence of the fact that the specific energy injection mechanism considered here has support in the UV: the energy shell is injected into the geometry near the boundary and then falls into the deeper regions of the AdS space, which represent IR modes of the boundary field theory. While this appears natural when viewed from the gravity side, it represents a radical deviation of the thermalization behavior of the dual field theory from that known in weakly coupled non-Abelian gauge theories, where thermalization occurs from the bottom up, independent of how the energy is injected into the field theory. For the reader unfamiliar with heavy ion collision phenomenology, we point out that the scale at which energy is deposited in the collision -the parton saturation scale Q s of the colliding nuclei -is higher by an inverse power of the coupling constant α s than the thermal scale T after thermalization. Thus, while the energy is not injected in the extreme UV in a nuclear collision, it is certainly not injected in the infrared. We also note that the difference in thermalization behavior is closely related to the observation that highly localized excitations in the strongly coupled gauge theory do not evolve into jets as is commonly found in weakly coupled, asymptotically free gauge theories that permit a perturbative treatment.
It is tempting to speculate to which extent our results may apply to the thermalization of QCD matter that is produced in relativistic heavy ion collisions. Phenomenologically, it is known that the thermalization has to occur fast, allowing the matter to expand according to the laws of nearly ideal hydrodynamics. If the relevant length scale for thermalization is given by the thermal scale ∼ /T , our results suggest that τ crit ∼ 0.5 /T if the matter is strongly coupled. For initial temperature value T ≈ 300 − 400 MeV at heavy ion collider energies, we obtain the estimate τ crit ∼ 0.3 fm/c, comfortably short enough to account for the experimental observations. We know that such a short thermalization time applies to all momentum scales, because jets are created in the nuclear collisions, albeit at suppressed levels. However, as already mentioned, the relevant momentum scale for energy deposition is the parton saturation scale Q s of the colliding nuclei, i.e. the transverse momentum below which most partons in the colliding nuclei are found. For achievable energies and large nuclei, Q s ≤ 3 GeV; in this range of momenta QCD may well be considered as a strongly coupled gauge theory.
Though this was not the emphasis in this paper, the rapid thermalization studied in this paper may also be of relevance for black hole physics. It has been suggested that black hole creation is the fastest possible form of thermalization that exists in nature [65] and it would therefore be of particular interest to more directly link the causality bound to this conjecture.
We also encountered phenomena that do not appear to have a simple explanation. For example, for long rectangular Wilson loops in the boundary gauge theory of AdS 4 , the area of the minimal surface exhibits a "swallow tail" behavior: the transition to the thermal limit proceeds via a discontinuous succession of shapes of the minimal surface. It is unclear why this phenomenon occurs for the rectangular Wilson loop but not for the two-point function, as the differential equation for the stationary surface is quite similar to that for the geodesic curve.
Finally, we briefly comment on the difference between dynamically collapsing shells and quasistatic shells that adiabatically approach the event horizon. As we show in the Appendix, the swallow tail behavior is ubiquitous for quasistatic shells close to the horizon, leading to significant differences between dynamical and quasistatic shells for large boundary separations.

a. Circular loops
Similar to the geodesic analysis, the four dimensional background can be divided into a black brane part (outside the shell) and a pure AdS part (inside the shell). However, we cannot solve the differential equations analytically and instead we rely on numerical methods. Outside the shell the solution satisfies the differential equation (68) in the black brane background (with d = 3), while inside the shell we must find a minimal surface in pure AdS (i.e. we set the mass M to zero in eq. (68). The minimal surface in pure AdS consists of a hemisphere. In the thin-shell limit, at the position of the shell z = z s , the minimal surface should satisfy the refraction condition This is derived analogously to the geodesic case (A2) in which case it agrees with numerical results obtained by minimizing the path of the geodesic. Note that in the dynamical case we do not use such a condition and solve the equations of motion in the whole space with mass function (4). However, here, taking advantage of this simple relation allows to use much less computationally intensive numerics in the limit of a thin shell (where the mass function is a step-function). We start by constructing a hemispherical string surface whose tip is given by z(0) = z * > z s and z (0) = 0, thus lying inside the shell. For this solution we can determine the intersection point ρ 1 with the shell and the gradient at that point analytically. We can then use the refraction condition to generate the rest of the solution in the black brane background, using the boundary condition z(ρ 1 ) = z s and z (ρ 1 ) = z out . Where this part of the solution intersects the boundary (cut-off at z 0 ) gives us the boundary radius R of the Wilson loop. We can also construct surfaces wholly outside the shell, which are just the thermal Wilson loops from section II B 1. Thus we get R(z * ) for any shell position z s . As with the geodesics case we can get multiple surfaces with the same boundary radius and shell position (such as in Fig. 22). An easy way to check whether these multiple solutions exist before calculating any areas (see also [28]) is to invert our numerical data appropriately so that we can plot it on a graph of z * against boundary time t 0 , the boundary time being obtained as a function of the shell position z s : The multiple solutions will show up as multiple possible values of z * for the same t 0 . We see this in Fig. 23(A) where we have plotted for three different boundary radii. The area of the minimal surface consists of three parts: we add the areas of the piece in empty AdS (inside the shell) and of the piece in the black brane background (outside the shell) and subtract the area of the solution in empty AdS with the same boundary radius, δA(R, t 0 ) = ρ1(t0) 0 dρ ρ z 2 1 + z 2 + R ρ1(t0) dρ ρ z 2 1 + In Fig. 23(B), we plot δA(R, t 0 ) as a function of R for a fixed shell position. The three branches appear in a very similar fashion to those for geodesics in AdS 3 shown in Fig. 22. Fig. 23(C) then shows δA(R, t 0 ) as a function of the boundary time t 0 for three different boundary radii. The appearance of the three branches gives a 'swallow tail' as found in [28] in the dynamical case for the infinite rectangular strip, and as we also find in Section III B 2. However, there are no such multiple solutions in the dynamical Wilson loop case as shown in [28] and Section III B 1. This can be taken as an indication of the limits of the quasistatic approximation in accurately capturing the behavior as the shell gets closer to the horizon (recall that to get the three branches in Figure 23(B) we have to wait till the shell approaches the would be horizon, here M 1/3 z s = 1/1.001). The appearance of three branches does not, however, signify that something is unphysical. As mentioned, they also appear in the dynamical case for the infinite strip.
What it does mean is that since we are using a saddle point approximation to calculate field theory quantities and our space time is Lorentzian, when there are multiple geodesics we should carefully follow a steepest descent procedure. The numerics are hard to pursue around the earliest time where the three branches occur in figure 23(C) and we see some noise around one corner of the 'swallow tail'. Here the tip of the sphere at z * , which is our initial condition, is very close to the shell, and very small variations in z * give rise to large variations in R and hence δA.

b. Strip
We can perform a similar analysis for the Wilson strip in a quasi-static background. The equation (76) is now the differential equation we solve numerically in the black brane background (with d = 3). Considering a similar boundary setup for the strip as in section II B 2, we follow a similar procedure to construct solutions as for the quasistatic Wilson loops. Also the refraction condition through the shell (in the thin shell approximation) is the same as for circular Wilson loops. Hence, we will not work this case out explicitly will immediately give the results. The string surface area consists of three different pieces: the part of the solution inside the shell (empty AdS) A AdS , the part of the solution outside the shell (black brane background) A BB , and a subtracting part corresponding to the full solution in empty AdS with the same boundary separation . The full area in the shell background is given by wherex represents the space-like separation of the part inside the shell. The area formula in empty AdS is given by the relation where R corresponds to the height of the Wilson strip at the boundary (see II B 2), and the black brane area is given by In Fig. 24(A), we plot the midpoint of z(x), namely z * , as a function of the boundary time t 0 . In the middle panel of Fig. 24, we plot the area δA( , t 0 ) (with R/π divided out) as a function the boundary separation for a fixed shell position z s (thus for a fixed boundary time t 0 ). On the right panel of Fig. 24 we plot the area δA( , t 0 ) (with R/π divided out) as a function of the boundary time t 0 for three fixed boundary separations. As in the dynamical case, we have three branches of solutions appearing, but here they show up already for smaller values of the radius.