How big are the smallest drops of quark-gluon plasma?

Using holographic duality, we present results for both head-on and off-center collisions of Gaussian shock waves in strongly coupled N=4\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \mathcal{N}=4 $$\end{document} supersymmetric Yang-Mills theory. The shock waves superficially resemble Lorentz contracted colliding protons. The collisions results in the formation of a plasma whose evolution is well described by viscous hydrodynamics. The size of the produced droplet is R ∼ 1/Teff where Teff is the effective temperature, which is the characteristic microscopic scale in strongly coupled plasma. These results demonstrate the applicability of hydrodynamics to microscopically small systems and bolster the notion that hydrodynamics can be applied to heavy-light ion collisions as well as some proton-proton collisions.


Introduction
Hydrodynamics is a long wavelength late time effective description of the transport of conserved degrees of freedom. Hydrodynamics lacks many of the excitations of the underlying quantum theory (i.e. those excitations which are not conserved). Other than their effect on transport coefficients, it is consistent to neglect such excitations in the macroscopic limit since they attenuate over microscopic scales.
A striking application of hydrodynamics is the modeling of heavy-light ion collisions and proton-proton collisions at RHIC and the LHC, where signatures of collective flow have been observed [1][2][3][4][5][6] which are consistent with hydrodynamic evolution of quarkgluon plasma [7][8][9][10][11]. The success of hydrodynamics is surprising since at experimentally accessible energies, microscopic scales such at the mean free path are probably not too different from the system size. This raises several interesting questions. For such tiny systems is it theoretically consistent to neglect nonhydrodynamic degrees of freedom? Does the presence of large gradients inherent to small systems excite nonhydrodynamic modes or otherwise spoil the hydrodynamic gradient expansion? How big are the smallest drops of matter which can meaningfully be identified as liquids? Are we deluding ourselves into thinking hydrodynamics can apply to systems as small as a proton?
Since microscopic length and time scales typically becomes smaller in the limit of strong coupling, it is natural to expect the domain of utility of hydrodynamics to be largest at strong coupling. However, strongly coupled dynamics in QCD are notoriously difficult to study. It is therefore of great utility to have a toy model of strongly coupled dynamics -which accounts for both hydrodynamic and nonhydrodynamic excitations -where one can answer the above questions in a controlled and systematic setting.
Holographic duality [12] maps the dynamics of certain strongly coupled non-Abelian gauge theories onto the dynamics of classical gravity in higher dimensions. The process of quark-gluon plasma formation maps onto the process of gravitational collapse and black hole formation, with the ring down of the black hole encoding the relaxation of the plasma JHEP03(2016)146 to a hydrodynamic description. The dual gravity calculation is microscopically complete, meaning it contains both hydrodynamic and nonhydrodynamic excitations. In other words, all stages of evolution -from far-from-equilibrium dynamics to hydrodynamics -are all encoded Einstein's equations. Since Einstein's equations can be solved numerically, holographic duality provides a unique arena where one can systematically test the domain of utility of hydrodynamics. We shall therefore use holographic duality to test hydrodynamics in extreme conditions.
The simplest theory with a dual gravity description is N = 4 supersymmetric Yang-Mills theory (SYM), which is dual to gravity in ten dimensional AdS 5 ×S 5 spacetime. A simple model of quark-gluon plasma production in SYM is the collision of gravitational shock waves [13][14][15][16][17][18][19], which can result in the formation of a black hole in AdS 5 ×S 5 . The shock's waveform, which is not fixed by Einstein's equations, determines the expectation value of the stress tensor in the dual field theory. Previously, in [18] we considered the collision of a localized shock wave with a planar shock wave, which in the dual field theory resembled a proton-nucleus collision. There it was found that the collision resulted in a microscopically small droplet plasma being produced, whose evolution was well described by viscous hydrodynamics. Here we shall extend our work to include the collisions which resemble proton-proton collisions. In particular, we shall consider the collision two shock waves moving in the ±z directions at the speed of light with expectation value of the SYM energy density where N c is the number of colors, µ is an energy scale, δ w is a smeared δ function, and x ⊥ = {x, y} are the coordinates transverse to the collision axis. Hence our shock waves are Gaussians with transverse width σ, impact parameter b, and total energy Superficially at least, the shocks resemble Lorentz contracted colliding protons. We shall refer to the shocks as "protons" with the quotes to emphasize to the reader than we are not colliding bound states in QCD but rather caricatures in SYM. The precollision bulk geometry contains a trapped surface and the collision results in the formation of a black hole. We numerically solve the Einstein equations for the geometry after the collision and report on the evolution of the expectation value of the SYM stress tensor T µν and test the validity of hydrodynamics. In strongly coupled SYM the hydrodynamic gradient expansion is an expansion in powers of 1/( T eff ) with the characteristic scale over which the stress tensor varies and T eff the effective temperature. Therefore, in strongly coupled SYM the scale 1/T eff plays the role of a mean free path. We focus on the low energy limit where T eff is small and 1/T eff is large. Without loss of generality, in the results presented below we measure all quantities in units of µ = 1. We choose shock parameters

