Holographic isotropization linearized

The holographic isotropization of a highly anisotropic, homogeneous, strongly coupled, non-Abelian plasma was simplified in arXiv:1202.0981 by linearizing Einstein's equations around the final, equilibrium state. This approximation reproduces the expectation value of the boundary stress tensor with a 20% accuracy. Here we elaborate on these results and extend them to observables that are directly sensitive to the bulk interior, focusing for simplicity on the entropy production on the event horizon. We also consider next-to-leading-order corrections and show that the leading terms alone provide a better description of the isotropization process for the states that are furthest from equilibrium.

Part of the motivation for these investigations comes from heavy ion collision (HIC) experiments at RHIC and LHC, in which the successful phenomenological description of data requires the applicability of viscous hydrodynamics with the quark-gluon plasma equation of state and a tiny shear viscosity η/s = O(1/4π) about 1 fm/c or less after the collision [12]. While this thermalization time is puzzlingly small from the viewpoint of bottom-up weak-coupling thermalization, it lies within the ball-park suggested by holographic calculations [1,[5][6][7][8][9][10][11]. Although the dynamics in a real HIC surely involves a combination of weakly and strongly coupled physics, the hope is that understanding the limit in which the physics is assumed to be strongly coupled at all length scales may help us bracket the real-world situation -see [13] for a review of applications of holography to the QGP.
In the gauge-gravity duality the gravitational representation of the thermal state of a holographic quantum field theory in the planar limit and at strong ('t Hooft) coupling is a bulk black brane [14]. Thus thermalization on the gravity side is the process whereby a bulk geometry approaches the form of (a patch of) a static black brane and the other bulk fields relax to their thermal form in the black brane background. Generically, this stage is preceded by one which one may dub 'hydrodynamization' whereby the system reaches a phase in which the system is still evolving but the dynamics is governed by hydrodynamic equations. On the gravity side this implies that the bulk geometry must have the form dictated by the so-called fluid-gravity duality [15] (see e.g. [16] for a review). As black brane metrics already depend on the bulk radial coordinate and we want to allow for at least time dependence, the studies of holographic thermalization processes require solving Einstein's equations in two or more variables, which in most cases requires the use of numerical relativity techniques.
Apart from toy-models based on the AdS-Vaidya geometry of infalling dust (see e.g. [17][18][19][20], but also [21]), the only existing approximation scheme to study holographic thermalization processes is the amplitude expansion, in which one linearizes Einstein's equations on top of the static black brane background. In this approximation the relaxation towards equilibrium is described by quasinormal modes with complex frequencies, whose imaginary parts lead to the damping of their amplitudes with time and hence to equilibration. These modes were thought so far to be appropriate for the description of only the late-time approach to equilibrium, when deviations from equilibrium are sufficiently small in amplitude [22].
An indication that this assumption might be too restrictive comes from black hole mergers in asymptotically flat four-dimensional spacetime. There, in the so-called closelimit approximation, the Einstein's equations linearized on top of the final black hole predict rather accurately the pattern of gravitational radiation at infinity provided the initial data have a single horizon surrounding the merging black holes [23,24]. This initial horizon, however, needs not to be a small perturbation of the final black hole for the close-limit approximation to work.
These features, together with the observation that the AdS analogue of gravitational radiation at infinity is the expectation value of the boundary stress tensor, motivated recent work [1] in which a linear approximation was applied to a simple example of far-fromequilibrium gravitational dynamics in AdS spacetime. Specifically, this reference considered the isotropization of a large number of initially homogeneous but anisotropic states in a Conformal Field Theory with a holographic description, or holographic Conformal Field Theory (hCFT) for short, with a large number of colors N c and strong 't Hooft coupling [25] λ = g 2 YM N c . The outcome of this analysis was that the linearized Einstein's equations indeed reproduce well the dynamics of the one-point function of the stress tensor. For small perturbations this (expectedly) provided excellent approximations to the time evolution of the stress tensor. Surprisingly, however, even for initial states which were not a small perturbation of the final state, the linearized evolution still predicted the stress-tensor with an accuracy of about 20%. The purpose of the present paper is to explore the applicability of the linearized Einstein's equations further, in particular for observables that are directly sensitive to the interior bulk dynamics such as the entropy production.
We start this article by reviewing the results of [1]. The key technical novelty of [1], compared to earlier works, was the ability to generate and follow the evolution of a large number (around 1000) of different non-equilibrium initial states. In the current article we provide a detailed explanation of how these initial states were obtained and evolved. We also comment on the relation between features of the initial states and the time it takes them to thermalize.
Subsequently, we review the linearized gravity approach to holographic isotropization and comment on the relation between the form of initial data and the accuracy of the linear approximation. We also study in detail the decomposition of the solutions into quasinormal modes. According to the common lore, quasinormal modes can only describe the late-time approach to equilibrium, for which only the lowest-lying modes are relevant. As already anticipated in [1], our finding is that, provided a sufficient number of modes is considered, the description in terms of quasinormal modes can actually be extended all the way back to early times for which the system is very far from equilibrium. The correct description of the system at all times then results from the heavy interference of the quasinormal modes with one another.
Finally, we discuss the corrections to the linearized Einstein's equations. For generic observables that are sensitive to the geometry deep in the bulk these corrections are crucial because they constitute the leading-order result. This is the case, for example, for the entropy production, on which we focus our attention.
Throughout our paper we will assume that the boundary stress tensor is diagonal. A systematic discussion of the amplitude expansion of Einstein's equations in the more general homogeneous but non-diagonal situation can be found in [26].

Holographic isotropization and the ansatz for the dual metric
The holographic isotropization that we will study is the dynamics of a hCFT, in which a strongly coupled (3 + 1)-dimensional plasma is homogeneous, but has a time-dependent pressure anisotropy leading to rich non-equilibrium physics [1,5]. In this setup all the local operators apart from the stress tensor have vanishing expectation values, which constitutes the universal sector of any hCFT [27].
This example of a holographic thermalization process is, most likely, the simplest one, as due to homogeneity it does not involve hydrodynamic modes and thermalization is achieved directly via relaxation to a thermal state without being preceded by a long hydrodynamic phase. The latter observation is the primary reason for considering it here, and previously in [1], in the context of the close-limit approximation for anti-de Sitter spacetimes.
As in [1], we will be interested in an unforced relaxation to the thermal state and we consider the situation in which there are no time-dependent sources pumping energy into our system, i.e. the boundary metric is flat. This is the crucial technical difference between our work and that of reference [5], which considered the holographic isotropization process in which non-equilibrium states were obtained by placing the boundary theory in a time-dependent metric for a finite period of time.
The form of the bulk metric is determined to a large extent by the ansatz for the boundary stress tensor. After imposing for simplicity rotational symmetry in two of the spacelike directions and taking into account the flatness of the boundary metric, the most general conserved and traceless stress tensor can be written as 1) where N c is the rank of the SU (N c ) gauge group and E is proportional to the energy density. The longitudinal P L (t) and transverse pressures P T (t) are then expressed in terms of a time-dependent anisotropy ∆P(t) through Note that the energy density in our setup is constant in time by virtue of the homogeneity assumption plus the conservation of the stress tensor. In this sense the energy density is part of the initial conditions. As the only possible static state with finite energy density is the isotropic and homogeneous plasma [28], the final state is known already from the start, without the need to solve any dynamical equation. This seems to be a rather non-generic feature of our setup, which we discuss at length in the conclusions section. The symmetries of the stress tensor -through the general near-boundary form of the metric in the Fefferman-Graham coordinates [29] -suggest the ansatz for the dual metric g ab with g tt , g LL and g TT components being dynamical, i.e. obtained by solving the equations of motion. 1 The freedom of defining what is meant by the bulk time t and the radial coordinate in AdS r has to be then (partly) fixed by imposing the form of the g rr and g tr components.
In our analysis it will be very convenient to follow [5] and express the dual metric g ab in the generalized ingoing Eddington-Finkelstein coordinates where A, Σ and B are functions of time t and the radial coordinate r. With this ansatz g rr = 0 and g tr = 1. In the ingoing Eddington-Finkelstein coordinates, constant time slices are null hypersurfaces as radial ingoing null geodesics, by construction, penetrate the bulk instantaneously.
The metric (2.4) has to solve the Einstein's equations with the negative cosmological constant R ab − 1 2 R g ab − 6 L 2 g ab = 0, (2.5) 1 It does not have to be necessarily the case and actually in the ADM formulation of general relativity it is the form of the lapse -the timelike warp factor (gtt) -and the shift (gta) that are fixed/imposed. For details in the context of AdS spacetimes see [9].
where in the following we set the radius of the AdS solution L to 1. The link between the form of the field theory stress tensor and the dual metric ansatz becomes clear after solving Einstein's equations in the near-boundary (large r) expansion and identifying normalizable modes with the components of the stress tensor using holographic renormalization [29]. For the specific case of SU (N c ) N = 4 super Yang-Mills theory this relation is The metric ansatz (2.4) leaves the residual diffeomorphisms freedom r → r + f (t). This freedom is fixed here by imposing the absence of the term proportional to r in the nearboundary expansion for the g tt metric component (2.6a). In (2.6) we suppressed the near-boundary expansion at relatively low order, but it is important to stress that the expansion has infinitely many terms carrying arbitrarily high derivatives of the pressure anisotropy. This inevitably leads to a general conclusion that a state given by the form of the geometry on a constant time slice is (partly) specified by infinitely many derivatives of the dual stress tensor, in our case the pressure anisotropy.

Thermalization criterium
Although E is constant in time, the physical temperature can only be assigned to the system once the equilibrium is reached. In this regime E = 3π 4 T 4 /4 and the metric takes the form and describes the AdS-Schwarzschild solution between the boundary and the future event horizon covering also the black brane interior. Although equilibration of a holographic system can be sampled with different field theory probes, including expectation values of local operators, two point functions, entanglement entropy and Wilson loops, in our studies we will primarily focus on tracing the evolution of the one-point function of the stress tensor. There are two reasons for this. In the first place this is the quantity of interest if one wants to make a phenomenological contact with the fast applicability of hydrodynamics at RHIC and LHC. Secondly, after the stress tensor eventually settles down to its thermal value, the geometry becomes a patch of the AdS-Schwarzschild black brane and from this moment on there is no need to evolve the Einstein's equations further.
We will hence define thermalization time as the isotropization time t iso , i.e. the time after which the stress tensor anisotropy ∆P(t) remains small compared to the energy density and eventually decays to zero. In our calculations, as in [1], we adopt the following criterium for t iso : but it is important to keep in mind that thermalization is never a clean-cut event and the threshold on the RHS of (2.10) can be always slightly raised or lowered without altering much the results.

The event horizon and its entropy
The event horizon is defined as the causal boundary of the black hole. It is a teleological object which can be located only after the black hole settles down to the state of permanent equilibrium.
For the purposes of the current paper, we will be interested in the event horizon's area as an example of one easy-to-compute bulk observable that is directly sensitive to the form of the geometry in the deep IR. Although no precise definition of the entropy density exists in a truly far-from-equilibrium situation, the change in the area density of the event horizon provides a crude measure of the total entropy produced in the thermalization process. For this reason we will loosely refer to the area density of the event horizon as 'entropy density', but we emphasize from the start that this terminology is rigorously applicable only near equilibrium. We also clarify that we chose to focus on the event horizon, as opposed to e.g. the apparent horizon, because in many cases there was no apparent horizon on our initial-time slice.
Due to the symmetries of our problem, the event horizon will be a hypersurface defined by r − r EH (t) = 0, (2.11) with the normal vector being null The latter is the geodesic equation for the outgoing light ray and needs to be supplemented with the following condition to be imposed in the asymptotic future In practical terms this condition implies that when the metric eventually approaches the form of the AdS-Schwarzschild black brane (2.8), r EH approaches the position of the event horizon of the static solution.
The area of the event horizon gives rise to the following expression for the entropy which is guaranteed to be a non-decreasing function of t. In (2.14) we implicitly assumed that the event horizon's area is mapped onto the boundary along ingoing null radial geodesics, i.e. along lines of constant t for the metric ansatz (2.4).
3 Far-from-equilibrium initial states and their nonlinear evolution 3.1 Solving Einstein's equations as an initial-value problem As originally noted in [5], the Einstein's equations for the metric ansatz (2.4) can be neatly written as denote respectively derivatives along the ingoing and outgoing radial null geodesics. In the following we will be interested in solving the initial-value problem, i.e. given the geometry on the initial-time slice we want to obtain the evolution of the dual stress tensor by computing the bulk spacetime outside the event horizon. Not all equations among (3.1) are evolution equations, i.e. specify the form of the metric on a neighboring time slice. Equations (3.1d) and (3.1e) are constraints in the sense that the remaining components of the Einstein's equations can be shown to guarantee that they are obeyed provided (3.1d) holds at the boundary and (3.1e) holds on the initial-time slice [5].
The null character of the coordinate system (2.4) leads to a nested algorithm for solving the initial-value problem in which one uses as evolution equations (3.1a)-(3.1c) and at each time step one only needs to solve linear ordinary differential equations in r. The precise scheme that we will follow is a slight modification of the one originally introduced in [5], and consists of the following steps: 1. we start with B as a function of r and the energy density E (constant in our setup); 2. the constraint equation (3.1e) allows us to solve for Σ as a function of r; 3. we then solve (3.1a) forΣ, with E being the integration constant; 4. having B, Σ andΣ, we solve (3.1b) forḂ; 5. with B, Σ,Ḃ andΣ at hand we can integrate (3.1c) for A; 6. knowingḂ and A and using (3.2) we get ∂ t B; 7. we proceed to the next time step using a finite difference scheme (for details see section 3.3).
In our setup the constraint (3.1d) is implemented in the near-boundary expansion, as it encodes the conservation of the stress tensor in the dual gauge theory. In order to monitor the accuracy of the numerical code we check the value of this constraint in the bulk when evaluated on the numerical solution (see also section 3.3). The algorithm outlined above needs to be supplemented with the initial conditions, the choice of which we discuss in the next subsection.

Specifying initial states
Gravity encodes dual initial states in the form of the geometry on a bulk initial-time slice. The conditions on the initial data arise from three sources: the constraint (3.1e), the nearboundary expansion (2.6) and bulk regularity. By the latter we mean that all possible singularities in the initial data must be hidden inside the event horizon.
One way to obtain a non-equilibrium state while automatically satisfying the conditions above is to start with vacuum AdS and perturb it by turning on a non-normalizable mode of the bulk metric or some other bulk field for a finite period of time [5]. The alternative approach, that we adopt here and which was used also in [8,9,30], is to start with nonequilibrium states defined as solutions of the constraints on the initial-time slice without invoking the way in which a particular state was created.
Equation (3.1e) imposes a constraint between the forms of B and Σ on the initialtime slice. Since B appears quadratically, we choose to specify the initial state through B and then use (3.1e) to solve for Σ. Note that this equation, together with the asymptotic behavior linear in r (2.6c), implies that Σ must necessarily be a convex function and hence that it must vanish for some r ≥ 0, with the inequality being saturated only for vacuum AdS and the AdS-Schwarzschild black brane. As our coordinate frame is spanned by the ingoing radial null geodesics and Σ measures the transverse area of the congruence, Σ = 0 implies reaching a caustic and hence the breakdown of our coordinate frame.
For the successful evolution of the initial data specified by some B we thus need to make sure that the locus where Σ vanishes is hidden behind the event horizon on the initial-time slice. As the event horizon is a teleological object, this cannot be verified a priori -we need to try to run a simulation and when it is successful we know that the initial state we started with was legitimate.
The contrary is not necessarily the case, as a caustic, a priori, is just a breakdown of a coordinate system. However, we verified numerically that in the neighborhood of a point where Σ vanishes we obtain very large curvatures. This suggests that this point must be hidden inside the event horizon. 2 We thus see there is an interesting interplay between the choice of B and the choice of the (initial) energy density E. Both quantities, a priori, seem to be very much independent when it comes to specifying the initial state. If, however, the point where Σ vanishes corresponds to a genuine curvature singularity, which is the case for the AdS-Schwarzschild black brane and which our numerical studies also indicate, then there must be a minimal energy density E for which this singularity is still covered by the event horizon on the initial-time slice. An indication that this might be the case comes from noting that for the AdS-Schwarzschild black brane the radial position of the event horizon is proportional to (the fourth root of) the energy density E. Similarly, in our setup we can always put the final event horizon at a smaller value of r than the point at which Σ vanished on the initial-time slice by decreasing the energy density. This discussion suggests that the states for which the initial position of the event horizon is very close to the point where Σ = 0 should be viewed as maximally far from equilibrium.
In our setup, we have a freedom of preparing arbitrary initial conditions, i.e. we can specify B as a function of r on the initial-time slice and E > 0, as long as B obeys the near-boundary expansion of the form (2.6b) and there are no naked singularities. We use this freedom to prepare and follow the evolution of states in which B has support very close to the boundary, very close to the horizon or spreads over a large range of the radial direction. In order to generate a large number of non-equilibrium initial states we followed the following procedure: 1. without loss of generality we choose units such that a 4 = −1, or equivalently E = 3 4 ; 2. we generate the initial B as a ratio of two 10th order polynomials in r with random coefficients in the range (0, 1); 3. we subtract from it a cubic expression so that the near-boundary expansion for B of the form (2.6b) is obeyed; 4. the whole expression is then normalized so that the maximal value of the B between the boundary and the position of the final event horizon (r = 1) is 1 2 ; 5. we then run a binary search algorithm to find the factor that B needs to be multiplied by such that the code is just stable, while storing successful runs. Typically, we repeat this step about 6 times per seed function generated in step 2.
In this way we can generate states which are as far from equilibrium as our numerical code allows. In the end this means there is some sensitivity to the number of grid points, since increasing the number of grid points would improve the stability.
Finally, it is interesting to note that a constraint of exactly the form (3.1e) also holds for metric ansatzes corresponding to a dual plasma expanding in one dimension [6,7]. This implies that our discussion about the specification of the initial states, including the maximally far from equilibrium ones, also applies in these other setups. However, if we relax the assumption of a homogeneity in the transverse plane, then Σ is no longer forced to be convex and there might be bulk states which do not lead to caustics/apparent singularities in the way described above.

Numerical implementation
In the numerical implementation instead of the variable r we used its inverse so that the AdS boundary is at z = 0. With the choice of units we made in our code, i.e. a 4 = −1, the horizon of the final black brane is then located at z = 1. Note however that, for definiteness, all dimensionful quantities that we will provide will be specified in terms of the energy density or, equivalently, the temperature of the final black brane, which is the only scale at the moment of thermalization. For discretization we use the pseudo-spectral collocation method for the spatial grid in the z direction (see [31,32] for accessible reviews of the spectral methods in the context of numerical GR). This allows us to maintain a very moderate grid, typically with 26 points, without a significant loss of accuracy. Each of the evolution equations (3.1a)-(3.1c) is written in terms of redefined variables in which the three terms of lowest order in the near-boundary expansion are subtracted, and solved as a set of linear equations.
The latter is a significant simplification, since it allows the use of standard procedures for solving linear equations. We implemented our numerical code in Mathematica and used the fourth order Adams-Bashforth method for evolution in time. The first five time steps were made with the use of low order Runge-Kutta methods. Typically, we used time steps of size 0.0025 ∼ (# grid points) −2 .
As a way of monitoring the accuracy of our code, we used the normalized constraint (3.1d) The convergence of our code is then illustrated by Fig. 1, which shows typical plots of the maximum value of the normalized constraint ξ(t).
The last feature that needs to be discussed is the choice of the position of the inner boundary of the computational grid. Note that the simulation is well-defined only if the grid covers the entire portion of the spacetime outside the event horizon. Initially this is hard to predict, since the position of the event horizon depends on the future evolution. Therefore one typically focuses attention on the presence of the apparent horizon because, if it can be found, it is guaranteed to lie inside the black hole. However, quite frequently in our case there was no apparent horizon on the initial-time slice and therefore we used the following procedure. We first tried to run simulations with the radial cut-off put at z = 1.01, which is right below the late-time position of the event horizon. This often worked, and when it did not we reran the simulation with z = 1.07 as a cut-off. The latter point turned out to almost always lie past the event horizon. In this way we could successfully evolve a large number of initial states.

Why the near-boundary expansion is not enough
The near-boundary expansion of the metric contains information about all time derivatives of the pressure anisotropy at the initial time. This allows the construction of a Taylor-series expansion for the pressure anisotropy, which at least for some time should be very close to the exact expression. Previous studies in the context of boost-invariant flow suggest that the convergence radius of a series of this type is too small to see thermalization [9]. We can convince ourselves that this is the case also here by considering the initial profile of the form According to the near-boundary expansion (2.6b) this profile at the initial time has a nonzero pressure anisotropy, but all its time derivatives vanish. The Taylor series expression predicts then a constant pressure anisotropy, which would lead to a singular spacetime [28]. Fig. 2 shows the numerical evolution of the profile (3.5) and proves that the radius of convergence of the Taylor expansion is indeed too small to see thermalization defined via (2.10). This result implies that the precise form of the metric g ab and the pressure anisotropy ∆P(t) at transient times need to be obtained via an explicit solution of the initial-value problem on the gravity side.

Evolution of a sample profile and expectations for thermalization times
To get intuition about how the dynamics proceeds on the gravity side and to get acquainted with the features following from the choice of a foliation by null constant-time slices, it will be instructive to discuss in detail the dynamics of the following initial state where z h = 2 1/2 3 1/4 E 1/4 . As B is supported at intermediate values of z, naive intuition from the physics of linear wave equations would suggest that the wave packet splits into two: one propagating inwards and the other propagating outwards. The one propagating outwards Figure 2. The left plot shows B(z, t) for the initial profile (3.5), which is shown as a thick red line at t = 0. By construction, all time derivatives of the pressure anisotropy vanish at t = 0. The thick blue curve at z = 0 shows the value of the gauge theory quantity ∆P(t)/E. One sees clearly that the Taylor expansion around the initial time (predicting constant pressure anisotropy) does not converge when thermalization is achieved. The right plot, in which we start with B(t = 0, z) = ( 4 3 E) 6 z 24 , shows similar behavior. The initial disturbance, which is localized in the IR part of the geometry, propagates to the boundary in a time limited by causality. This creates the pressure anisotropy, which quickly relaxes back to zero.
is expected to eventually reach the boundary, bounce back and fall into the bulk. Both wave packets will be eventually absorbed by the event horizon (which is guaranteed to be present given that E = 0) leading to the increase in its area.
These expectations are confirmed by the outcome of the numerical simulation, as illustrated by Fig. 3, which depicts the bulk anisotropy (left plot) and the square of the Riemann tensor, the Kretschmann scalar (right plot). We can clearly see the rise in the curvature due to the outgoing wave packet as it approaches the boundary of AdS. Closer inspection reveals also the presence of a wave packet resulting from the bouncing off the boundary of the outgoing packet. This wave packet, due to the null nature of our coordinate frame, propagates towards from the boundary to the horizon along lines of constant Eddington-Finkelstein time. Note also that this signal falls through the black brane event horizon without significant scattering. This feature persisted for other choices of initial states and seems to be related to the high degree of symmetry of our problem.
It is interesting to note that the initial ingoing part of the wave packet seems to be mostly taken care of by the solution of the constraints. Indeed, although B is supported only over some small range of z centered around z h /3, the metric functions A and Σ deviate from their vacuum values all the way from this point to the horizon, as required by causality. In contrast, the curvature outside the outgoing wave packet is very close to the curvature of the static black brane.
These observations suggest that the states which take the longest time to thermalize are those that are initially localized close to the horizon on the initial-time slice. An example is provided by B(t = 0, z) ∼ z 24 , whose evolution is shown in Fig. 2 (right). The reason is that the outgoing wave packet needs to escape the neighborhood of the horizon and travel all the way to the boundary to bounce off and finally fall into the black brane  .6), which is shown as a thick red curve at t = 0 and which is initially localized near z = 1 3 z h . The blue curve at the boundary (z = 0) depicts the pressure anisotropy as a function of time in the gauge theory. The right figure shows the Kretschmann scalar (with the value for an equilibrium black brane with the same energy density subtracted) as a function of time and radial coordinate for the same initial profile. One clearly sees on this plot the wave bouncing off the boundary and falling into the black brane. In the adopted generalized ingoing Eddington-Finkelstein coordinates this happens instantaneously.
horizon. By localizing the initial profiles close to the horizon, the longest isotropization times that we were able to obtain with our numerics, which used rather moderate grids, were about 1.1/T − 1.2/T , with T the final equilibrium temperature (see Fig. 5).

Leading order correction to the pressure anisotropy
Linearizing Einstein's equations in the setup of holographic isotropization can be formally phrased as an expansion in the amplitude of perturbations on top of the AdS-Schwarzschild black brane. We thus write where α is a formal parameter counting the order in the amplitude expansion. The smallness of the initial data can be physically quantified by either measuring the total entropy production on the event horizon or by following the amplitude of the pressure anisotropy during the evolution process and comparing it to the energy density. It is important to restress that we want to use the linearized approximation without necessarily restricting to the initial data being small perturbations of the AdS-Schwarzschild black brane, precisely in the spirit of the original close-limit approximation [23,24] but now in the context of AdS spacetime.
The initial data for the full nonlinear Einstein's equations are given by specifying the energy density E and the form of B as a function of the radial coordinate on the initial-time slice. As anticipated earlier, the main motivation for choosing B over Σ in specifying the initial data was that the former appears quadratically in the constraint (3.1e). This feature persists also with the other components of the Einstein's equations apart from the equation (3.1b), which immediately leads to  The energy density E, which is constant in our setup and is the remaining part of the initial state specification, is already included in the background we linearize on top of. In full detail, the equation for δB (1) This is a simple, first order in time, partial differential equation which can be solved once supplemented with the Dirichlet boundary condition at z = 0. To improve the stability in our computations we actually worked with and required that δB Equation (4.4) captures the leading order dynamics of the whole stress tensor, as the stress tensor is completely specified in terms of the energy density E (encoded in the background) and the pressure anisotropy ∆P(t) (encoded in δB (1) ).
In our previous work [1] we compared the predictions for ∆P(t) following from the nonlinear Einstein's equations with the ones obtained using (4.4) and we found that the results agreed with a 20% accuracy for the vast majority of the profiles we considered. The match of the qualitative features can be seen in Fig. 4. The histogram in Fig. 5 shows the quantitative match for the isotropization time for more than 800 non-equilibrium initial states which are all large perturbations of the AdS-Schwarzschild black brane. We see that the linearized Einstein's equations predict the isotropization time with a 20% accuracy in the vast majority of cases.
An interesting feature visible in Fig. 4, as well as earlier in Fig. 3, is that the pressure anisotropy obtained from the linearized Einstein's equations always follows closely the full result at early times and, if at all, the curves start to differ only at some transient time. We understand this feature as a natural consequence of the fact that, due to the fixed asymptotics, the dynamics is approximately linear near the boundary of AdS, which is the bulk region causally responsible for the early-times dynamics. The pressure anisotropy curves start to differ only after the boundary is reached by a signal propagating from the interior of the geometry, as seen in Fig. 3, and this is the moment at which in the full analysis nonlinear effects become most visible. Our main result can be phrased as the The total number of initial states is more than 800. We see both that holographic isotropization proceeds quickly, at most over a time scale set by the inverse temperature, and that the linearized Einstein's equations correctly reproduce the isotropization time with a 20% accuracy in most cases. Note that the histogram is based on a different sample of initial states than those originally considered in [1]. In particular, we incorporated the binary search algorithm absent in [1] and were stricter about the maximum violation of the constraint that we allowed. (Botom) Close inspection of one of the few profiles for which the linearized approximation seemingly fails by more than 20% (∆t iso /t iso = −0.5) shows that it is the imperfect isotropization criterium which leads to the mismatch rather than the failure of the linear approximation. Indeed, the left plot shows that, on the scale of the initial anisotropy, the linear result yields a good approximation. However, the isotropization criterium makes no reference to this scale, and results in a 50% difference in the isotropization times, indicated by the arrows on the right plot. See [9] for a related discussion of subtleties involved in defining the thermalization (or more accurately hydrodynamization) time in a similar setup.
surprising statement that these linearities did not result in a large effect on the boundary stress tensor for any of the initial states that we analyzed.

Connection with quasinormal modes
Equation (4.4) can be solved either as an evolution equation given some initial profile for δB (1) , as discussed in the previous section, or by decomposing δB (1) as a superposition of modes with factorized time dependence: These modes are known as quasinormal modes, and they are characterized by the requirements that they be normalizable near the boundary (z = 0) and that they obey ingoing boundary conditions at the event horizon (z = 1). 3 The latter condition makes the frequencies ω j complex with imaginary parts responsible for the exponential decay in time. The quasinormal modes (4.6) appear in pairs, as taking the complex conjugate of the equation (4.4) for the quasinormal mode with frequency ω j leads to the equation for the quasinormal mode with frequency −ω * j . This feature can be seen in Fig. 6. Im Ω z h Figure 6. The plot on the left shows the frequencies of the ten lowest quasinormal modes including their complex conjugates. The mode with the smallest negative imaginary part will by the dominant mode at late times. Notice that the spacing between the modes is approximately constant (it differs by about 0.1%). The plot on the right displays the lowest three quasinormal modes as a function of the radial coordinate z, where blue and red denote their real and imaginary parts. The normalization we use is such that the real part at the horizon (z/z h = 1) is equal to unity, whereas the imaginary part vanishes there. One clearly sees that higher modes (which decay faster) are more dominant near the boundary.
In the context of gravitational collapse, the lowest quasinormal modes are known to govern the late-time decay of black hole perturbations (see e.g. [34]) and this is also expected in the current setup. On the other hand, the results from [1], reviewed in the previous section, suggest that the equation (4.4) predicts rather well the full time dependence of the large-z behavior of the warp factor B. Hence a natural question is what is the quasinormal mode content of the perturbations that we considered.
In order to answer this question we followed the prescription of [22] and computed the lowest 10 quasinormal modes (4.6) by solving equation (4.4) for the ansatz (4.6) in the near-horizon expansion and evaluating the resulting expression at the boundary to find ω j 's leading to normalizable modes. The (somewhat arbitrary) normalization of our modes is fixed by demanding that at the horizon (z = 1) b j (1) = 1. (4.7) On Fig. 6 we plot the obtained frequencies ω j of the lowest 10 quasinormal modes, as well as bulk profiles for the real and imaginary parts of b 1 (z), b 2 (z) and b 3 (z) normalized according to (4.7).
The idea now is to use the quasinormal modes to decompose solutions of (4.4), i.e. to write a solution of (4.4) in the form where we truncated the expansion at some N QNM , although formally we could set N QNM = ∞.
In our calculations we used N QNM = 10. One can view (4.8) as a further simplification as compared with solving numerically (4.4), which approximates very well (within 20%) the full Einstein's equations when it comes to predicting the form of the dual stress tensor. The reason for this extra simplification is that now the solution is specified by providing a few complex numbers 4 (say 10 complex coefficients c j 's) which due to the linearity of the problem can be fitted on the initial-time slice to B(t = 0, z).
As a way of generating coefficients c j 's we minimized  by using the least squares method on a discrete sample of the radial position z. Naturally, one needs far more points than the number of quasinormal modes included in (4.8).
The subtlety in using (4.9) lies in the choice of the mutiplicative factor under the integral, which we set to be 1/z 3 . We checked that both 1/z and 1/z 4 do not work well, as the first one does not take sufficiently into account and the other overcounts the nearboundary behavior of B(t = 0, z). On the other hand, 1/z 2 seems to work equally well as 1/z 3 , but for definiteness we focused here on the latter. Fig. 7 displays the difference between B(t = 0, z) and δB (1) QNM (t = 0, z) as a function of the number of quasinormal modes in two representative examples. Clearly, if a good fit is possible, then the profile (4.8) will solve the linearized Einstein's equations nicely since each quasinormal mode solves them individually. 4 One may construct exceptional initial profiles, which are for instance very close to the boundary, or very rapidly oscillating. Including more quasinormal modes (taking NQNM in (4.8) somewhat bigger than 10) would allow us to treat these cases more accurately.
In Fig. 8 we compare the linearized evolution obtained from a direct solution of (4.4) and from a solution based on a decomposition into quasinormal modes. One can see that the contribution from each individual quasinormal mode can be large, but that the final sum approximates the linearized evolution very well. Finally, in Fig. 9 we plot three representative examples, where the profile with B(t = 0, z) having support mostly in the IR displays this interference phenomenon particularly nicely. The plot in the middle shows the error for the same profile as a function of the bulk coordinate z while using the 10 lowest quasinormal modes. The right plot displays the error for B(t = 0, z) = z 25 and shows clearly that a profile which is dominated in the IR is much harder to fit by the quasinormal modes. This causes oscillations in the evolution, as can be seen in Fig. 9. Figure 8. On the left one sees the pressure anisotropy ∆P/E as predicted by the linearized evolution, or indistinguishably by the sum of the lowest 10 quasinormal modes as a thick blue line. One can also see there the sum of the first 1 (blue), 2 (green), 3 (orange) or 4 (red) quasinormal modes. As becomes apparent, the late time dynamics is well approximated already by keeping only the lowest quasinormal mode, but if one uses more the fit starts matching earlier. Note that the coefficients are computed such that the sum of the 10 fits the initial state best. On the right we plot the individual quasinormal modes with the same coloring. One clearly sees that each of them carries very large anisotropy, but that their interference matches the linearized solution. Figure 9. In this figure we illustrate anisotropies of the full (blue), linearized (green) and quasinormal mode (red) evolution of three representative initial profiles, located in the IR, spread-out and located in the UV respectively. Clearly, the initial profile located in the IR takes some time before exciting the anisotropy at the boundary, which also explains the late thermalization. The UV profile can have a very large anisotropy, but isotropises very fast. These features are nicely described when looking at the quasinormal modes coefficients c j 's (below, real (blue) and imaginary (red) part). For the IR profile each individual contribution is very large, but they interfere in such a way to give only moderate anisotropies. In this way it is also possible to reach isotropization as late as 6 times the lowest QNM e-folding time. We also see that one would need to compute more quasinormal modes to accurately fit this profile.

