Holographic approach to thermalization in general anisotropic theories

We employ the holographic approach to study the thermalization in the quenched strongly-coupled field theories with very general anisotropic scalings involving Lifshitz and hyperscaling violating fixed points. The holographic dual is the Vaidya-like time-dependent geometries where the asymptotic metric has general anisotropic scaling isometry. We find the Ryu-Takanayagi extremal surface with which to calculate the time-dependent entanglement entropy between a strip region with width $2R$ and its outside region. In the special case with the isotropic metric, we explore the entanglement entropy also in the spherical region of radius $R$. The growth of the entanglement entropy characterizes the thermalization rate after the quench. We study the thermalization process at the early times and late times in both large $R$ and small $R$ limits. The allowed scaling parameter regions are constrained by the null energy conditions as well as the requirement of the existence of the extremal surface solutions. This generalizes the previous works on this subject. All obtained results can be compared with experiments and other methods of probing thermalization.


I. INTRODUCTION
The holographic duality provides a unique method to investigate the dynamics of strongly coupled field theories, which links the geometrical quantities to quantum observables. The use of this interplay to study thermal phases of strongly coupled field theories by its holographically dual AdS black hole/brane geometries has attracted lots of attention since its first proposal [1]. It was then suggested that the small perturbations of the AdS metric are dual to the hydrodynamics or linear response regions of boundary fields [2][3][4]. Also, by probed strings in AdS background, the diffusion behavior of Brownian particles in boundary fields can be derived(see [5]for reviews). In particular, the time evolution of entanglement entropy between Brownian particles and boundary fields had been studied by means of the probed string method [6].
Along this line of thoughts, it is very interesting to consider holographic analysis of strongly coupled field theories far from equilibrium. For example, a system, which starts from a highly excited state after a quench, is expected to evolve towards a stationary state at thermal equilibrium. The holographic dual of this far from equilibrium problem is a process of gravitational collapse ending in the formation of a black hole/brane, described by time-dependent Vaidya geometries. To understand the dynamics of thermalization we use the time-dependent entanglement entropy in the holographic setup. According to Ryu and Takayanagi [7], the entanglement entropy between a spatial region, Σ of dimension d − 1 and its outside region in the d-dimensional boundary theory is dual to the area of the extremal surface Γ in the d + 1-dimensional bulk geometry, which is holomorphic to Σ and has the same boundary ∂Σ as Σ. The prescription for the time-dependent holographic entanglement was proposed in [9]. The use of this prescription to study the thermalization following the quenches in various Vaidya-like backgrounds can be found in . In these studies, the region Σ in the boundary theories is bounded by ∂Σ, which is either a d − 2-dimensional sphere of radius R or two planes separated by distance 2R. Computing the time dependent entanglement entropy as a function of spatial scale thus provides a probe of scale-dependent thermalization. The bulk Vaidya geometry describes the falling of the d−1-dimensional thin shell along the light cone from the boundary at the time t = 0. If we consider the boundary to be the plane, the black brane eventually forms when the thin shell falls within the horizon distance, y h . In the boundary theory, this process is dual to the input of energy at t = 0 (quench) driving the system far from equilibrium, and the subsequent thermalization process of the system.
In this paper we extend some of these works of [25] and [27], and consider the Vaidyalike geometries that describes the formation the black brane with the general asymptotic anisotropic Lifshitz metric including hyperscaling violation. The holographic models in this background have been used to study the diffusion of heavy quarks in anisotropy plasma when it is slightly out-of-equilibrium [31], and also the properties of dissipation and fluctuations of Brownian particles [32]. Here, in the cases of far-from-equilibrium, we consider the time evolution of entanglement entropy between a strip region of width 2R and its outside region.
To have the analytical expressions we focus on both large R and small R limits as compared to the horizon scale. In both cases, the early time and late time entanglement growth and its dependence on the scaling parameters in this model will be explored. Notice that in this study, the anisotropy effect is encoded in the effective bulk spacial dimension,d as will be defined later, which is just d in the isotropic theories. In addition, the entangled region we propose is to probe the anisotropic spatial coordinate with the scaling parameter, say a 1 relative to the scaling parameter a y in the bulk direction where there exists a free parameter ∆ = a 1 + a y − 4, which is zero in [25] and [27]. Apart from the strip case, we study the entanglement region bounded by a sphere of radius R when the backgrounds are isotropic with ∆ = 0. We thus mainly focus on the contributions from the nonzero ∆ to the dynamics. Additionally, the constraints on these parameters from the null energy condition together with the constraints from the existences of the extremal surface are derived so that the allowed scaling behaviors of holographic entanglement entropy can be used to justify (or falsify) the holographic method from experiment tests and other methods on strongly coupled problems.
Layout of the paper is as follows. In Sec. II we introduce Lifshitz-like anisotropic hyperscaling violation theories, and corresponding static black hole metric. The null energy condition is developed to impose the constraints on the scaling parameters of theories. The dynamics of the extremal surface and its corresponding equations of motion are studied in Sec. III. In Sec. IV, we study the entanglement entropy at thermal equilibrium. The early time entanglement entropy growth and its late time saturation will be discussed in Sec. V and VI respectively. In Sec. VII, we choose the Einstein-Axion-Dilaton theory as an example to realize the allowed scaling parameter regions given by the constraint from the null energy condition as well as the requirement of the existence of the extremal surface solutions. Finally, Sec. VIII concludes the work. In Appendix A, we provide the detailed analysis about the regimes of the scaling parameters for the system undergoing continuous and discontinuous saturation in the strip case.