JHEP03(2016)146
We shall refer to the b = 0 and b = 3x collisions simply as "head-on" and "off-center," respectively. Both collisions result in the formation of a microscopically small droplet plasma of size R ∼ 1/T eff , whose evolution is well described by viscous hydrodynamics. These results demonstrate that hydrodynamics can work in microscopically small systems and bolster the notion that the collisional debris in heavy-light ion collisions and protonproton collisions can be modeled using hydrodynamics. An outline of the rest of the paper is as follows. In section 2 we construct a simple test to see whether the evolution of T µν is governed by hydrodynamics. In section 3 we outline the gravitational formulation of the problem. In section 4 we present our results for the evolution of T µν and in section 5 we discuss our results and generalizations to confining theories. Readers not interesting the details of the gravitational calculation can skip section 3.

A litmus test for hydrodynamic evolution of T µν
With the gravitational calculation presented below we shall obtain the exact expectation value of the SYM stress tensor T µν . At time and length scales than microscopic scales the evolution of T µν must be governed by hydrodynamics. This means that where T µν hydro ( , u µ ) is given in terms of the proper energy density and the fluid velocity u µ via the hydrodynamic constitutive relations. When (2.1) is satisfied we shall say that the system has hydrodynamized. In this section our goal is develop a simple test, which requires only T µν , to see whether the system has hydrodynamized. We shall not give a through review of hydrodynamics here and instead shall only state the salient features of hydrodynamics needed for our analysis. For more detailed discussions of relativistic hydrodynamics we refer the reader to [20,21].
We use the mostly positive Minkowski space metric η µν = diag(−1, 1, 1, 1) and define the fluid velocity to be the timelike (u 0 > 0) normalized (u µ u µ = −1) eigenvector of T µν hydro with − the associated eigenvalue, With these conventions, at second order in gradients the hydrodynamic constitutive relations in strongly coupled SYM read [20,22] T µν where p is the pressure given by the conformal equation of state and T µν (1) and T µν (2) are the first and second order gradient corrections, respectively. Explicitly,