Full leading order correction to the bulk
As reported in [1] and briefly reviewed in section 4.1, linearizing the Einstein's equations leads to a very significant simplification in our setup when it comes to predicting the dynamics of the dual stress tensor, which is the only excited local operator. A natural question arising is whether also the dynamics of nonlocal observables, such as Wilson loops or the entanglement entropy, can be determined in a satisfactory way by using the linearized Einstein's equations alone. As nonlocal observables probing sufficiently large domain in a dual field theory obtain direct contributions from the interior of bulk spacetime, this question translates to which extend the linearized gravity reproduces the full bulk.
Already at the superficial level it is easy to see that in order to answer this question we need to extend the analysis from [1] and include quadratic corrections to linearized Einstein's equations. A simple way to understand it is to notice that at the linear level both A and Σ retain their background form and hence the equation (4.4) does not capture the phenomenon of the entropy production. The latter appears only when A or Σ are different from their background values, as follows from equations (2.12) and (2.14).
The lack of entropy production at first order in the amplitude expansion does not matter for initial states for which the initial condition B(t = 0, z) is tiny. In this case the quadratic corrections will be even tinier and as negligible as the actual entropy production in this particular process. However, we already saw that the linearized Einstein's equations do a good job in predicting ∆P(t) also in the cases when the perturbation is not small by any means. In these situations though, they would do a miserable job in predicting the entropy production.
This argument implies that actually the full leading order correction to the bulk spacetime in the current setup comprises δB (1) , which captures the dynamics of ∆P(t) in the leading order, and δA (2) and δΣ (2) . The latter quantities are directly responsible for the leading order entropy production in the close-limit approximation in our setup.
From a more general perspective, we will treat the entropy production as an example of a quantity that is highly sensitive to the IR part of the bulk geometry and we expect that the corrections δA (2) and δΣ (2) will be crucial for computing other observables of that sort, including two-point functions, entanglement entropy, Wilson loops, etc.
By looking closer at the full set of nonlinear Einstein's equations (3.1) one can easily anticipate that δΣ (2) can be obtained by solving the constraint (3.1e) expanded up to the quadratic order in the amplitude This is a second order linear ordinary differential equation in z, which, having obtained δB (1) (t, z), can be solved using very basic numerical methods on each time slice. The boundary condition for δΣ (2) is that it vanishes at z = 0, as the boundary behavior is always incorporated in the background.
Having both B (1) (t, z) and Σ (2) (t, z) we can use equation (3.1c) expanded to the quadratic order in amplitude to solve for δA (2) (1 − z 4 )(δΣ (2) + z∂ z δΣ (2) ) + 12∂ t δΣ (2) . (4.11) Again, one obtains a second order, linear, ordinary differential equation in z, which can be integrated with the use of basic numerical techniques. The choice of integration constants is dictated by the near-boundary expansion (2.6), which at low orders is the same for the linearized and the full Einstein's equations. At the second order of the amplitude expansion the correction to B vanishes. The reason is that equation (3.1b) at this order is homogeneous and first order in time, and we already used the initial profile for B as the initial condition for δB (1) so that B (2) (t = 0, z) = 0. Finally, it is worth noting that in order to improve the accuracy of our procedures, in numerical calculations we actually redefined δΣ (2) and δA (2) to be δΣ (2) reg = z −5 δΣ (2) and δA (2) reg = z −4 δA (2) . (4.12) Although formally we refer to our expansion scheme as 'the amplitude expansion', in the spirit of the close limit approximation we are actually interested in evaluating the whole leading order contribution to the geometry (δB (1) , δΣ (2) and δA (2) ) for initial data of arbitrary amplitude and in comparing it with the full answer.

