Holographic thermalization in a top-down confining model

It is interesting to ask how a confinement scale affects the thermalization of strongly coupled gauge theories with gravity duals. We study this question for the AdS soliton model, which underlies top-down holographic models for Yang-Mills theory and QCD. Injecting energy via a homogeneous massless scalar source that is briefly turned on, our fully backreacted numerical analysis finds two regimes. Either a black brane forms, possibly after one or more bounces, after which the pressure components relax according to the lowest quasinormal mode. Or the scalar shell keeps scattering, in which case the pressure components oscillate and undergo modulation on time scales independent of the (small) shell amplitude. We show analytically that the scattering shell cannot relax to a homogeneous equilibrium state, and explain the modulation as due to a near-resonance between a normal mode frequency of the metric and the frequency with which the scalar shell oscillates.


Introduction
What happens to a strongly coupled field theory when it is brought far from equilibrium? This question is important in many areas of physics, including the formation of a quark gluon plasma in ultrarelativistic heavy ion collisions and quantum quenches in cold atom systems. It is a difficult question, however, because conventional techniques fail in the strongly coupled, far-from-equilibrium regime. In recent years, progress has been made using the gauge/gravity duality, also known as "holography". The simplest AdS/CFT models describe conformal field theories, and the original studies of holographic thermalization focused on those. Interestingly, when extrapolating to heavy ion collisions, one typically finds thermalization times that are short enough to be compatible with experiment [1,2,3,4,5,6,7,8,9,10]. Another noteworthy result is that short-wavelength modes thermalize first in the simplest holographic models [11,7]. An obvious question is whether there are interesting new effects for non-conformal models, in particular for confining ones.
The study of holographic thermalization in confining models was initiated in [12,13], where the hard wall model was considered, first in a weak field approximation [12] and then using fully backreacted numerical simulations [13]. Following [14] (see also [3]), starting from the ground state, energy was injected by turning on a homogeneous scalar source of amplitude for a brief time interval δt. In the bulk, this leads to a planar shell falling towards the interior of AdS, which in the non-confining context of [14] always led to the formation of a black brane, corresponding to thermalization in field theory. In the hard wall model, however, two regimes were found [12,13]. Certain shells collapsed into large black branes, while others kept scattering between the hard wall and the AdS boundary. In the scattering phase, for certain boundary conditions at the hard wall, the oscillating scalar expectation values underwent interesting modulation on time scales scaling like the inverse amplitude squared, due to resonant transfer of energy [13] similar to that discovered in [15] for collapse in global AdS.
The hard wall model is simple, but is sometimes criticized for being crude and ad hoc (an interior region of AdS being artificially cut away by the hard wall), and for being rather different from large-N Yang-Mills theory in certain respects (see, for instance, Section 1.1 of [12] for a brief discussion). We therefore decided to re-examine these issues in the context of the top-down AdS soliton model [16], which underlies top-down holographic models for Yang-Mills theory [16] and QCD [17]. In the bulk spacetime corresponding to the ground state, the radial direction and a circle combine into a cigar-shaped geometry, causing the radial direction to cap off smoothly at a radius that sets the confinement scale. While the starting point of the construction in [16] was the AdS 7 soliton (which after compactification on two circles left the radial direction and 3+1 large field theory dimensions), we will consider AdS solitons in 4, 5, 6 and 7 dimensions.
Starting from the AdS soliton spacetime, we will again inject energy using a minimally coupled massless scalar field with amplitude of order . An important difference with the hard wall model is that now the metric itself contains dynamics (in the sense that it is not completely determined by constraints), because there is no isotropy between the circle and the other spatial field theory dimensions. If and when a black brane forms, isotropy is restored in the metric components. Given a bulk solution, holographic renormalizaton can be used to extract field theory quantities such as the expectation values of the energy and of the pressure components.
We perform a fully backreacted numerical analysis, and identify a regime in which the infalling shell collapses into a black brane, possibly after one or more bounces, as well as a regime in which the infalling shell keeps scattering between the tip of the cigar (which we will henceforth refer to as the IR) and the AdS boundary. In the former case, we find that the pressure components relax to their (isotropic) equilibrium values according to their lowest quasinormal mode. In particular, the difference in pressure components along the noncompact and compact spatial dimensions of the AdS boundary relaxes to zero as an oscillating exponential. In the scattering phase, when the injected energy is sufficiently small, we show analytically that the shell cannot relax to a homogeneous equilibrium state, and we find numerically that the pressure anisotropy oscillates and undergoes modulation on a time scale that is -independent in the limit of small amplitude and short injection time. This modulation time scale is very different from the 1/ 2 time scale found for the hard wall model, and indeed the physical mechanism is different as well: the oscillations are due to an almost resonance between the oscillation frequency of the scalar shell and the lowest normal mode frequency of the dynamical metric component. Just above the threshhold for black hole formation, we also find solutions that bounce a few times against the AdS boundary before collapsing into a black brane, similar to solutions found in global AdS [15].
Recently, several other papers have studied holographic thermalization in non-conformal models. Based on a quasinormal mode analysis of top-down non-conformal (but gapless) models, it was conjectured in [18] that, as soon as a horizon is formed in the bulk, deviations from conformality do not significantly affect thermalization times. Similarly, the numerical analysis in [19] found that the equilibration dynamics of N = 4 SYM theory does not change much when chemical potentials or magnetic fields are added. In [20], quasinormal modes were computed for bottom-up models mimicking the equation of state of QCD, and a nontrivial dependence on the equation of state was found. A confining bottom-up model for QCD was studied nonlinearly in [21], but only for initial conditions that already contain a small black hole; in this case, good agreement with a quasinormal mode analysis was found. In the present paper, we study top-down confining models at the nonlinear level, in regimes with and without horizons. If and when a horizon forms, the subsequent dynamics is welldescribed by a quasinormal mode analysis, and the confinement scale does not play much of a role. In the parameter regime in which no horizon forms, the dynamics is dramatically different.
The remainder of this paper is organized as follows. In Section 2 we describe our setup, including our metric ansatz, the equations of motion and how to extract field theory quantities from the bulk solutions. In Section 3 we briefly discuss the numerical methods we have used. Section 4 starts introducing the reader to our numerical results and to the two different phases. In Section 5 we discuss the black hole phase and in Section 6 the scattering phase. We conclude in Section 7. A number of technical details are contained in three appendices.