JHEP03(2016)146
where D ≡ u · ∂, and the shear σ µν and vorticity Ω µν tensors are with the transverse projector P µν ≡ g µν + u µ u ν . (2.7) Note u µ P µν = 0 since u µ u µ = −1. For any tensor A µν the bracketed tensor A µν is defined to be the symmetric (A µν = A νµ ), traceless (η µν A µν = 0) and transverse (u µ A µν = 0) component of A µν , meaning In strongly coupled SYM the shear viscosity η, second order transport coefficients τ Π , λ 1 , λ 2 and effective temperature T eff read [20,22,23] The hydrodynamic equations of motion for and u µ are given by the energy-momentum conservation equation ∂ µ T µν hydro = 0. One option to test for hydrodynamization is to construct initial data for and u µ and then to evolve and u µ forward in time via the hydrodynamic equations of motion. One can then reconstruct T µν hydro and compare it to the full stress T µν to test the applicability of hydrodynamics.
However, given that we will have the full expectation value T µν , solving the hydrodynamic equations of motion is unnecessary. Consider the eigenvalues p (λ) and associated eigenvectors e µ (λ) of T µ ν T µ ν e ν (λ) = p (λ) e µ (λ) , (2. 10) with no sum over λ implied. If the system has hydrodynamized, then T µ ν should have one timelike eigenvector e µ (0) , which according to (2.1) and (2.2), is just the fluid velocity. We therefore define the fluid velocity and proper energy to be With the complete spacetime dependence of the flow field u µ and energy density determined from the exact stress, we can construct T µν hydro using eqs. (2.3)-(2.9). It will be useful to quantify the degree in which the system has hydrodynamized. To this end we define the dimensionless residual measure where ||A µν || denotes the L 2 norm of A µν , which for symmetric A µν is the magnitude of the largest magnitude eigenvalue of A µν . Why do we take the max in (2.12)? Note that in general T µν − T µν hydro is nonmonotonic in time due to nonhydrodynamic modes,

JHEP03(2016)146
such as quasinormal modes, oscillating in time while decaying. Taking the max in (2.12) ameliorates the effect of the oscillations and aids in avoiding the false identification of hydrodynamic evolution if T µν − T µν hydro is momentarily small. If ∆(t, x) 1, then T µν is well described by the hydrodynamic constitutive relations at x at time t and all future times.
We therefore have achieved our goal of quantifying hydrodynamization purely in terms of T µν . To reiterate the strategy, one first extracts the proper energy density and fluid velocity u µ from T µν . One then uses the hydrodynamic constitutive relations to compute T µν hydro . Finally, one then constructs the hydrodynamic residual ∆ and looks for regions of spacetime where ∆ 1 to identify hydrodynamic behavior.

