On thermalization of a boost-invariant non Abelian plasma

Using a holographic method, we further investigate the relaxation towards the hydrodynamic regime of a boost-invariant non-Abelian plasma taken out-of-equilibrium. In the dual description, the system is driven out-of-equilibrium by boundary sourcing, a deformation of the boundary metric, as proposed by Chesler and Yaffe. The effects of several deformation profiles on the bulk geometry are investigated by the analysis of the corresponding solutions of the Einstein equations. The time of restoration of the hydrodynamic regime is investigated: setting the effective temperature of the system at the end of the boundary quenching to $T_{eff}(\tau^*)=500$ MeV, the hydrodynamic regime is reached after a lapse of time of ${\cal O}$(1 fm/c).


Introduction
The possibility of employing the gauge/gravity duality correspondence [1][2][3] to the analysis of far-from-equilibrium processes in strongly coupled systems 1 has enlarged the realm of nonperturbative phenomena to which this theoretical method can be applied in a quite controllable way, while traditional approaches are less effective. Among the different systems, of particular interest is the one produced in ultrarelativistic heavy ion collisions (HI), as those taking place at RHIC and at LHC. The features of this system, for time scales larger than about 1 fm/c, seem to be well reproduced in a hydrodynamic setup involving a strongly coupled/low viscosity fluid [5,6], a framework allowing to correctly describe experimental observables such as the light hadron transverse momentum spectra. In the hydrodynamic scheme of crucial importance are, after the pre-equilibrium regime following the collisions, the conditions at the time when the hydrodynamic behavior sets up. In particular, at this time the system stress-energy tensor T µ ν is important. We write it as in terms of the system energy density and of p ⊥ , p the pressures along one of the two transverse directions (with respect to the collision axis) and in the longitudinal direction, respectively 2 . For example, in approaches based on the idea of an initial state described by color glass condensate with a saturation scale Q s of the order of a few GeV, the initial stress-energy tensor comes from classical chromo-electric and chromo-magnetic fields and has the form T µ ν ∝ diag(− , , , − ). However, such an initial condition has been shown to be unstable, as soon as one-loop corrections in the strong coupling constant are switched on [7].
In the analysis presented in our study, a prominent role is therefore played by the stress-energy tensor, which is needed to understand the features of the equilibrium regime, the time needed to reach equilibrium, and the properties of the dynamics driving the system towards equilibrium. Questions to be answered are, for a system driven out-of-the-equilibrium, whether the late-time state is described by hydrodynamics, and whether the elapsed time needed to reach equilibrium is related to the dynamics during the out-of-equilibrium state. These issues concerning the transient process from a farfrom-equilibrium to a hydrodynamic regime can be faced in the holographic approach [8], using as a guideline the conjectured correspondence between a supersymmetric, conformal field theory defined in a four dimensional (4D) space-time and a gravity theory in a higher dimensional space-time, a five dimensional Anti de-Sitter space (AdS 5 ) times a five dimensional sphere S 5 [1]. The gauge theory is defined on the 4D boundary of AdS 5 , and the connection has been extended to the case of nonlinear fluid dynamics [9]. Non-equilibrium can be studied by solving in the bulk the 5D Einstein equations for the metric subject to suitable time-dependent boundary conditions. Information about the late-time regime can be obtained computing various invariants. For example, in Ref. [10] the square of the Riemann tensor 2 (the Kretschmann scalar) has been studied as an expansion in inverse powers of the proper time up to second order, finding that, in the asymptotic τ → ∞ regime, the scalar 2 is free of singularities in the perfect fluid case. The boundary theory stress-energy tensor results from the solutions of the Einstein equations in the bulk. 1 For a review see [4]. 2 Along the paper, we refer to energy density and pressures without considering the factor Analyses in the case of boost-invariant fluid [11,12], for some choices of the profile distortion, have shown the formation of a horizon in the bulk, and have given access to the components of the boundary theory stress-energy tensor. The investigations have been extended to samples of initial states characterized by different values of the components of the stress-energy tensor, without reference to the mechanism producing each state [13][14][15][16][17][18].
As formulated in [11,12], the initial off-equilibrium states can be thought as being produced by external sources distorting the boundary metric for short time intervals. We focus on this issue, with a study of external distortion profiles, step-shaped or shortduration sequences repeating themselves with different intensities, with the aim of describing phenomena where a small number of collisions takes place before the system starts evolving to thermal equilibrium. This study in the boost-invariant framework is useful for more involved cases where a lower symmetry is assumed for the system, for example in shock-wave collisions [19][20][21][22].
There are several reasons to study complex distortion shapes in the boundary metric. The bulk Einstein equations are nonlinear, and different distorsions produce unrelated responses. The resulting effective temperature and entropy density follow different profiles as a consequence of the boundary sourcing. Moreover, it is interesting to investigate whether the onset of the hydrodynamic regime is related to the distortion curve. Finally, these analyses have to face problems concerning the stability of the solution of GR equations with different boundary conditions.
The plan of the present paper, in which a few of the above issues are examined, is the following. In Section 2 we describe how a system is taken out-of-equilibrium through the distortion of the boundary metric, as done in [11,12]. We discuss a suitable procedure for solving the resulting 5D Einstein equations, providing a few details of the algorithm and the issue of the stability of the solutions. In Section 3 we describe the profiles of the boundary geometry distortion we have chosen to study, with the motivations. In Section 4 we discuss our findings, in particular concerning the equilibration time. The conclusions are presented in the last Section. Details dealing with the calculation of the boundary stress-energy tensor T µ ν are collected in appendix A, while in appendix B we discuss some aspects of the numerical algorithm.