II. THE HOLOGRAPHY BACKGROUND
To make our analysis of thermalization as general as possible, we consider the gravitational collapse that eventually forms the following black brane with the metric where near the boundary y = 0, the emblackening factor is assumed to be h(0) = 1. Then, in the boundary, the metric has scaling symmetry, We also assume that there is a simple zero for h(y) at y = y h , corresponding to the position of the horizon, in which near the boundary, h(y) has the leading term, The temperature of the black brane is given by T = |h (y h )| 4πy a 0 2 + ay 2 −2 h . In order to describe the formation of the black brane of (1), we introduce Eddington-Finkelstein coordinates, Then the metric in (1) becomes In this work, we consider the gravitational collapse described by the following Vaidya-type metric where f (v, y) = 1 − Θ(v)g(y), and g(y) = 1 − h(y). The metric describes the formation of the black brane of (1) by the infalling of a d − 1-dimensional delta-function shell along the trajectory v = 0. The region v > 0 outside the shell is with metric as in (1)  gives a quench on the system at t = 0, and subsequently the system evolves into the thermal equilibrium states with finite temperature.
Nevertheless, the null energy condition (NEC) constrains the parameters introduced in the metric [33], which is In Einstein gravity, this condition is equivalent to R µν µ ν ≥ 0 where R µν is the curvature tensor obtained from (7), and then gives the following constraint equations, a i a 1 + 1, and ∆ 0 = a 0 + a y − 4. In the case of the isotropic background with all a i 's equal, the NEC reduces to the one in the hyperscaling violating Lifshitz theory [27].

III. DYNAMICS OF THE EXTREMAL SURFACE
In this section, we derive the equations of motion for the extremal surfaces in the cases that entanglement region Σ in the boundary theory is a strip of width 2R in x 1 direction and is a region bounded by a sphere of radius R. In this work, we consider the d + 1-dimensional bulk geometry. Thus the boundary theory living at y = 0 is in the d-dimensional spacetime.
The entanglement region Σ in the boundary theory is bounded by d − 2-dimensional surface ∂Σ. We consider that the in-falling planar shell propagates from the boundary at the boundary time t = 0. Starting from this section, the time t means the boundary time and is different from the coordinate time t in the previous section. Then the d − 1-dimensional extremal surface Γ in the bulk, if exists, is uniquely fixed by the boundary conditions that it should match ∂Σ at the boundary y = 0 and at the boundary time t. Once we find the time-dependent extremal surface Γ and its area A Γ , the entanglement entropy between Σ and its outside region is given by the Ryu-Takayanagi formula The approach we adopt in this paper mainly follow the works of [25] and [27]. Here we highlight the key equations and showing relevant solutions in the following sections, which are straightforward generalization of their results to our model.

