Thermalization in a Holographic Confining Gauge Theory

Time dependent perturbations of states in the holographic dual of a 3+1 dimensional confining theory are considered. The perturbations are induced by varying the coupling to the theory's most relevant operator. The dual gravitational theory belongs to a class of Einstein-dilaton theories which exhibit a mass gap at zero temperature and a first order deconfining phase transition at finite temperature. The perturbation is realized in various thermal bulk solutions by specifying time dependent boundary conditions on the scalar, and we solve the fully backreacted Einstein-dilaton equations of motion subject to these boundary conditions. We compute the characteristic time scale of many thermalization processes, noting that in every case we examine, this time scale is determined by the imaginary part of the lowest lying quasi-normal mode of the final state black brane. We quantify the dependence of this final state on parameters of the quench, and construct a dynamical phase diagram. Further support for a universal scaling regime in the abrupt quench limit is provided.


Contents
1 Overview

Holographic Thermalization
To date, holographic theories of strongly interacting matter have been deformed, probed, perturbed and otherwise studied in an enormous number of theoretical experiments. The results of these investigations have provided many important lessons about the nature of strongly coupled field theories, even in applications that fall outside the well established examples of AdS/CFT duality. While many earlier studies of holographic matter focused on the system's linear response to the insertion of various operators, quite a bit of attention has recently turned towards understanding the full non-linear dynamics of these holographic field theories. Broadly, the aim of these applications is to uncover results about the processes by which strongly coupled gauge theories respond to arbitrary time dependent perturbations. Often, as in the case of the ground-breaking early examples of numerical holography [1][2][3][4], the gauge theory responds by "thermalizing" the perturbation. This means that after some characteristic time scale, typically set by the Hawking temperature of a black hole in the dual gravity solution, the field theory arrives at a final static state in thermodynamic equilibrium. In this equilibrium state, all correlation functions assume their thermal values, which is to say the trace is taken with respect to the thermal probability distribution appropriate to the ensemble. In other examples [5][6][7], perturbations of the gauge theory have been found which apparently never thermalize. These "islands of stability" correspond to initial gravitational data that never leads to gravitational collapse and horizon formation.
A conceptually simple means of perturbing a gauge theory is to turn on a time dependent source for some operator and evaluate the system's response. Beyond the linear regime, such a source can be used to continuously drive the system, or to "quench" it. By the latter, one typically means that a coupling in the gauge theory is varied over a compact timescaleτ , whereτ is sometimes chosen to be small compared to other scales in the theory.
When quenched in this way, it is clear that the system's response will depend on the properties of the operator to which the source couples as well as the parameters which specify the quench. In the weak field approach pioneered in [8], perturbative results have been obtained for the quench of a marginal operator in the case where the perturbation's amplitude and characteristic timescale are both sufficiently small. This analysis allowed the authors to study the formation of black holes and black branes in asymptotically Anti-de Sitter spacetimes.
Specifically, they constructed the limiting behavior of the dynamical phase diagram for the outcome of massless scalar collapse in global and Poincare patch AdS. Notably, they found that arbitrarily small perturbations of the Poincare patch solution always resulted in the formation of a black brane, whereas finite volume effects in the global solution prevented horizon formation under certain conditions. Interestingly, they also observe a crossover between final states with small and large black holes, as well as Choptuik behavior in the vicinity of the transition between final states with a horizon and those without. This scaling behavior strongly suggests the presence of a second order line separating these phases. More recently, the weak field approach has been extended to provide interesting insights into thermalization in non-isotropic quenches [9], in finite density states of a gauge theory [10], and in a simple model of strongly coupled matter with a mass gap [11].
The holographic toolkit makes it conceptually straightforward to move beyond the weak field perturbative scheme as well as to generalize the quench to perturbations by relevant operators. Along the first direction, calculating beyond the reach of perturbative methods clearly implies grappling with the full non-linearity of the Einstein equations. In general, this is an exercise in "numerical holography", loosely defined as the set of all holographic calculations that require numerical techniques more involved than a call to Mathematica's NDSOLVE. Instead, one typically appeals to any one of an assortment of numerical methods specifically tailored to the problem of gravitational in-fall.
A recent example of holographic thermalization which connects perturbative weak field calculations to the full numerical evolution of a marginal perturbation appears in [7,11]. In this setup, the authors consider the response of a confining gauge theory to time dependent perturbations by studying the dynamical evolution of a massless scalar quench in the AdS hardwall geometry. This gravitational theory is characterized by a length scale z 0 = 1/Λ at which the geometry is artificially terminated. In turn, this length scale gives rise to a mass gap ∼ Λ in the dual gauge theory, as well as a temperature above which a large black brane can appear in the bulk. Accordingly, the AdS hardwall serves as a crude model for a confining gauge theory with a deconfined phase dual to the large black brane solution.
Weak field calculations in the hardwall background indicated that the dynamical phase diagram ought to include a transition between perturbations which result in black brane formation and those that remain in a horizon-less scattering state. This should be contrasted with the situation in [8] in which such a distinction between final states was only possible when the dual field theory is defined on a sphere. Unfortunately, the weak field approach was unable to answer questions about the late time fate of the scattering solutions, which is of particular importance for the field theory interpretation. To resolve the late time behavior of the scattering solutions, a numerical method was introduced to solve the Einstein equations and evolve the system arbitrarily far forward in time.
A particularly noteworthy result of the numerical investigation was the observation that there exist perturbations such that the subsequent scattering solutions never undergo gravitational collapse. Since no horizon is formed in such a process, the implication is that the holographically dual strongly coupled matter does not achieve thermodynamic equilibrium at late times. This finding differs in important ways from the islands of stability identified in scalar collapse in global AdS. In particular, it demonstrates that a perturbation initiated by explicitly sourcing a marginal operator need not necessarily thermalize in the infinite volume gauge theory. To rephrase this result in a way that will be more directly applicable to our present work, the introduction of a mass gap in the dual gauge theory can strongly effect the thermalization time for a certain class of perturbations.
A related, but logically distinct line of research involves the deformation of strongly coupled matter by time dependent perturbations of a relevant scalar operator. Numerical investigations into this situation appear in both systems with a finite density of charge carriers (for example [12]) and without [13][14][15]. We will focus on the latter scenario, in which an uncharged black brane solution is perturbed by varying the nonnormalizable boundary mode of a massive bulk scalar in time. Holographically, the dual picture is that an initial state of the gauge theory in thermodynamic equilibrium is deformed by turning on a relevant operator over some timescaleτ . The strongly coupled matter generically responds by passing through a non-linear regime before settling once more into thermodynamic equilibrium, albeit in a different thermodynamic macrostate.
One of the most interesting results to emerge from the numerics of [13][14][15] has been the appearance of a "universal fast quench regime" in which the change in energy density E after the quench scaled as a power law in the quench widthτ , (1.1) This scaling is argued to manifest when the quench width is small compared to any other scales in the theory, such as the characteristic amplitude of the perturbationδ. Moreover, the scaling was shown to be robust for many choices of the relevant operator's conformal dimension ∆, and its appearance was further elucidated analytically in a recent series of papers [16][17][18]. Given the non-trivial dynamical phase structure and strong dependence of the thermalization time on perturbations in the hardwall model, as well as the appearance of universal regimes in the space of quench parameters, it is natural to wonder how these properties might manifest in a holographic dual to a gauge theory that is more similar to SU (3) Yang-Mills. In a sense we make more precise in the following subsection, this means we are interested in finding a natural way to "soften" the hardwall while retaining the mass gap ∼ Λ that was responsible for the interesting dynamical features of the quench. The construction and subsequent perturbation of such a holographic model will be the focus of this work.

