Real-Time Dynamics of Plasma Balls from Holography

Hans Bantilan,1, 2, ∗ Pau Figueras,1, † and David Mateos3, 4, ‡ 1School of Mathematical Sciences, Queen Mary University of London, Mile End Road, London E1 4NS, United Kingdom 2Department of Applied Mathematics and Theoretical Physics (DAMTP), Centre for Mathematical Sciences, University of Cambridge, Wilberforce Road, Cambridge CB3 0WA, United Kingdom 3Departament de Fı́sica Quàntica i Astrofı́sica and Institut de Ciències del Cosmos (ICC), Universitat de Barcelona, Martı́ i Franquès 1, ES-08028, Barcelona, Spain 4Institució Catalana de Recerca i Estudis Avançats (ICREA), Passeig Lluı́s Companys 23, ES-08010, Barcelona, Spain

The possibility of modeling a confining gauge theory with a gravity dual was initiated by [32], and the AdS soliton geometry [33] provides a concrete realization of such a gravity dual. This geometry consists of AdS compactified on a circle that shrinks to zero size at a region called the infrared (IR) bottom, smoothly capping off the spacetime there. The dual gauge theory is also compactified on a circle whose inverse size sets the confinement scale Λ. The AdS soliton is the preferred homogeneous phase at temperatures below the deconfinement temperature T c = Λ. Above T c , the preferred homogeneous phase is a translationally invariant black brane dual to an infinitely extended deconfined plasma. Homogeneous gravitational collapse in the AdS soliton that leads to black brane configurations was studied in [34].
It was conjectured, in [35], that localized, finite-energy black hole solutions also exist in the AdS soliton background and that they should be dual to stable, mixed-phase configurations in the gauge theory known as plasma balls: finitesize droplets of deconfined plasma surrounded by the confining vacuum. Subsequent work, see, for example, [36][37][38][39][40][41][42][43], confirmed this picture. In particular, [42] numerically con-structed static solutions of finite-energy black holes localized at the IR bottom of the AdS soliton. It was found that small black holes of this type are well described by spherically symmetric Schwarzschild solutions that have T T c , whereas large black holes resemble "pancakes" with T T c and are well approximated by the AdS black brane solution at points away from the interface with the surrounding vacuum.
In a set of articles, we are initiating a new program to study the far-from-equilibrium physics of confining gauge theories with gravity duals. In these theories, any collision with sufficiently high but finite energy is expected to result in the formation of an excited plasma ball. Such collisions were analyzed in the approximation of linearized gravity in [44]. Although this captures the expected set of asymptotic states, the linear approximation is unable to describe the formation of a black hole horizon and its subsequent time evolution.
In this first Letter we present the first fully nonlinear evolution of finite-energy black holes in a confining geometry. Specifically, we consider gravity in five spacetime dimensions with a negative cosmological constant coupled to a massless scalar field and solve for spacetimes with AdS soliton asymptotics. This gravity model is firmly embedded in string theory, since it is dual to a confining gauge theory [32] in four dimensions obtained by compactifying N = 4, SU(N c ) super Yang-Mills (SYM) theory on a circle with supersymmetry-breaking boundary conditions. Throughout the Letter we relegate several technical details to the Supplemental Material.
Numerical Scheme.-The results presented in this Letter are obtained with a numerical code that solves the generalized harmonic form of the Einstein field equations. We obtain solutions g µν =ĝ µν + h µν , where the deformation h µν is not small and the background metricĝ µν is that of the AdS soli-arXiv:2001.05476v2 [hep-th] 12 May 2020 4 . We have set the AdS radius and the location of the AdS boundary to unity, so all coordinates are effectively dimensionless. However, we reinstate dimensions whenever we quote specific values below. In these coordinates the IR bottom is at ρ = 0 and the boundary directions −∞ < t, x 1 , x 2 < ∞ are noncompact. There are only three of these boundary directions, reflecting the fact that four-dimensional SYM compactified on a circle behaves effectively as a three-dimensional theory at distances longer than the size of the circle.
In deforming the full metric away from the AdS soliton, we take the general form of the metric to be In other words, we assume no symmetries except for the U(1) of the circle parametrized by θ which, furthermore, is assumed to be hypersurface orthogonal. The period ∆θ of this circle sets both the confinement scale Λ and the deconfinement temperature T c = 1/∆θ = Λ. Light rays along the ρ-direction take a boundary time t bounce ≈ 0.8Λ −1 to travel from the IR bottom to the AdS boundary and back. From the near-boundary falloff of the metric we extract the gauge theory stress tensor, which we multiply by a factor of 2π 2 /N 2 c so that the rescaled quantities are finite in the large-N c limit. We choose a renormalization scheme in which the AdS soliton has a vanishing stress tensor, which implies that a translationally invariant black brane has energy density and pressures We couple gravity to a massless real scalar field ϕ. In the gauge theory this means that we consider the coupled dynamics of the stress tensor and a scalar operator dual to ϕ. For concreteness in this Letter we will only display the expectation values of the stress tensor. Since we do not turn on a source for this field, the equation of state of the dual gauge theory remains the same as in a CFT, namely T µ µ = 0. We construct time-symmetric data on the initial-time slice by solving the Hamiltonian constraint subject to a freely chosen initial scalar field profile. We take the latter to be a superposition of Gaussian lumps of the form where and (ρ c , x 1c , x 2c ) is the center of the lump. Results.-We dynamically construct localized black holes via gravitational collapse of a superposition of lumps of the form (4). We have performed several simulations, including some that do not result in the formation of a black hole. For concreteness, in this Letter we focus on a simulation with A 0 = 1.0, δ = 0.2, and varying ρ c , x 1c , x 2c for a superposition of three Gaussian lumps, but we have verified that the results that we describe here do not depend on the number of lumps, or the particular values of A 0 , δ , ρ c , x 1c , x 2c . In order to introduce an initial anisotropy in the x 1 -x 2 plane, we choose the three lumps to be centered at ρ c = x 1c = 0 but at different locations along x 2 given by Λx 2c = {−0.1, 0, 0.1}. This results in the formation of a single black hole at the IR bottom whose horizon is initially more extended along x 2 than along x 1 , with a total mass M ≈ 0.226Λ. The solution is invariant under x 1 → −x 1 and/or x 2 → −x 2 , so figures will only show half of each axis. Scalar and gravitational waves propagate outward in all spatial directions as the black hole relaxes to equilibrium. In this context, the timelike AdS boundary and the IR bottom play the role of a waveguide: as outgoing waves propagate along the boundary directions, these waves bounce back and forth between the boundary at ρ = 1 and the IR bottom at ρ = 0 at intervals given by t bounce . Figure 1 illustrates these dynamics with several snapshots of the scalar field variable.
In the vicinity of the black hole in the x 1 -x 2 plane, there is a long-wavelength feature extending in the ρ direction that is left behind by waves without enough momentum to escape in the boundary directions, but whose wavelengths are too long to be efficiently absorbed by the black hole. The outgoing wavefronts are evident as the snapshots progress in time, propagating in the boundary directions. The initial anisotropy in the x 1 -x 2 plane is imprinted in the outgoing waves as differences in the spatial profiles along the x 1 and x 2 directions, as well as in the shape of the apparent horizon. k/Λ = 0 k/Λ = 2 k/Λ = 4 k/Λ = 6 n = 1 9.6, 10.4, 10.5 (10.7) 9.8, 10.7, 10.9 (11.1) 10.9, 11.7, 11.8 (12.1) 11.7, 13.1, 13.6 (13.7) n = 2 16.8, 17. TABLE I. Frequencies ω n (k)/Λ of some low-lying modes with k 1 = k 2 = k extracted from sinusoidal fits to the projections φ n, k (t) shown in Fig 2. The first three entries in each cell correspond to fits within the three time windows tΛ ∈ (0, 1.6), (0.8, 2.4), (1.6, 3.2), respectively. The values inside the parentheses correspond to the normal frequencies of the AdS soliton.
Projections φ n, k (t) of the bulk scalar field ϕ onto the normal modes of the AdS soliton according to (6), with k 1 = k 2 ≡ k. Times indicated on the horizontal axis are in multiples of t bounce ≈ 0.8Λ −1 .
As outgoing waves disperse away into the asymptotic region x 1 , x 2 → ±∞, the scalar field is expected to be described by a superposition of the discrete and gapped set of normal modes of the pure AdS soliton background. Each mode is characterized by a radial quantum number n with a corresponding radial profile Ψ n (ρ), a momentum k a in the x a direction and an energy ω n (k) 2 = k 2 + m 2 n . In the gauge theory, these are precisely the asymptotic states of masses m n of the confining vacuum. The projections into these modes are where µ(ρ) = 2ρ/(1 − ρ 2 ) is the measure factor that ensures unit normalization of the Ψ n . At asymptotically late times we expect the different modes to decouple from one another and, hence, that φ n, k (t) → c n, k exp [−iω n (k)t], where the coefficients c n, k are time independent and each mode oscillates with well-defined frequency ω n (k). Figure 2 shows the result of applying (6) to project the scalar field ϕ(t, ρ, x a ) from the fully nonlinear evolution into some of the lowestlying modes. Table I lists the frequencies ω n (k) extracted from sinusoidal fits to these projections within different time windows. On the one hand, we see that these frequencies approach the known normal frequencies of a scalar field in the AdS soliton [45] as time progresses, showing that the latetime behavior can be approximately interpreted in terms of the expected asymptotic states. On the other hand, the small discrepancies in the frequencies at the latest time window, and the fact that the projections shown in Fig 2 still exhibit amplitude modulations, suggest that nonlinear mode-mixing is still taking place for the evolution times considered here. In gauge theory language, complete freeze-out has not yet been reached.
It is illustrative to see the time evolution from the view- Snapshot of the gauge theory energy density at t = 1.6Λ −1 ≈ 2t bounce . point of the boundary stress tensor. Figure 3 shows a snapshot of the energy density, whose bulk dual is the finiteenergy black hole depicted in the second panel of Fig 1. Dispersing waves in the bulk give rise to the ripples visible in Fig 3, leaving behind a black hole that is imprinted on the boundary as a central lump of localized energy density that oscillates around a nonzero spatial profile, as it settles down to equilibrium. We have checked that the late-time energydensity profile of this dynamical plasma ball is approaching the profile of an equilibrium plasma ball of [42] with a temperature T equil ≈ 2.1 T c (see below). This means that, at late times, our plasma balls are better approximated by small and nearly spherically symmetric Schwarzschild black holes than by pancakelike configurations. The spacelike eigenvalues P 1 , P 2 , P θ of the boundary stress tensor correspond to the pressures in the local rest frame. Figure 4 depicts the time dependence of all the eigenvalues of the stress tensor at the center of the plasma ball, x 1 = x 2 = 0. Each curve in Fig 4  is periodically disturbed, first at t = t bounce /2 when waves from the IR bottom first reach the boundary and, thereafter, in intervals of ∆t = t bounce . To emphasize the fact that these bounces are associated with the AdS soliton geometry in Figure 4, we also show the energy density for a simulation that does not lead to the formation of a black hole.
Spatial profiles of the gauge theory energy density are shown in Fig 5 at several representative times. As waves from the IR bottom reach the boundary, which happens around t = t bounce /2, the shorter-wavelength spatial inhomogeneities in the IR become more evidently imprinted in the boundary one-point functions. These short-wavelength features are subsequently either (i) carried away to the asymptotic region x 1 , x 2 → ±∞ by outgoing waves with momentum in the boundary directions, or (ii) carried by radially propagating waves back to the vicinity of the black hole, where they are efficiently absorbed. Both mechanisms for shedding short-wavelength features are consistent with the bulk snapshots in Figure 1. For comparison, Fig 5 also shows the spatial profile of an equilibrium plasma ball with T = 2.1 T c [42]. We see that this is close to, but not exactly on top of, the profile of the latest-time dynamical plasma ball. This small discrepancy is expected because the dynamical ball has not yet equilibrated at the time shown. Moreover, the discrepancy is larger away from the center of the ball, since there the oscillations are more pronounced than at the center.
Discussion.-We have presented the first study of the realtime evolution for finite-energy black holes in the AdS soliton background. The black holes are dual to plasma balls: localized droplets of deconfined matter surrounded by the confining vacuum. A collision with sufficiently high but finite energy in a confining theory with a gravity dual will generically produce an excited plasma ball. Our results suggest that small excited plasma balls relax to classically nonlinearly stable configurations with long-lived excitations. Qualitatively, the radiation emitted during the relaxation process can be understood as consisting of two components. The first one disperses away to infinity in the noncompact directions, and at late times it has an interpretation in terms of individual particles. The entropy carried away by these modes is O(1) and, hence, subleading with respect to the O(N 2 c ) entropy that remains in the plasma ball. In other words, in our classical approximation the radiated field is a coherent field with effectively zero entropy. In this approximation, the black hole and the dual plasma ball eventually settle down to equilibrium, whereas in a finite-N c theory like QCD the black hole would Hawking evaporate and the plasma ball would hadronize. The second component of the radiation is associated with long-lived waves bouncing back and forth between the IR bottom and the AdS boundary at intervals t bounce ≈ 0.8Λ −1 . These waves cause periodic disturbances of the black hole and lead to an equilibration time longer than what would be naively expected. For example, an estimate based on the final temperature would give Λt equil ≈ Λ/T equil ≈ 0.5 for the plasma ball studied here.
Our analysis suggests that these long-lived disturbances are a robust feature of small, finite-energy black holes in a confining background, or equivalently, of small plasma balls in large-N c gauge theories with a gravity dual. This feature can be understood as follows. Modes with wavelengths much longer than the size of the black hole interact weakly with it and easily disperse away in the noncompact directions, and modes with wavelengths much shorter than the size of the black hole are efficiently absorbed by it. Meanwhile, modes with wavelengths comparable to the size of the black hole interact strongly enough to be attracted to it and, hence, do not disperse away easily, but at the same time, are not efficiently absorbed by it. These waves can bounce back and forth between the IR bottom and the AdS boundary.
In the gauge theory, these periodic disturbances correspond to a periodic transfer of energy between IR and ultraviolet (UV) modes inside the plasma ball, coupled to oscillations in the shape of the ball. We have seen that they have an important effect on the dynamics of the plasma ball that we have considered. Since its mass is M < Λ, this plasma ball is not a good model for the fireball created in a HIC. For this reason we are currently investigating whether or not the periodic disturbances persist for large plasma balls. The obvious next step in our program is the study of finite-energy black hole collisions, which is a natural mechanism by which these large plasma balls can be formed. These will offer new holographic insights on confinement and finite-size effects in the equilibration of large-N c , strongly coupled gauge theories. Some aspects of these theories will certainly differ from those in finite-N c theories. For example, the final state in these theories at asymptotically late times will contain an equilibrated plasma ball that will not hadronize. However, some features in the evolution may be similar between finite and infinite N c . Assessing which ones will require further analyses beyond this first study.