A. sphere
We first consider the case that ∂Σ is a sphere with radius R. Thus, to accommodate the d − 2-dimensional rotation symmetry, we assume the metric (7) with such symmetry by setting all a i 's to be the same, say a i = a 1 , and in the spherical coordinates, it can be written as where ∆ = a 1 + a y − 4. The embedding of d − 1-dimensional extremal surface Γ in (13) is described by two functions v(ρ) and y(ρ) together with an extension in all Ω d−2 directions.
The area of Γ is then given by where 2 ) R d−2 is the area of d − 2-sphere with radius R. The functions v(ρ) and y(ρ) are determined by minimizing the area in (14), giving the following equations of motion, √ Qy The boundary conditions arė Again, here and later, the time t will label the boundary time. Since f (y, v) = 1 for v < 0 = 0 in the both v < 0 and v > 0 regions. From (16) there exists a constant of motion E given by Solving (19) forv givesv . Usingv in (20) can further rewrite (17) as Here we in particular define ∆ = a 1 + a y − 4 to highlight the nonzero ∆ effects as compared to [27] with ∆ = 0 although the sphere case has the same spatially rotational symmetry as in [27]. Solutions of y(ρ) in v > 0 and v < 0 are matched at ρ c with v(ρ c ) = 0. Integrating both (16) and (17) across ρ c , the matching conditions are thatv(ρ) is continuous across ρ c andẏ where y c = y(ρ c ) and subscripts + and − denote the solution in the region of v > 0 and v < 0 separately. Once the solution y(ρ) is found, further integratingv in (20) with the initial conditions in (18) for v > 0 gives the boundary time In the end, the integral formula for the area (25) can be calculated from the contributions of v < 0 and v > 0 respectively as where we have used the fact that for v > 0,Q = h+ẏ 2 y ∆ h+B 2 E 2 by plugging (20) in (15), and for v < 0, Q = 1 +ẏ 2 y ∆ with h = 1 and E = 0. Translating ρ c to the boundary time t using the relation in (23), we are able to find the time-dependent area, and also the corresponding entanglement entropy.

B. strip
We now consider the entanglement region Σ of a strip extending along say In what follows, we denote x 1 as x for simplicity. In this case, the area of the extremal surface Γ is expressed as with , which is the effective bulk spacial dimension as also defined in [31].
Varying A Γ with respect to the functions v(x) and y(x) leads to the following equations of motion, The translational symmetry in (25) of A Γ in the variable x gives the conserved quantity The value of J can be obtained by evaluating it at x = 0 in the tip of the extremal surface and v > 0 regions, there exists another conserved quantity E given by (27). Together with (29), we have In the vacuum region with v < 0, the value of the x-independent E can be determined at a particular x = 0 where E = 0 given by the initial conditions. Also, f = 1 in the vacuum regime leads to the relation betweenv andẏ at arbitrary Substituting all above relations into (29), for v < 0 and x > 0 with no loss of generality, However, requiring dx/dy is finite as y → 0 gives the constraints Integrating (27) and (28) across the null shell allows us to find the matching conditions, which are the same forms of (22) by replacing ρ with x. Thus, the matching conditions in the strip case can determine the conserved quantity E of (30) for v > 0 in the black brane region given by Then, via (30), the relation betweenv andẏ in the black brane region becomeṡ The counterpart of the equation for y(x) in the black brane region with v > 0 as in (32) can be found from (35) with the relation in (29) aṡ With (35) and the square root of (36), we obtain for In the end, from (32) and (36), we can write down the strip R in terms of y c and y t as Due to (37) with the boundary condition v(R) = t, the boundary time t thus can also be expressed as a function of y c , Moreover, substituting (32), (35), (31) and (36) into (26), we summarize the integral formula for the area in (25) from the contributions of v < 0 and v > 0 respectively as Then, through (38) and (39), the area are expressed in terms of the boundary time t and the strip width R. Notice that with the additional constant J in the strip case, the area in (25) and the corresponding the entanglement entropy can be computed even without knowing the solution of y(x).