Confining Gauge Theories and Holography
Holographic models of confining gauge theories generically involve a gravitational theory whose solutions include those which explicitly break scale invariance in the bulk. This condition is necessary but not sufficient to allow for non-trivial temperature dependence of states in the dual gauge theory, as well as a mass gap at zero temperature. These models also include solutions whose metric breaks bulk conformal invariance, the most relevant of which will be branches of black hole solutions which are holographically dual to thermal states with non-zero stress-energy. The deconfinement transition will be a thermodynamic phase transition from solutions without a black hole, to solutions with one.
In what follows, we will focus on a particular class of holographic models for confining gauge theories which are in the Einstein-Dilaton family. As this title implies, our model will involve a non-trivial scalar whose profile explicitly breaks conformal invariance in the radial direction, and can be tuned to produce gravitational solutions dual to states on either side of a deconfinement transition. Within this rather simple class of models there exists a remarkable wealth of phenomenological possibilities, which include sharp deconfining transitions of first or second order, crossover behavior between phases, and theories driven from their UV fixed point by relevant operators of arbitrary dimension.
When the gravitational theory is dimensionally reduced to five dimensions, all such bulk solutions as described above (with the exception of flows to a non-trivial IR CFT which are not relevant for holographic models of Yang-Mills) will have a naked singularity in the IR [19,20]. The presence of this IR singularity is a signal that something is missing in the effective holographic description. In known examples, typically the missing ingredients are an assortment of other bulk fields that will necessarily obtain vevs in the ground state solution. If such fields are reintroduced, the singularity is expected to be resolved. In [21], Gubser provides a simple criterion for a resolvable singularity which states that for a singularity to be resolvable it should be possible to construct a solution with an infinitesimal regular horizon surrounding the singularity. Such singularities are called "good" or resolvable singularities. Even in such "good" solutions, holographic calculations of correlators may be illdefined in the presence of the naked IR singularity. This happens when in the neighborhood of the singularity the Sturm-Liouville problem that governs the bulk fluctuation has two normalizable solutions. In such a case an extra boundary condition is needed to determine the solution, and this is the signal that the correct solution is sensitive to the details of the singularity's resolution. Such cases are called "holographically unreliable" because without knowing the details of the singularity's resolution holographic computations are inherently ambiguous. There are, however, numerous cases where the Sturm-Liouville problem near the naked singularity has only a single normalizable solution. In that case the solution is unique and does not depend on the resolution of the singularity. When this occurs, the singularity is said to be repulsive, [19,20] and boundary theory correlation functions can be reliably calculated from holography in such backgrounds.
In confining states of a gauge theory, the expectation value of Wilson loops exhibit area law scaling. Holographically, the Wilson loop has a natural definition as a minimal surface in the bulk affixed to the edges of a Wilson loop on the boundary. In [22], it was shown that holographic Wilson loops will give area law behavior providing that the string frame metric scale factor has a minimum, and that the scale factor at this minimum is non-zero. In Einstein-Dilaton theories, this requirement can be reformulated as a constraint on the IR behavior of the dilaton potential, as we will discuss in more detail below.
In our model, the zero temperature, zero entropy 1 ground state of the gauge theory will be dual to a bulk solution with a singularity in the far IR, but no event horizon. This singularity allows the bulk scalar to diverge in the IR, as well as the spatial metric scale factor to vanish. 2 The latter is the feature responsible for vanishing entropy density, and is a familiar aspect of many candidate holographic ground states. As discussed above, the IR singularity can be thought of as a rather innocuous consequence of working in a limit of non-critical string theory, providing the singularity satisfies several criteria. Our singularity is of this "good" type.
As the temperature increases in the gauge theory, no new solutions are immediately available in the dual gravitational theory. Instead, the dual bulk solutions are given by the zero-temperature (confining) solution with its time coordinate compactified into the thermal circle. The temperature T associated to this solution is governed by the size of the thermal circle in the familiar way, and such a state is often called the "thermal gas" saddle point. As it lacks a black hole horizon, it is obvious that the thermal gas is characterized by thermodynamics subleading in the large N c (where N c 1 is roughly the rank of the dual gauge theory) expansion inherent to our holographic setup.
Increasing the temperature further, one eventually arrives at a special temperature T = T 0 which marks the appearance of a new branch of solutions in the bulk theory. These solutions are of the black brane type, with a non-compact horizon of planar topology. As seen in the left plot of figure 3 for any T > T 0 , there are two black hole solutions. One is in the left branch of the curve while the other is in the right branch of the curve. The black holes on the left branch are called large black holes as their horizon size is larger than any of those in the right-branch. They have positive specific heat and are therefore locally thermodynamically stable. On the contrary, the small black hole branch on the right has negative specific heat and the associated black holes are thermodynamically unstable. Qualitatively the diagram is similar to that of global AdS space, although here the space is flat and the volume infinite. In the case of global AdS, small black holes asymptote to Schwarzschild black holes in flat space in the limit of small horizon size (with T → ∞). Here also, small black holes in the limit of small horizon size have T → ∞ but they are very different from Schwarzschild black holes, [20,23,24]. Finally, the unique black hole that is at T = T 0 is very special, as here the specific heat diverges.
The thermodynamics of such black holes is governed by the standard black-brane relations, which implies they have an entropy density proportional to L 3 /κ 2 where L is the characteristic length scale of the metric (for example the AdS radius) and κ is the five-dimensional gravitational constant. A standard application of the holographic dictionary for known five-dimensional gravitational duals gives L 3 /κ 2 ∝ N 2 c . The proportionality constant depends on the details of the duality under consideration. In holographic applications that do not descend from a known string theory solution, the proportionality constant is an undetermined parameter of the theory, but the N 2 c scaling of various thermodynamic quantities is expected to be robust.
As mentioned above, by studying the susceptibilities of the small black hole branch, for example the specific heat, it is straightforward to show that these black brane so- Figure 1. The phase diagram of our model. The thermal gas solution (grey) is the only solution below T = T 0 . Above T 0 the black hole branches appear (black plane), but the thermal gas is thermodynamically preferred. At T = T c there is a first order deconfinement transition from the thermal gas to large black brane solutions (large black plane). Cartoon from [26].
lutions are thermodynamically unstable. This is a local statement, independent of the global thermodynamic (in)stability of these solutions. In light of the Correlated Stability Conjecture (CSC) [25], one might worry that these solutions are also dynamically unstable, which is to say their fluctuation spectrum may contain a mode in the upperhalf complex frequency plane. Such an instability obviously leads to an exponential growth in time, and the conjecture posits that the existence of such a mode indicates a Gregory-Laflamme type instability towards a "lumpy" horizon. This dynamical instability will play no role in the present work as we will study s-wave perturbations of spatially homogeneous solutions.
Above T 0 , the thermal gas still provides the thermodynamically-dominant solution, as it has lower free energy than the black brane solutions at the same temperature. This remains true until T = T c , at which point the model predicts a first order phase transition from the thermal gas to the black-brane solutions. Above T c , the system is described by the large black brane solutions that are thermodynamically stable. This first order phase transition is a holographic realization of a deconfining transition to states with an entropy density that scales like N 2 c and perimeter law behavior of the Wilson loop. The full phase structure of our model is illustrated in figure 1.

An Einstein-Dilaton Model of Thermalization in a Confining Theory
Our goal in this work is to continue the investigation of the effects of confinement on thermalization processes in holographic gauge theories. We accomplish this by perturbing a state of a confining gauge theory with a relevant operator that we loosely associate with TrF 2 , and then studying the system's response. From the bulk perspective, we are choosing a solution of our model and deforming it briefly through time-dependent boundary conditions on the dilaton. The deformation will generically throw the system out of equilibrium, which will then undergo non-linear evolution governed by the Einstein-Dilaton equations.
An important question is whether or not an arbitrary perturbation of a confined state in the strongly interacting gauge theory will thermalize. In this context, "to thermalize" means that the late time behavior of the state is in thermodynamic equilibrium, and all correlation functions are time independent and assume their thermal values. The gravitational formulation of this question is whether or not an arbitrary perturbation of the horizon-less thermal gas solutions of our model necessarily results in black brane formation. Put another way, is the exponentially diverging dilaton potential of our model sufficiently similar to the hard wall of [7,11] to encourage scattering solutions and prohibit horizon formation under certain conditions?
At present, the answers to these questions remain beyond our grasp. In the absence of a horizon in the initial state, the diverging scalar in the IR leads to a number of challenging numerical obstacles. To address these issues, one should regulate the bulk solutions in the IR through the introduction of a numerical cut-off, and study the effects of the position of this cut-off on the computation's results. A related approach, which we will adopt, forfeits the ability to comment on black hole formation while hopefully retaining qualitative effects of the confining dilaton potential. This scheme replaces the hard numerical cut-off with a small black brane horizon, which then acts as a natural IR regulator.
In our model, this necessarily implies that our initial bulk solution will not be holographically dual to the ground state of the gauge theory. At first glance it actually looks much worse, since these small black brane solutions aren't the thermodynamically preferred state at sufficiently low temperature, nor are they necessarily dynamically stable. Nonetheless, we argue that they provide a sensible starting point for our calculation, as in the limit of vanishing horizon size these black holes coincide with the ground state of the system.
The small black-brane solutions asymptotically approach the zero-temperature solution in the limit where the horizon recedes infinitely far from the boundary into the IR. Thus, provided we accept the concession that a horizon has already been formed in the bulk, for sufficiently large separation between the UV boundary and the position of the horizon the gravitational solution shares many important features with the ground state solution. The most important of these are a scalar profile which is rapidly growing towards the IR, and a rapidly vanishing warp factor in the Einstein frame metric. We expect that these features will be responsible for most of the interesting dynamics in our system, such as the timescale with which thermalization takes place, and accordingly that the small black holes can provide us with a qualitative understanding of the consequences of these features.
Moreover, the thermodynamic instability of these black branes and the possible presence of a dynamical instability in their fluctuation spectrum posited by the CSC is not expected to have much relevance to our discussion. Gregory-Laflamme type instabilities are related to a "clumping" of mass or charge at the horizon, and should thus involve finite wave-number perturbations in the directions transverse to the brane. We study translationally invariant perturbations in the non-critical (five dimensional) theory, and hence do not anticipate the possibility of exciting such instabilities. Indeed, we have not yet seen any direct evidence for a dynamical instability in the small black hole branch of our model.
The main results of our study are a portion of the dynamical phase diagram, and the scaling properties of various equilibration processes in our model. The former is realized as a map between perturbation parameters and the final state achieved after equilibration. Our perturbation will be controlled by two parameters corresponding to the amplitude of the perturbation and its duration. Thus, the dynamical phase diagram is a two dimensional plot with a curve marking the boundary between perturbations that result in a small black hole and those that thermalize into large black holes.
By focusing on the details of the final state, we are also able to make claims about the dependence of the characteristic thermalization time in the boundary gauge theory on the quench parameters. Not surprisingly, this timescale turns out to be quantified by the quasi-normal mode spectrum of the final state solution. Moreover, by quantifying the relationship between the features of the quench and the change in energy density induced by the perturbation, we demonstrate explicitly that the simple scaling anticipated in (1.1) appears in our model as well.
The remainder of this work is devoted to detailing the setup, implementation, and interpretation of time dependent quenches of a holographic confining gauge theory. In section 2 we introduce the specifics of the bulk theory we are interested in perturbing, as well as the equations of motion that govern its behavior. The behavior of these solutions near the conformal boundary are of particular importance for the application of holographic methods, and thus we provide these near boundary solutions in some detail.
The numerical construction and thermodynamic properties of the static, equilibrium solutions of our model are discussed in section 3. This section primarily serves as an orientation to our holographic gauge theory. In it one finds the dependence of the free energy and entropy density of our strongly coupled matter on temperature, as well as a comparison of the temperature dependence of the speed of sound in our model to that of SU (3) Yang-Mills. It is complemented by the results of section 4 which further characterize the static states in terms of their thermodynamic one-point functions.
Section 5 details the computational approach we adopt to evolve our gravitational system in time. This approach is based on an assortment of well known numerical techniques that we carefully tune to accommodate the specifics of our model. Most notably, we explain how we handle the copious logarithmic fall-offs in the near boundary behavior of bulk fields. These fall offs are generic to gravitational theories in odd bulk dimensions, and present obstacles to the accurate determination of the normalizable and non-normalizable UV coefficients of the various bulk fields. This section will be primarily interesting to those who would like to numerically study dynamical quenches in related models.
The output from our numerical method and the interpretation of this output is contained in sections 6 and 7. These sections constitute the primary results of our study. They include examples of the response our model to various classes of quench by a relevant scalar operator, as well as a dynamical phase diagram for the outcome of these quenches when performed in a particular initial state. We examine the dependence of the final equilibrium state on the quench parameters, and comment on the appearance of an anticipated universal scaling regime in the fast quench limit. The connection and applicability of these results to similar processes in other theories is discussed in section 8. Specifically, we comment on the implications of our calculations for the thermalization of probes in the strongly coupled matter produced in heavy ion collisions, as well as future directions one might wish to pursue.