Holographic setup
The model we will consider is the Einstein-Hilbert action with a minimally coupled massless scalar field, where the cosmological constant Λ is related to the AdS radius L by Λ = −d(d − 1)/2L 2 . The boundary term is the Gibbons-Hawking-York term, which is necessary to render the variational principle well defined [22], but it will not play a role in this work. We will start with the AdS soliton as an initial condition, and then inject energy into the system by perturbing the scalar field. The AdS soliton background corresponds to a Euclidean black brane with an extra time direction. An explicit metric can be written as We will henceforth work in units with L = 1. The AdS boundary is located at z = 0 and the confinement scale is set by z 0 (which would correspond to the horizon if we Wick rotated back the above to a black brane). Note that the θ coordinate is compact in order to avoid a conical singularity, and this metric breaks rotational invariance between the x coordinates and the θ coordinate. This will have the implication that in order to solve for the time-dependence of the metric, we will need the second order dynamical Einstein equations. This should be contrasted with rotationally invariant metrics, for which the metric can be determined using first order constraint equations alone. The massless scalar field will be dual to a marginal operator in the dual field theory. To quench the system we will use a source J coupled to this operator and turn it on for a short period of time. While we imagine the source to vanish outside a finite time interval, in our numerical computations we choose for simplicity a Gaussian profile of the form which is indeed negligibly for |t| δt. The total injected energy will scale like E ∼ 2 /δt d for small and small δt [14]. In the dual gravitational description, this source term corresponds to the value of φ at the AdS boundary. After the source has been turned off, the system's energy will have increased and the gravitational bulk solution will have a nontrivial time dependence, governed by the action (1). The main question is how this time-dependent solution behaves, in particular if it collapses into a black brane solution or not. To avoid an extra scalar field, we will also consider quenching the metric, and compare this to the case of perturbing the scalar. We can inject energy by turning on a short time dependence for η θθ /η x j x j , where η is the boundary metric, which breaks the isotropy of the boundary metric between the θ and x coordinates (see Section 2.1, in particular (5)). This can also be interpreted as quenching the size of the compactified dimension. In this case, only the gravitational mode will be turned on, and the dynamics will be qualitatively different from the case where both the scalar and the metric mode are turned on. Although we will just very briefly consider such quenches in the metric, it is important to remember that it is possible in this setup to inject energy via the metric without a scalar field and without breaking additional rotational symmetries.

Ansatz and equations of motion
To solve the Einstein equations we need to choose specific coordinates. We can constrain the form of the metric by using the symmetries of the problem and by suitable gauge transformations (diffeomorphisms), but otherwise the metric will be completely general. In particular, our metric will only depend on time and on the radial bulk coordinate. Also, note in particular that due to parity invariance in the θ and x coordinates, we can set all off-diagonal terms involving these coordinates to zero. We may then use our gauge transformation to bring the metric to a diagonal form with three free functions. Note that the absence of rotational symmetry between the θ and x coordinates in the AdS soliton background forces us to choose a more general ansatz than in setups with rotational symmetry in all spatial coordinates [14,23,12,13,21], in which case there are usually only two free functions in the metric, which are completely determined by the constraint equations (there are no propagating degrees of freedom in the metric). We have found that the following ansatz is useful: with the initial conditions that f = 1, h = 1 and b = 0 before the injection. We will refer to the boundary at z = 0 as the UV and the point z = z 0 as the IR. The coordinate θ is periodic with period 4πz 0 /d to avoid a conical singularity. This form of the metric has a remaining gauge symmetry corresponding to rescaling of all coordinates. In the numerics, we will use this to set z 0 = 1, but for now we will keep z 0 explicit. The coordinates in (4) are very badly behaved at z = z 0 , however, so for numerics we will use a different ansatz, see Section 3 and Appendix A. To inject energy via the metric, we do this via the function b(z, t), namely we can assume a boundary profile on the form and turn off the scalar. It turns out to be convenient to define the following variables where means derivative with respect to the z coordinate. 1 Introducing 1 We warn the reader that will have a different meaning in Section 6.1 and Appendix A.
the equations of motion following from the ansatz (4) arė Note that in this gauge, only derivatives of b appear in the equations of motion, so we do not need to integrate at every time step to obtain b. This is the reason for the particular parametrization in (4), which makes the equations of motion decouple nicely. A similar ansatz is used in [3]. Evaluating equation (8) at the point z = z 0 , and using (6), we obtaiṅ so that for some constant C. Since initially f = 1 and b = 0, we have that C = 1. Thus we can state this result as f e − d−2 which will be crucial for the analysis in Section 6.1. This condition actually is the statement that the regularity (absence of a conical singularity) at z = z 0 is preserved in time. This can be easily seen in the ansatz (28), which is used in the numerical analysis, and we refer the reader to Section 3 for further discussion.