Supplemental Material Preliminaries
The AdS soliton [33] in d + 1 spacetime dimensions is a double Wick rotated black brane, with an R 1,d−2 × S 1 boundary. In this paper, we specialize to d = 4. Starting with the metric g 5 of the black brane on the Poincaré patch of AdS 5 the double Wick rotation t = iθ , with period θ ∼ θ + ∆θ , and x 3 = it gives the metricĝ of the AdS soliton The period ∆θ of the circle parametrized by θ is fixed to ∆θ = (4π/d)z 0 = πz 0 in order to ensure that the spacetime z > z 0 terminates smoothly at z = z 0 . This period ∆θ sets the confinement scale Λ = 1/∆θ = 1/(πz 0 ) of the boundary gauge theory, so that the confinement/deconfinement temperature is T c = Λ. In this paper we consider dynamical (i.e., out-of-equilibrium) finiteenergy black holes sitting at z = z 0 and localized in the x 1 -x 2 plane. We can bring the boundary to a finite coordinate location by introducing a coordinate ρ through so that the boundary is located at ρ = and the IR bottom is located at ρ = 0. The AdS soliton metric in terms of this compactified coordinate iŝ 4 and we have chosen z 0 = 1 and = 1 without loss of generality.

Evolved Variables
The metric is evolved in terms of variablesḡ µν , which are constructed out of the full spacetime metric g µν and the AdS soliton metricĝ µν : where ρ is the compactified holographic coordinate defined in the previous section. Similarly, the generalized harmonic (GH) source functions is expressed in terms of evolved vari-ablesH µ that are constructed out of the full spacetime GH source functions H µ and the valuesĤ µ they evaluate to in the AdS soliton: Finally, the scalar field is also evolved in terms of a variablē ϕ: In the GH formalism, a gauge choice is made by specifying GH source functions H µ . Here, we make this choice by specifying the evolved variables (S8) using the following specific form: are the initial values of the source functions, and F µ is the target gauge that we calculated using the procedure outlined in Ref. [9] F t ≡ 2 f 1ḡtρ , (S10) Here, .