The Action
We will be interested in Einstein-Dilaton theories that are tuned to qualitatively reproduce some important properties of QCD. These theories can be described by an action of the form where V (ϕ) is the dilaton potential which will govern the dynamics of the system, and K is the Gibbons-Hawking-York term necessary to help render the variational problem well defined on the boundary. The model we study can be classified as a "bottom-up" description of a confining gauge theory, in the sense that V (ϕ) is not derived from a known supergravity theory. Instead, this potential is manufactured to satisfy certain criteria in the IR and/or UV limits of the gravitational theory so as to induce desirable features in the dual boundary theory. For example, these asymptotics determine whether or not there is a mass gap, while other asymptotics determine the scaling dimension of the dual scalar operator. These criteria have been elucidated and categorized in a series of papers beginning with [19,20]. Importantly, the asymptotic behaviors of the scalar potential that we choose for the current work are present in known solutions of gauged supergravity. An example of this which is closely related to the model investigated in this work is described in [27]. As many qualitative features of holographic models are dominated by these asymptotic properties, it is reasonable to expect that predictions from our model may apply to a broad class of holographic systems, including those with a more esteemed string theory lineage. While there are a variety of scalar potentials whose distinct IR asymptotics lead to confining theories, it was found in [19,20] that a particular IR behavior also produces linear radial trajectories for glueballs, as well as the structure of the Yang-Mills phase diagram just above the first order deconfinement transition. This behavior requires which, apart from the subleading √ ϕ, is the non-critical string theory dilaton potential in five dimensions in the Einstein frame.
At high energies, the UV properties of the gauge theory depend on the nature of the scalar operator which perturbs it. In [19,20] a marginal operator was used, and the potential chosen such that the UV fixed point was at ϕ → −∞. The potential is of the form so that λ → 0 where λ in the UV is identified with the 't Hooft coupling constant. This fixed point is extremely shallow as all finite derivatives of the potential vanish at the fixed point. This potential is tuned to match the thermodynamics of SU (3) Yang-Mills, and the resulting theory, dubbed Improved Holographic QCD (IHQCD), is capable of quantitatively reproducing many features of QCD across all energy scales. Clearly it would be desirable to study dynamical processes in IHQCD, beyond the linearized level. However, the presence of a marginally relevant operator in the UV creates serious computational difficulties. One consequence of the marginal deformation inherent to IHQCD is that the scalar λ vanishes logarithmically near the UV boundary. This mild falloff is very difficult to handle numerically, as it requires retaining a very high level of numerical precision to correctly identify the leading and subleading coefficients governing the near boundary behavior of the solution. To circumvent this issue, we will sacrifice the quantitative match to Yang-Mills theory provided by IHQCD in favor of computational convenience.
One can do this without departing drastically from the coarse features of Yang-Mills by following the approach advocated in [28]. In this approach, the marginal scalar operator dual to λ is traded for a relevant operator with dimension not far from 4, which one hopes to roughly identify with a boundary operator of the form TrF 2 where F is the Yang-Mills field strength. The fact that this operator is no longer marginal is meant to capture the anomalous dimension that the operator acquires after running some ways towards the IR, reminiscent of what happens in Yang-Mills.
We will therefore assume that there is a regular UV fixed point at ϕ = 0, without loss of generality, and the potential near this fixed point takes the standard form The various coefficients V (i) are fixed by the symmetries of the gauge theory and were shown in [19,20,29] to be in one to one correspondence to the coefficients of the holographic β-function. The AdS scale L of the UV AdS space is determined by L 2 V (0) = 12. The relationship between the mass of the scalar and the conformal dimension of the dual gauge theory operator ∆ can be made precise. In the standard quantization ∆ is defined to be the larger of the two roots of An example of a potential for which these properties are realized first appeared in [28]: and a parameter set which produces a bulk theory dual to a confining gauge theory deformed by a dimension ∆ = 3 operator is (a, b) = (1/500, 10009/1500). This theory is decidedly not real world YM. Unlike the more refined IHQCD models of [30], it fails to produce a very good match to the thermodynamics of real world SU (3) Yang-Mills (see section 3.1). Nevertheless, this holographic theory does exhibit a mass gap, the thermodynamic properties determined by its potential are qualitatively reminiscent of SU (3) glue, and it renders the bulk system relatively amenable to numerical investigation. Accordingly, in what follows we will primarily be concerned with solutions to the equations of motion derived from (2.1) with the potential written in (2.6).

Solutions
Among the solutions to this model's equations of motion are the planar geometries (with and without a horizon) of the form where v is an ingoing null-coordinate, z is the radial direction, and x describe the planar R 3 . The various metric functions appearing in this ansatz, as well as the scalar, are taken to be functions of both v and z. The solutions satisfy the equations of motion where the prime ( ) denotes differentiation with respect to z, and d + is a modified derivative defined by which is the derivative along the outgoing null vector.

Boundary Analysis
Near the UV boundary, the various metric functions and a bulk scalar with m 2 = −3/L 2 can be expanded like and we adopt coordinates in which the asymptotically AdS 5 boundary has a 0 (v) = s 0 (v) = 1 and α 0m (v) = σ 0m (v) = 0. In general, logarithmic terms are expected in even boundary theory dimensions, d, and when the scalar operator has dimension ∆ such that ∆ − d/2 is an integer. The scalar's source, f 0 (v), remains a boundary condition to be implemented. Inserting these expansions into the Einstein equations and expanding near the boundary at z = 0 provides the asymptotic behaviors of the various fields. In what follows, we leave the potential largely unspecified, requiring only that it gives rise to a dual conformal gauge theory which is deformed by a parity invariant relevant operator of dimension ∆ = 3. This constrains a few terms in an expansion of the dilaton potential about the AdS fixed point like for positive integers, n. Accordingly, the fields behave like this: where The functions ζ, a 4 , f 0 and f 2 all depend on time, and are not determined by the asymptotic series expansion. The function ζ is completely unfixed, and reflects residual reparametrization invariance in z. The coefficient a 4 is related to the background energy density, and f 0 and f 2 can be related to the source for and response of the dual scalar operator, respectively. The AdS radius L has been (and will continue to be) set to one, and can be reinstated in any formula by dimensional analysis. In order for the full set of Einstein equations to be consistently solved, the time derivative of a 4 is subject to a constraint, This constraint is related to the Ward identity governing the divergence of the stress tensor, which we revisit in section 4. From these expressions it is simple to work out the asymptotic behavior of the modified derivatives d + Σ and d + ϕ. They are where C Σ is the exhausting constant For later convenience, we write these asymptotic expansions in the form Solutions with non-zero f 0 correspond to adding to the boundary Lagrangian a term like for mass scale Λ. For the model of interest (2.6), we have ∆ = 3 and thus the source for the relevant perturbation has dimensions of energy. Regularity of the solutions in the IR implies one constraint between the UV coefficients f 0 , f 2 and a 4 , and thus the UV deformations are specified by a single dimensionless quantity which we can take to be T tt /f 4 0 . In other words, we anticipate a one parameter family of black hole solutions characterized by their energy density in units of f 0 . Often times it will be more instructive to parametrize various solutions by the value their scalar obtains in the IR, a method introduced in the following section.