Boundary expansion and holographic renormalization
To compute field theory observables, one resorts to a process called "holographic renormalization" [24,25], which requires adding counterterms to the action to cancel divergences from the near-boundary region. These counterterms, which in odd dimensions give finite contributions to the various one-point functions, must be evaluated explicitly for every dimension and quickly become quite involved for increasing dimension. In addition, these contributions make the one-point functions scheme-dependent. However, when the source is turned off, the first non-trivial term in the boundary expansion is of order z d and no counterterms are needed. Since we will be interested in the evolution of the one-point functions after the source has been turned off, we will therefore be able to ignore the counterterms.
In even dimensions the counterterms do not give finite contributions to the one-point functions even when the source is nonzero and thus in this case the counterterms can always be ignored [24,26]. For d = 3 we provide the full asymptotic boundary expansion in Appendix C, which will be used in some of the figures. The asymptotic behaviour of the various fields after the source has been turned off is given by where the z d coefficients of f and of h have been related by the equations of motion. We will see later that E will be the total injected energy, while the coefficient φ d is related to the vacuum expectation value of the dual operator. We also have thatĖ = 0, which follows from the equations of motion or from holographic Ward identities.
To identify the stress energy components at the boundary, we want to write the metric in the Fefferman-Graham gauge Doing this asymptotically, we can identify , which gives us the metric Now it is easy to read off the non-zero stress energy components of the boundary field theory [24,26]: from which we see that E is indeed the total injected energy (up to a factor of 1/16πG N ), and − 1 is the initial AdS soliton energy density. Note also that T µ µ = 0. 2 The vacuum expectation value of the scalar is Note that taking the difference T θθ − T xx cancels the total injected energy E and isolates the dynamical mode b, which is why we will prefer to plot this quantity instead of the individual pressure components.

Temperature of black brane solutions
As seen in (24), the energy density will be positive for energies E > 1/z d 0 , and we expect that black branes will form. A black brane can be written as the metric Note in particular that in the case of dynamically evolving from the AdS soliton background into the black brane (26), isotropy between the x and θ coordinates must then be restored. From (24) this means that we must have b d = 1 and this is indeed verified numerically. The temperature of such a black brane, as obtained by the standard procedure of requiring the absence of a conical singularity for the Euclidean version of (26), is given by T = d/4πξ h . Asymptotically, the radial coordinates ξ and ζ are related by ξ = ζ −ζ d+1 /2dξ d h , from which, comparing with (23), we can obtain the temperature of the black brane,

Numerical methods
In this section we will list some important tricks that we had to employ to achieve stable numerical evolution. We used a fourth order finite difference method to discretize the radial direction, and then we used the ordinary differential equation solver scipy.integrate.ode from the Python library scipy [27] to evolve the resulting system of ordinary differential equations in time. We have as initial conditions f = 1, h = 1, b = 0 and φ = 0, corresponding to the AdS soliton geometry. The boundary conditions we impose in the UV are f (0, t) = 1 and h(0, t) = 1 as well as φ(0, t) = J(t) and b(0, t) = 0 (b(0, t) = J(t) and φ ≡ 0 if we quench the metric instead of the scalar), and the source is always taken as a gaussian J(t) = e −t 2 /δt 2 . In the IR, we do not have to impose any boundary conditions, since regularity already follows from the equations of motion. However, there are some potential sources for numerical instability and inaccuracy, the coordinate singularity in the IR and the AdS boundary being two examples.