IV. ENTANGLEMENT ENTROPY IN THERMAL EQUILIBRIUM STATES
In this section we consider the entanglement entropy for a final equilibrium state dual in the bulk to a black brane in (6). We take the large R (R y ) limits respectively, where the former limit corresponds to the case that the tip of the black brane, y b to be defined later, approaches the black brane horizon, y h and the latter one is to assume y b y h so that the relevant metric we consider is that of the pure hyperscaling violating anisotropic Lifshitz spacetime of the form (1) in the case of h → 1.
The corresponding A Γ in both spherical and strip entanglement regions will be computed accordingly.

A. sphere
In the case that ∂Σ is a sphere, the extremal surface in the black brane background (6) by setting all a i 's equal to a 1 is denoted as Σ BH . In this case, the boundary conditionṡ v(0) =ẏ(0) = 0 give E = 0 in the equation of motion (21). We also denote the solution of (21) with E = 0, f = h(y) and satisfying the boundary condition (18) by y BH (ρ). And the area of the extremal surface can then be obtained from (24) by setting ρ c = R and E = 0 given by We first consider the extremal surface in the large R limit (R y ∆ 2 +1 h ). In this case, it is anticipated that the tip of the Σ BH , denoted as y BH (0) ≡ y b , is very close to the horizon y h , namely where 1. The solution of y BH near the horizon is in this form with the first order perturbation y 1 (ρ) given by where I µ (x) is the modified Bessel function of the first kind. Let us denote the inverse function of y BH (ρ) as ρ BH (y). We also have the expansion of ρ BH (y) near the boundary where To find the relation between and R, the existence of the matching regimes to do the matching of above two solutions is crucial [34]. We extend the solution (44) ( near horizon) to the region of relatively large ρ but still satisfying ρ R − O(R 0 ), and then the solution (46) (near boundary) to the region with 1−y/y h 1 but still satisfying R−ρ O(R 1 ).
We find that with such extensions both solutions (43) and (45) have the similar structure in the matching region where the relation of and R can be read off by comparing their leading order behaviour given by . Apparently, the large R limit drives the value of to a small value. The nonzero ∆ does not change the relation (47) found in AdS spacetime in [25].
The area of the extremal surface in the large R limit can be computed from the solutions (43) and (45) by splitting the integral (41) into the areas near horizon (IR) and near boundary (UV). Then the UV divergence part of A Γeq is where y U V is the UV cutoff in the y integration, and the finite part is obtained as Similar results are also obtained in [27] in the limits of ∆ = 0 of our system.
In the small R limit, namely (R y ∆ 2 +1 h ), the area of the extremal surface is determined by (41) in the limit of h → 1. In particular, in the case of a y = 2, we can find the exact solution of (21), where the relation of the tip of the extremal surface y and R, an essential information to analytically find the finite part of the area, can be read off as Substituting (50) into (41) (h → 1), and again dividing the area into divergent and finite parts, we have where (54) Here we restrict ourselves to the case d ≥ 3 so that the dimension of ∂Σ is larger than 2 for considering the different geometries of the extremal surfaces such as the sphere and strip cases. The above results can reproduce the extremal surface in AdS spacetime by choosing the appropriate values of the scaling parameters in [35]. Although a y = 2 is chosen, our results with ∆ = 0 and the obtained finite part of the extremal area generalize the result of [27]. Through the Ryu-Takayanagi formula (12), the corresponding entanglement entropy can also be obtained where its finite part is important to make a comparison with the field theory results.