The Initial State
We eventually wish to thermalize the scalar perturbation in an initially static state of the gauge theory dual to the bulk theory described by (2.1) and (2.6), at finite temperature. These states are constructed by restricting the background functions in the ansatz to vary only radially, that is removing v-dependent terms in (2.8-2.12), and then integrating the resulting system of ordinary differential equations to obtain solutions for A(z), Σ(z) and ϕ(z). These bulk solutions then provide the starting point for the subsequent deformation and evolution.
To construct numerical solutions to the equations of motion (2.8-2.12) describing the initial state, we exploit the fact that we desire a regular horizon at z = z H and thus expect that the background functions can be expanded around the horizon like where F is any of {A, Σ, ϕ}. The coefficients F i can be partially fixed by appealing to symmetries of the underlying gravity theory. For example, A 0 = 0 by assumption of a regular horizon, and Σ 0 = 1 can be enforced by rescalings of the spatial coordinates. As mentioned previously, the metric ansatz (2.7) does not fully fix the gauge freedoms of the gravitational theory, and this is reflected in horizon data as the ability to choose A 1 arbitrarily. In practice, we have found it convenient to fix this residual reparametrization invariance by requiring that the boundary coefficient ζ of the previous section vanishes. The value of the scalar at the horizon, ϕ 0 ≡ ϕ H thus parametrizes the available solutions in this model. All the higher order coefficients are fixed in terms of these first few coefficients, and accordingly one can use the expansion to sufficiently high order to provide IR data at a radial location near the horizon with which to seed a numerical routine.

Thermodynamics
The solutions to the equations of motion are asymptotically AdS 5 near the boundary at z = 0, with a scalar that vanishes linearly. However, as our integration strategy relies on fixing the values of various parameters at the horizon, it will often happen that the numerical solutions obtained are characterized by different energy scales, f 0 and appear in different coordinates. More explicitly, static solutions near the boundary generically behave as for constants Σ F andf 0 . As the aim is to compare states of different temperature in the same theory, we should arrange that the solutions take a canonical form at the boundary, with identical metric and energy scale. This can be accomplished by performing a coordinate transformation to a new radial coordinate in whichf 0 = 1, as well as simple rescalings of the other coordinates. To wit, the transformatioñ ensures a near boundary solution of the form which is then suitable for comparing thermodynamic properties between solutions. Following the standard prescriptions for black brane thermodynamics, temperatures and entropies can be readily extracted from the near horizon geometry. In the "canonical" coordinates of (3.4-3.6), the temperature is easily computed from the surface gravityκ at the horizon since and ξ is a unit Killing vector, timelike outside the Killing horizon atz =z H . Meanwhile, the entropy density is proportional to the area of the horizon, so that 10) Figure 3. Plots of the temperature scaled by the critical temperature as a function of λ H /λ c (left) and the entropy density scaled by the third power of the temperature as a function of T /T c (right). The rightmost plot becomes "dotted" as one passes through the phase transition by lowering the temperature from above. This is meant to indicate that in the field theory, this low temperature phase is governed by the thermal gas solutions, whose entropy is subleading in the number of colors N c .
where s is the entropy density andγ is the determinant of the spatial part of the metric at the horizon. For the solutions we consider, inserting the near horizon form of the canonical metric into (3.10) gives As different solutions are fully characterized by the value of their scalar at the horizon, it is often useful to think of their thermodynamic properties as functions of λ H = exp ϕ H . In this spirit, once the temperature and entropy of a given state can be reliably computed, one can measure the free energy of the solution from the integrated first law: For the model considered here, the free energy density is plotted as a function of T in figure 2. From these plots, one learns that there exists a critical λ H ≡ λ c , located where the free energy changes sign, at which there is a first order phase transition from the thermal gas to the black brane phase. One also notes the presence of two black hole branches in the plot of F as a function of T . Since the specific heat is given by wherein T is manifestly positive, one finds that the "upper" branch of black hole solutions, with positive curvature, have C v < 0 and are thus thermodynamically unstable. Accordingly we will anoint this the "small black hole" branch in analogy to similar solutions in global AdS. In figure 3 the temperature is plotted as a function of λ H , and the entropy density is plotted against the temperature. From the former one notes the existence of a minimum temperature, which we denote T 0 , while from the latter one finds that this model is characterized by a large jump in the entropy density at the first order phase transition. Moreover, as the squared speed of sound is just 14) this model attains the high-T conformal value of c 2 s = 1/3 more rapidly than anticipated by SU (3) glue on the lattice. The results are shown in figure 4.
In what follows, we will focus our attention on finite temperature initial states in thermal equilibrium. This is a subset of thermalization processes which excludes the possibility of commenting on a variety of interesting questions related to the conditions required for the formation of a black hole in this theory. Included in this subset, however, are processes which probe the non-trivial gravitational phase structure of the bulk theory. Generically, one expects that this phase structure will result in an interplay between physics on the large and small black hole branches, which in turn might better inform our understanding of thermalization in confining gauge theories.

More Boundary Theory Observables
The dual gauge theory data is encoded in correlation functions of the boundary theory operators. In this work, we will primarily be interested in the one-point functions of the stress-energy tensor T ij and the dimension three operator O. Holographically, the values of these one point functions are related to the boundary coefficients of the normalizable modes of the metric and bulk scalar, a 4 and f 2 respectively.

Renormalized One-Point Functions
The precise relationship, however, requires a careful analysis of the near boundary onshell bulk action, which generically has both power law and logarithmic divergences at z = 0. These divergences can be regularized and renormalized following the standard dogma of Holographic Renormalization [32,33], the end result being a (possibly scheme dependent) identification of bulk falloffs with boundary theory correlation functions. This procedure is by now a familiar aspect of many holographic calculations, and accordingly the details will be left to appendix A.
The main results of this analysis are the set of local counterterms required to regulate the on shell action, and the one point functions of the stress energy tensor and the scalar operator as functions of the bulk boundary data. In the present system, with flat boundary metric and ϕ dual to a dimension three operator, the relevant counterterms turn out to be where γ is the pull back of the metric to the radial cutoff at , A is a set of finite counterterms which define the renormalization scheme, and F 4 = 16/27 − V (4) /24. Absent supersymmetry, or another guiding principle with which to fix the coefficients of the finite counterterms, we shall henceforth adopt a holographic minimal subtraction scheme and simply ignore them. The corresponding one point functions are then given by where we have used (A.25) to express the correlation functions in terms of the coefficients in (2.20). As expected, these correlation functions are constrained by the presence of Ward identities. In terms of these near boundary expansion coefficients, the conformal Ward identity reads while the divergence of the stress tensor yields which illustrates the non-conservation of the system's energy in the presence of time dependent sources. To better interpret the values of these various one-point functions, it is convenient to define subtracted correlators which effectively measure deviations of the energy, pressure, and/or scalar expectation value from the static zero temperature solution obtained in the limit ϕ H → ∞. One way to do this is to modify directly the renormalized on shell action by addition of appropriate finite local counterterms. For example, a non-zero A = c ϕ ϕ 4 introduces a parameter c ϕ which can be tuned such that the energy density of the zero temperature solution vanishes. Since the counterterms are properties of the theory and not a specific solution, this subtraction will be manifest in all one-point functions computed in the boundary theory.
A related, and perhaps more practical prescription is to simply define new correlators with the limiting ϕ H → ∞ values explicitly removed. Accordingly, in the examples that follow we often construct "hatted" one-point functions defined such that Ω ≡ Ω − Ω ϕ H =∞ (4.7) for any boundary theory operator Ω. In figure 5 we plot the hatted correlators as functions of the value the scalar obtains in the IR for future reference.
The numerical procedure we adopt is initialized by the boundary coefficients a 4 , f 2 and f 0 , which are in turn extracted from numerically generated initial states as described in the previous section. For this reason, the evolution's stability depends crucially on the accuracy of these values. It is thus reassuring that we find excellent agreement between the free energy density computed from (3.12), and from boundary data using F = −p ≡ − T xx .
We also compute the dependence of the energy density T tt on the system's entropy density s. The former is calculated from the UV boundary coefficients extracted from our numerical solutions, while the latter is computed from horizon data. The result appears in figure 6, in units off 0 and κ = 1. At large energy density, our model correctly reproduces the expectation for a conformal theory, T tt ∝ s 4/3 . As the energy density decreases, the dimensionfull source's explicit breaking of the UV theory's conformal invariance becomes increasingly important. In the extremal limit, the small black hole branch is characterized by a logarithmic dependence of the energy density on entropy density, T tt ∝ s √ − ln s. From the point of view of the boundary gauge theory, the latter behavior does not manifest as the thermodynamically preferred solution is the thermal gas whose entropy density is subleading in N c .