Gravitational description
According to holographic duality, the evolution of the expectation value of the stress tensor in strongly coupled SYM is encoded in the gravitational field in asymptotically AdS 5 ×S 5 spacetime. Einstein's equations are consistent with no dynamics on the S 5 . We shall make this assumption for our numerical analysis and revisit it in the Discussion section below. Hence we shall focus on gravitational dynamics in asymptotically AdS 5 spacetime.
Gravitational states dual to shock waves in SYM moving in the ±z directions can be constructed by looking for steady state solutions to Einstein's equations which only depend on time through the combination z ∓ ≡ z ∓ t [13,24]. Consider the Fefferman-Graham coordinate system ansatz for the metric, where x ≡ {x, y, z}, x ⊥ ≡ {x, y} and r is the AdS radial coordinate with r = ∞ the AdS boundary. Note that here and in what follows we set the AdS radius to unity. The anzatz (3.1) satisfies Einstein's equations provided h ± satisfy the linear partial differential equation The solution to (3.2) satisfying vanishing boundary conditions at the AdS boundary is where H ± is an arbitrary function. The metric (3.1) with h ± given by (3.3) consists of a gravitational wave moving in the ±z direction -parallel to the AdS boundary -at the speed of light. The fact that H ± is arbitrary reflects the fact that Einstein's equations do not fix the waveform of gravitational waves. Likewise, in the dual field theory nothing fixes the waveform of the shock waves. Indeed, SYM is conformal and contains no bound states or scales which would otherwise determine the structure of shock waves. We shall exploit the arbitrariness of H ± to construct states in SYM which superficially at least resemble localized and Lorentz JHEP03(2016)146 contracted protons. To do so we note that the Fourier transform of H ± , H ± , determines the expectation value of the SYM stress tensor via [25] with all other components vanishing. Therefore, once we specify the SYM energy density T 00 for a single shock, we specify the dual gravitational wave. We shall employ the single shock energy density given in eq. (1.1) in the Introduction. This means For the smeared δ function we choose a Gaussian with width µw = 0.375. Note that evolution inside the future light cone of planar shock collisions with width µw = 0.375 well approximates that of the δ function limit [17]. The shock parameters we employ are given in eq. (1.3) in the Introduction. At early times, t −w, the profiles H ± have negligible overlap and the precollision geometry can be constructed from (3.1) by replacing the last term with the sum of corresponding terms from left and right moving shocks. The resulting metric satisfies Einstein's equations, at early times, up to exponentially small errors.
To evolve the precollision geometry forward in time we use the characteristic formulation of gravitational dynamics in asymptotically AdS spacetimes. Here we simply outline the salient details of the characteristic formulation in asymptotically AdS spacetime. For details we refer the reader to [26]. Our metric ansatz reads with Greek indices denoting spacetime boundary coordinates, x µ = (t, x, y, z). Near the boundary, g µν = η µν + g (4) µν /r 4 + O(1/r 5 ). The subleading coefficients g (4) µν determine the SYM stress tensor [25] 00 . (3.8) An important practical matter in numerical relativity is determining the computational domain and excising singularities behind the event horizon. To excise singularities, we look for the position of an apparent horizon, which if it exists must lie behind an event horizon, and stop numerical evolution at its location. Note that the metric ansatz (3.7) is invariant under the residual diffeomorphism for arbitrary λ. If the apparent horizon has planar topology and lies at say r = r h (t, x), then we may use this residual diffeomorphism invariance to set the apparent horizon to be JHEP03(2016)146 at r = 1 by choosing λ = −1 + r h . Let ∇ be the covariant derivative under the spatial metric g ij . Fixing the apparent horizon to be at constant r, at the horizon the metric must satisfy [26] We note that eq. (3.10) is a constraint on initial data: if (3.10) is satisfied at any fixed time and appropriate boundary conditions are enforced, then Einstein's equations guarantee (3.10) remains satisfied at future times. See [26] for more details. Hence, we only need to enforce (3.10) on initial data. For given initial data we compute the l.h.s. of (3.10) at r = 1 and search for a radial shift λ such that (3.10) is satisfied to some chosen accuracy. The shift λ is computed using Newton's method. In our numerical simulations we demand that the norm of the l.h.s. of (3.10) is 5 × 10 −3 or smaller.
To generate initial data for our characteristic evolution, we numerically transform the precollision metric in Fefferman-Graham coordinates to the metric ansatz (3.7). We begin time evolution at t = −1.5. For both head-on and off-center initial data we find an apparent horizon at time t = −1.5before the collision on the boundary. Since the apparent horizon must lie inside an event horizon, this means that there is an event horizon even before the collision takes place on the boundary. Why must an event horizon exist before the collision? Empty AdS contains a cosmological horizon, which is simply the Poincare horizon. It takes a light ray an infinite amount of (boundary) time to travel from the Poincare horizon to any bulk point in the geometry. This means that if a planar topology event horizon exists at any time -say after the collision -then it must have also existed at all times in the past and coincided with the Poincare horizon at t = −∞. Similar conclusions were reached within the context of a quench studied in [27].
For numerical evolution we periodically compactly spatial directions with transverse size L x = L y = 18 and longitudinal size L z = 9 and and evolve from t = −1.5 to t = 2.5. For the head-on collision, evolution is performed using a spectral grid of size N x = N y = 39, N z = 155 and N r = 48 while for the off-center collision we increase the number of transverse points to N x = N y = 45. 1 To monitor the convergence of our numerical solution, we monitor violations of the horizon fixing condition (3.10). Since in the continuum limit (3.10) is an initial value constraint, violations of (3.10) after the initialization time t = −1.5 provide a simple measure of the degree in which our numerics are converging to the continuum limit. For both head-on and off-center collisions violations of (3.10) are 5 × 10 −3 or smaller at all times, which is the same size of the allowed violations of (3.10) at the initialization time. The fact that (3.10) remains satisfied after the initialization time provides a nontrivial check of the convergence of our numerics.