Leading order correction to entropy production
Having obtained the geometry at leading order, i.e. having solved equations (4.4), (4.10) and (4.11), we use (2.12) and (2.14) to calculate the position of the event horizon as well as its entropy.
Contrary to the nonlinear case, we are not guaranteed that the entropy will be a nondecreasing function of time. An easy way to understand it is by noting that the solution of linearized equations corresponds to the solution of full equations with some 'matter' stress tensor in the bulk. As this matter is not guaranteed to obey the energy conditions, the approximated entropy may in principle decrease. Moreover, the linearized formula for the area of the event horizon following from e.g. (2.14) is not guaranteed to lead to a nonnegative result.
It turns out, however, that most often these two caveats do not apply, and even if they do, they do not change the results significantly when compared with the full answer, as can be seen in Figs. 10 and 11, which are the entropy analogues of the anisotropy Figs. 4 and 5. The results depicted in Figs. 10 and 11 show that the leading order correction consisting of δB (1) , δΣ (2) and δA (2) reproduces a bulk IR-sensitive observable such as the entropy production with a 20% accuracy for the vast majority of states. The numerical techniques needed for obtaining these leading order approximation are much simpler than needed to obtain the nonlinear results. In particular, the built-in Mathematica command NDSolve is sufficient. This is already a significant simplification, as the actual procedures used to solve the nonlinear problem required a fully-fledged numerical implementation.