Numerical Strategy
The primary goal of the present work is to discover how an initial state described by a solution from section 3.1 responds to time dependent perturbation by a dimension three operator. Holographically, we implement this perturbation by varying the UV boundary condition on the scalar in time. For example, a source whose profile is given by corresponds to perturbing the initial state (with coupling to O fixed byf 0 ≡ f 0 (−∞)) by a gaussian deformation with amplitude −δf 0 and variance τ 2 centered at v = 0. Accordingly, the boundary theory experiment we have in mind is dialing down the coupling to the relevant operator by an amount δf 0 and for a time τ .
In the remainder of this section, a recipe for computing the time evolution of an initial state in the presence of time dependent scalar source f 0 (v) is explained from both generic and practical viewpoints.

Integration Routine
The equations of motion (2.8-2.12) have been arranged into a "nested" structure convenient for numerical study. We employ and consequently review the method detailed in [34], which exploits this nested structure to reduce the partial differential equations to a sequence of ordinary differential equations that can be solved time step by time step.
Before addressing the task of solving the Einstein equations describing this system, it is important to understand the boundary information needed to fully specify a given solution. For this, one may turn to the near boundary behavior of the homogeneous Einstein equations. For example, in the UV equation (2.10) is which has the general solution Thus, near the boundary linearly independent solutions behave like z −1 and z 0 . Comparing to (2.19) it is clear that Σ can be fully specified by the first two terms in the UV expansion-the AdS boundary condition and the radial reparametrization artifact ζ. Similarly, from the unsourced (2.9) one finds near the boundary which has the solution d + Σ ∼ z 2 . From (2.26) it is clear that one must specify the value of C Σ in order to determine d + Σ, which in turn depends on f 0 , f 2 , and a 4 . Finally, from the unsourced (2.8) in the UV one solves with d + ϕ ∼ z 3/2 . From (2.27), evidently the desired solution is the one in which the coefficient of z 3/2 vanishes. In principle, the same sort of analysis can be applied to (2.11) to explore the boundary behavior of A. The homogeneous equation is readily solved by A ∼ z −1 + z 0 , which (2.18) shows are fixed by ζ andζ on the given time slice. Since ζ parameterizes a residual gauge degree of freedom, one may adopt a gauge whereζ = 0. Alternatively, it may be desirable to determine ζ(v) dynamically. In this case one may use the value of A on the apparent horizon instead ofζ as an integration constant. In practice this can be accomplished by immobilizing the location of the horizon (e.g. requiring z H = 1 always). The condition for the location of an apparent horizon in these backgrounds (see Appendix C) is simply and the horizon stationarity equation thus requires The upshot of the preceding analysis is that solving the Einstein equations requires knowledge of a 4 , ζ and ϕ(v 0 , z)-from which f 2 can be extracted-on the initial time slice, together with a choice of the forcing function f 0 (v).
After obtaining this initial data, one is well poised to evolve the system. The procedure is as follows: 1. From ϕ(v 0 , z), the AdS boundary condition s 0 (v 0 ) = 1, and the value of ζ(v 0 ), equation (2.10) can be integrated to obtain Σ(v 0 , z).

With this knowledge, and the value of
3. Then equation (2.8) can be solved with the boundary condition that d + ϕ has no fall off like z 3/2 at the boundary, resulting in the knowledge of d + ϕ(v 0 , z).
4. Finally, we turn to the second order ODE given by equation (2.11). As indicated above, this equation can be solved given the values of ζ(v 0 ) and eitherζ(v 0 ) or Already this routine is enough to evolve the fields Σ and ϕ, as well as ζ in time. More precisely, for any field F , knowledge of d + F permits one to writė In practice, it is not necessary to evolve Σ explicitly, as its profile on the next time slice will be constructed when step 1 repeats. The value of f 2 (v 0 + ∆v) can be computed by extractingḟ 2 (v 0 ) from the boundary behavior ofφ(v 0 , z) and integrating, or it can be extracted directly from the boundary behavior of ϕ(v 0 + ∆v). In practice we have found the latter numerically favorable. At this juncture, all that is needed to update the routine on the next time step iṡ a 4 (v 0 ). There are options for computing this. One method is to solve equation (2.12), rearrange the d + derivative to obtain ∂ v d + Σ at v 0 , and then study the UV behavior to extractȧ 4 (v 0 ). A better way is to use (2.25) and directly evolve a 4 (v 0 ) to the next time slice. In this way one arrives at time v = v 0 + ∆v with updated values of f 2 , a 4 , ζ, and ϕ(v 0 + ∆v, z), which is all the information required to begin the integration routine anew.

A practical method
The radial integration of the equations of motion can in principle be performed by any of a number of well worn techniques, including pseudo spectral and finite difference methods. In practice, we found that the logarithmic terms in the near boundary behaviors of the fields rendered pseudo spectral methods unstable unless sufficient care was exercised to explicitly remove these terms (an observation also made by the authors of [15]). While one can typically ameliorate this issue through suitable field redefinitions, we have found it more convenient to employ a finite difference discretization method and integrate from the boundary to the IR. This method is similar in spirit to the one that appears in [35].
The response of the boundary gauge theory to these time dependent perturbations is encoded in the time evolution of f 2 and a 4 , which show up in higher orders of the boundary series solutions (2.18-2.20). Therefore, to determine them accurately, we work with redefined fields where f 2 and a 4 appear in the leading order at the boundary expansions by subtracting source contributions in lower orders. To surmount some of the challenges posed by the logarithmic behavior near the boundary, we also subtract the first few logarithmic terms. Moreover, we rescale by appropriate factors of z so that the redefined fields do not vanish at the boundary. We thus work with the subtracted fields defined by z 2Ā = A − 3 n=0 a n z n−2 − 6 n=4 α n1 z n−2 log z − α 62 z 4 log 2 z, (5.9) (Dσ) n1 z n−2 log z − (Dσ) 62 z 4 log 2 z, (5.12) In our choice, the redefined fields are C 2 over the radial domain. Substituting these expressions into the equations of motion, we obtain evolution equations of the redefined fields, which are solved with the strategy described in section 5.1. The equations of motion are discretized by finite difference grids in the z and v coordinates whose sizes are ∆z and ∆v, and we arrange our numerical computation so that second order accuracy is supposed to be achieved. We use a fourth-order Runge-Kutta method in z-integration, and integrate from the boundary toward the interior. The z-derivatives of the fields which are not evolved in each differential equation are computed by a central finite difference scheme. We adopt a gauge where ζ(v) = 0 is fixed for all times. This generically implies that the location of the black hole horizon can vary throughout the evolution. We monitor the location of the apparent horizon on each time step and terminate the integration slightly inside it, which is also inside the event horizon. To updateφ to the next time step, we use an upwind difference scheme for the advection term in (2.13). The grid sizes are chosen such that the Courant-Friedrichs-Lewy condition is satisfied, and we use ∆v ≤ ∆z. For v-evolution, we use a modified Euler's method, where to update a 4 we integrate the constraint (2.25) with a fourth order Adams-Bashforth formula, 3 ∆v. (5.14) Once the final static configuration is reached, when the apparent and event horizons coincide, we solve (C.8) backward in time to compute the evolution of the event horizon z EH (v). We can then calculate the area of the event horizon on constant v-slices as The monotonic increase of the event horizon is utilized as a consistency check of our numerical computation. The area of the apparent horizon is computed analogously, as

