A Equations of motion suitable for numerics

# Holographic thermalization in a top-down confining model

Holographic thermalization in a top-down confining model

B. Craps, E. J. Lindgren, A. Taliotis

Theoretische Natuurkunde, Vrije Universiteit Brussel, and

International Solvay Institutes, Pleinlaan 2, B-1050 Brussels, Belgium

Physique Théorique et Mathématique, Université Libre de Bruxelles,

Campus Plaine C.P. 231, B-1050 Bruxelles, Belgium

Ben.Craps@vub.ac.be, Jonathan.Lindgren@vub.ac.be, Anastasios.Taliotis@vub.ac.be

ABSTRACT

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.

## 1 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 . 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- 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 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 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 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 non-trivial 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 well-described 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.

## 2 Holographic setup

The model we will consider is the Einstein-Hilbert action with a minimally coupled massless scalar field,

 S=12κ2∫dd+1x√−g(R−2Λ−12(∂ϕ2))−1κ2∫∂ddx√−γK, (1)

where the cosmological constant is related to the AdS radius by . 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

 ds2=L2z2⎛⎜ ⎜⎝−dt2+dz21−zdzd0+(1−zdzd0)dθ2+d→x2d−2⎞⎟ ⎟⎠. (2)

We will henceforth work in units with . The AdS boundary is located at and the confinement scale is set by (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 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 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

 J(t)=ϵe−t2δt2, (3)

which is indeed negligibly for . The total injected energy will scale like for small and small [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 , where is the boundary metric, which breaks the isotropy of the boundary metric between the and 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.

### 2.1 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 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 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:

 ds2=L2z2⎛⎜ ⎜⎝−h(z,t)2dt2+f(z,t)21−zdzd0dz2+(1−zdzd0)e(d−2)b(z,t)dθ2+e−b(z,t)d→x2d−2⎞⎟ ⎟⎠, (4)

with the initial conditions that , and before the injection. We will refer to the boundary at as the UV and the point as the IR. The coordinate is periodic with period 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 , but for now we will keep explicit. The coordinates in (4) are very badly behaved at , 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 , namely we can assume a boundary profile on the form

 b(0,t)=ϵe−t2δt2, (5)

and turn off the scalar. It turns out to be convenient to define the following variables

 Π=˙ϕfh,P=˙bfh,Φ=ϕ′,B=b′, (6)

where means derivative with respect to the coordinate.1 Introducing

 G(z)=1−zdzd0, (7)

the equations of motion following from the ansatz (4) are

 ˙f= z1−2d(d−2)(d−1)Gh2(Gz−2(d−1))′(PBz+2P)+z2−2dΦΠGh(Gz−2(d−1))′+d−22Ph, (8) h′h= 1z2(d−1)(Gz−2(d−1))′((d−1)(d−2)2(P2+GB2+4GBz)+GΦ2+Π2)+ +(d−2)B−f′f, (9) h′h= 2d(d−1)(f2−1)z2d(Gz−2(d−1))′+f′f, (10) ˙P= 1d−1((Ge(d−1)b)′he−(d−1)bfzd−1)′zd−1, (11) ˙Π= (hGΦfzd−1)′zd−1, (12) ˙B= (Phf)′, (13) ˙Φ= (Πhf)′. (14)

Note that in this gauge, only derivatives of appear in the equations of motion, so we do not need to integrate at every time step to obtain . 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 , and using (6), we obtain

 ˙fz=z0=(d−22Ph)z=z0=(d−22˙bf)z=z0, (15)

so that

 fz=z0=Ced−22b|r=0 (16)

for some constant . Since initially and , we have that . Thus we can state this result as

 (fe−d−22b)z=z0=1, (17)

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 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.

### 2.2 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 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 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

 f= 1+E2(d−1)zd+…, (18) h= 1−E2(d−1)zd+…, (19) b= bdzd+…, (20) ϕ= ϕdzd+…, (21)

where the coefficients of and of have been related by the equations of motion. We will see later that will be the total injected energy, while the coefficient is related to the vacuum expectation value of the dual operator. We also have that , 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

 ds2=dζ2ζ2+1ζ2gαβdxαdxβ. (22)

Doing this asymptotically, we can identify , which gives us the metric

 ds2=dζ2ζ2−1ζ2{ (1−(E−1zd0)ζdd+O(ζd+1))dt2 +(1+(Ed−1+1−dzd0)ζdd+(d−2)bdζd+O(ζd+1))dθ2 +(1+(Ed−1+1zd0)ζdd−bdζd+O(ζd+1))d→x2}. (23)

Now it is easy to read off the non-zero stress energy components of the boundary field theory [24, 26]:

 ⟨Ttt⟩= 116πGN(E−1zd0), (24a) ⟨Tθθ⟩= 116πGN(Ed−1+1−dzd0+d(d−2)bd), (24b) ⟨Txx⟩= 116πGN(Ed−1+1zd0−dbd), (24c)

from which we see that is indeed the total injected energy (up to a factor of ), and is the initial AdS soliton energy density. Note also that .2 The vacuum expectation value of the scalar is

 ⟨O⟩=116πGNϕd. (25)

Note that taking the difference cancels the total injected energy and isolates the dynamical mode , 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 , and we expect that black branes will form. A black brane can be written as the metric

 ds2=1ξ2⎛⎜ ⎜ ⎜⎝−dt2(1−ξdξdh)+dξ21−ξdξdh+dθ2+d→x2d−2⎞⎟ ⎟ ⎟⎠. (26)

Note in particular that in the case of dynamically evolving from the AdS soliton background into the black brane (26), isotropy between the and coordinates must then be restored. From (24) this means that we must have 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 . Asymptotically, the radial coordinates and are related by , from which, comparing with (23), we can obtain the temperature of the black brane,

 T=d4πξh=d4π⎡⎢ ⎢⎣E−1zd0d−1⎤⎥ ⎥⎦1/d. (27)

## 3 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 , , and , corresponding to the AdS soliton geometry. The boundary conditions we impose in the UV are and as well as and ( and if we quench the metric instead of the scalar), and the source is always taken as a gaussian . 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.

### 3.1 The coordinate singularity

The 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 plane. However, the relation to the radial coordinate in this locally flat space is , and thus . This means that a small grid spacing in will be mapped to a very large grid spacing in (which is the natural coordinate around the point ), so a linearly spaced discretization in the coordinate will become incredibly bad at this point. We thus found it convenient to instead work with the coordinate , and use the metric ansatz

 ds2=1s(r)2(−h(r,t)2dt2+4f(r,t)2dg(r)dr2+r2g(r)e(d−2)b(r,t)d~θ2+e−b(r,t)d→x2d−2), (28)

where and . The advantage of this parametrization of the metric is that now is a finite slowly varying non-zero function and . The new periodic coordinate now has period . 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 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 , , , and are evolved in time, while the function is solved for at each time step using equation (A), 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 we do not have a boundary, and we can imagine continuing past to negative values and thus it is possible to still use central difference schemes when computing derivatives close to . An equivalent way of reaching the same result is to use the fact that all functions must be even functions of when computing the derivatives.

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 satisfying the free wave equation . Defining and we obtain

 ˙V=P′,˙P=V′, (29)

which should be compared to the equations of motion for the scalar field and the metric component . Now, if we discretize the coordinate by , and consider the derivative approximation , this will compute an approximation to the derivative at the point . We thus see that it might be convenient to put and on two different grids, one on , and one one where , 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 and 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 on one grid, and , and 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.

### 3.3 Extracting boundary data

To extract the boundary data, we will have to compute quantities like when . This becomes increasingly difficult when the dimension increases, since we are taking the ratio of two very small numbers. In particular for , 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.

## 4 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 . For small we have the relation , which is expected since the injected energy (which is the only parameter associated to the shell in the thin shell limit) goes like [14]. The shapes of the phase diagrams resemble those found for the hard wall model in [13]. In particular, for large we find numerically the relation , which is the same as in the hard wall model with Neumann boundary conditions.

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 . 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 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 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 .

## 5 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 black brane. This in particular means that isotropy between the coordinate and the coordinates will be restored, which in particular means from equation (24) that , and this is indeed verified numerically. Thus the pressure anisotropy will dynamically evolve from 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 . 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.

### 5.1 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 Fig. 5, a few examples of the deviation of the pressure difference from 0 are shown. We see that, as expected, the decay time is set by the lowest quasinormal mode, since the decay constants (for ), (for ) and (for ) are in good agreement with the values for the lowest quasinormal mode frequencies of AdS Schwarzschild black branes obtained in [29], namely , and , respectively.

## 6 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 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.

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 [30], where turbulent transfer of energy to short wavelengths was interpreted as an instability of AdS. In Fig. 9 we repeat our analysis for a smaller-amplitude 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.

### 6.1 Static solutions and non-thermalization

From equation (24) we see that we are not able to form a black brane if . 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 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 (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 and are trivially excluded. The only solutions for which 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 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

 ds2=1^zd(−^h2dt2+d^z2^f2+dθ2e(d−2)^b+d→x2d−2e−^b). (30)

This can be obtained by the coordinate transformations and field redefinitions given by

 b=^b−logGd−1,z=^zG12(d−1),^f2=f2G(1−zG′2(d−1)G)2,^h2=h2G−1d−1, (31)

where and . The coordinate now ranges from to . The (static) equations of motion for such an ansatz are

 ^h′^h=^f′^f−d^z(^f2−1), (32)
 ^h′^h=−^f′^f−d−24^z^B2−12(d−1)^z^Φ2, (33)
 (^B^h^f^zd−1)′=0, (34)
 (^Φ^h^f^zd−1)′=0, (35)

where , and prime now denotes derivative with respect to . We can integrate (34) and (35) to obtain

 ^Φ=Cϕ^f^h^zd−1, (36)
 ^B=Cb^f^h^zd−1. (37)

where and are integration constants, tuning the UV behavior. From (32) we have , so that we obtain the following formulas for and

 ^Φ=Cϕe∫^z0d^z′(^f2−1)^zd−1, (38)
 ^B=Cbe∫^z0d^z′(^f2−1)^zd−1. (39)

By eliminating from (32) and (33) and substituting the expressions in (39) and (38) for and we obtain that must satisfy

 2^f′^f−d^z(^f2−1)=−d−24C2b^z2d−1e∫^z02d^z′(^f2−1)d^z′−12(d−1)C2ϕ^z2d−1e∫^z02d^z′(^f2−1)d^z′. (40)

With the boundary expansion of being (see (18)), we obtain the boundary expansion of as . We expect that black branes will form when the total energy density is positive, which from (24) corresponds to . Here we will now consider the case (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 in the IR when or equivalently . 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 derivative in (40) is negative if . Since close to the boundary, we obtain that for all .

It is also easy to see that cannot become negative (because (40) implies that if then around that point for an so can only go to zero asymptotically). Also, can not asymptote to any other constant than zero. This can be seen by assuming that , and then (40) implies that for some which is inconsistent with (unless if , in which case for , which is also inconsistent, or if ). Since is strictly decreasing, it thus follows that we must have when . When we thus have that

 e∫^z0d^z′(^f2−1)d^z′≈C′^z−2d (41)

for some constant . For simplicity of notation we can thus redefine and . The asymptotic behaviour of is then

 ^b≈Cb,IRlog^z. (42)

Equation (40) now becomes in the IR

 2^f′^f≈−(d−24C2b,IR+12(d−2)C2ϕ,IR+d)^z−1, (43)

from which we can obtain the asymptotic behaviour

 ^f∼^z−d−28C2b,IR−14(d−2)C2ϕ,IR−d2. (44)

We must now compute the asymptotic relations between (, ) and (, ), by using the expressions in (31). We have that , which directly implies the asymptotic relations

 b≈^b+2log^zz0 (45)

and

 f≈d2(d−1)^f(^zz0)d−1, (46)

which imply

 fe−d−22b≈d2(d−1)^zz0^fe−d−22^b. (47)

From the above relations and (42) and (44) we now obtain the asymptotic behaviour

 b≈(Cb,IR+2)log^z, (48)
 f2∼^z(2(d−1)−d−24C2b,IR−12(d−2)C2ϕ,IR−d), (49)

so that we finally obtain

 fe−d−22b∼^z−d−22(Cb,IR2+1)2−14(d−2)C2ϕ,IR. (50)

We thus see that will generically vanish in the IR, trivially violating the condition that . The only way to have approach a constant in the IR, is when the power in (50) vanishes. This only happens when and , which in particular implies that the scalar identically vanishes. Only for these particular IR parameters will 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 and and then integrate to the IR, which should give a unique solution. Specializing to 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 , can be given by the metric (4) with and with replaced by . After transforming to the metric (30), by using the transformations in (31), we can obtain the asymptotic behavior for and as and . We can now easily obtain from (47) that . Thus, the only possible solution that can be obtained is which corresponds to our initial condition, and which corresponds in the UV to .

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.

### 6.2 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 , for reasons we will explain below.) The time scale can be seen to be independent of and as long as both parameters are sufficiently small. This is different from the modulation time scale observed in [13] for the 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 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 . Thus the metric is of order and the next order corrections to the scalar are . Working to order , 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 frequency3 . In the limit of small (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 the bounce frequency becomes independent of .

However, the metric also has some intrinsic frequencies (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 for some , we would expect a resonance, such that the amplitude of the ’th normal mode will increase linearly with time. But if , such that we are close to resonance, it would be natural to expect another time scale showing up, namely . 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 in its normal modes, where the decomposition is of the form . The functions (corresponding to the th normal mode with frequency ) satisfy the equation (75) with , which constitutes a Sturm-Liouville problem (which makes this decomposition possible), and they are normalized with respect to the inner product . 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 (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 , where is some localized profile and is the width. Thus when extracting the coefficient when computing the boundary observables, this will be proportional to , while the derivatives of the normal modes are of . 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 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 . 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.

#### 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 . (We denote the frequency of the kicks by .) The equation of motion is

 ¨x+ω2x=∑n≥0δ(t−nT), (51)

subject to the initial condition that should vanish for . To compare with our gravitational setup, is the analogue of (the metric backreaction), the delta functions are analogous to the stress energy tensor for the scalar which sources the metric, the frequency is analogous to the oscillation frequency of