Distorting the boundary
The main idea, as in [11,12], is to study the effects of a time-dependent deformation of the metric of the 4D boundary. As a consequence of the deformation, gravitational radiation is produced and propagates in the 5D bulk, and a black hole is formed together with its horizon. In this way, one can trade the study of the reach of equilibrium in the 4D boundary theory for the analysis of the time evolution of the geometry in the 5D bulk.
We set the 4D coordinates x µ = (x 0 , x 1 , x 2 , x 3 ), identifying x 3 = x with the axis in the direction of collisions and along which the plasma expands. Boost-invariance along that axis is imposed, together with translation invariance and O(2) rotation invariance in the orthogonal plane x ⊥ = {x 1 , x 2 }. In terms of the proper time τ and rapidity y, given by x 0 = τ cosh y, x = τ sinh y, the 4D Minkowski line element is given by ds 2 = −dτ 2 + dx 2 ⊥ + τ 2 dy 2 . A distortion of the boundary metric, leaving the spatial three-volume unchanged and respecting the imposed symmetries, can be obtained through a function γ(τ ) which encodes information about a deformation: The space-time with the metric (2) is considered as the boundary of the 5D space-time in which the gravity dual is defined. We adopt 5D Eddington-Finkelstein coordinates, with r the extra-dimension coordinate and the boundary reached for r → ∞. In general, the 5D bulk metric can be chosen as with the functions A, Σ and B only depending on r and τ to respect the adopted symmetries. The behavior of these functions against the deformation allows us to describe how the bulk metric changes as a consequence of the distortion of the boundary. With the metric (3), at fixed x ⊥ and rapidity y, infalling radial null geodesics correspond to constant τ , while outgoing radial null geodesics are obtained from dr dτ = A(r, τ ) 2 . Therefore, for a generic function ξ(r, τ ) the derivatives represent directional derivatives along the infalling radial null geodesics and the outgoing radial null geodesics, respectively. The 5D metric (3) is the solution of Einstein's equations with negative cosmological constant. In terms of A(r, τ ), Σ(r, τ ) and B(r, τ ) the equations can be rephrased as [12]: The five equations (6)-(10) can be considered as three dynamical and two constraint equations. The condition that the metric (3) produces the 4D metric (2) at the boundary, for r → ∞, constrains the large r behavior of A(r, τ ), Σ(r, τ ) and B(r, τ ): The invariance of (3) under the diffeomorphism r → r + λ(τ ) has been exploited to impose the large r condition for A. We switch on the distortion of the boundary metric at the initial time τ = τ i , starting from the AdS 5 bulk metric: As discussed in [14], taking the limit r → ∞ and τ → 0 with the Eddington-Finkelstein coordinates, Eq. (14), is ambiguous. For this reason, we set τ i > 0. The initial conditions for A, Σ and B are therefore: For r → ∞, the Einstein's equations can be solved starting from the relations (11)- (13). As suggested in [11,12], for A(r, τ ), Σ(r, τ ) and B(r, τ ) the large-r expansions can be written as The conditions (11)-(13) fix the n = 0 coefficients in (18)- (20), as well as those for n = 1 in A asy . The other coefficients can be determined imposing that Eqs. (6)-(8) are satisfied by the functions (18)- (20), but this leaves two coefficients undetermined, a 4 (τ ) and b 4 (τ ). From Eq. (10) a relation follows between a 4 and b 4 : with the primes denoting derivatives with respect to τ , and the function G(τ ) expressed in terms of γ (n) (τ ) = d n γ(τ ) dτ n : In our procedure, a 4 (τ ) and b 4 (τ ) are determined by Eq. (21) and matching the computed functions Σ and B with their asymptotic expressions (19) and (20) for large r. Using a 4 (τ ) and b 4 (τ ), important observables can be obtained, such as the components of the boundary stress-energy tensor T µ ν , as discussed in next sections. The details of our method for solving Eqs. (6)-(10) are as follows. In the set of functions Σ,Σ, B,Ḃ and A, we treat a given function and its derivative (5) as independent. At τ = τ i , Σ(r, τ i ), B(r, τ i ) and A(r, τ i ) are known, Eqs. (15)- (17). To complete the set of initial conditions, we computeΣ(r, τ i ) andḂ(r, τ i ) solving Eqs. (6) and (7), respectively.
At τ > τ i the algorithm is organized in steps.
1. Using the definition (5), we determine ∂ τ Σ(r, τ i ) and ∂ τ B(r, τ i ), which allow us to obtain Σ and B at the subsequent time to be valid.
2. Once Σ(r, τ i + dτ ) and B(r, τ i + dτ ) are determined, we compute a 4 (τ i + dτ ) and b 4 (τ i + dτ ). Indeed, a 4 (τ i + dτ ) can be obtained fitting the computed functions Σ(r, τ i + dτ ) and B(r, τ i + dτ ) with their asymptotic expressions (19) and (20) in the range r as r r max (later on in this section, it will be shown how to set r as ).
With the functions known at τ i + dτ , the cycle is repeated until a chosen final value of τ is reached. Typically, with the γ functions considered in our analysis we use dτ = 10 −2 .
At each τ step, our solution algorithm requires three values of the bulk coordinate r: r max and r min , which are the maximum and minimum value of r used in the numerical calculation, and r as , a value employed to determine a 4 (τ ) from the asymptotic expressions (18)- (20) used in the range [r as , r max ]. Stability against variation of such parameters is used as a criterium to check the numerical results. Particularly relevant is the minimum value r min , required to perform an excision in the bulk coordinates. This value is set in the following way. At fixed τ , for each one of Eqs. (6)-(10) we construct the function R i (r, τ ) defined as the absolute value of the l.h.s. of the equation divided by the sum of the absolute value of its addendi (see appendix B for definitions). At τ = τ i , we choose a small initial value of r min ; then, at each τ , r min is increased by steps of size typically 0.005, until the differential equations (6)-(8) and one constraint equation (10) are fulfilled, namely by requiring the condition R i (r min , τ ) < 0.03 for i = 6, 7, 8, 10. As a result, the bulk excision is always beyond the apparent horizon. In r min , the function A(r min , τ ) becomes negative but not large, hence possible singularities beyond the horizon are avoided. A similar procedure is carried out for r as . We initially set r as , and then we gradually increase it until the condition R i (r as , τ ) < 0.01 for i = 6, 7, 8, 10 is satisfied. Therefore, at this stage, we monitor the accuracy of the three differential equations and one constraint equation in r min and r as . At the end of the procedure, we check that the condition R i ≤ 0.01 is fulfilled for all Eqs. (6)-(10), i.e. i = 6 − 10, and in the whole (r, τ ) range. We find that this is indeed the case except for few tiny regions. In the largest part of the (r, τ ) plane deviations from zero are smaller than O(10 −4 ), as discussed in appendix B. During the numerical evaluation, no corrections are required to account for the finite value of the time step. In spite of the intricacy of the boundary conditions, stable solutions are found. 3 3 Models for the distortion of the boundary geometry We are interested in investigating deformations of the boundary metric characterized by different duration, intensity and structure. We set the generic form for γ(τ ): with and This expression generalizes the one used in [12], and represents sequences of N "pulses", each one of intensity c j and duration proportional to ∆ j , with the possibility of superimposing a smooth deformation, of intensity w, which asymptotically changes the scales of the coordinates. After their generation, pulses travel in the longitudinal direction x , driving the boundary state out-of-equilibrium. The response of the 5D bulk geometry to this deformation describes how equilibrium is reached at late times. We choose three different models, A, B and C, setting the various parameters in Eqs. (23)- (25). Model A represents two pulses in the boundary metric. The parameters are set to w = 0 and N = 2, c 1 = 1, ∆ 1,2 = 1, τ 0,1 = 5 4 ∆ 1 , c 2 = 2. Moreover, two different values of τ 0,2 : τ 0,2 = 17 4 ∆ 2 and τ 0,2 = 9 4 ∆ 2 are considered, changing the time interval between the pulses. The final time of the distortion is τ A f = 5. 25 and τ A f = 3.25 for the two cases, respectively, two quantities important for our discussion.
Model C combines both the previous ones, and is obtained using 3 Other methods for the numerical solution of analogous General Relativity equations in the case of initial value problems are described in [23] and in references therein.
For this model the last short pulse ends at τ C f = 9.45. As mentioned, to avoid in the Eddington-Finkelstein coordinates the ambiguity in the limits r → ∞ and τ → 0, in all cases the initial time is τ i > 0. The three profiles γ(τ ) are depicted in Figs. 2, 5 and 8.
The distortion profiles are chosen with the purpose of studying different situations. Model A is the repetition, in the boundary metric, of two short pulses of different intensity, Fig. 2. In particular, we consider two cases with different interval of time between the pulses, ranging from the condition of distant perturbations to the case of overlapping pulses. This model can provide us with information on the horizon formation, in comparison with the case of isolated pulse studied in [12], and on the effects related to the time elapsed between the distortions. Model B has the purpose of studying a combination of two effects with very different time scales, a short pulse and a slow continuum asymptotically producing a rescaling of the boundary coordinates. Model C, a combination of the previous cases, accounts for the effects of several pulses with various intensities, and is closer to physical processes driving systems out-of equilibrium in relativistic heavy ion collisions. In all cases, the horizon formation and behavior are investigated, together with the various components the boundary stress-energy tensor, looking at the time when the hydrodynamic regime sets in.