Thermalization Examples
In this section we provide examples of the thermalization processes produced by our numerical computation. As the scalar source of the static states is nonzero, quenches are not symmetric between positive and negative δf 0 , and we consider quenches with δf 0 > 0 for simplicity. Several tests suggest that broad stroke, qualitative features of the resulting time evolution seem not to be significantly altered by the choice of this sign. 4 Thought this section and the next we will work in units off 0 and henceforth set κ 2 = 1. The stability and accuracy of initial data is discussed in Appendix B. An important question in the interpretation of these results is "how far" our initial states are along a given black hole branch. A practical metric for quantifying this distance is the area of the black hole horizon. The zero temperature solution, which can be thought of as the λ H → ∞ limit of the small black hole branch, has vanishing horizon area, and the horizon area grows monotonically as the scalar vanishes at λ H → 1. This further justifies our "large" and "small" naming conventions for the black hole branches. In the studies which follow, the smallest black hole we perturb has λ H /λ c = 3.23. In figure 7, the area of the black hole horizons as a function of λ is given, and this smallest black hole is indicated by a red dot. Evidently, the horizon area of λ H /λ c = 3.23 in this case is about 1/5 of that at the first order phase transition.
Generically, a time dependent perturbation of a black hole will take a static initial state through a non-linear regime controlled by the details of the quench profile, followed by a linear regime governed by the "ring-down" to the final steady state configuration (there may also be late time power law tails in some situations, but we will not be concerned with these here). The ring-down is fully determined by the quasi-normal modes of the final state black hole, and is dominated by the mode closest to the real axis, ω 1 . In turn, the quasi-normal modes characterize the linear response of the final state to small perturbations in any of several available channels. In the present case, where the perturbations preserve the homogeneity of the spatial R 3 , the gauge invariant perturbations organize themselves into representations of SO(3) [36,37] transforming as the transverse-traceless (spin-2), vector, or scalar.
In figure 8 the late time ring-down of one particular time dependent perturbation is shown on a logarithmic scale. The figure clearly indicates the presence of an excited mode of the form for some real constants Z 1 , ω * , and Γ > 0. Although figure 8 illustrates the general late time features of any perturbation in our study, it is important to note that the Figure 9. The temperature dependence of the decay width Γ for the lowest lying scalar quasinormal mode in several states of our theory. The blue circles are large black branes whose temperature is an integer multiple of T c . The orange squares correspond to the minimum temperature black brane (top) and the smallest black hole we perturb in our study (bottom). The ratio Γ/πT approaches 1.75953 (the dashed line) at high temperatures, which coincides with the expected value for perturbations of AdS 5 Schwarzschild by a dimension 3 scalar operator [38].
quasi-normal mode spectrum, and hence the particular value of ω 1 , is a property of the final state achieved after the quench. In other words, one should keep in mind that ω 1 varies as one traverses the static state solution space so that ω 1 = ω 1 (λ H ). In figure 9 we quantify this by plotting Γ as a function of temperature for several thermal states of our holographic theory. As expected, linear scaling of Γ with temperature appears in the conformal limit, and deviations from this behavior are already evident at the first order phase transition where T = T c . The qualitative properties of this plot are anticipated by the thermodynamic features of our model shown in figure 4 and the right plot of figure 3 which display a similar approach to conformal behavior above T c . In practice, we find that our quenches are well described by a rapid transition from the regime driven by the quench to the ring down characterized by linear response. We will thus be particularly interested in the decay width Γ, as it provides a natural time scale for our equilibration processes. The gaussian perturbations that we define in 5.1 are controlled by two dimensionless parameters,δ ≡ δf 0 /f 0 andτ ≡ τf 0 . Thus, we can define a thermalization time which immediately suggests two natural questions we would like to address: 1. For a given initial state, what is the dependence of the thermalization time onδ andτ ? Figure 10. Large amplitude quench. The blue and purple lines correspond toδ = 0.5 and 1, respectively. In the case of the larger amplitude quench, it is interesting to note that the energy density appears to be driven below the ground state energy density (i.e. negative) in the first moments of the quench.

2.
How does the answer to the previous question depend on the choice of initial state?
Our answers to these questions motivate the remainder of this section.

Thermalization Between Branches
The form of the ward identity (4.6) encourages the intuitive expectation that a large, fast quench will result in the most significant increase in the system's energy. In figure 10, examples of the evolution of the boundary operators in such a scenario are shown. These computations are performed with a widthτ = 0.168 and amplitudesδ = 0.5 and 1. These perturbations basically dominate the dynamics, and can generically be tuned to result in a final configuration characterized by a large black hole independent of the particulars of the initial state. Indeed, we have performed analogous quenches in many different initial states, all of which show qualitatively similar time evolution during and after the quench. In line with our previous discussion, the plots of Ô and T ii show late time behavior that is well fit by an exponential decay related to the ring-down of the black hole. As expected, we find that the lowest lying quasinormal mode for these large black hole final states retreats from the real axis as the final state energy density increases. Accordingly, we observe that these large amplitude perturbations imply rapid thermalization in the dual field theory.
Once the entire evolution between static initial and final states is known, the time evolution of the apparent and event horizons can be computed. This provides a useful view of the equilibration process, as well as a consistency check of the time evolution. The area of the event horizon monotonically increases as shown in figure 11, and its location should always lie outside (closer to the boundary than) that of the apparent horizon. We verify explicitly that both of these statements hold for all perturbations we study.

Thermalization Along a Branch
When the quench amplitude is very small compared to its width, the Ward identity (4.6) implies that the system's energy will be very nearly conserved, and thus we anticipate final configurations that are very close to the initial state. In particular, an initial state characterized by a small black hole can remain within the small black hole branch.
An example of this kind of quench is shown in figure 12. There the perturbation width is taken to be the same as in the example of previous section,τ = 0.168, but now the amplitude is chosen to beδ = 0.02.
With these parameters, the area of the event horizon in the final state is only about 2.5% larger than that of the initial state. Pictorially, on the scale of figure 7 this quench results in the red dot which marks the location of the initial state moving imperceptibly to the left after the perturbation. At late times, the ring-down of O and T xx is far more pronounced in this case, with easily identifiable oscillations shown in figure 12.   similar behavior. In line with expectations, unlike the large amplitude quenches this class of perturbation is sensitive to the particulars of the initial state. For example, for the quench shown in figure 12, we find that T therm ∼ 2 in units off 0 -but the same perturbation in an initial state with about 13.5 times the energy density has a thermalization time T therm ∼ 0.3, more than 6.6 times faster.
Of course there are a variety of perturbations one can perform which result in a final state thermodynamically similar to the initial state. As mentioned above, these all share the feature that the perturbation amplitude is relatively small compared to its width. We will see in detail below how to quantify relatively small. At present, we turn our attention to a similar example of thermalization along a branch. This is the large amplitude, slow quench shown in figure 13. In this case,δ = 0.5 andτ = 1.19, so that the ratio of the amplitude to width here is about 3.5 times the previous example. Accordingly, we anticipate a greater change in energy density and that this perturbation will thermalize more quickly. From the late time behavior we find that this is indeed the case, with T therm ∼ 0.86.
Our code is able to evolve even slower perturbations, including those whose widths are many times larger than their amplitude. Predictably, in this case the various one point functions respond by being gradually deformed from (and then nearly returned to) their initial vales. For these adiabatic perturbations, the "bumpy" features in the late time behavior shown in figure 13 disappear, and the exponential decay at late times can not be seen. This can be viewed as a consequence of the fact that for such slow perturbations, the quench width is larger than the time scale defined by the final state's lowest lying quasi-normal mode.

Parameter Dependence and Scaling Regimes
To better understand the features of different quenches, we look in more detail at the dependence of the final state on the perturbation parameters. In figure 14, we show several results for the dependence of the final state energy density onδ for fixedτ and vice-versa. In figure 14(a), we compare perturbations with many different amplitudes but one of three fixed quench widths:τ = 0.168, 0.433 and 0.838. In figure 14(b), we fix the perturbation amplitude atδ = 0.05, 0.5 or 1 and consider quenches with many different widths.
The fixed τ scenario is studied closely in figure 15. For sufficiently smallδ and τ , the change in T tt across the quench is well approximated by a quadratic in the amplitude of the perturbation. In fact, this behavior manifests over a wide range of quench amplitudes providing that the quench is very fast, so thatτ is very small (or τ f 0 ). As the quench width is increased, however, deviations from this simple scaling become more pronounced. In figure 15(b), in whichτ is no longer small, the quadratic amplitude dependence of the final state energy density is only present forδ 0.1. Outside of this small amplitude region we find that in this case T tt ∼ δf 0 .
In figure 16 we analyze the dependence of the final state energy density on the quench speed with fixedδ. Again, there is a pronounced scaling regime evident when the quench is fast. In this case we find that T tt asymptotically behaves as 1/τ 2 even when δf 0 is comparable in size tof 0 . The scaling behavior degrades asτ becomes large, as is obvious from the lower left corners of the plots in figure 16. This suggests that it is roughly the duration of the quench as compared to the characteristic scale of the theoryf 0 that determines whether we are in this scaling regime, as opposed to  Figure 16. Analysis of two fixed τ cases. In both cases, the small τ asymptotic behavior is given by T tt ∼ 1/τ 2 . In 16(b), there is a distinct intermediate region near τ 1 where T tt ∼ 1/τ 2 once again before transitioning to T tt ∼ 1/τ in the adiabatic limit (inset). the size of the width relative to the quench amplitude. In figure 16(b) we consider perturbations with fixed small amplitude. As shown in the inset, in this case there appears to be an intermediate region aroundτ 1 where the T tt ∼ 1/τ 2 scaling momentarily returns. However, as the quench width is further reduced we find that this scaling is not maintained, and the functional form appears to approach T tt ∼ 1/τ in the adiabatic limit. This adiabatic behavior at smallδ is consistent with what was found in the absence of a confining potential in [14].
Combing the lessons of this subsection's results, we find that T tt ∝ (δf 0 /τ ) 2 when the quench width is small compared to all other scales. This implies that a constant T tt line on the quench parameter plane is given by δf 0 ∼ τ . The scaling appears to be approximately respected even for small amplitude perturbations, which is consistent with the linear behavior found in the bottom left corner of the dynamical phase diagram in figure 17. Perhaps not surprisingly, this behavior is also consistent with a result first observed numerically under fairly different circumstances [14]. It was subsequently demonstrated [16][17][18] that such a scaling is universal in the sense that it only depends on the near boundary features of the bulk theory, which is asymptotically AdS in many applications relevant for holography. Roughly, the idea is that very abrupt processes in the UV do not probe very deeply into the bulk, as the quench concludes before the disturbance has time to propagate far from the boundary. Since all the interesting features of our model appear deep in the IR, they play no role in the dynamics of this fast quench regime and it is therefore natural for our model to exhibit the same universal scaling behavior.