The coordinate singularity
The z coordinate ansatz (4) is very inconvenient in the IR. The reason is that at this point the geometry looks locally like Minkowski space in cylinder coordinates, with rotational invariance in the (z, θ) plane. However, the relation to the radial coordinate in this locally flat space is z = z 0 (1−r 2 ), and thus dz = −2rz 0 dr. This means that a small grid spacing in z will be mapped to a very large grid spacing in r (which is the natural coordinate around the point z = z 0 ), so a linearly spaced discretization in the z coordinate will become incredibly bad at this point. We thus found it convenient to instead work with the coordinate r = 1 − z/z 0 , and use the metric ansatz where s(r) = z 0 (1 − r 2 ) and g(r) The advantage of this parametrization of the metric is that now g(r) is a finite slowly varying non-zero function and g(0) = 1/z 2 0 . The new periodic coordinateθ now has period 4πz 2 0 /d 3/2 . While the coordinate system (4) is convenient to derive analytic results and to extract the boundary field theory observables, the coordinate system (28) will be used for the numerical evolution. It is also clear now in these coordinates, that the regularity condition (absence of conical singularity), means that f e (d−2)b/2 remains constant in time, which is exactly the statement (17). The equations of motion for the ansatz (28) can be found in Appendix A. The functions f , P , Π, Φ and B are evolved in time, while the function h is solved for at each time step using equation (60), and equation (62) is checked for consistency during the time evolution.
There is also a convenient trick that can be employed to compute derivatives close to an origin of a polar coordinate grid. Usually if one were to employ finite differences close to a boundary, one would have to resort to non-symmetric stencils which can induce instabilities or numerical inaccuracies. However, at r = 0 we do not have a boundary, and we can imagine continuing r past r = 0 to negative values and thus it is possible to still use central difference schemes when computing derivatives close to r = 0. An equivalent way of reaching the same result is to use the fact that all functions must be even functions of r when computing the derivatives.

Radial discretization
We have found that high order finite difference discretization has worked well. However, to avoid high frequency spurious oscillations, we have found that it is convenient to put different functions on two different grids. To motivate this, consider first a function v(t, z) satisfying the free wave equationv = v . Defining V = v and P =v we obtaiṅ which should be compared to the equations of motion for the scalar field and the metric component b. Now, if we discretize the z coordinate by {z j } n j=0 , and consider the derivative approximation (z j+1 − z j )/∆z, this will compute an approximation to the derivative at the point (z j + z j+1 )/2. We thus see that it might be convenient to put V and P on two different grids, one on {z j } n j=0 , and one one {x j } n j=0 where x j = (z j +z j+1 )/2, to improve the accuracy of the derivative approximations. If one were to use a central difference scheme, we find that it typically induces high frequency noise. This high frequency noise is still present when using higher order central difference schemes, but disappears when putting V and P on different grids (also when we use a higher order finite difference scheme).
In our more complicated setup, the same reasoning holds for the free wave equation in AdS, and we have found it very useful to employ the same trick even when including backreaction. Thus we have put Φ and B on one grid, and Π, P and f on the other. Function values are then interpolated to the other grid when necessary. This proves to result in very stable evolution and the high frequency noise that is present when using central difference schemes with all functions on the same grid disappears.

Extracting boundary data
To extract the boundary data, we will have to compute quantities like (f (z) − f (0))/z d when z → 0. This becomes increasingly difficult when the dimension increases, since we are taking the ratio of two very small numbers. In particular for d = 6, there is high frequency noise which makes it difficult to extract the observables. For the simulations of black hole formation (Fig. 5), we therefore found it appropriate to use a Savitzky-Golay [28] filter to get rid of this noise and to make the boundary observables more smooth in time.

Phase diagram
When injecting energy into (the Poincaré patch of) vacuum AdS, we always form a black brane. However, since the energy density of the AdS soliton is negative, and any black brane has positive energy density, there should be a threshold for black hole formation when injecting energy into the AdS soliton. The obvious question is then, what solution do we obtain below the threshold? In the probe limit, the scalar field will just bounce forever between the boundary and the IR, so one could ask if this behaviour will still remain when turning on backreaction, or if the system will equilibrate into some static solution after a long time. In Section 6.1, we prove that the system cannot equilibrate into any static solution. We will thus refer to these solutions as the "scattering phase", and the solutions that thermalize into black holes as the "black hole phase". In Fig. 1, we show the separation between the two different phases, in terms of the parameters and δt. For small δt we have the relation ∼ δt d/2 , which is expected since the injected energy (which is the only parameter associated to the shell in the thin shell limit) goes like E ∼ 2 /δt d [14]. The shapes of the phase diagrams resemble those found for the hard wall model in [13]. In particular, for large δt we find numerically the relation ∼ δt, which is the same as in the hard wall model with Neumann boundary conditions.  Figure 1: The separation between the black hole phase and the scattering phase. For small δt, we see that ∼ δt d/2 , which is expected since the total injected energy goes like E ∼ 2 /δt d [14]. For large δt, we find the relation ∼ δt.
Another interesting question is if we can find scattering solutions above the energy threshold. Intuitively, right above the threshhold, a wave packet should bounce before collapsing into a black brane due to the finite width of the wave packet, and this is indeed what we find: Right above the threshhold when black brane formation is possible (the energy density is positive), there is a region where solutions reflect many times against the boundary without collapsing (although we are not able to say whether they eventually collapse, due to numerical difficulties in following the solutions for a long time). We have also found solutions that bounce a few times and then collapse into a black hole, similar to what was found in global AdS [15]. In Fig. 3 we plot the number of scatterings before collapse, as a function of amplitude for fixed δt = 0.24z 0 . We see that when decreasing the number of reflections against the boundary before collapse varies between 0 and 3, and then for smaller there is a large region where the solutions do not seem to collapse.
In Fig. 2, we show the vacuum expectation value of the scalar operator and min{f /h} as a function of time, for a solution that bounces twice before collapsing into a black hole. After two reflections (identified by the sharp peaks in the vacuum expectation value) we see that min{f /h} approaches zero, which indicates the formation of an apparent horizon. If the wave packet is very close to collapsing to a black hole while it scatters in the IR, the wave packet usually becomes very squeezed and comes out almost like a shock wave, resulting in the very sharp peaks in the expectation value O . Note that for smaller than 0.06304, there is a parameter region where the injected energy density is above the black brane threshold, but nevertheless the solutions seem to scatter for as long as we have been able to follow them. This region is relatively large, since the threshhold where the energy density becomes negative is ≈ 0.0607.