JHEP03(2016)146 4 Results
Define the rescaled stress In figure 1 we plot the rescaled energy density T 00 and the rescaled momentum density T 0i in the plane y = 0 for head-on and off-center collisions at several values of time. The color scaling in the plots of T 0i denotes the magnitude of the momentum density | T 0i |, while the streamlines indicate the direction of T 0i . At time t = −1.125 the systems consist of two separated "protons" centered on z = ±1.125, x = 0 for the head-on collision and z = ±1.125, x = ∓ b 2 for the off-center collision. In both collisions the "protons" move towards each other at the speed of light and collide at t = z = 0. Note that at t = 0 the momentum density of the head-on collision nearly cancels out. Indeed, at t = 0 the energy and momentum densities for both collisions are at the order 10% level just the linear superposition of that of the incoming "protons." This indicates that nonlinear interactions haven't yet modified the stress. Similar observations were made within the context of planar shock collisions in [16]. Nevertheless, at subsequent times we see that the collisions dramatically alter the outgoing state. For both collisions the amplitude of the energy and momentum densities near the forward light cone have decreased by order 50% relative to the past light cone, with the lost energy being deposited inside the forward light cone. Moreover, after t = 0 we see the presence of transverse flow for both impact parameters. As we demonstrate below, the evolution of the postcollision debris for both impact parameters is well described by viscous hydrodynamics.
In figure 2 we plot T xx and T zz and the hydrodynamic approximation T xx hydro and T zz hydro at the origin x = 0 for both impact parameters. At this point the stress tensor is diagonal, the flow velocity u = 0, and T xx and T zz are simply the pressures in the x and z directions. As shown in the figure, the pressures start off at zero before the collision. Note that the pressures are remarkably similar for both collisions. Evidently, the finite impact parameter employed here has little effect on the pressures at the origin. During the collisions the pressures increase dramatically, reflecting a system which is highly anisotropic and far from equilibrium. Nevertheless, after time the system has hydrodynamized at the origin and the evolution of the pressures for both impact parameters is well described by the hydrodynamic constitutive relations. Remarkably, at this time T xx is order 10 times greater than T zz . In figure 3 we plot the hydrodynamic residual ∆ at x = 0 for both impact parameters. We see from the figure that ∆ is large at times t < t hydro . However, as t → t hydro the residual ∆ rapidly decays to ∆ 0.2. After t = t hydro the residual continues to decay, but at a much slower rate than before t = t hydro . The rapid decay of ∆ prior to t = t hydro indicates the rapid decay of nonhydrodynamic modes.
Over what region of space has the system hydrodynamized? To study the spatial domain of applicability of hydrodynamics, in figure 4 Figure 3. The hydrodynamic residual ∆ at at the spatial origin, x = y = z = 0, as a function of time for head-on and off-center collisions. Note that ∆ need not be differentiable due to the max and matrix norm in its definition (2.12). Regions with ∆ 1 have hydrodynamized. Prior to t = t hydro ≈ 1.2, the residual is large and hydrodynamics is not a good description of the evolution. However, as t → t hydro the residual dramatically decreases. Thereafter, ∆ continues to decrease but at a much slower rate. The rapid decay of ∆ prior to t = t hydro reflects the rapid decay of nonhydrodynamic modes in the far-from-equilibrium state created by the collisions.  for both head-on and off-center collisions. Regions with ∆ 1 have hydrodynamized. Note that we have restricted the plot to ∆ < 1 in order to highlight the hydrodynamic behavior. The black curve in the plots is the surface ∆ = 0.2. For both impact parameters there is a crisply defined region -whose boundary is well approximated by the ∆ = 0.2 surface -where ∆ 1. We identify the matter in the interior of the ∆ = 0.2 surface as a droplet of liquid. Outside the ∆ = 0.2 surface ∆ rapidly increases, indicating the presence of nonhydrodynamic modes on the surface of the droplet. Note the irregularity in the off-center collision droplet shape in the x−z plane is due to nonhydrodynamic modes and not fluid rotation in the x−z plane. For the head-on collision the ∆ = 0.2 surface is circular in the x−y plane. In contrast, for the off-center collision the ∆ = 0.2 surface is elliptical in the x−y plane, with the the short axis of the ellipse oriented in the same direction as the impact parameter b = 3x. Nevertheless, for both collisions the transverse radius of the ∆ = 0.2 surface is roughly the same and equal to R ∼ 3, which is just the radius σ of our "protons." Note the irregularity in the off-center collision droplet shape in the x−z plane is due to nonhydrodynamic modes and not fluid rotation in the x−z plane. Indeed, in the interior of the droplet the vorticity is small with ||Ω µν || ∼ 0.1||σ µν ||. Additionally, note the difference in droplet shape in the x−y plane for head-on and off-center collisions. For the headon collision droplet's surface is circular in the x−y plane. In contrast, for the off-center collision the droplet's surface is elliptical in the x−y plane, with the the short axis of the ellipse oriented in the same direction as the impact parameter b = 3x. Nevertheless, for both collisions the transverse radius of the droplet is roughly the same and equal to R ∼ 3, which is just the radius σ of our "protons" employed in (1.3).