Dynamical Phase Structure
The coarse features of the quench dynamics can be efficiently summarized in a dynamical phase diagram, which we present in figure 17. To construct this diagram, we begin in the smallest black hole with which we can reliably evolve arbitrary perturbations, and scan the (τ ,δ) parameter plane. For each of the many perturbations, we follow the evolution until a static final state solution is obtained, and then use figure 5 to decide if that solution describes a large or small black hole. The resulting phase diagram has several noteworthy features. First, the thick dashed line separating the perturbations which result in a small black hole from those that end up in a large black hole likely marks the location of a crossover, as we have not been able to find any non-analytic behavior across this transition. To look for the sort of non-analytic behavior we have in mind, one can consider the set of susceptibilities that describe the system's dynamical response to various perturbations. For example, derivatives of the form where ∆E is the change in the system's energy density and the subscript is an instruction to take derivatives along the direction where that parameter is held fixed. In our quenches, these susceptibilities appear to be continuous and smooth as one approaches the transition line. This can be inferred from the plots shown in figure 18 which shows an example typical of what we find when we differentiate the final state  energy density with respect to the quench parameters at fixedδ orτ . In general we expect that most dynamical observables in our model such as thermalization times and entropy production depend smoothly on ∆E ≡ T tt FINAL − T tt INIT , and thus we do not anticipate critical behavior to arise in other observables in the vicinity of this line. This is analogous to what happens in weak field perturbations of global AdS, in which the line separating small and large black hole final states was also consistent with a crossover [8].
Finally, the functional form of this dashed curve appears to be suitably described by a simple power law in the "small-fast" regime. Here the data can be well fit to a function of the formδ ∼τ , which is in line with general expectations as we discussed above. Independent of its precise functional form, we find that the crossover line monotonically increases inδ asτ increases for the parameter space covered by our study.

Discussion
The line of research initiated in the present work provides another step forward in the holographic study of gauge theories that are qualitatively similar to SU (3) Yang-Mills in 3+1 dimensions. Our gravitational model is characterized by a non-trivial (and carefully tuned) dilaton potential which introduces an additional level of complexity to the gravitational infall calculations dual to boundary theory thermalization. Specifically, the gravitational solutions of our model all involve a running dilaton in the bulk whose profile encodes the breaking of conformal invariance in the dual gauge theory.
The fact that these "initial state" static solutions have nontrivial dilaton profiles requires a careful treatment of their near-boundary behavior in order to achieve stable evolution. We have applied a numerical technique based on finite difference integration to successfully evolve a broad class of perturbations, solving the fully backreacted Einstein-dilaton equations for the duration of the perturbation/equilibration process.
By studying the dependence of the final state on the parameters of the quench, we have collected several interesting lessons which we expect to prove insensitive to the precise details of the dilaton potential. Among these are the observation of a rapid transition from the perturbation (initial condition) dominated regime to the quasinormal mode dominated regime of the quench dynamics. This behavior appears to be ubiquitous in holographic realizations of thermal quenches, and likely reflects the fact that for solutions with sufficiently large black hole horizons the thermalization time is governed by the Hawking temperature T therm ∼ 1/T (which is now the dominant length scale). The extent to which this continues to be true in the extremal limit of our model is largely an open question. From figure 9 it is clear that the decay width of the lowest lying quasi-normal mode has already departed from the simple conformal temperature scaling for final states with temperature near T c . However, as we illustrate in section 6.2, our perturbations which remain on the small black hole branch also manifest this quick transition to the linear regime. An interesting question is then where might one expect to see deviations from this "ubiquitous" behavior?
As we move farther along the small black hole branch, the length scale introduced by the presence of the dimensionful source becomes increasingly important. This suggests that a sensible guess for where deviations from the rapid transition to the linear regime appear is for those perturbations whose final state has an energy density which is small compared to the source-in other words T tt /f 0 4 1. The smallest black hole which we perturb in the present work has T tt /f 0 4 ∼ O(1), and thus it is perhaps not surprising that we see no evidence for a new non-linear regime sensitive to the presence of the confining potential. Pushing the reach of our numerics closer to the extremal solution is an ongoing direction of our research which we hope to report on in the future.
Ultimately, we would like to go beyond the small black hole limit we have studied in this work and consider perturbations of the horizon-less zero temperature solution directly. The most compelling motivation for this is to determine whether or not the diverging dilaton potential characteristic of a broad class of holographic models of QCD is sufficiently repulsive to give rise to scattering type solutions which never result in black brane formation. If this is indeed the case, our dynamical phase diagram would gain a new line dividing those perturbations which result in black brane formation from those which do not. Finding the associated scattering solutions would have interesting implications for the dual gauge theory, suggesting a class of perturbations that the strongly coupled matter can not thermalize.
Moreover, the boundary of these "unthermalized" perturbations in the dynamical phase diagram would be interesting in its own right. By analogy with more familiar examples in asymptotically flat space, one might hope to find critical behavior akin to Choptuik phenomena [39], in which a bulk scaling solution appears on the boundary of the perturbations that do and do not form a small black hole. Solutions of this sort are typically accompanied by various power law scalings characteristic of a second order phase transition. In our model, for example, this could manifest as the final state energy density assuming the form ∆E|τ ∼ (δ −δ c ) γ whereδ c is the amplitude of the critical perturbation and γ is a critical exponent quantifying the universality class of our model. The exponent may be different from that of Choptuik, as the small black holes in our theories (unlike AdS) depart importantly from flat space black holes.
The appearance of a universal scaling regime in the fast quench limit is another noteworthy output from our calculations. As previously discussed, this scaling regime was anticipated on very general grounds, and its manifestation in our abrupt quench data is in some ways an encouraging check on our numerics. It is interesting to wonder how this scaling might be effected by the potential barrier inherent to the zero temperature solution. In so far as the universality of this scaling depends only on the UV features of the bulk solution, it is likely that quenches of the extremal solution whose width is much smaller than the scalar source will again result in a final state whose energy density satisfies (1.1). However, because the arguments for this scaling behavior are closely tied to the early time dynamics of the quench, it remains unclear whether the increase in energy density will manifest in black brane formation or a non-thermal scattering solution.
Investigating the properties of other probes which can be used to characterize our perturbations is another interesting future direction. In [15,40] a variety of non-local probes were used to measure the approach to thermal equilibrium. These include various two point correlation functions, Wilson loops, and the system's entanglement entropy. The advantage of these non-local probes is that they are capable of moving beyond the binary thermal/non-thermal characterization of the perturbation, as they are sensitive to the scale dependence of the thermalization process. In the gravitational picture, they accomplish this by sampling the geometry away from the UV boundary.
To understand why this is true, it is instructive to consider the equal time two-point correlator for some gauge invariant boundary operator with large conformal dimension. In the semi-classical (optical) limit, this correlation function can be computed holographically by calculating the length of the bulk geodesic that connects two spatially separated points on the boundary. As the separation distance between these points increases, the bulk geodesic droops increasingly deeper into the IR. By measuring the deviation of the length of this bulk geodesic from the value one obtains in the final thermal static state as a function of boundary separation and time, one arrives at a picture of the approach to thermal equilibrium at different length scales. Constructing this picture in our model is currently in progress.
Finally an important question is the implications of our results for heavy-ion collisions. The fast thermalization of sufficient high-density initial states should be comparable to the processes studied here upon translating appropriately the quench parameters. Of course it should be kept in mind that the heavy ion collisions are anisotropic in space, but there are good reason to believe that this is not very important for initial thermalization but more important for the subsequent evolution of the plasma. In this respect our numerical results on thermalization time should provide reliable estimates for the analogous heavy-ion thermalization process. Most importantly, the part of the physics that is not descried here, namely the boundary between thermalization and non-thermalization will also provide important clues for thermalization in the less understood pp collisions, where recently CMS reported the first ever evidence for collective effects.
NOTE ADDED: As this work was being prepared for submission, several papers [41][42][43] focusing on quasi-normal modes in nonconformal theories/backgrounds simultaneously appeared on the arXiv. One of them [43] considered bottom-up Einstein-Dilaton theories which were not confining, and quasi-normal modes for a probe massless scalar were computed. These works have an overlap with our results on thermalization times in a non-conformal field theory, and appear to broadly agree with our conclusions.