B. strip
To find the area of the extremal surface of the strip, we substitute the relation between R and y b found in (38) to (40) by setting y c = y t = y b and E = 0, we have where H(y) is defined in (36). In the large R limit, the tip of the extremal surface y b is assumed to be close to y h in terms of the expansion of (43) for small where again can be related to R by (47). The straightforward calculations show that the UV divergent part of the area is and the finite part is Note that A Γdiv vanishes when a 1 (d − 1) < 2 + ∆ as y U V → 0. The contribution from the scaling parameter a 1 to the extremal surface plays the same role as in the sphere case, leading to the same entanglement entropy for the sphere and strip cases in our setting. Nevertheless, it will be seen that the departure from the thermal equilibrium and the subsequent time evolution for two cases are very different, in particular during the late-time thermalization processes.
In the small R limit, from (38) with y t → y (0) b and y c → 0, we obtain the relation between R and y 1) ) a dimensionless constant, and for a sensible result. The area of (55) for a 1 (d − 1) = 2 + ∆ becomes with A Γdiv in (56). However, for a 1 (d − 1) = 2 + ∆, The critical value determined by a 1 (d − 1) = 2 + ∆ with the logarithmic divergence rather than the power-law ones generalizes the result in [27] in that the isotropic background is considered. In particular, when a y = 2 and all a i 's equal to a 1 , the area of the extremal surface still has quite different behavior from that of the sphere case in the small R limit.
Starting from the next section, we will focus on the nonequilibrium aspect of thermalization processes encoded in the time-dependent entanglement entropy with respect to the spatial scale R with the gravity dual of the probe of the geometry of a collapsing black brane. Here the black brane horizon y h presumably sets the time scales for the nonequilibrium system to reach the "local equilibrium" as to cease the production of thermodynamical entropy. In the large R limit with R y , it is anticipated that the tip of the extremal surface y t in the end is still far away from the horizon y h , namely y t y h . It means that after the initial growth, the entropy of the system will directly reach its saturation value as the time t reaches the saturation time t s determined by the size of the system R, which will be studied later. For the quantitative comparison of the entanglement entropy between the small and large R limits, we just consider the timedependent entanglement entropy during the early growth and the final saturation dynamics to be studied in the following sections.

V. THE EARLY TIME ENTANGLEMENT ENTROPY GROWTH
In this section, we study the entanglement growth at the early times when the infalling sheet meets the extremal surface at y c when y ∆ 2 +1 c R for both large and small R limits.
It means that the in-falling shell is very close to the boundary y = 0 so that the extremal surface Γ is mostly in the v < 0 region with the pure hyperscaling violating anisotropic Lifshitz metric. We expect that the same growth rate will be found for large and small R limits in either the sphere or strip case. Also, during such early times, the infalling sheet does not have much enough time to probe the whole geometry in a way that the same growth rates will be expectedly achieved for both the sphere and strip cases.
A. sphere Here we start from considering the zeroth order solution of Γ that satisfies the equation of motion for y in (21)  For such a small y c , the relevant zeroth order solution during the early times is near the boundary for small y obtained from (45) and (46) as the blackening factor can safely be approximated by h(y) → 1 for small y. The existence of the extremal surface holomorphic to Σ requires that C > 0 and the finite dρ/dy as y → 0 giving the further constraints on the scaling parameters Then, we rewrite the area integral (14) in terms of the variable ρ(y), which is the inverse function of y(ρ), as where the prime means the derivative with respect to y and y t is the tip of the extremal surface, that is, y t = y(0). For the early times with small t, let us consider the perturbations around the zeroth order solution where the time-dependent A Γ (t) just slightly departs from The partial derivatives are evaluated at y = y (0) or ρ = ρ (0) and f = 1 where y t = y (0) (0) and v = v(y (0) ) = v (0) . The first term vanishes due to the vanishing of the area at the tip y = y (0) . The second and the third terms vanish due to the equations of motion for ρ and v. For the last term, we have δf = h(y) for y < y c , δf = 0 for y > y c , where y c is the position where the infalling sheet crosses the extremal surface, namely . Note that is required so that the boundary time t increases in y c . Thus, in terms of the boundary time, the early time entanglement entropy growth for both small and large R is given by where to ensure that the entanglement entropy increases in time. The extra ∆ dependence of the growth rate of the entanglement entropy generalizes the results in [27] with ∆ = 0. In particular, for the positive (negative) value of ∆, the power of t increases (decreases) with ∆ (|∆|) that speeds up (slows down) the growth rates as compared with the case of ∆ = 0.