Results and discussions
In the numerical calculation for the models previously described, the solutions of the bulk Einstein equations allow to obtain the apparent horizon and the event horizon, which define the system temperature T ef f . The apparent horizon is determined as the locus of points where the conditionΣ(r h , τ ) = 0 is fulfilled. 4 The event horizon, which separates causally disconnected regions, is found drawing outgoing radial null geodesic curves and looking for the critical curve which separates the geodesics escaping towards the boundary from the ones plunging back in the bulk. From these curves, the function r h (τ ) (r h (τ ) being the position of the horizon at the proper time τ ) is determined, which defines the temperature T ef f (τ ). Asymptotically in proper time, the apparent horizon and the event horizon coincide. 5 From the solution, Σ(r h , τ ) 3 , the area of each horizon per unit of rapidity, can be computed.
The three models share common features, which we discuss together with the differences. Before presenting this analysis, let us recall a few results concerning the boundary stress-energy tensor T µ ν . As discussed in the Bjorken's seminal paper [27], under the assumptions of homogeneity, boost-invariance and invariance under rotations in the transverse plane with respect to the fluid velocity, the fluid stress-energy tensor T ν µ has well-defined properties, and its components only depend on the proper time τ , a manifestly boost-invariant condition. Imposing T ν µ conservation and traceless conditions, the tensor components can 4 The apparent horizon is the outermost trapped null surface. To compute it, a foliation of spacetime at fixed times is considered. In each spacelike hypersurface, the apparent horizon is a closed surface such that [24]: with γ ab the induced metric and K ab the extrinsic curvature of the hypersurface, K = g ab K ab , and s a the outward-pointing spacelike unit vector, normal to the apparent horizon and tangent to the hypersurface. For the spacetime described by Eq. (3) this corresponds to the conditionΣ = 0. 5 Studies of the behavior of the horizons in the gravity dual of a boost-invariant flow can be found in [25,26].
be expressed in terms of a single function f (τ ), that can be chosen to be −T 0 0 : For a perfect fluid, the equation of state = 3p, with p = p = p ⊥ , fixes the τ dependence: (τ ) = const τ 4/3 , stemming from the relations Deviations from the ideal behavior can be included as corrections taking viscous effects into account, in a gradient expansion. On the gauge theory side, the gradient expansion corresponds to a late time expansion, with subleading terms identified with the dissipative ones in the hydrodynamical theory. As shown in [10], the function f (τ ) deviates from f (τ ) = const/τ 4/3 if the assumption of perfect fluid behavior is removed. For example, the case (τ ) ∼ 1 τ s corresponds to p ∼ (s − 1) τ s and p ⊥ ∼ (2 − s) 2τ s . 6 Subleading corrections to the perfect fluid behavior relate the energy density to the definition of an effective temperature T ef f (τ ). The calculation in N = 4 SYM gives [14,28,29]: This leads to and to the pressure ratio and anisotropy: with c 1 = 1 3π and c 2 = 1 + 2 log 2 18π 2 . Λ is a parameter to be determined for each one of the considered models.
Understanding the results of the calculation of the solutions of the 5D Einstein equations, hence, requires the computation of the boundary stress-energy tensor T µ ν in Eq. (1) and the comparison with the proper time dependence following the previous equations. The various components of the stress-energy tensor of the boundary theory can be determined using the holographic renormalization procedure developed in [30,31] and shortly described in appendix A. The functions a 4 (τ ) and b 4 (τ ) in the expansions (18) and (20), computed through the Einstein equations, are required, and the energy density and the longitudinal and parallel pressures are obtained in terms of such functions: The functions˜ γ (τ ),p ⊥,γ (τ ) andp ,γ (τ ) are related to the distortion profile γ(τ ) in the boundary metric, and are reported in appendix A.
We can now discuss the three considered models, focusing on the effects of the boundary distortion and on the transient from the far-from-equilibrium state to the hydrodynamic regime.
Model A.
In the case of two distant pulses, the horizon is formed as a two-step process which follows the sequence of the boundary distortion, as one can infer considering the computed function A(r, τ )/r 2 which is depicted in Fig. 1. The outgoing radial null geodesics are clearly separated by a critical geodesic, the event horizon, which starts plunging back in the bulk immediately after each pulse in the boundary. This can be even more clearly seen in the plot of the τ dependence of the effective temperature T ef f (τ ) and in the plot of the area of the apparent horizon per unit of rapidity, Σ(r h , τ ) 3 , collected in Fig. 2. The temperature has two decreasing regimes, the first one after the first pulse, with the system starting relaxation immediately after the end of the boundary distortion. Time intervals of about twice the duration of the pulse are very long with respect to the relaxation time. The same phenomenon is observed considering the function Σ(r h , τ ) 3 , which shows two time ranges in which it reaches a constant value, soon after the pulses. However, on such time scales, the components of the stress energy tensor have their own (computed) behavior, and no isotropization is observed after the first pulse.
On the other hand, for a distortion represented by two close (nearly overlapping) pulses, no particular structure is observed before the end of the boundary perturbation. The horizon area per unit of rapidity, Σ(r h , τ ) 3 , has a sharp increase during the distortion, and reaches the asymptotic large τ value Σ(r h , τ ) 3 = π 3 Λ 2 soon after the end of the distortion.
The hydrodynamic τ −dependence of T ef f (τ ) and of the energy density (τ ), Eqs. (30) and (31), is reached immediately after the end of the distortion of the boundary τ A f , with values of the parameter Λ which are different for the two cases of different time intervals between the pulses. The values of Λ obtained from T ef f (τ ) are collected in Table 1. The gray lines are radial null outgoing geodesics, the dashed dark blue line is the apparent horizon, and the continuous cyan line is the event horizon obtained as the critical geodesic. The excision in the low-r region, used in the calculation, is shown. In the right panels, the 3d representation of A(r, τ )/r 2 is displayed together with the apparent horizon (continuous line).
The temperature can be obtained from the conditionΣ(r h , τ ) = 0, or from the energy density Eq.(29), with numerical results differing at the level of 2%. For a quantitative discussion, a time τ * can be fixed imposing that the energy density (τ * ) differs from the hydrodynamic value H (τ * ) in (31) by less than 5%: In both cases of distant or overlapping pulses, τ * essentially coincides with the end of the boundary metric distortion. From the plot of T µ ν in Fig. 3, we observe that at τ * the pressure anisotropy is still sizeable, and the values of ∆p/ and p || /p ⊥ are still much larger than the ones expected by viscous hydrodynamics. Therefore, we define the "thermalization time" τ p from the condition where (p || (τ p )/p ⊥ (τ p )) H is given by (34). Quantitatively, from the condition that p /p ⊥ differs by less than 5% from the asymptotic NNLO expression, we can set τ p = 6.8 for two distant pulses. The difference τ p − τ * can be expressed in physical units if one scale in the system is set. Imposing T ef f (τ * ) = 500 MeV corresponds to τ p − τ * 0.60 fm/c, which indicates the elapsed time between the end of the pulse and the restoration of the hydrodynamical regime. For two overlapping pulses, the time difference is τ p − τ * 1.03 fm/c. The numerical values are collected in Table 1. Notice that τ * and τ p differ in general, although Eqs. (27) and (28) hold. Indeed, at a generic τ the relation between energy density and its asymptotic limit can be written as: with ξ(τ * ) such that (39) is verified. For, e.g., the longitudinal pressure, using (27) and (41) one has p (τ * ) = ξ(τ * ) p , H (τ * ) − τ * ξ (τ * ) H (τ * ). Therefore, due to the last term, the value τ p verifying (40) does not necessarily coincide with τ * .
Model B.
The purpose of this model is to investigate a boundary distortion characterized by two largely different time scales. The horizon is formed immediately after the boundary deformation is switched on, as depicted in Fig. 4, then starts plugging down in the bulk until the short pulse becomes active. This can be seen clearly in Fig. 5, where the effective temperature and Σ 3 display a decreasing and constant regime, respectively, in the proper time interval τ i < τ < τ 0,1 , after which a violent perturbation takes the system out of equilibrium. The time τ * can be fixed using (31). For τ > τ * τ B f the determination of Λ from T ef f and the energy density gives Λ = 1.12, with a variation of 3% around this value. Thermalization is reached at the time τ p obtained by the condition (40), as one can infer considering the observables depicted in Figs. 5 and 6. Setting the scale so that the temperature is T ef f (τ * ) = 500 MeV, we find τ p − τ * 0.42 fm/c.