Including higher order corrections
Having computed the full leading order correction, a natural question is whether including higher order terms in the amplitude expansion will get us closer to the full nonlinear result. A direct calculation shows that the first subleading correction to the pressure anisotropy comes from δB (3) , whereas the first subleading correction to the entropy production comes from the combined effect of including δA (4) and δΣ (4) .
We started with some initial B(t = 0, z) and, after each run, we increased its amplitude by fixed amount. For definiteness, in Figs Figure 10. Comparison between the leading order prediction for the entropy production coming from the Einstein's equations expanded to quadratic order (red curves) and the exact entropy production the obtained from the full nonlinear analysis (blue curves), for the same 15 initial states as those in Fig. 4 with the same ordering. Note that none of the initial conditions can be regarded as a small perturbation either in the sense of the amplitude or the total entropy produced during the equilibration process. On Fig. 11 it can be seen that the leading order prediction for the entropy production from linearized Einstein's equations matches remarkably well, within 20%, with the full nonlinear result for the vast majority of the initial states that we analyzed. by increasing the amplitude by 1 3 E at each run. As expected, the leading order approximation predicts well both the pressure anisotropy and the entropy production for all amplitudes, including the largest ones that we could simulate. Also unsurprisingly, for small enough amplitudes the subleading correction improves the result as compared to the lead- Figure 11. Histogram for the error in the entropy production analogous to that in Fig. 5 for the isotropization time. Both histograms are based on the same sample of more than 800 initial states. We define the error in the entropy production as ∆Entropy = ∆s (2) −∆s ∆s , where ∆s = s(t = ∞) − s(t = 0) and ∆s (2) = s (2) (t = ∞) − s (2) (t = 0) are the exact and the leading order entropy productions, respectively. We see that the leading order approximation reproduces the correct result with a 20% accuracy for the vast majority of states.
ing order result. However, for sufficiently large amplitudes the expansion does not seem to converge. In the case of the pressure anisotropy, including the subleading correction overshoots the prediction for the initial states with the largest amplitudes. In the case of the entropy density, including the subleading correction for large-amplitude profiles leads to rather unphysical results since the entropy density decreases over certain periods of time. A further indication that adding subleading corrections does not improve the result for generic states comes from Fig. 14, where we see that the histogram for ∆t iso /t iso does not become narrower along the ∆t iso axis after including the third order correction for ∆P(t).
With the prospect of studying other more complicated setups, all these results point to the conclusion that it is safe, and generically better, to keep only the leading order terms.