B. strip
In the strip case, the area of the extremal surface can be obtained from (40) by setting where again the prime means the derivative with respect to y. As in the case of the sphere, the non-vanishing term of variation of A Γ (t) evaluated at the zeroth order solutions y (0) and v (0) = v(y (0) ), which are obtained from (32) and (5), is With δf in (67), straightforward calculations find the same entanglement growth at early times as in the sphere case in (69) by replacing d withd. As compared with [27] with ∆ = 0, the nontrivial dependence of ∆ in the early time growth rates as stated in the sphere case deserves the tests by the quantum field theories.

VI. THE LATE TIME SATURATION
Although the entanglement entropies in the early times exhibit the same growth rates between the geometries of the sphere and strip, the late time thermalization process will find the different time dependent behavior for two geometries, as the entanglement growth can be realized as an " entanglement tsunami" led by a sharp front moving inward from the boundary Σ [25]. Also, for a given geometry, we will work on the large and small R limits separatively. It is anticipated that the time dependent saturation behavior will be the same but the saturation time scales might be very different.

A. sphere
Let us start with the sphere case in the large R limit. Near the saturation, the extremal surface Γ is mostly in the v > 0 black hole region, and will become very close to the extremal surface in the purely black hole background (6). Thus, the in-falling shell is very near to the tip of Γ near ρ = 0, which can be parametrized as with δ 1. Also, near ρ = 0 for v < 0, the extremal surface y(ρ) satisfies (21) with E = 0 and f = 1 where the leading order solution y(ρ) can be approximated by From (73) and the definition y c = y(ρ c ), the relation between the tip of the extremal surface y t and the infalling sheet y c is obtained as Also, from the matching condition (22) at ρ * ≈ ρ c and the approximation solution (74), the variables y(ρ) and v(ρ) at the matching point are given respectively bẏ Plugging these into (19), the constant of motion E for v > 0 in the black hole region can be found to be Next, we can also do the expansion of y(ρ) around y BH (ρ) in the v > 0 region, in particular in the small ρ approximation as In the sphere case, d > 3 is considered. Given y(ρ → 0) = y t and y BH (ρ → 0) = y b in (78), the relation between y t and y b can be given from (75) as − 1 where the boundary time t depends on the infalling sheet y c through (23). Thus, near saturation, the boundary time t almost reaches the saturation time t s to be defined later.
We then can calculate the time dependent entanglement entropy near the saturation after the quench. The area of the extremal surface can be divided into v < 0 and v > 0 parts, Here we are interested in the time-dependent part of the δ-dependence of the extremal surface. Near saturation, the solution around ρ = 0 in (74) for v < 0 and (78) for v > 0 will be relevant. We denote ∆A Γ = A Γv<0 + A Γv>0 − A Γeq , by subtracting the extremal surface in the case of thermal equilibrium due to the pure black hole in (41), and it becomes From (23), the boundary time can be written as t = t 1 + t 2 where We then define the saturation time t s below and evaluate it in the large R limit from (20)with E = 0, f = h and the conditions v BH (R) = t and v BH (y b ) = 0 at t = t s , which is given by In that, T is the black brane temperature defined above, and the integral is dominated in the near horizon region. In terms of T , c E ∝ T together with (80), arriving at Note that d ≥ 4 is required. Based upon the scaling parameters constraints in (4) from the property of the black brane, (9),(10), and (11) due to null energy conditions, and (63) and (64) from the solution of the extremal surface as well as (68) and (70) for the sensible entanglement entropy, a 1 is always positive so that the boundary time t approaches to t s from its below. Thus, the saturation can continuously reach. It will be compared with the strip case where the continuous saturation will occur only for some parameter regions to be discussed later. The nonzero ∆ in the general model here does not affect the saturation behaviour of the entanglement entropy in the sense that the power law saturation depends on the dimension d only. As expected, the nonzero value of ∆ makes contribution to determine the saturation time scale t s .
Near the saturation in the small R limit, the analysis is different from the large R case as follows. In the large R limit, the tip the extremal surface in the pure black hole geometry y b is exponentially close to the horizon y h in a way that h(y b ) ∝ e −2γR . However, in the small R limit, y b y h giving h(y b ) 1. Moreover, we will have the same relation between y b and y c as in the large R limit (79) with the different d 1 to account for the fact that h(y b ) 1.
Also, the small δ expansion of ∆A Γ has the same powers of δ in the leading order of δ as in (80). In the end, the same saturation behavior in time for the extremal surface as in (84) is found in the small R. Moreover for small R, since a 1 is always positive, the saturation time t s becomes larger (smaller) as a 1 increases (decreases).