A Holographic Renormalization and Boundary Theory Correlators
The near boundary analysis of section 2.2.1 is convenient for developing an algorithm for solving the Einstein equations, but the in-going null coordinates are ill suited to constructing the generating functional of the dual gauge theory. Following [32], the analysis of the on shell action near the boundary is most directly performed in Fefferman-Graham coordinates, in which the radial direction is orthogonal to the boundary directions.
In Fefferman-Graham coordinates, the metric takes the form and the metric and scalar permit the following expansions: Substituting these expansions into the Einstein equations and solving them order by order in ρ allows one to determine many of the coefficients in the expansions algebraically. The primary exceptions are the leading coefficients of the normalizable modes, g (4) and ϕ (2) which can only be determined given the full radial profile in the bulk. Nevertheless, the near boundary analysis does constrain these undetermined coefficients, a fact realized in the boundary gauge theory by the existence of Ward identities. The regularized on shell action (2.1) can be written where Einstein's equations have been used to eliminate the Ricci scalar. This action exhibits the following divergences in the limit → 0: To construct the appropriate counter terms, the regulated action must be expressed in terms of the fields living on the surface ρ = . These fields are the pullback of the metric, γ(t, ), and the scalar ϕ(t, ). Performing this inversion yields where A contains O( 0 ) finite counter terms, and the coefficient F 4 depends on the details of the higher order terms in the scalar potential. For the potential in (2.6), it reads In the absence of an organizing principle such as supersymmetry, the finite counter terms are left unfixed and thus lead to scheme dependent ambiguities in the correlation functions.
Once the counterterms have been identified, one can form the subtracted action like and the renormalized correlation functions are then computed as follows: where T ij [γ] is the stress tensor of the theory at ρ = . This boundary stress tensor is generically the sum of the contribution from the regularized action, and the contribution due to the presence of the counterterms. The regularized action gives with K ij the extrinsic curvature of the regulating surface and K = γ ij K ij its trace. To obtain the contribution from the counterterms, it is convenient to first catalogue the metric variations of a boundary action of the form where A, B, C, D and E are arbitrary scalar functionals independent of the metric. Obviously this action contains (A.5) as a special case. With a bit of effort, one can show that the metric variations of S B yield All contractions and curvatures in T B ij implicitly refer to the metric γ, and derivatives differentiate everything to their right. From the terms in (A.12) with coefficients determined by comparison to (A.5), it is straightforward to obtain the contribution to the boundary stress tensor resulting from the counterterms.
Performing this maneuver, summing the result with (A.10) and then inserting the on-shell near boundary expansions from (A.2) into (A.8, A.9), one obtains the renormalized one point functions: (A.14) (A.15) The schematic notation ∂A refers to the contributions coming from the set of finite counter terms contained in A. These contributions are scheme dependent, and will henceforth be neglected for simplicity. Note that in deriving these expressions the boundary metric has been assumed to be flat, g (0) = η ij . Written in terms of the near boundary expansion coefficients, it is straightforward to demonstrate that these one-point functions respect the anticipated Ward identities. For example (in units of κ 2 ), demonstrates the breaking of conformal symmetry in the presence of a dimensionful source in terms of the classical result (first term) and terms due to the matter anomaly in four dimensions. Similarly describes the change in the system's energy in terms of the work done on it by a time dependent source. The derivation of these identities requires additional constraints which are easily obtained from the equations of motion. They relate, for example, time derivatives of the undetermined near boundary coefficients. Because the numerical computations directly access the boundary expansion coefficients given in (2.20), it is convenient to relate these coefficients to those appearing in (A.2). The coordinate change is given by which implies two particularly useful equations: These equations can be solved perturbatively in ρ to obtain z(t, ρ) and v(t, ρ). One straightforward method to this end is to expand the ingoing null coordinates in powers of ρ, like z(t, ρ) = ρ + n=2 ρ n s n (t) +s n (t) log ρ , and upon inserting these expansions into (2.18) and (2.20) and regrouping terms, one directly obtains the independent coefficients ϕ (0) , ϕ (2) and g (4)tt in terms of f 0 , f 2 , and a 4 . The result is

B Initial Data and Convergence
The initial data we wish to perturb and evolve in time are the solutions to the static equations of motion given by setting v-derivatives to zero in (2.8-2.12). These solutions are constructed by integrating the static equations from the horizon to the boundary and matching to the boundary behavior of solutions in the gauge ζ(v) = 0. This computation can be performed at very high precisions without much effort. The results are then exported on the discretized grids desired for time evolution. Since the numerical codes used to produce and evolve the static solutions are distinct, it is an important check on our numerical package that the unperturbed initial state can be stably evolved. This procedure also helps us determine suitable grid sizes for satisfactorily suppressing numerical errors. When the energy density of the initial state background is large, the dilaton is small throughout the bulk and it is fairly easy to achieve robust evolution. These static solutions can be evolved with only moderately small grid sizes, and numerical errors are very small. The scaling of the discretization error is second order in ∆z, as desired.
As the background scalar becomes larger, however, we find that even simply maintaining the static initial state becomes a challenge. In particular, the difficulty drastically increases as the black hole size decreases as it heads towards the small black hole branch. In figure 19, time evolution of a static solution with λ H /λ c = 3.23 is shown. This corresponds to almost the smallest black brane in which we can reliably perform many distinct perturbations. We divide the domain between the boundary and the initial black hole horizon z H | stat into N intervals, and hence the grid size is ∆z = z H | stat /N . In the left panel is the time evolution of Ô . When N = 200 and 400, there are some numerical oscillations which persist on the order of ∆z, which is  Figure 19. Convergence of static time evolution by changing the grid size for an initial data set with λ H /λ c = 3.23. In 19(b), the magnitude of the difference of T tt (v) from the input value is shown, δ T tt ≡ | T tt (v) − T tt (0) |/| T tt (0) |. The blue, purple and green lines correspond to N = 200, 400, 800, respectively. In the late time, the green lines do not oscillate and converge to the static value exponentially.
larger than the desired numerical error. When the grid size is decreased by one half, although oscillations remain, the decrease of their magnitude is demonstrably second order in ∆z. This implies that the presence of these oscillations is not indicating any real numerical instability, but is instead noise. We observe that this noise is largely insensitive to the size of ∆v. For the initial data shown in figure 19, we find that the noise is heavily suppressed and the late time behaviors of both panels converge to constant values when N = 800. These results suggest that for large background scalars, taking a very small grid size is necessary in order to rid the computation of unwanted noise. We have also made several tests with larger λ H /λ c , and in those cases the required smallness of the grid size for suppressing the noise continues to quickly increase. The static evolution of T xx behaves similarly to Ô .
These oscillations are far less pronounced in T tt , as indicated by (4.6). In the right panel of figure 19, we plot the magnitude of the difference of T tt for the unperturbed solution evolved in time from the input value. Notice the difference of scale compared to the left panel. In this logarithmic plot, quadratic convergence against the grid size of the finite discretization is evident. We also checked the bulk constraint equation (2.12), and verified that this too converges quadratically.
One consequence of evolving our unperturbed initial state solutions in time is that any numerical irregularities in the solution are damped by the presence of the horizon as the numerical solution "rings down" to a numerically stable configuration. This procedure thus cleans the numerical data for static initial states from the discrepancy between using different numerical codes for static and dynamical computations. In practice, before we perform any sort of quench, we typically allow the numerical data for the initial state to evolve unperturbed for a bit to clean the noise, and then compute the evolution of the perturbed state with suitably small (quench dependent) grid sizes.
C Geometrical Aside

C.1 The Apparent Horizon
The formulation of our gravitational problem relies crucially on the notion of an apparent horizon. This horizon is primarily important because unlike its more familiar sibling the event horizon, it is not teleological. This is to say that the location of the apparent horizon can be determined on each time slice, whereas the event horizon can only be deduced once the final state of the geometry is known.
An operational definition of an apparent horizon is a spacelike surface on which an outgoing null congruence normal to the surface has zero expansion. Following [44] one may study this expansion θ via a null vector tangent to an outgoing null geodesic, k. In the basis provided by the coordinates of (2.7), such a vector is given by This vector is not affinely parametrized, which means that it satisfies the geodesic equation k µ ∇ µ k ν = κ k ν (C.2) for non-zero κ. It is easy to show that in the present case, where the prime denotes differentiation with respect to z. Because κ is non-zero, the expansion equation is modified as follows: In this expression, the exponential pre-factor is necessary to convert to an affine parametrization. As it is manifestly positive and non-zero, it will play no role in the present discussion. From (2.7), (C.1) and (C.3) a short calculation reveals that in this background ansatz 5) and combining this result with the definition of the outwards directed derivative d + from (2.13), and the requirement that the expansion vanish at the location of the apparent horizon z H leads to this section's main result: Apparent horizons have several interesting properties that make them particularly well suited to the studying of dynamical processes in gravity. Specifically, for static spacetimes with black holes in the bulk, the apparent, event, and Killing horizons all coincide. Moreover, it is possible to show on general grounds that in the non-static case an apparent horizon will always lie within the event horizon. This suggests that an apparent horizon provides an IR cutoff for the numeric problem that is both natural and practical, as well as available on each time slice.

C.2 The Event Horizon
When the dynamics is such that the gravitational system settles back into a static state, as is the case in the holographic thermalization processes we study here, it is also straightforward to compute the location and area of the event horizon over time.
For this one needs only the observation (mentioned above) that the event and apparent horizons coincide in static spacetimes, and that the event horizon necessarily travels along a geodesic in the bulk. The line element traversed by a light ray in the ingoing null coordinates satisfies 0 = −A dv 2 − 2 z 2 dvdz, (C.7) which implies the following geodesic equation for the location of the event horizon, z EH : to be solved subject to the boundary condition z EH (v → ∞) = z AH (v → ∞), where z AH is the location of the apparent horizon. In practice, equation (C.8) can be integrated backwards in time along the geodesic using the numerically determined metric function A at each step.