Black hole phase
In the black hole phase, the space-time will collapse into a black brane, and a horizon will form. The resulting solution will be an AdS d+1 black brane. This in particular means that isotropy between the θ coordinate and the x coordinates will be restored, which in particular means from equation (24) , and this is indeed verified numerically. Thus the pressure anisotropy T θθ − T xx will dynamically evolve from −d/16z d 0 πG N to 0. A relevant question is what this isotropization process looks like. In Fig. 4, we show a typical evolution of the pressure anisotropy for d = 3. We see that the system quite rapidly enters a regime where it is isotropic up to some small fluctuations. In Section 5.1, we will compare these small fluctuations with the quasinormal modes of the resulting black brane. Our numerics will not allow us to follow the evolution for very long times after a black hole has formed, but long enough to see the quasinormal mode behaviour.

Quasinormal modes and late time behaviour
For late times we can view the solution as being composed of a black brane background with small fluctuations. The late-time relaxation is thus expected to be governed by the lowest lying quasinormal modes for this black brane. A standard way to illustrate this behaviour is to plot the logarithm of the absolute value of the deviation of some observable from its final value. In

Scattering phase
In the scattering phase, the scalar field wave packet that falls from the boundary, will bounce in the deep IR and return to the boundary. When it reaches the boundary there will typically be some excitation of the boundary observables. The wave packet then reflects from the boundary and the scattering repeats. There will be a similar quasiperiodic behavior in the metric, since due to the broken rotational symmetry between the x and θ coordinates the metric has dynamical degrees of freedom of its own. For all figures we have varied the grid spacing to make sure that the results are not numerical artifacts.
In Fig. 6, we show a typical scattering solution. As we can see, every time the scalar field wave packet reaches the boundary, there is a bump in the expectation value, and this oscillation goes on forever as far as we know. We can also see that the dynamical degrees of freedom in the metric are excited, as expected, leading to non-trivial behavior in the boundary pressure components. One interesting feature is that the interpretation of the scattering solution as a localized wave packet persists for very long times, even for solutions where the non-linearities play a significant role. This is not obvious; one could have imagined that the wave packet would broaden and that at late times we would have seemingly random fluctuations, but instead we see that the wave packet remains approximately localized for long times. However, the shape of the wave packet can change with time due to the non-trivial dynamics of the full Einstein equations, as is reflected in Fig. 7 for a scattering solution close to the black hole threshold. ). There is a clear distortion of the wave packet which is due to the non-trivial dynamics in the full Einstein equations and can not be seen in the probe limit, although the wavepacket remains fairly localized. This example is for d = 3, with parameters = 0.01 and δt = 0.1z 0 , and time is in units of z 0 . Note that this is already quite far into the non-linear regime since black hole formation occurs around ≈ 0.016.
In Fig. 8 we show a typical long time scattering solution when the metric is quenched according to (5). We notice that the pressure anisotropy develops increasingly sharp features after long times, which suggests transfer of energy to high frequency modes. At first sight, this might seem reminiscent of what happens to small-amplitude spherical scalar perturbations in global AdS 3 [30], where turbulent transfer of energy to short wavelengths was interpreted as an instability of AdS 3 . In Fig. 9 we repeat our analysis for a smalleramplitude source. While a Fourier analysis of the early and late time behavior in Fig. 8 confirmed the transfer of energy to higher frequencies, a similar analysis for Fig. 9 showed no significant transfer to higher frequencies in the time range studied: in the latter case, the spectrum is dominated by normal mode frequencies, with roughly the same strength at early and late times. Decreasing the amplitude of the source further would simply rescale the vertical axis in Fig. 9, showing that for small amplitude the dynamics we see happens on timescales independent of the amplitude. The limited time range of our numerical simulation does not allow us to exclude transfer of energy on longer time scales (e.g., scaling as the inverse amplitude). However, based on the absence of resonances in the normal mode spectrum in our setup, we expect no such energy transfer for small amplitudes. If so, the energy transfer observed in Fig. 8 is the result of strong nonlinearity and quite different from that of [30].
One important question is whether or not these scattering solutions will go on forever, or whether the system will approach some static solution. In section 6.1 we will show that, if the injected energy density is below the black brane threshold, the system must keep scattering forever.  , when is decreased further. For these small amplitudes, we find no significant transfer to high frequency modes.