B. strip
In the case of the strip, let us consider that the relation between x c and y c near saturation is the same as in (73) by replacing ρ c by x c . From (32), we are able to find the solution of IR behaviour for v < 0 Due to the relation between x c and y c , and the definition of y c = y(x c ), from (86) we have Also, from (87) and (32), the conserved quantity E (34) can be approximated in terms of δ Recall that in large (small) R limit, h(y b ) ∆ h e −2γR (h(y b ) 1). We then assume the relation between y b and y t to be with a small parameter ξ near saturation that will be linked to another small parameter δ later. To find their relation, we rewrite (38) as where and the function R is defined by Also, notice that R can also be expressed at thermal equilibrium by R = R(y b ). The integration (91) can be expanded in small δ, and to the leading order they are R 1 R 2 y 1+ ∆ 2 t δ, and R 3 2y . With their expansions for small ξ and δ, from (90) and (92) the relation between ξ and δ is found to be ξ We then can divide the boundary time t (39) into four parts where and t s is the saturation time, obtained from (39) by setting y c = y b and E = 0. For the large R limit, since the relation between y b and R is identical to that in the sphere case in (42) and (47), the leading behaviour of t s is also the same as (82). For a positive (negative) ∆, the saturation time scale t s for fixed R and T becomes larger (smaller) than the case of ∆ = 0. For the small R limit, the relation between the tip of the extremal surface y b and the size of the boundary R in (58) allows to write down the approximate saturation time t s t s 2 ∆ 0 + 2 with R 0 defined in (58). Again, the contribution of ∆ = 0 to the powers of R can be checked from the saturation time t s in the field theories. For a positive (negative) ∆, in the small R limit, the saturation time t s becomes larger (smaller) than the case of ∆ = 0 for fixed R and T .
Expanding (94) for small δ gives where the function I is defined by . (96) After collecting the results of t 1 , t 2 and t 3 in their expansions for small δ and the relation between δ and ξ above, we arrive at where For continuous saturation in the case that t s is always great than t, u 1 should be smaller than zero. Nevertheless, when u 1 > 0 that has been discussed in [25], y c is still far away from y t near saturation so that y t might jump to y b at t = t s , causing discontinuous saturation.
Specially, the AdS case undergoes discontinuous saturation for both large R and small R limits. In the case of the general anisotropic model, the constraints on the scaling parameters due to u 1 < 0, giving continuous saturation are for large R, and for small R whereĨ 0 andR 0 are defined in (114) and (115) (see the details in Appendix).
Recall that the scaling parameters are constrained by (4) from the the black brane, (9),(10), and (11) due to null energy conditions, and (33) and (59)  To deal with the integration of area, it is convenient to reformulate (40) as where and the function A is defined by Again, A(y b ) = A Γeq A ∂Σ . Next we do the small δ expansion of A Γv<0 in (40) and A 1 , . Note that the leading order term of A 1 is the same as A Γv<0 in both h(y b ) → 1 (small R) or h(y b ) → 0 (large R). In the end, we conclude that in both large and small R limits, the area reaches its saturated value in the way of Following [25], straightforward calculations also shows that in the large and small R limits is equal to 1, which leads to A Γ → A Γeq quadratically in t − t s . However, in general R, the leading order terms of A 1 and A Γv<0 in the small δ expansion will be different, leading to the non-zero coefficient of the linear δ term in (104) In general, the powers of t − t s in A Γ − A Γeq are independent of the scaling parameters as well as the spatial dimensions [27].