Boundary Conditions
We set Dirichlet boundary conditions at spatial infinity for the metric, source functions, and scalar field: (S11) At the IR bottom, ρ = 0, we impose regularity conditions on all evolved variables according to their even or odd character there: (S12) Local flatness near the IR bottom ρ = 0 imposes an additional relation among the metric variables, ensuring that no conical singularities arise there. In terms of our regularized metric variables defined in (S7), the absence of a conical singularity at ρ = 0 implies We impose this condition onḡ θ θ at ρ = 0 instead of imposing the regularity condition for this variable as in (S12).
In the noncompact directions x 1 , x 2 , we impose Dirichlet boundary conditions for the metric, source functions, and scalar field at an outer boundary x 1 , x 2 = ±x max : (S14) We ensure that x max is chosen to be sufficiently large in units of Λ −1 to avoid spurious boundary effects throughout the course of the evolution times considered. For the near boundary gauge choice (S10) to be consistent with the ρ-component of the GH constraints C µ = H µ − x µ , we haveḡ tt −ḡ ρρ /4 −ḡ x 1 x 1 −ḡ x 2 x 2 −ḡ θ θ = 0, to leading order in the approach to ρ = 1. We thus ensure that on the initial slice, which is preserved by evolution that satisfies the GH constraints.

Time evolution
We monitor the evolution of black holes by keeping track of trapped surfaces. We excise a portion of the interior of any apparent horizon (AH) that forms to remove any singularities from the computational domain. No boundary conditions are imposed on the excision surface; instead, the Einstein equations are solved there using one-sided stencils. We use Kreiss-Oliger dissipation [46] to damp unphysical highfrequency modes that can arise at grid boundaries, with a typical dissipation parameter of 0.35.
We numerically solve the Hamiltonian constraint for initial data at t = 0 and the Einstein equations for subsequent time slices t > 0 using the PAMR/AMRD libraries [47], and discretize the equations using second-order finite differences. The evolution equations for the metric and scalar field are integrated in time using an iterative Newton-Gauss-Seidel relaxation procedure. The numerical grid is in (t, ρ, where we typically use t max = (10/π)Λ −1 and x max = 3Λ −1 . A typical base grid has N ρ = 33, N x 1 = N x 2 = 385 grid points with an addition 4 levels of refinement, and with equal grid spacings ∆ρ = ∆x 1 = ∆x 2 . We use a typical Courant factor of λ ≡ ∆t/∆ρ = 0.1. The results presented here were obtained with fixed mesh refinement centered around the plasma ball, although the code has adaptive mesh refinement capabilities.

Numerical tests
To check that our numerical solutions are converging to a solution of the Einstein equations, we compute an independent residual by taking the numerical solution and substituting it into a discretized version of G µν + Λg µν − 8πT µν . Since the numerical solution was found solving the generalized harmonic form of the Einstein equations, we expect the independent residual to only be due to numerical truncation error and thus converge to zero at the quadratic rate determined by our second-order accurate discretization. We can compute a convergence factor for the independent residual by Here, f h denotes a component of G µν + Λg µν − 8πT µν , and a 0 h, a 1 h are the grid spacings of two different resolutions. Given our second-order accurate finite difference stencils, we expect Q to approach Q = 2 as h → 0. FIG. S1. Convergence factors for the independent residual from the same simulation. The L 2 norm of the convergence factors is taken over the entire grid.