Model C.
This model is constructed with the aim of representing a situation close to the physical system during the scattering process, with two different time scales and several pulses of various number and intensity. Immediately after switching on the boundary deformation, the horizon is formed in the bulk geometry and follows the profile of the distortion. Even in short time intervals between the pulses, the horizon starts to plung in the bulk, with a decreasing temperature and a saturation of Σ 3 . This can be inferred considering the computed function A/r 2 depicted in Fig. 7 and, in details, studying the τ dependence of the effective temperature and horizon area per unit rapidity in Fig. 8. Σ 3 (r h , τ ) displays a step-like behavior in proper time, closely following the distortion γ(τ ), and reaches a constant value from τ = τ C f on. On the contrary, the various components of the stressenergy tensor have structures which become regular only for τ > τ C f , after the last pulse. As one can observe analyzing the results in Fig. 9, the energy density reaches the NNLO hydrodynamic behavior at τ = τ * τ C f , while pressure anisotropy persists for longer times. The parameter Λ, obtained from different observables, takes the value Λ = 1.59, with a variation similar to the one in models A and B. The pressure anisotropy ∆p/ and ratio p /p ⊥ set the value τ p of proper time, as one can infer considering the results in Fig. 9; setting T ef f (τ * ) = 500 MeV, we obtain τ p − τ * 0.2 fm/c. of the stress-energy tensor, pressure anisotropy ∆p/ = (p ⊥ − p )/ and ratio p /p ⊥ , computed for τ > τ A f , for the boundary distortion with two distant (left) and two overlapping pulses (right panels). The short and long dashed lines correspond to the hydrodynamic result and to the NNLO result in the 1/τ expansion.