JHEP03(2016)146
head-on o↵-center (1) ||/p ||T µ⌫ (2) ||/p Figure 5. The effective temperature T eff and the first and second order gradient norms ||T µν (1) || and ||T µν (2) || in units of the average pressure p for head-on and off-center collisions. All plots are shown at time t = 1.25 t hydro = 1.5 and are restricted to the domain inside the ∆ = 0.2 surface shown in figure 4, which is also shown here as the solid curve. In this region hydrodynamics is a good approximation to the evolution of the stress. For both collisions the average effective temperature in the displayed region is T eff ∼ 0.25. Inside the ∆ = 0.2 surface ||T µν (2) || ||T µν (1) ||, meaning second order gradient corrections are negligible. Note however, that near the surface of the droplet, where nonhydrodynamic modes are excited, ||T µν (2) || begins to become comparable to ||T µν (1) ||.
In figures 5-7 we restrict our attention to time t = 1.25 t hydro = 1.5 and to the region inside the droplet of fluid. We explain the coloring of the different plots below. In the left column of figure 5 we show the effective temperature T eff , defined in eq. (2.9). From the figure it is evident that the average value of the temperature inside the droplet is T eff ∼ 0.25 for both collisions. We therefore obtain the dimensionless measure of the transverse size of the droplet, Similar conclusions were reached in [18], where a holographic model of proton-nucleus collisions was studied. Additionally, note indicating rapid hydrodynamization. Similar hydrodynamization times were observed in [16,18,28,29]. We therefore conclude that hydrodynamic evolution applies even when both the system size and time after the collision are on the order of microscopic scale 1/T eff .

JHEP03(2016)146
In the middle and right columns of figures 5 we plot the size of first and second order gradient corrections to the hydrodynamic constitutive relations, ||T µν (1) || and ||T µν (2) ||, normalized by the average pressure p. In the interior of the droplet ||T µν (2) || is nearly an order of magnitude smaller than ||T µν (1) || for both impact parameters, meaning that second order gradient corrections are negligible. 2