VII. AN EXAMPLE FROM THE EINSTEIN-AXION-DILATON THEORY
Here we study the scaling parameter regions given by all sort of the constraints in the cases of the sphere and strip separately. As an example, we consider the anisotropy background in the Einstein-Axion-Dilaton theory [36]. The background metric obtained there is where C R and C z are constants, and f (r) . The temperature of the black brane is T = |d+1−θ| 4πr h . According to our notations, the scaling parameters read Thus, ∆ h = d + 1−θ z , ∆ = 2 z − 2, ∆ 0 = 0. Since ∂ v Θ(v) = 0 for v = 0, the null energy conditions in (9), (10), and (11) reduce to two constraints In the sphere case where a 1 = a 2 = ... = a d−2 , this has z = 1 and ∆ h = d + 1 − θ, ∆ = 0.
While z = 1, the null energy condition (110) gives the allowed ranges of the parameter θ as Also, for ∆ = 0, (63) leads to ∆ h > 3 consistent with the requirement of (4), giving with which, the constraint (64) holds for d ≥ 4. The constraint (68) is satisfied since ∆ 0 = 0 in this model. Together with (70), We summarize the scaling parameter θ in the Einstein-Axion-Dilaton theory for the sphere case by choosing different dimension d in Fig.(1). As for the strip case, collecting all constraints from (33) Fig.(2) for illustration. Also, due to (99) for large R and (100) for small R, within the allowed parameter regions, the model undergoes discontinuous saturation as in AdS space [25], although it involves more general scaling parameters. Note that the constraint in (68) is satisfied since ∆ 0 = 0 in the Einstein-Axion-Dilaton theory.
Together with (126) or (130), we find that the parameters in the blue regime always lead to discontinuous saturation for the system.

VIII. CONCLUSIONS
In this paper we employ the holographic method to study the thermalization of the In section (VI B), near saturation we assume y c is close to y t , and derive (97). However, in the case of u 1 > 0, y c is still far away from y b near saturation, and y t will jump to y b at t = t s . In this appendix, we will find the condition of u 1 < 0 in the large and small R limits.
To analyze R (y t ) and I(y b ) in (98) where β = y b y h . In (114) and (115), we have used the binomial identity and the Euler integral of the first kind. From (114) and R(y b ) = R, it is straightforward to obtain In the large R limit, the top of extremal surface y b is very close to horizon y h of β → 1.
In the small R limit, the top of extremal surface y b is very close to the boundary y = 0 of β → 0.
The leading order behaviour of I(y b ) in the large R limit relies on the asymptotic approximation in (115) at j → ∞ given bỹ Then in the large R limit (115) can be approximated by Apparently when β → 1, the summation above diverges. This divergence can be translated into the singular behavior of I(y b ) in the case of y b = y h (1 − ), when → 0, whereas the most singular behaviour is given by (119). Moreover, the last expression is obtained using the relation ln( ) = −2Rγ + O(R 0 ). Similarly, to discover what R (y b ) behaves in the large R limit, the associated large j behaviour in (116) is obtained as Then the leading terms when β → 1 becomes As a result, due to (120) and (124), the behavior of u 1 (98) in the large R limit is that Note that we have used h(y b ) ∆ h e −2γR . Since we assume ∆ h > 0, and from (33) having a 1 > 0, we find that In the small R limit, we can approximate (115) and (116) as and Then we find (98) in the small R limit to be By (33), (115) and (116) whereĨ 0 > 0,R 0 > 0 and a 1 > 0, thus 4a 1 (d−1) (2+∆)Ĩ 0R0 > 0. Finally, from (129), we conclude (2 + ∆)Ĩ 0R0 < 1 for the small R limit (130)