Conclusions
We have investigated the effects of different types of distortions of the boundary metric, in a boost-invariant setup. The quenches are introduced as a way to take the boundary theory out-of-equilibrium, and the distortion profiles are used to describe processes with different time scales and intensities. As a common feature, we observe a rapid formation of the horizon in the bulk metric, with the possibility of defining an effective τ −dependent temperature. In coincidence with the end of each main distortion, the relaxation starts with the horizon plunging in the bulk. We have seen how sequences of distortions break the relaxation process. We find that for all the considered distortion profiles, T ef f (τ ) starts to follow a viscous hydrodynamic expression as soon as the quench is switched off. At this time the pressures are different, and evolve towards a common value which  Table 1: Numerical values of the relevant quantities for models A, B and C of boundary distortion.
is reached at a later time. Setting T ef f (τ * ) = 500 MeV, the elapsed time before the restoration of the hydrodynamic regime is always a fraction of fm/c. Although the system described by the holographic approach is in several respects different from the real QCD system, and boundary sourcing a quite abstract representation of the heavy ion collision process, the obtained results allow us to argue what can be expected in realistic situations. in terms of the induced metric h µν , of the lapse function N and of the shift vector field N µ . The stress tensor for the space foliation M r obtained at constant r, is given bỹ with h ≡ det (h µν ) and S the gravitational action. The boundary stress energy tensor is obtained from the limit: The action S in (43) includes as a counterterm a local functional S ct of h µν , whose contribution toT µν regularizes the divergences at r → ∞ [33], [31]. As a consequence, the stress-energy tensor is written as with K µν the extrinsic curvature, (4) G µν the boundary Einstein tensor with curvature tensors defined with respect to h µν . The counterterms −3h µν + 1 2 (4) G µν cancel the powers from r 2 to r −1 from the first two terms in (45), and introduce terms of order r −2 which contribute to the final result (45) [33]. The last term proportional to σ µν , with in our case σ µν = diag − 3 2 α 4 (τ ), e γ(τ ) − 1 2 α 4 (τ ) + 2β 4 (τ ) , e γ(τ ) − 1 2 α 4 (τ ) + 2β 4 (τ ) , e −2γ(τ ) τ 2 − 1 2 α 4 (τ ) − 4β 4 (τ ) and α 4 (τ ) and β 4 (τ ) in (15)- (17), cancels the 1 r 2 log r contributions in (45). The results