Discussion
Given that 1/T eff is the salient microscopic scale in strongly coupled plasma -akin to the mean free path at weak coupling -it is remarkable that hydrodynamics can describe the evolution of systems as small at 1/T eff . By setting the "proton" radius equal to the actual proton radius, σ ∼ 1 fm, and using N c = 3 colors, the single shock energy (1.2) is ∼ 20 GeV. Likewise, the effective temperature of the produced plasma is T eff ∼ 200 MeV. Given that collisions at RHIC and the LHC have higher energies and temperatures, it is natural to expect RT eff to be larger in RHIC and LHC collisions than the simulated collision presented here. We therefore conclude that droplets of liquid the size of the proton need not be thought of as unnaturally small. Evidently, there are theories -namely strongly coupled SYM -which enjoy hydrodynamic evolution in even smaller systems.
While hydrodynamics is a good description of the evolution of our collision, it should be noted that the produced hydrodynamic flow has an extreme character to it. In particular, when the system size is on the order of the microscopic scale 1/T eff , gradients must be large. To quantify the size of gradients we introduce the anisotropy function where the pressures p (i) are determined by the eigenvalue equation (2.10) and again p is the average pressure, p = 3 = avg(p (i) ). In ideal hydrodynamics, where all the pressures p (i) are equal, the anisotropy vanishes. It therefore follows that after the system has hydrodynamized, any nonzero anisotropy must be due to gradient corrections in the hydrodynamic constitutive relations (2.3). In the left column of figure 6 we plot A at time t = 1.25 t hydro = 1.5 in the interior of the droplet. For both impact parameters A ∼ 1, with larger anisotropies near the droplet's surface. Evidently, gradients are large and ideal hydrodynamics is not a good approximation. This is also evident in the pressures displayed in figure 2, where the transverse and longitudinal pressures differ by a factor of 5 to 10. Moreover, as shown in figure 5, ||T µν (1) || ∼ p, so the first order gradient correction T µν (1) is comparable to the ideal hydrodynamic stress. However, despite the presence of large gradients, the above observed large anisotropy must be almost entirely due to the first 2 We note however, that inside the ∆ = 0.2 surface shown in figure 4 we have || T µν − T µν hydro || ∼ ||T µν (2) ||. Given that T µν hydro is computed to second order in gradients, this could mean that third order gradient corrections to T µν hydro are comparable to the second order gradient corrections. Alternatively, and in our opinion more likely, this could reflect the fact that the spectrum of the discretized Einstein equations we solved is different from that of the continuum limit used to obtain the gradient expansion of T µν hydro . This difference must manifest itself at suitably high order in gradients and can be ameliorated by using a finer discretization scheme.

JHEP03(2016)146
head-on o↵-center  figure 4, which is also shown here as the solid curve. In this region hydrodynamics is a good approximation of the evolution of the stress. In the limit of ideal hydrodynamics the anisotropy vanishes. The order 1 anisotropy observed here for both collisions reflects the fact that there are large gradients in the system. Likewise, for both collisions we see that v x and v z are similar in magnitude. The large transverse velocity is also a signature of large gradients.
order gradient correction T µν (1) alone, since ||T µν (2) || is nearly an order of magnitude smaller than ||T µν (1) || inside the droplet of fluid. This means that viscous hydrodynamics alone is sufficient to capture the observed large anisotropies and evolution of the system.
It is remarkable that the second order gradient correction T µν (2) is negligible even when gradients are large enough to produce an order 1 anisotropy. It is also remarkable that the presence of large gradients does not substantially excite nonhydrodynamic modes (e.g. nonhydrodynamic quasinormal modes). In an infinite box of liquid, nonhydrodynamic modes decay over a timescale of order 1/T eff . Evidently, in the interior of the small droplet studied here, the rate of decay of nonhydrodynamic modes is stronger than the rate in which they are excited.
We now turn to the nature of the produced hydrodynamic flow. Also shown in figure  driven by transverse gradients, the rapid development of transverse flow observed here must reflect the presence of large transverse gradients. In figure 7 we plot the 3-velocity in the x−y plane at time t = 1.25 t hydro in the interior of the droplet produced in the off-center collision. The coloring in the plot denotes the magnitude of the fluid 3-velocity while the streamlines denote its direction. The fluid velocity is not radial, which is natural given the finite impact parameter and the elliptical shape of the droplet in the x−y plane. However, the flow of momentum is nearly radial. To quantify this, define the Fourier coefficient c n ≡ d 3 x T 0i ρ i cos nφ whereρ = cos φx + sin φŷ is the radial unit vector with φ the azimuthal angle. Deviations in radial flow are encoded in c 2 . However, at time t = 1.25 t hydro , the coefficient c 2 is only 1% of c 0 , indicating small deviations from radial momentum flow. It would be interesting to see how c 2 grows as time progresses. It would also be interesting to see how c 2 is affected by different choices of shock profiles. Indeed, off-center collisions of Gaussian shock profiles with large transverse width σ generate purely radial flow [19], since the overlap of two off-center Gaussians is still a Gaussian and thus rotationally invariant about the z axis.
Let us now ask how the droplet of fluid produced in the collision further evolves at later times. In conformal SYM the droplet must expand forever with size R ∼ t. This is very different from a confining theory, where the droplet cannot expand indefinitely. Likewise, as the droplet expands it must cool. Its easy to reason from energy conservation that the temperature must decrease slower than 1/t, meaning that RT eff → ∞ as t → ∞, so hydrodynamics becomes a better and better description of the evolution. We expect the surface of the droplet to become smoother and smoother as time progresses. Likewise, we expect nonhydrodynamic modes near the surface of the droplet to continue to decay with the fraction of energy behaving hydrodynamically approaching 100% as t → ∞. Indeed, such behavior was observed within the context of planar shock collisions [17], where nonhydrodynamic modes near the light cone decayed, with the lost energy transported inside the light cone where it hydrodynamizes.