Static solutions and non-thermalization
From equation (24) we see that we are not able to form a black brane if E < 1/z d 0 . However, one could imagine that there are other static solutions that the system can end up in. We will in this section show that if E < 1/z d 0 there are no static solutions that can be obtained through time evolution. To summarize the argument, the key information we get from the dynamical equations is the relation (17). This condition is essentially the requirement that the spacetime should be regular at z = z 0 (such that a conical singularity can not be formed at this point during time evolution). We will then consider static solutions, by looking at the static equations of motion, and show that any possible static solution is incompatible with (17). Actually, most of the solutions have a completely different asymptotic behavior at z = z 0 and are trivially excluded. The only solutions for which f e − d−2 2 b goes to a constant, turns out to be the AdS soliton solutions. However, as we will see, all AdS solitons except our initial condition soliton will have f e − d−2 2 b approaching a different constant than 1, violating (17), so if the injected energy is non-zero, no static solutions can form. This reasoning is reminiscent of the argument for non-thermalization in the hard wall model, given in [13], where we also had a relation similar to (17).
To investigate this we will start by considering a different coordinate system, such that the metric takes the form This can be obtained by the coordinate transformations and field redefinitions given by where Theẑ coordinate now ranges from 0 to ∞. The (static) equations of motion for such an ansatz arê whereB =b ,Φ = φ and prime now denotes derivative with respect toẑ. We can integrate (34) and (35) to obtainΦ where C b and C φ are integration constants, tuning the UV behavior. From (32) we havê f /ĥ = e ẑ 0 d z (f 2 −1) , so that we obtain the following formulas forB andΦ By eliminatingĥ /ĥ from (32) and (33) and substituting the expressions in (39) and (38) forB andΦ we obtain thatf must satisfy With the boundary expansion of f being f = 1 + E 2(d−1) z d + . . . (see (18)), we obtain the boundary expansion off asf = 1 + . .. We expect that black branes will form when the total energy density is positive, which from (24) corresponds to E − 1/z d 0 > 0. Here we will now consider the case E < 1/z d 0 (negative energy density) and show that any possible static solutions with negative energy density cannot be obtained dynamically. We emphasize that some solutions of (40) might have singular behaviours and should be excluded as relevant solutions by other arguments, but we will not care about such issues here, and just directly show that any static solutions, which must satisfy (40), cannot be formed dynamically with the AdS soliton as initial condition. Recall the relation (17), which says that we must have f e − d−2 2 b = 1 in the IR whenẑ → ∞ or equivalently z rightarrowz 0 . The idea is now to show that all solutions obtained by solving (40) are inconsistent with this requirement. We will use the notation ≈ to mean that two quantities are equal asymptotically, while ∼ means that they are equal asymptotically up to an overall constant.
To show this we will have to compute the IR asymptotic behaviour of the solutions (40). We first note that the derivativef in (40) is negative if 0 <f < 1. Sincef = 1 + (E − 1/z d 0 )ẑ d /2(d − 1) + . . . < 1 close to the boundary, we obtain thatf < 1 for allẑ. It is also easy to see thatf cannot become negative (because (40) implies that iff = 0 then around that pointf ∼ẑ −α for an α > 0 sof can only go to zero asymptotically). Also, f can not asymptote to any other constant than zero. This can be seen by assuming that f → c > 0, and then (40) implies thatf ∼ e −αẑ β for some α, β > 0 which is inconsistent withf → c (unless if C φ = C b = 0, in which casef ∼ẑ −α for α > 0, which is also inconsistent, or iff ≡ 1). Sincef is strictly decreasing, it thus follows that we must havê f → 0 whenẑ → ∞. Whenẑ → ∞ we thus have that for some constant C . For simplicity of notation we can thus redefine C b C = C b,IR and C φ C = C φ,IR . The asymptotic behaviour ofb is then Equation (40) now becomes in the IR from which we can obtain the asymptotic behaviour We must now compute the asymptotic relations between (f , b) and (f ,b), by using the expressions in (31). We have thatẑ ≈ z 0 G −1/(2(d−1)) , which directly implies the asymptotic relations b ≈b + 2 logẑ z 0 and which imply From the above relations and (42) and (44) we now obtain the asymptotic behaviour so that we finally obtain We thus see that f e − d−2 2 b will generically vanish in the IR, trivially violating the condition that f e − d−2 2 b → 1. The only way to have f e − d−2 2 b approach a constant in the IR, is when the power in (50) vanishes. This only happens when C φ,IR = 0 and C b,IR = −2, which in particular implies that the scalar identically vanishes. Only for these particular IR parameters will f e − d−2 2 b go to a constant. We will now show, however, that it will go to a constant different from 1.
To specify a solution in the bulk, it would be customary to specify the UV behavior, meaning that we specify E and C b and then integrate to the IR, which should give a unique solution. Specializing to C b,IR = −2 should leave us a one parameter family of static solutions. Below we will construct this one parameter family of solutions, which turns out to be all the AdS solitons.
An AdS soliton solution with a general confinement scale z 1 , can be given by the metric (4) with b = 0 and f = h = 1 with z 0 replaced by z 1 . After transforming to the metric (30), by using the transformations in (31), we can obtain the asymptotic behavior forf andb asb ≈ −2 logẑ z 1 andf ≈ 2(d−1) d (z 1 /ẑ) d−1 . We can now easily obtain from (47) that f e − d−2 2 b → z 1 /z 0 . Thus, the only possible solution that can be obtained is z 1 = z 0 which corresponds to our initial condition, and which corresponds in the UV to E = 0.
To conclude, when the total energy density lies between that of the AdS soliton and zero (the threshold for black brane formation), no static solutions exists.