Summary and open problems
In this article we presented a thorough study of the simplest example of holographic thermalization process, the holographic isotropization. The primary motivation for these investigations was understanding to which extent the full process of holographic isotropization can be simplified by describing it using the linearized Einstein's equations alone. Such an  Figure 12. Pressure anisotropy as a function of time as predicted by the full Einstein's equations (blue curves), by the first order approximation (red curves) and by the first+third order approximation (green curves). In all cases the initial state is of the form (4.13) with an amplitude that increases by 1 3 E from one plot to the next (from left to right and from top to bottom). We see that for small (large) amplitudes adding the third order result improves (worsens) the agreement with the exact result.  Figure 13. Entropy of the event horizon as a function of time as predicted by the full Einstein's equations (blue curves), by the second order approximation (red curves) and by the second+fourth order approximation (green curves). In all cases the initial state is of the form (4.13) with an amplitude that increases by 1 3 E from one plot to the next (from left to right and from top to bottom). We see that the second order prediction works better for all cases shown. approximation is an AdS generalization of the closed limit approximation for asymptotically flat spacetimes, in which in certain black hole mergers one linearizes the Einstein's equations on top of the final black hole predicting quite accurately the gravitational wave pattern at infinity. Our previous studies in [1] showed that such generalization of the close limit approximation works remarkably well, i.e. with a 20% accuracy, when it comes to predicting the time dependence of the expectation value of the boundary stress tensor. Figure 14. Same histogram as in Fig. 5 except that the first order approximation for the pressure anisotropy is replaced by the first+third order approximation. Comparing with Fig. 5 we see that the histogram does not become narrower along the ∆t iso axis after including the third order correction.
The main phenomenological motivation for studying holographic thermalization is learning possible lessons about the way the thermalization (or rather hydrodynamization) process proceeds in relativistic heavy ion collisions at RHIC and LHC. By replacing QCD by a theory with a gravity dual one only expects to obtain either qualitative insights or quantitative ball-park estimates [35]. In this sense a 20% accuracy is more than what is needed in order to understand the phenomena we are interested in, and at the same time may allow to address otherwise technically hard-to-tackle questions. Two examples of such problems in the relativistic heavy ion collisions context are the pre-equilibrium phase of the elliptic flow and the equilibration of transverse-plane inhomogeneities following from event-by-event fluctuations. Solving their holographic analogues in full generality will require complicated simulations of AdS-black hole spacetimes depending on all five bulk coordinates and our hope is that a suitably developed linear approximation may allow us to obtain results with a reasonable accuracy at a much smaller cost.
In addition to elaborating in detail on our previous studies, the chief aim of the current paper was to establish to what extent our approximation reproduces the full thermalization process, i.e. not only the behavior of fields in the vicinity of the boundary (responsible for the expectation values of gauge theory operators), but also the entire bulk spacetime (which carries direct information about e.g. nonlocal probes in the gauge theory).
In our previous studies we linearized Einstein's equations around the final-state black brane solution and we kept only first order terms in the amplitude expansion. We then solved these linear equations supplemented with appropriate AdS boundary conditions and with the same initial conditions as required to solve the full Einstein's equations. In the spirit of the closed-limit approximation the initial conditions were not necessarily small perturbations of the final black brane. This procedure correctly captured the evolution of the expectation value of the gauge theory stress tensor with a 20% accuracy, but it predicted no entropy production despite the fact that this was often substantial.
From this perspective, our previous studies in [1] provided only a partial contribution to the full leading order approximation to the bulk spacetime, which also includes terms quadratic in the amplitude expansion. These quadratic terms do not contribute to the expectation value of the stress tensor, but they do modify the bulk geometry. In order to probe this effect, in this paper we considered the contribution of the quadratic terms to the entropy production and concluded that this is also reproduced with a 20% accuracy.
We also analyzed the effect of higher order corrections in the amplitude expansion. Not surprisingly, we found that for initial states that start off sufficiently close to equilibrium the inclusion of these terms improves the agreement with the nonlinear evolution. However, for large deviations from equilibrium (large amplitudes) the inclusion of these terms actually spoils the 20% agreement between the leading order approximation and the exact result. This suggests that in future cases in which a full nonlinear simulation might not be feasible, it might be enough to work with the leading order approximation.
An interesting byproduct of our analysis is that the quasinormal modes, which so far were thought to be responsible only for the late time approach to equilibrium in the expectation values of dual operators, might actually predict quite accurately the full time dependence of the dual stress tensor provided a sufficient number of them is included.
The most important open problem raised by our previous work [1] and the current paper is how to generalize the close-limit approximation to more phenomenologically interesting setups, for example to expanding plasmas such as those considered in [6][7][8]. Those setups are qualitatively different form the one we have considered because of two interrelated features. First, the late time physics will be described by a solution of hydrodynamic equations which is not known a priori in terms of the initial conditions without solving for the evolution of the system. This contrasts with the current case, in which the final state is a homogeneous thermal state whose energy density is known from the start. Second, the late time state will not be static, so one does not expect to be able to describe it by linearizing around the AdS-Schwarzschild black brane. Instead, since the energy density of an expanding plasma gets more and more diluted as time progresses, one might naively hope to be able to linearize around the Poincare patch of vacuum AdS. An indication that this is not the right procedure comes from Ref. [21], whose authors analyzed (mainly) the gravitational collapse of a massless scalar field triggered by turning on its non-normalizable mode with a small amplitude. The naive small amplitude expansion on top of vacuum AdS had to be resummed so that the background became the AdS-Vaidya solution instead. The latter geometry has a macroscopic horizon which governs the dissipation of perturbations propagating on top of it. This suggests that, in cases of expanding plasmas, the correct procedure is to linearize around some other expanding configuration. We will report elsewhere on studies in this direction in the context of the simplest expanding plasma system, the boost-invariant flow [36].