B Testing the numerical algorithm
To quantify the accuracy of our numerical algorithm, based on the Runge-Kutta method for the solution of the differential equations in the variable r, we monitor the ratios of the stress-energy tensor, pressure anisotropy ∆p/ = (p ⊥ − p )/ and ratio p /p ⊥ , computed for τ > τ B f . The curves have the same meaning as in Fig. 6.
in the (r, τ ) domain behind the excision r > r min . The final results are shown in Fig. 10 in the case of model B. The resulting ratios R 6 , R 7 , R 8 deviate from zero at a level smaller than O(10 −4 ), but for a tiny region close to the excision in the case of R 6 and  R 8 , and for a few spots in the range 3 r 5 (where the source function γ(τ ) is peaked) in the case of R 7 . A similar result is obtained for R 9 and R 10 , with larger deviations from zero: this is a consequence, for R 10 , of the smallness of the two addendi in (10), and for R 9 of the computation ofΣ = ∂ τΣ + 1 2 A (Σ) through a discretized time derivative. To monitor the time slicing, we compute R 9 for two different time steps in the region 1.2 r 3.0 at τ = 4, where there is the largest deviation from zero. Comparing the results for dτ = 10 −2 and dτ = 0.5 · 10 −2 , as shown in Fig. 11, one observes the decreasing of R 9 by doubling the grid points.