Long time amplitude modulation
For small-amplitude scattering solutions (small ), we observe an amplitude modulation in the pressure anisotropy on a long time scale, see Fig. 10. (The relevant timescale is actually hard to see for d = 6, for reasons we will explain below.) The time scale can be seen to be independent of and δt as long as both parameters are sufficiently small. This is different from the 1/ 2 modulation time scale observed in [13] for the d = 3 hard wall model with Neumann boundary conditions. As we will now explain, in the present case the modulation is due a near-resonance between a metric mode and the bouncing scalar shell. (In the d = 3 hard wall model with Neumann boundary conditions, no dynamical metric modes were excited, and the modulation was due to a resonant spectrum of scalar field normal mode frequencies.) In the small regime, the scalar field is O( ). Thus the metric is of order O( 2 ) and the next order corrections to the scalar are O( 3 ). Working to order O( 2 ), we can therefore consider φ as a probe scalar acting as an external source on the metric. Since the scalar φ bounces back and forth between the IR and the boundary, the source for the metric backreaction can be characterized by a frequency 3 f φ . In the limit of small δt (the thin shell limit), this frequency will be the same as for a lightlike particle (following lightlike geodesics) that bounces between the boundary and the interior. In particular, for small δt the bounce frequency becomes independent of δt.
However, the metric also has some intrinsic frequencies f i (the normal mode frequencies, see Appendix B). Every time the scalar crosses the space-time it kicks the metric. It is useful to decompose the metric fluctuation into its normal modes. If f φ = f j for some j, we would expect a resonance, such that the amplitude of the j'th normal mode will increase linearly with time. But if f φ ≈ f j , such that we are close to resonance, it would be natural to expect another time scale showing up, namely T = 1/|f φ −f j |. The results are summarized in table 1, and can be compared with the numerical results in Fig. 10 and Fig. 11. The latter figure shows the decomposition of b in its normal modes, where the decomposition is of the form b(z, t) = n≥0 a n (t)Q n (z). The functions Q n (corresponding to the nth normal mode with frequency ω n ) satisfy the equation (75) with ω = ω n , which constitutes a Sturm-Liouville problem (which makes this decomposition possible), and they are normalized with respect to the inner product z 0 0 Q n (z)Q m (z)z −(d−1) dz = δ mn . Note that replacing the frequency of the scalar wave packet by that of a light-like thin shell still gives decent result, but using the true frequency is required to get accurate results, especially for d = 5 (where the system is very close to resonance). To summarize, the modulation can be traced back to a near resonance between the lowest normal mode frequency of metric perturbations and the frequency of a bouncing scalar shell. As expected from this picture, this type of modulation does not show up when we quench the metric instead of the scalar field, as can be seen in Fig. 9.
One can see from the numerics, however, that the metric perturbation consists not only of a few normal modes: it has a slowly moving normal mode part and a rapidly moving wave packet part. The wave packet part is in general smaller than the normal mode part. However, close to the boundary, the wave packet part can still give large contributions to the boundary observables. The intuitive explanation is as follows. Close to the boundary a wave packet looks typically like ψ((z − t)/δt), where ψ is some localized profile and δt is the width. Thus when extracting the z d coefficient when computing the boundary observables, this will be proportional to ∂ d z ψ((z − t)/δt) ∼ 1/δt d , while the derivatives of the normal modes are of O(1). We thus see that for larger dimensions, the wave packet part is expected to become more important. These are exactly the sharp peaks one can see in Fig. 10, and indeed they become larger for larger dimansions. For d = 6 they completely dominate and this is the reason why we cannot see the modulation due to the first normal mode in the vacuum expectation value in d = 6. However, it can still be seen in the normal mode decomposition in Fig. 11, since here the contribution from the wave packet part is still small.
16πG N z 6 0 T θθ − T xx +6 Figure 10: The pressure anisotropy after a weak scalar perturbation with small and δt = 0.1z 0 . Time is measured in units of z 0 , and the vertical axis has been rescaled by 2 /δt d , which is the expected dependence of the total energy of the system for small , δt. For d = 3, 4, 5 we see that the amplitude undergoes an amplitude modulation on a much longer time scale which is in excellent agreement with the result in table 1. For d = 6, the modulation due to the first normal mode is hidden by the peaks from the bouncing wave packet part; it is clearly visible, however, in the normal mode decomposition in Fig. 11.  Figure 11: The metric function b decomposed in normal modes, after a weak scalar perturbation with small and δt = 0.1z 0 . Time is measured in units of z 0 . We see that, as expected, the lowest mode is more excited than higher modes, and undergoes an amplitude modulation which agrees with the result in Table 1.