JHEP03(2016)146
It would be interesting to push our analysis further. How big are the smallest drops of liquid? Clearly RT eff cannot be made arbitrarily small since the hydrodynamic gradient expansion (2.3) breaks down when RT eff 1. Moreover, in the gravitational description there likely exists a critical energy E c (or alternatively a critical impact parameter b c [30]) below (above) which no black hole is formed and which for E = E c + 0 + (or b = b c + 0 − ) critical gravitational collapse occurs [31]. The absence of a black hole when E < E c (or b > b c ) means that in the dual field theory the collisional debris will not evolve hydrodynamically at any future time. Clearly it would be interesting to look for criticality and the associated signatures of hydrodynamics turning off.
It would also be of great interest to study collisions in confining theories. Aside from being more realistic, collisions in confining theories have the added feature that the produced droplet -a plasma ball -cannot expand indefinitely. As the system cools it must eventually reach the confinement/deconfinement transition and freeze into a gas of hadrons. In the large N c limit freezeout is suppressed and plasma balls are metastable with a lifetime of order N 2 c [32]. At N c = ∞, plasma balls must eventually equilibrate and become static due to internal frictional forces. Therefore, in a large N c confining theory the question of how big are the smallest drops of quark-gluon plasma is tantamount to how big are the smallest plasma balls. An illuminating warmup problem is that of SYM on a three sphere of radius R. In the large N c limit the theory becomes confining [33][34][35]. Therefore, instead of studying plasma balls produced in a collision, one can study plasma confined to the three sphere and analyze the plasma's behavior as the radius of the sphere is changed.
The dual gravitational description of SYM on a three sphere is that of gravity in global AdS 5 ×S 5 . At strong coupling, and in the canonical ensemble, the deconfinement temperature T c is RT c = 3 2π ≈ 0.48, (5.2) with the deconfinement transition manifesting itself in the gravity description as the Hawking-Page phase transition [33]. Above the deconfinement transition the dual geometry contains a black hole, the boundary energy density is order N 2 c , and the boundary state is that of a quark-gluon plasma with hydrodynamic evolution. Below the transition the dual geometry is just thermal AdS and the boundary state has order N 0 c energy and does not admit a hydrodynamic description. Hence in the canonical ensemble of SYM living on a three sphere, the answer to the question of how big are the smallest drops of quark-gluon plasma is given by (5.2).
The situation is less clear in the microcannonical ensemble. In this case the dual black hole geometry experiences a Gregory-Laflamme instability [36] on the S 5 [37][38][39]. The instability sets in at RT eff ≈ 0.50. (5.3) It is unknown what the final state of the Gregory-Laflamme instability is. One possibility is that there do not exist stable black hole solutions below (5.3). In this case the instability could go on indefinitely, perhaps generating structures on the order of the Planck scale where quantum effects will be important. Another possibility is that there exists JHEP03(2016)146 equilibrium black holes with broken rotational invariance on the S 5 . In this case the final state of the Gregory-Laflamme instability would asymptote to this solution. In either case, it seems likely that the onset of the Gregory-Laflamme instability destroys hydrodynamic evolution in the dual field theory. It therefore seems that in either ensemble, the smallest drops of quark-gluon plasma in the confining phase of SYM have size RT eff ≈ 1/2.