Harmonic oscillator toy model
To develop a better understanding of the modulations we have just described, we now study a sourced harmonic oscillator which is conceptually similar to our gravitational setup (in the small-amplitude scattering phase) and which experiences a similar amplitude modulation phenomenon. Consider a harmonic oscillator with angular frequency ω, sourced by a sequence of local kicks (modelled by delta functions) with period T . (We denote the frequency of the kicks by f = 1/T .) The equation of motion is subject to the initial condition that x(t) should vanish for t < 0. To compare with our gravitational setup, x is the analogue of B (the metric backreaction), the delta functions are analogous to the stress energy tensor for the scalar φ which sources the metric, the frequency f = 1/T is analogous to the oscillation frequency f φ of φ, and ω is analogous to  Table 1: The lowest normal mode frequency f 0 , the oscillation frequency of the scalar wave packet f for δt = 0.1z 0 , the oscillation frequency of a lightlike thin shell f and the expected modulation times T = 1/|f 0 − f | and T = 1/|f 0 − f | using f respectively f . Note that the frequency of a thin shell is extremely close to the frequency of the bouncing scalar field. However, note also that in AdS 6 (d = 5) we are extremely close to resonance, and to get an accurate modulation time we must use the true frequency of the scalar wave packet (compare with Fig. 10 and Fig. 11). The estimated error comes from reading off the oscillation frequencies of the wave packet from the numerical simulations, while the errors of f 0 and f are negligible.
the normal mode frequencies of the metric perturbations. We can solve (51) by performing a Laplace transform. For the Laplace transformed field X we have It is now easy to do the inverse Laplace transform, to obtain where θ(t) is the Heaviside step function. By letting N = t/T (the largest integer less than or equal to t/T ), we can write this as sin(ω(t − nT )) = Im e iωt N n=0 e −iωnT = Im e iωt 1 − e −iω(N +1) under the assumption that T ω ≡ 0 (mod 2π). Extracting the imaginary part and using some trigonometric identities, we obtain . (55) If the system is almost at resonance, T ω ≈ 2π, the middle factor in (55) will give rise to fast oscillations, while the last factor gives rise to slow amplitude modulations. To see this, we write f − ω/2π = f (so is the difference between the source frequency and the oscillator frequency). The third factor in (55) now becomes sin ω N + 1 2 T = sin [π(1 − T )(N + 1)] = ± sin [π T (N + 1)] where in the last step we approximated N = t/T ≈ t/T , and we see that we indeed obtain an overall amplitude modulation with period 1/ . An example with ω = 1 and = 0.05/2π is shown in Fig. 12. Further, we note that it is the small denominator in (55) that causes a near-resonant normal mode to dominate the other normal modes. If T ω = 2πk, for some non-zero integer k, the summation of the geometric series in (54) yields instead which is a sine function with a (step-wise) increasing amplitude, as expected when we are at resonance. This can also be obtained as a limit T ω → 2πk in (55). x(t) Figure 12: Evolution of the sourced harmonic oscillator given by equation (51), close to resonance with parameters ω = 1 and T = 2π/1.05. We see that the result is a rapidly oscillating solution with a long time modulation with period 2π/0.05 ≈ 125.7.

Conclusions and outlook
A common feature of the hard wall model and the AdS soliton is the energy gap between the ground state and the lightest black hole. This feature underlies the dynamical phase structure uncovered in [12,13] (for the hard wall model) and in the present work (for the AdS soliton).
In the black hole phase, we have found that the late-time decay is governed by the lowest quasinormal mode, as could have been expected. It was nice, however, to find numerical solutions in this model that collapse after one or more bounces, which is reminiscent of spherical collapse in global AdS [15] -it would be interesting to study the similarities and differences in more detail. A good starting point could be the observation that the AdS 3 soliton is identical to global AdS 3 .
In the scattering phase, we have established that, if the injected energy density is below the black hole threshold, relaxation to a static solution that is translationally invariant in the spatial directions along the boundary is impossible. An interesting question for future work is whether there are instabilities towards the formation of inhomogeneities in the x or θ directions. 4 In addition to the mere existence of scattering solutions, we have found that the pressure components exhibit clear amplitude modulation, and we have explained this as due to a near-resonance between a scalar shell bouncing frequency and a metric normal mode. It would be interesting to see whether scattering solutions can also be established for more phenomenological holographic models, and if so, whether they exhibit similar features. radial function and a harmonically oscillating function of time. The normal mode spectrum is expected to be discrete in confined geometries, and we saw in Section 6.2 that these normal mode frequencies explain the amplitude modulation of the pressure anisotropy.
To find such solutions, we assume that b = O(µ), for some small parameter µ, and that the scalar field vanishes. We then solve the equations to linear order in µ. So letting f = 1 + µf 1 + . . ., h = 1 + µh 1 + . . ., B = 1 + µB 1 + . . . and P = 1 + µP 1 + . . . we obtaiṅ We also obtain the Ward identity When going to Fefferman-Graham gauge as in Section 2.2, the intermediate z 2 terms will not affect the z 3 terms. Moreover, since in even dimensional AdS spaces the counterterms do not affect the field theory observables, the boundary observables are given by the same formulas (24) and (25), even when the sources are turned on. The Ward identity then takes the form