A Holographic Path to the Turbulent Side of Gravity

We study the dynamics of a 2+1 dimensional relativistic viscous conformal fluid in Minkowski spacetime. Such fluid solutions arise as duals, under the"gravity/fluid correspondence", to 3+1 dimensional asymptotically anti-de Sitter (AAdS) black brane solutions to the Einstein equation. We examine stability properties of shear flows, which correspond to hydrodynamic quasinormal modes of the black brane. We find that, for sufficiently high Reynolds number, the solution undergoes an inverse turbulent cascade to long wavelength modes. We then map this fluid solution, via the gravity/fluid duality, into a bulk metric. This suggests a new and interesting feature of the behavior of perturbed AAdS black holes and black branes, which is not readily captured by a standard quasinormal mode analysis. Namely, for sufficiently large perturbed black objects (with long-lived quasinormal modes), nonlinear effects transfer energy from short to long wavelength modes via a turbulent cascade within the metric perturbation. As long wavelength modes have slower decay, this lengthens the overall lifetime of the perturbation. We also discuss various implications of this behavior, including expectations for higher dimensions, and the possibility of predicting turbulence in more general gravitational scenarios.


I. INTRODUCTION
The AdS/CFT correspondence [1,2] proposes a remarkable connection between quantum gravity in d+1 dimensions and quantum field theory in d dimensions. In a certain classical limit, this correspondence can be utilized to link the behavior of perturbed asymptotically anti-de Sitter (AAdS) black branes in general relativity to that of viscous conformal fluids on the AdS boundary, provided the perturbations are of sufficiently long wavelength [3][4][5][6]. This limit of the AdS/CFT correspondence is known as the gravity/fluid correspondence.
The gravity/fluid correspondence can also be derived on its own in a purely classical manner without any appeal to AdS/CFT, as a derivative expansion within general relativity (see, e.g., Ref. [4]). This derivation provides an explicit perturbative mapping between solutions, which can be exploited to relate gravitational and fluid behavior. In particular, interesting known phenomena on one side of the duality should have counterparts on the other, which can lead to new predictions, or to new methods of analysis. This has been used as a means to frame fluid dynamics questions in terms of gravitational physics. For example, Ref. [7] explored the relation between the Penrose inequalities -which predict the onset of naked singularities in general relativity -and finite-time blowup of solutions in hydrodynamics. In Ref. [8], it was suggested that a gravity dual could be utilized to understand the complex phenomenon of fluid turbulence. In the present work, following the analysis presented in Ref. [9], we follow the opposite route, namely to study the implications that turbulent phenomena -that can arise in fluid dynamics -have for our understanding of general relativity.
Turbulence is a ubiquitous property of fluid flows observed in nature at sufficiently high Reynolds number, R. Such behavior has recently been shown to also arise in inviscid conformal relativistic hydrodynamics [9,10]. While the actual fluid dual to an AAdS black brane has nonzero shear viscosity, this viscosity is subleading in the black hole temperature, so the inviscid approximation is valid at sufficiently high temperature [9]. This suggests that there should be a corresponding regime where long-wavelength black hole perturbations in asymptotically AdS spacetimes behave in a turbulent manner. Furthermore, in 2 + 1 dimensions such inviscid conformal fluids display an inverse cascade of energy to large scales [9], in accord with intuition from Navier-Stokes fluids [11]. This ensures that provided the initial condition falls within the regime of applicability of the gravity/fluid correspondence (i.e., sufficiently long wavelength perturbations), so should its time evolution, and therefore there should exist a black brane that behaves in a dual manner. The intuition described here has been borne out in very recent ground-breaking work [12], which demonstrated the development of turbulence in gravitational perturbations in 3 + 1 spacetime dimensions. Thus, gravitational behavior in this regime is effectively captured by a hydrodynamic analysis.
Drawing again on intuition from fluids, one should also expect on the gravity side, behavior akin both to turbulent and laminar flow. It is important to emphasize that these two phenomena can arise on the fluid side irrespective of the velocity of the background flow alone. Rather, the behavior depends on the value of the Reynolds number, where ρ, v, L, and η are the characteristic energy density, velocity fluctuation, distance scale, and shear viscosity, describing the flow, respectively 1 . For high values of R, turbulence occurs, whereas for small values the flow is laminar. These results and observations concerning the turbulent nature of perturbed AAdS black branes appear to be in tension with the standard expectation that such perturbations decay exponentially via quasinormal modes [13,14]. Indeed, for small-amplitude gravitational perturbations (which are dual to fluid flows with small velocity fluctuations v), a linear analysis should be valid. Due to the symmetries of the black brane, a mode decomposition is then possible. Such an analysis indicates that the modes decay in time as radiation is absorbed by the black brane (the only place where energy is lost, as the AdS boundary acts a mirror). Quasinormal modes of black branes in AdS can be grouped into three "channels" based upon their transformation properties under rotations: the sound, shear, and scalar channels. The longest lived families of modes within the sound and shear channels, in turn, are known as the "hydrodynamic" quasinormal modes of the black brane [15]. The dual fluid captures the behavior associated with these quasinormal modes. (Of course, the fluid satisfies nonlinear equations, whereas quasinormal modes are solutions to linear equations.) Higher quasinormal modes decay more rapidly, and are neglected in the gravity/fluid correspondence.
Of course, any behavior dual to quasinormal mode decay is surely absent in analyses of fluids with vanishing viscosity. In this paper, in order to examine this issue more closely, we extend the previous analysis of Ref. [9] to include viscosity, which captures the role of the black brane as a sink of energy. We numerically study turbulent (and laminar) solutions for the viscous relativistic conformal fluid (in d = 2 + 1) which arises in the gravity/fluid correspondence. We then contrast our results with the expectations we have laid out for the gravity dual, and we draw conclusions about the regime of applicability of linear perturbation theory about black holes.
In the following section we review the gravity/fluid correspondence in more detail. We sketch the perturbative derivation from general relativity, and we write down the relevant equations for our work. The dissipative relativistic hydrodynamic equations are closely related to those of Israel and Stewart [16][17][18], and are thus suitable for numerical implementation [3]. In Sec. III we proceed to describe our numerical setup, as well as the initial data. We work in Minkowski spacetime on R × T 2 , which is dual to a (periodically identified) black brane 2 in a Poincaré patch of AdS. Our initial data consists of a shear flow, which corresponds on the gravity side to a hydrodynamic shear quasinormal mode of the black brane. Due to the presence of viscosity, the shear flow is expected to decay exponentially in the absence of turbulence, until the fluid reaches an equilibrium state, corresponding to a (uniformly boosted) black brane.
We present our results in Sec. IV. Our simulations confirm that turbulent behavior and the inverse cascade continue to manifest beyond a critical Reynolds number R c , which we determine numerically. For R > R c , the background decaying shear flow is linearly unstable to perturbations. Such perturbations can grow until they reach the amplitude of the background shear flow, at which point fully developed turbulence is attained. As in the inviscid case [9], an inverse cascade of energy is observed, eventually leaving two large counter-rotating vortices. On the other hand, for R < R c , the shear flow is stable to perturbations, and it decays exponentially.
Finally, in Sec. V (as well as Appendices A and B) we use the gravity/fluid correspondence to relate our results for the fluid to the AAdS black brane. The case of laminar shear flow corresponds directly to the ordinary decay of the hydrodynamic shear quasinormal mode of the black brane. However, for R > R c , the instability of the fluid flow corresponds to an instability of the quasinormal mode. (We stress that this does not imply an instability of the black brane, since an overall decay continues to occur.) Once the growing mode becomes of order the original quasinormal mode, the original decay is interrupted and the overall behavior is strongly modified by a fully developed turbulent behavior. In the 4-dimensional bulk, the energy cascades to the longest wavelength that fits within our torus 3 .
We conclude that ordinary perturbation theory about the uniform black brane background is not the most suitable method of analysis for capturing such effects analytically. In fact, the instability is only clearly apparent if one linearizes the Einstein equation about the decaying quasinormal mode solution itself. (Perturbation theory about the uniform black brane would have to be implemented to higher orders before the exponential growth could be recognized.) Physically, the reason for this behavior is that, for high-temperature black branes in AdS, the lowest lying quasinormal modes become very long lived. Thus, for a given perturbation, as the temperature is increased, the linear viscous damping term becomes small compared with nonlinear terms in the Einstein equation. The regime of applicability of linear perturbation theory is thus pushed to very small metric perturbations. On the fluid side, such properties are conveniently captured by the Reynolds number (although, as noted above, a relativistic generalization is desirable for relativistic fluids). Thus it would be very interesting to obtain a geometrical realization of the Reynolds number, in order to predict the onset of turbulence in gravity [19].
More generally, the unstable nature of certain long-lived quasinormal modes suggests that the decay of a sufficiently perturbed black brane can deviate from the picture suggested by ordinary perturbation theory. Rather than being describable by quasinormal decay, the black brane can undergo a turbulent cascade with a power law decay. Only after the energy cascades to long wavelengths will a quasinormal mode decay take hold.
In this work, we follow all notation and sign conventions of [20]. We use lower case Greek letters (µ, ν = 0, 1, . . . , d−1) for indices of boundary quantities, and we upper case Latin letters (M, N = 0, 1, . . . , d) in the bulk. Boundary indices are raised and lowered with the boundary Minkowski metric η µν .

II. GRAVITY/FLUID CORRESPONDENCE
In this section we review the basic results of the gravity/fluid correspondence. We sketch the derivation from Einstein's equation in the bulk. We also discuss issues concerning the well-posedness of viscous relativistic fluids, and we write down suitable equations of motion that will be used in our simulations [3]. The derivation which we review below follows that of Bhattacharyya et al [4].
As noted in the introduction, we restrict to boundary fluids in Minkowski spacetime, which are dual to perturbed AAdS black branes. Our simulations adopt d = 3, but in this section we keep d arbitrary. We also take the boundary manifold to be R × T 2 ; that is, we impose periodic boundary conditions along boundary spatial directions. Results in d = 3 were derived in Ref. [6], while the arbitrary d case whose equations we write down was analyzed in Refs. [5,21].
The starting point for the derivation of the gravity/fluid correspondence is a uniform boosted black brane spacetime, which, written in ingoing Eddington-Finkelstein coordinates, reads, (2.1) Here, the fields b and u µ (satisfying u µ u µ = −1) are constants. This is a solution to the bulk Einstein equation, with the cosmological constant Λ = −d(d − 1)/2. The boosted black brane is related to the static black brane by a coordinate transformation. The coordinates x µ = (t, x, y) are to be thought of as "boundary" coordinates, while the coordinate r is the "bulk" radial coordinate.
To each asymptotically AdS bulk solution, is associated a metric and conserved stress-energy tensor on the timelike boundary of the spacetime at r → ∞ (see, e.g., Ref. [22] or Appendix A for the precise definition). The boundary metric, in the case of (2.1) is η µν , while the boundary stress-energy tensor is This is a fluid stress-energy tensor, so one may read off the energy density and pressure, (2.5) The stress-energy tensor is traceless, with equation of state as required by conformal invariance. Imposing the first law of thermodynamics, dρ = T ds, as well as the relation ρ + P = sT , gives the entropy density s and fluid temperature T , Here, A is a constant of integration. This is fixed to A ≡ (4π) d /(16πG d+1 d d−1 ) by equating T with the Hawking temperature T Hawking = d/(4πb) of the black brane 4 . To move beyond the uniform fluid, b and u µ are promoted to functions of the boundary coordinates x µ , which are slowly varying; that is, if L is the length scale of variation of these fields, then L b. At this point, the metric (2.1) is no longer a solution to Einstein's equation. However, due to the fact that the fields are slowly varying, it is possible to systematically correct the metric order by order in a derivative expansion, so that Einstein's equation is solved to any given order in derivatives. One can then compute the boundary stress-energy tensor corresponding to the metric at each order, and take this as defining the boundary fluid.
In this setup, the boundary metric is fixed to η µν throughout. This can be thought of as a "Dirichlet condition" on the boundary. The Einstein equation reduces to a set of "constraints" along the timelike boundary of AdS, as well as evolution equations into the bulk. The "momentum constraint" gives rise to conservation of boundary stress-energy, while the "Hamiltonian constraint" ensures tracelessness. The "evolution equations" reduce to ordinary differential equations along r. Regularity at the future black brane horizon is imposed as one of the boundary conditions for these ODEs. This corresponds to the imposition of an ingoing boundary condition, and it is responsible for the breaking of time-reversal symmetry inherent in the fluid viscosity term. The "Landau frame" gauge condition, is also imposed. After a rather long, but direct, calculation, the resulting boundary stress-energy tensor -to second order in derivatives -is found to be The shear and vorticity tensors are given by The angled brackets denote the symmetric traceless part of the projection orthogonal to u µ , while P µν is the spatial projector orthogonal to u µ , It may be verified that Π µν is symmetric and satisfies The transport coefficients {η, τ Π , λ i } have been worked out explicitly in various dimensions [5,6,21], Thus, the boundary fluid has a nonzero shear viscosity η, but the bulk viscosity vanishes, so that the stress-energy tensor remains traceless. Conformal invariance can also be used to deduce the presence of the particular nonzero second order transport coefficients directly [3]. For completeness, we note that conservation of the boundary stress-energy tensor leads to the equations of motion for ρ and u µ , It is easy to see that in this derivative expansion, Π µν is subleading in b/L, as compared with the perfect fluid stress-energy tensor T (0) µν . Thus, given a fixed L, the viscous part may be neglected for small b, or, equivalently, large T . This is the limit which was taken in Ref. [9]. In this work, however, we wish to move beyond the T → ∞ limit, so viscosity must be included in our simulations. As explained in the introduction, this corresponds on the gravity side to the effects of energy losses through the horizon.
At this point, one may wonder why we have bothered to include terms to second order in derivatives, since the shear viscosity appears at first order. The reason is that relativistic viscous fluid formulations which are first order in derivatives, as originally laid out by Eckart [23], lead to acausal propagation, and are generally ill-posed [16]. It turns out to be possible to resolve these issues and to produce a hyperbolic system by including second order terms; in particular, the term involving τ Π [16][17][18]24]. The second order terms which appear in Eq. (2.11) resolve the problems introduced by the viscosity, but they bring about analogous issues at higher order. To fully resolve these difficulties, it is necessary to promote Π µν to an independent field. Then, one reduces the order of the system of equations by substituting −2ησ µν → Π µν on the second line of (2.11). This substitution is consistent to the order to which we are working in the derivative expansion. Furthermore, this assumption will remain valid in 2 + 1 dimensions under time evolution by virtue of the expected inverse cascade [9]. Therefore, following Baier et al [3], we obtain The formulation we have described above also includes additional second order terms with coefficients {λ i }. We have decided to include these in the interest of completeness, although we find that they have no effect on our results. Indeed, as discussed by Geroch [25,26], all hyperbolic relativistic theories of fluids with viscosity should be physically equivalent. By this, one means that any additional terms in the equations of motion, when evaluated within the domain of applicability of the theory, should be small, as compared with the lower order terms. That is, if the higher order terms became important, then there would be no justification in not including even higher order terms, and the perturbative expansion would break down. This also means that the specific value of τ Π is unimportant, so long as it is sufficiently large that the theory is causal. (We will use this fact later to increase its value in order to speed up our numerical simulations.) To summarize, the system of interest is described by Eqs. (2.23), (2.24), and (2.25). These equations, however, require further manipulation prior to numerical implementation. For comparison, we recall that in the context of inviscid hydrodynamics, it is covenient to express the hydrodynamic equations in conservation form (see, e.g., [27]). That is, t-derivatives of energy and momentum density are equated with x i -derivatives of fluxes. Such a form of the equations is particularly advantageous when studying solutions that can develop sharp gradients or discontinuities, and was employed in Ref. [9]. However, a simple extension of this approach is not possible in the presence of viscosity since the evolution equation (2.25) for Π µν is not of the desired form. In particular, since P µν projects orthogonally to u µ rather than to ∂ µ t , this equation contains a t-derivative of u µ in addition to that of Π µν . However, since the presence of viscocity prevents the development of steep gradients and discontinuities in our solutions, adopting a conservative form is not necessary.
We therefore follow an approach similar to that employed in [28,29] within the context of heavy ion collisions. To begin, we note that in d = 3 the conditions u µ u µ = −1, u ν Π µν = 0, and Π µ µ = 0 reduce the number of dynamical variables to five. We take these to be U ≡ (ρ, u x , u y , Π xx , Π xy ). Equations (2.23), (2.24), and (2.25) are quite complicated when expressed in terms of U, and in order to evolve the equations numerically, one must solve for the t-derivatives of the fields. Using computational algebra software, we can write our equations in the desired form, and these are the equations we implement in our code.

III. SIMULATIONS
In this section we describe our choice of initial data and details of our numerical setup.
A. Initial data As described in the introduction, we choose initial data corresponding to a shear hydrodynamic quasinormal mode of the black brane. Our studies concentrate on nonlinear phenomena described by the system; however, for future reference, in this subsection we analyse the evolution under the linearized equations of motion.
Consider perturbations about a uniform fluid solution, which is dual to a non-boosted uniform black brane in the bulk. Solutions to the linearized equations of motion whose only nonzero perturbed fields are u x (1) = u x (1) (t, y) and Π (1) xy (t, y) describe shear flow. That is, fluid flow orthogonal to the velocity gradient [3]. The linearized equations of motion for shear flow reduce to This has two solutions for small k, both of which describe pure exponential decay. The first solution corresponds, in the bulk, to the hydrodynamic shear quasinormal mode of the black brane [3,4,15]. (The second solution shows that τ Π is the decay timescale for Π µν to approach −2ησ µν .) Thus, since we are interested in understanding the corresponding black brane quasinormal mode in a nonlinear context, we choose initial data with and all other fields zero 5 . We vary the background energy density ρ 0 , velocity amplitude v 0 , the number of modes n, and the torus size [0, D] 2 . The reason n and D are varied separately (rather than as the wavelength λ = D/n) is that effects due to the finite size of the box can come into play for small n.
Moving from the linear to nonlinear level, we expect the pure decay of this shear flow to persist, at least for small velocities. That this is the case will be verified in Sec. IV. We also keep the velocities small in most of our simulations in order to match to the linear predictions.
Our main reason for setting up this flow, of course, is to study its stability. In order to do so, we also initially seed u x with a very small random perturbation 6 . By studying the effects of this perturbation, one can learn about the robustness (or lack thereof) of pure quasinormal mode decay.
We note that the presence of the short viscous timescale τ Π [see Eq. (3.7) in the previous subsection] imposes a harsh constraint on the timestep length for an explicit integration method, ∆t ∝ τ Π . (3.10) For our simulations, this is a much stronger constraint than that arising from the finite propagation speeds of the solution (i.e., the CFL condition). However, as we discussed in Sec. II, the precise value of τ Π should not have physical significance, as long as the equations of motion remain hyperbolic. Therefore, to allow for a more efficient numerical integration, we increased τ Π by a factor of 100 for many of our runs. We verified that this had no significant effects on any of the physical properties we measured.

IV. RESULTS AND ANALYSIS
In this section we present and analyze results. We first define the Reynolds number for the shear flow, which we find to accurately predict the onset of instability. We then describe the three observed "phases" of a fully developed turbulent flow: initial growth of instabilities, inverse turbulent energy cascade, and final exponential decay (see Fig. 1 for a preview). Our results are largely consistent with expectations drawn from solutions to the Navier-Stokes equations in 2 + 1 dimensions.

A. Reynolds number
For steady flows, the Reynolds number is a useful dimensionless quantity which can be used to predict stability (see, e.g., [32]). It is generally true that for sufficiently low Reynolds numbers, the flow is stable with respect to small perturbations, in which case it is said to be laminar. In contrast, for large Reynolds numbers the flow is unstable, which eventually leads to turbulence. The critical Reynolds number which separates these two regimes depends upon the particular flow under consideration.
The shear flow that we study is not steady, causing the Reynolds number to change with time.
[The shear viscosity causes it to decay according to Eq. (3.7).] Nevertheless, it is a useful quantity to consider, as the decay can be treated in a quasi-stationary manner 7 . In this case, the stability properties of the flow depend only upon the instantaneous value of the Reynolds number.
For our flow, we define the Reynolds number to be 8 Substituting for η and λ, Thus initially, and with time, R decays with max(u x ). (ρ, D and n are all either constants, or nearly constant.) Later in this section, we will verify that there exists a critical Reynolds number R c , and we will determine its value. Strictly speaking, flows which have different values of n are not geometrically similar, meaning that they are not related by a universal scaling of distances. This is due to the presence of two independent associated length scales (the wavelength λ, and box size D). So, one should exercise caution when comparing the Reynolds numbers of two such flows. However, for large values of n, the finite box size should not play an important role in governing stability, and our definition (4.1) makes sense. (We will address effects at small n later in this section.) Our definition of R, and our decision to compare flows at different n, is motivated both by simplicity, as well as the desire to address the infinite brane limit (n, D → ∞ while holding λ fixed).
We note that, for fluids dual to black holes with compact spatial sections, one should also be careful when using a value of R c for high angular quantum number flows, to predict stability of low-l flows. In particular, as we will discuss further, we would expect l = 2 shear modes to be stable for any value of R.

B. Stability of shear flow
In this subsection we analyze the early stages of the flow, where the overall properties are governed by the shear decay in u x .
Recall that, in addition to the shearing configuration, the initial data is seeded with a small-amplitude random perturbation. This has the potential to grow or decay, depending on the Reynolds number of the flow. To track the presence of growing instabilities, we monitored the u y field. In the linear analysis of Sec. III A, u y remains exactly zero, so its growth reflects the growing unstable mode. (In addition, even for stable flows, u y becomes nonzero due to nonlinearities; but this remains small.) As expected, as we varied the initial data we found that the various solutions could be categorized into several groups, based on the growth of u y 2 .
At one extreme, the typical turbulent run is illustrated in Fig. 1. This shows the vorticity field at several times. In Fig. 2, we show the corresponding L 2 norms of the vorticity and the components of u i . We see that initially, the vorticity and u x decay exponentially in the manner expected for shear flow. However, during that time, u y 2 undergoes exponential growth, until it reaches the same amplitude as u x 2 . This brings the solution into an equipartition of energy between u x and u y . At this point the initial decay has been disrupted, and turbulence sets in, as we will describe in more detail in the following subsections. The onset of this behavior is known as the Kelvin-Helmholtz instability in fluid dynamics.
Visual inspection of the u y field (see Fig. 3) indicates that the growing mode itself is also very roughly a shear mode, typically with n/2 < m < n. As we have alluded to earlier, the finite box size plays a role at small n, and it comes into play here. For example, we find that the n = 1 case is stable, even for large R, because the box does not admit a mode with m < 1. This is the expected behavior as there is no room for an inverse cascade, given an n = 1 initial configuration. For large n, there is no obstacle in fitting the growing mode into the box, and the box size D plays no role. Thus, fixing λ and extrapolating to the infinite box (D → ∞), the instability should be present for sufficiently large R.
At the other extreme, at low values of R, the flow is laminar. This is illustrated in Fig. 4. In contrast to the turbulent case, u y 2 remains small throughout the run. (u y is being continuously driven nonlinearly by u x , so its amplitude decays with u x .) With small initial velocities, the linear analysis of Sec. III A is applicable, and the measured decay rate should be consistent with the prediction (3.7). Fig. 4a shows that this is indeed the case.
Between these two extremes, we found that there were certain intermediate flows which provided important physical information. Such an example is illustrated in Fig. 5. In this case, the flow begins at high Reynolds number but it decays before turbulence can fully develop. The plot of u y 2 shows initially small perturbations growing nearly exponentially for some time -as in the turbulent case -before peaking, then decaying exponentially.
The time dependence of R in such runs is clearly evident, as it is directly proportional to max(u x ) (also plotted). Initially, R > R c , but as time progresses, R decreases. Since the background flow is slowly varying, we assume that the instantaneous growth rate of u y 2 depends only on the background value of max(u x ). Thus, at the peak of u y 2 , R = R c , while for R < R c , u y 2 decays. This allows us to extract R c . Indeed in Fig. 5b we plot the growth rate of u y 2 versus max(u x ). Where the curve crosses through zero corresponds to the peak of u y 2 in Fig. 5a, and the value of R c may be read off. Thus, such intermediate runs provide detailed information on the stability of the background flow.
By searching for such critical runs (by adjusting v 0 ) for different values of (ρ 0 , D, n), we found that the critical Reynolds number of the shear flow is   Fig. 1). In (a) there is an initial exponential decay, followed by a power law during the inverse cascade, followed by a slower exponential decay. In (b) uy grows exponentially until it is of similar amplitude to ux.
FIG. 3. uy field at t = 300 for the same run as Fig. 1.
The dependence upon n (see also Fig. 6) is due to the finite size of the box, as discussed above. We measured this within the range 6 ≤ n ≤ 20. (We verified, by varying also ρ 0 and D, that R c is independent of these quantities.) For large n, we expect R c to approach a limit, as the finite box size should play no role. Fig. 5b also provides some information about the growth rate of u y away from its zero value at R = R c . Fixing n and ρ 0 , this shows that for R > R c , higher R [higher max(u x )] gives faster growth. For R R c , the growth rate is linear in max(u x ). In contrast, as R is lowered, there is a bound on the decay rate. This can be understood as a competition between a driving effect from the background shear flow in u x , and a viscous decay (3.7) associated with the shear mode (4.4) in u y . Once the driving term drops to the point of irrelevancy, all that is left is the decay term, which gives a fixed decay rate, as it is a property of the mode in u y . Indeed, the asymptotic value of the growth rate of u y 2 in Fig. 5b is -0.0013. This matches very nicely the prediction from Eq. (3.7) using the observed m = 6.
We also found that for large R, the growth rate of u y 2 is inversely proportional to λ, at fixed max(u x ). Together with the above, this points to the contribution of the driving term to the growth rate as being proportional to ∂ y u x .

C. Turbulent regime
We now turn our attention back to the turbulent case, at the point where equipartition of energy is reached between u y and u x (i.e., just subsequent to Fig. 1b). As seen in Fig. 1c, the overall flow is completely disrupted, and the vorticity field displays a number of turbulent eddies. The exponential decay of vorticity in Fig. 2a ceases, and is replaced with a power law decay. We typically observed ω 2 ∝ t −α , with α 1.2 ± 0.2. There have been several previous studies of unforced turbulent decay of Navier-Stokes fluids in d = 3, which have also found power law decays (e.g., [33][34][35]).
During the power law decay, the eddies merge into vortices, which continue to merge into increasingly large vortices (Figs. 1d and 1e). This inverse energy cascade can be attributed to the conservation of enstrophy [11]. Co-rotating vortices merge, while counter-rotating vorticies repel. Thus, one is finally left with two counter-rotating vortices, as  in the inviscid case [9].

D. Late time decay
The vortices which form are the relativistic analog of the Oseen vortex, which is an attractor solution to the Navier-Stokes equation [36]. Its functional form is where the parameter C 2 (t) = √ 4νt, and the kinematic viscosity ν = η/ρ. A fit for the parameters C 1 and C 2 is shown in Fig. 7. We see that the vortex is close to, but does not match exactly, the Oseen profile. We attribute this to differences betwen Navier-Stokes fluids, and the relativistic compressible fluids we study 9 . Fitting also in time, it is possible to extract ν from these profiles. As with the profile fit, we found a value of ν which, while not quite the predicted value, was within a factor of two.
At late times, the two vortices continue to decay in amplitude, until the fluid finally becomes linear. The solution is then the sum of long-wavelength shear and sound modes, which decay exponentially. We measured the decay rates at late times. The decay was slightly faster than the linear prediction, although the measured difference decreases at later times.
As a result of the increased wavelength, the decay rate at late times is lower than the decay rate of the initial flow. Thus, the presence of the turbulent cascade drastically lengthens the time period before the fluid settles down from its initial state to a uniform flow. In contrast, for d > 3, where we expect a direct cascade, turbulence causes a more rapid decay than the linear behavior. This is due to higher modes decaying faster and strong dissipation at the viscous scale.

V. DISCUSSION: BLACK BRANES AND TURBULENCE
In the previous section, we have established conditions for the onset of turbulent phenomena in conformal, viscous relativistic fluids in 2 + 1 dimensions, as well as the subsequent behavior once this develops. Having studied solutions to the dual hydrodynamic theory, we turn in this section to general relativity, and to the behavior of perturbed AAdS black holes and black branes.

A. Decay of perturbations and turbulence in the bulk
The decay properties of the shear fluid flow that we have analyzed in Sec. IV carry over directly to the shear hydrodynamic quasinormal mode of the black brane. Indeed, as described in the introduction, the gravity/fluid correspondence naturally captures the behavior of the lowest lying shear and sound families of quasinormal modes. (As illustrated in Ref. [12], the higher order quasinormal modes typically decay very rapidly, and the metric produced via the duality is a very good approximation to a solution of Einstein's equation.) Translating to the black brane language, our results imply that, for R > R c , hydrodynamic shear quasinormal modes are unstable to small perturbations (the "instability" refers to the quasinormal mode; not to the black brane). More precisely, certain deviations from the pure quasinormal mode undergo exponential growth until either they reach the amplitude of the quasinormal mode (fully developed turbulence) or the Reynolds number -which decays in time -becomes smaller than R c and an exponential decay ensues. Once turbulence sets in, for 4-dimensional AAdS black branes, energy is transferred to longer wavelength modes and a polynomial decay is induced. Eventually, when metric deviations about the uniform black brane become small enough, exponential decay resumes. In both cases, the final decay is at a slower rate than the original decay, as the perturbation is of longer wavelength.
On the other hand, for R < R c , the quasinormal mode is stable, so it exhibits the usual clean exponential decay. We stress that in all cases described, the global norm of the solution decays in time.
More generally, one is interested in the behavior of a generic black brane perturbation, containing many modes of small but comparable magnitude. The standard picture states that if the amplitude of the perturbation is small enough, then at sufficiently late times it asymptotically approaches 10 a sum of quasinormal modes, which evolve in time independently. However, the question is how small the amplitude must be for this to be realized; this is determined by the Reynolds number. At high Reynolds number, certain quasinormal modes are unstable; and an unstable mode will never be realized in a decay. The new picture that emerges is that of both laminar and turbulent phases. The laminar phase corresponds to the standard quasinormal mode picture. At high Reynolds number, however, the decaying perturbation immediately enters the turbulent phase, which has been uncovered through the gravity/fluid correspondence.
This new turbulent phase displays a far richer phenomenology. Our results indicate that turbulence, when it develops, induces eddies. Eddies with vorticity of the same sign merge, leading to increasingly large vortices as time proceeds [39]. The form of the bulk metric [to leading order, Eq. (2.1)], indicates that the boundary fluid structure is carried unperturbed along ingoing null geodesic "tubes" [5], connecting the boundary at spatial infinity to the black brane event horizon. Thus the Oseen-like vortices present at late times in the turbulent fluid solution describe rather compact distributions of gravitational radiation connecting the asymptotic and black hole regions. They can be regarded as natural realizations of "extended geons", which extend through the bulk as "gravitational tornadoes" or "funnels" 11 . The structure that we describe is illustrated in Fig. 8, where we plot the curvature invariants I 1 and I 2 (see Appendix B). As discussed in Appendix A, a possible way to understand this behavior is related to the fact that a map to relativistic hydrodynamics is also possible away from the boundary, at suitably defined timelike surfaces in the bulk. Thus, the solution can also be analyzed at these surfaces to show, in particular, that a conserved quantity related to the enstrophy can be defined away from the boundary. Such a quantity is a key element needed to argue that an inverse energy cascade occurs.
Naturally, as already pointed out in [19], it would be very useful to develop a spacetime definition of the Reynolds number. This would provide an intrinsic way to predict the onset of turbulence in gravity and could thus be applied in broader contexts. Using the gravity/fluid correspondence, this would also lead to a Reynolds number suitable for relativistic hydrodynamics. Based upon the form of the bulk metric (2.1), and the fluid Reynolds number (4.1), the form of the Reynolds number for black hole perturbations is, roughly, where we have substituted (certain components of) the metric perturbation h AB ≡ g AB − g AB for the velocity fluctuation; and L is the characteristic length scale of the perturbation. Of course, whether or not this is applicable in more general contexts would require further investigation. In particular, a suitable definition of R should be gauge invariant.
As a final comment, we point out an important application of the fact that the inverse cascade guarantees the system stays within the domain of validity of the gravity/fluid correspondence. The relativistic hydrodynamic equations in 2 + 1 dimensions are dual to (long-wavelength) perturbed black branes in the bulk. Thus, a Newtonian limit in the bulk -where time derivatives are taken to be one-order subleading to space derivatives and velocities are taken to also be small 12 -corresponds to a Navier-Stokes limit on the boundary. Since we know that the Navier-Stokes equation admits global, well-behaved solutions [43], one can surmise that general relativity is similarly well-behaved in the bulk. 11 These should not be confused with "black funnels" [40], which are bulk black holes with a horizon that connects to the conformal boundary of the spacetime. Evslin [19] has conjectured these black funnels to correspond to singular fluid vortices. 12 For some alternative discussions, see, e.g., [41,42].

B. Connection to ordinary perturbation theory
Linear perturbation theory predicts that small-amplitude metric perturbations can be decomposed (asymptotically) into independent modes which undergo simple exponential decay. On the other hand, the picture described in the previous section indicates the presence of a qualitatively distinct, turbulent, behavior for sufficiently high Reynolds numbers, regardless of the perturbation amplitude. How are these two notions reconciled? The short answer is that at very high black hole temperatures (in AdS), the regime of applicability of linear perturbation theory is very small. To study this further, we take a closer look at perturbation theory.
In ordinary perturbation theory, the full metric is expanded as where g (0) AB is taken to be the background metric; in our case, the uniform AdS black brane. The first order metric perturbation h (1) AB satisfies the homogeneous partial differential equation, , h (1) ) is the linearized Einstein tensor. For the black brane, the symmetry properties admit a mode decomposition, and by solving Eq. (5.3), the quasinormal mode spectrum is determined. All of the modes decay for the AdS black brane [13,14].
At second order in perturbation theory, where the second order Einstein tensor on the right hand side is quadratic in the first order perturbation h AB . Since the homogeneous part of this equation is unchanged from the first order case, the quasinormal mode spectrum of h (2) AB is also unchanged. These decaying modes are excited by the inhomogeneous source term. At any finite order in perturbation theory, the same applies. Thus, the growth we describe can be only captured by carrying out the analysis to sufficiently high orders in perturbation theory to recognize the underlying exponential behavior.
If, rather than taking the background metric g (0) AB to be the uniform AdS black brane, it is instead taken to be the AdS black brane plus the shear hydrodynamic quasinormal mode, then the growth is easily seen to be possible at the linearized level (5.3). This is best illustrated through a simple toy model of ordinary differential equations (inspired by a local analysis of the Navier-Stokes equations) which exhibits similar mode-coupling properties, namely dx dt + αx = 0, (5.5) dy dt + βy − γxy = 0, (5.6) with α, β, and γ all positive constants. The variable x(t) in this system corresponds to u y in the black brane system, which initially describes a shear mode. On the other hand, y(t) corresponds to the initially zero u y . Both x and y are subject to dissipation (in the black brane case, α ≈ β), and we have included a mode-coupling γxy in (5.6). Let us solve the system (5.5)-(5.6) perturbatively in two different ways 13 . We expand both x and y as The exact solution to (5.5) is x ∝ e −αt . In a manner analogous to perturbation theory about the uniform black brane, we take the "background solution" to be x (0) = y (0) = 0, while the "linearized solution" corresponds to the quasinormal mode, i.e., x (1) = a 1 e −αt . Then, the equation for y (1) is with solution y (1) = b 1 e −βt . At second order, The solution to this is y (2) = −a 1 b 1 γe −(α+β)t /α + b 2 e −βt ; higher order corrections can be computed in a similar manner. At each order in perturbation theory the solution is a sum of decaying exponentials. Suppose instead that we take x (0) = a 0 e −αt and y (0) = 0. This is analogous to taking the black brane perturbed by a quasinormal mode as the "background solution". Then, at linear order we have dy (1) dt which has solution y (1) = b 1 exp [−βt + a 0 γ (1 − e −αt ) /α]. In the limit of small α, this becomes Thus, if the "background" is long-lived, the growth rate of y depends upon a competition between the driving term γxy and the dissipative term βy. For finite α > 0, the driving term decreases with time relative to the dissipative term. In this case, y (1) is eventually dominated by the exponential decay (cf. Fig. 5). One can define a "Reynolds number" of the flow x(t) in the toy model, by taking the ratio of the mode-coupling term γxy, to the linear dissipative term, βy, in (5.6). If this "Reynolds number", γa 0 e −αt /β, is large, then the nonlinear term γxy should be kept in any perturbative analysis, and one should solve Eq. (5.11) -rather than Eq. (5.9) -to determine y.
A similar story holds in for black brane perturbations in general relativity. The problem with trying to analyze high Reynolds number perturbations by performing ordinary perturbation theory about the uniform black brane background, is that this drops certain large nonlinear terms while keeping small linear terms. For the black brane, the dissipation rate (the analog of α ≈ β) depends inversely on the temperature. So, increasing the temperature, while holding fixed the amplitude and wavelength of a perturbation (i.e., increasing the Reynolds number), the linear term evantually becomes small relative to nonlinear terms. Thus, the regime of validity of ordinary perturbation theory is reduced as the Reynolds number is increased.

C. Beyond 4-d AAdS black branes
Here we discuss some possible implications and extensions of this work.

Higher dimensions
Based upon our results for 4-dimensional bulk spacetimes, together with the gravity/fluid correspondence established in arbitrary dimensions, as well as numerical results confirming the expectation of direct energy cascades for inviscid conformal relativistic fluids [10], it is possible to anticipate properties of 5 (and higher) dimensional spacetimes (for an early discussion, see [6]). Three immediate consequences are: • First, as with the 4-dimensional case, at sufficiently high Reynolds number, the quasinormal mode description fails to accurately describe the decay of black brane perturbations.
• Second, in contrast to the 4-dimensional case, turbulence is characterized by a direct energy cascade to short wavelengths. Since shorter wavelengths have a more rapid decay rate, thermalization in 5 dimensions will be attained in a shorter timescale than in an analog 4-dimensional case. Notice that a potential concern here is that the cascade to short wavelengths may cause the solution to exit the regime of validity of the gravity/fluid correspondence. Two comments are relevant here. (i) Even if this were the case, perturbations initially satisfying LT 1 (as required by the correspondence) still undergo turbulent dynamics which induce structure at shorter wavelengths. Thus, at the moment of "exiting" the regime of validity, many modes would be present and their subsequent behavior should be studied within general relativity. (ii) It is also possible for the cascade to occur completely within the regime of validity of the correspondence if the viscous scale (again defined as in the non-relativistic case) L η = (η 3 / ) 1/4 (where is the rate of energy dissipation by viscosity) satisfies L η T 1. Notice that L η grows with temperature, so at sufficiently high T , this condition is satisfied. In this case energy would be expected to cascade down to the viscous scale, and then dissipate.
As a consequence of the turbulent behavior on the hydrodynamical side -which displays a self-similar behavior -one would expect a similar (fractal) structure for black hole perturbations on the gravitational side 14 . This structure is expected to smooth out in time yielding a slightly hotter, uniform black brane. If this is the case, then in this high temperature limit [case (ii)], a global solution to the dual relativistic hydrodynamics problem would seem to be a natural consequence. As in the d = 2 + 1 fluid case, this would have obvious implications for global solutions to the Navier-Stokes equation. Establishing this decay result rigorously on the gravity side amounts to determining nonlinear stability 15 of large black holes in AdS.
• Finally, a word of caution with regard to nonlinear numerical studies of gravitational perturbations in 5 dimensions: Due to the high computational cost of such simulations, symmetries are usually imposed on the solution, effectively reducing the number of dimensions which are actually simulated. However, this may restrict or affect the development of turbulence. In particular, 5-dimensional spacetimes which are dimensionally reduced to 4 dimensions will give rise to an inverse energy cascade -rather than a direct one -while 3-dimensional treatments will eliminate turbulence altogether. As one is often interested in using the AdS/CFT correspondence to describe high temperature CFTs through a gravitational analysis, it is important to bear in mind that the imposition of computation-saving symmetries can impact the extracted physics.

Black holes
We have investigated the decay of black branes in a Poincaré patch of AdS, with torus topology. However, it is also of interest to study black holes in global AdS. Based on earlier work in the inviscid case [9], we expect a qualitatively similar turbulent behavior. One primary difference, though, is the final number of vortices which remain (e.g., fluids dual to Kerr-Schwarzschild settle down to two clockwise and two counter-clockwise rotating vortices, while Kerr-AdS settles to just one of each sign). Other particular details, such as the power law decay rate during turbulence, and the critical Reynolds number, may also differ. The power law behavior, ω 2 ∝ t −α , can also be estimated in the inviscid case. An examination of results presented in [9] indicates that for black holes, the decaying behavior is realized with a similar exponent to black branes, in the range 0.5 < α < 1.5.

Beyond AdS
At a speculative but certainly tantalizing level, that gravity displays turbulent behavior in AAdS spacetimes suggests that more general asymptotic conditions should also be investigated. Is this a special property of AdS, or could it arise in the asymptotically flat or de Sitter cases? What about Dirichlet boundary conditions, but without a cosmological constant? There are two elements to consider: The boundary conditions imposed on the solutions to the Einstein equation, and the presence of the cosmological constant in the equation of motion. From a partial differential equations point of view the cosmological constant introduces a lower order term in the equations, which does not affect local propagation.
On the other hand, it is well known that linearized perturbations of AAdS black hole spacetimes have a family of very slowly decaying modes (the hydrodynamic modes), which (as we have shown) play a key role in terms of being unstable to perturbations. Such modes are also present in certain vacuum solutions bounded by an accelerating mirror, which have been shown to be dual to Navier-Stokes fluids [42]. Interestingly, at least some hydrodynamic modes can also be connected to QNMs of asymptotically flat spacetimes [47]. However, in this regime, the modes decay rapidly, and therefore do not govern the long term behavior of the system. However, they could play a role in channeling energy in the transient stages. Furthermore, it is well known that massive fields introduce effective boundaries which could induce the long-lived mode behavior, and thus allow for turbulent-like phenomena. As discussed earlier, if a Reynolds number can be suitably defined for gravity, it would help predict the onset of turbulence in these varied scenarios.

D. Final words
As we have stressed throughout this work, the gravity/fluid correspondence translates intuition of fluid behavior into the realm of gravity. This has allowed us to identify key features of the behavior of perturbed AAdS black holes, including a new dynamical phase where fully developed "turbulence" gives rise to a polynomial decay. Furthermore, turbulent behavior also indicates that perturbed black holes can behave in a strongly dimensionally-dependent way, both qualitatively and quantitatively. Establishing that gravity can behave in a turbulent manner opens new doors to searching for other situations where this can take place. For instance, it provides motivation to look for scenarios where slowly decaying perturbations might give rise to interesting non-linear interactions. Finally, insights from turbulence may shed new light on particular systems known to exhibit related behavior, such as the chaotic behavior of spacelike singularities in early-universe mixmaster dynamics [19].
where the exponent contains an indefinite integral, and we assumed the integration can be performed for a given equation of state P (ρ).
Notice that Ω µν , is strongly conserved by the flow by virtue of Cartan's identity. That is, Motivated by this observation, we can construct a current of the form J µ Z ≡ λΩ 2 u µ , where we have denoted Ω αβ Ω αβ ≡ Ω 2 , and λ a scalar field to be fixed by the conservation requirement. Computing the divergence of J µ one obtains, We have used the definition of the Lie derivative and (A19) on the second line, and the third line follows from, an identity which is valid in two spatial dimensions [9]. Clearly, the divergence in Eq. (A20) vanishes if λ(ρ) ∝ exp − dρ ρ +P (ρ ) .
Notice also that Ω 2 is related to the square vorticity (ω 2 ≡ ω µν ω µν ) by Ω 2 = 4ρ 2 ω 2 . Finally, by integrating the current J µ Z over a constant-t hypersurface Σ t , the expression for the enstrophy takes the form It is straightforward to check that this expression reduces to the one reported previously in [9] for the conformal fluid on the boundary, when the equation of state is given by P (ρ) = 1 2 ρ.
We are now in a position to evaluate the enstrophy on an arbitrary r = constant timelike hypersurface in terms of boundary variables, while for the square vorticity we find that,ω 2 = 1 α 2 ω 2 .
As a final remark, let us note that to zeroth order in the gradient expansion, the transformation from one surface element to the other is given byû Consequently, the enstrophy calculated at any given radius r, not only is conserved, but is also equivalent (up to an unimportant constant factor which can be trivially absorbed into the definition) to the enstrophy of the boundary fluid. That is, The fact that it is possible to define an enstrophy for the fluid on each slice of the bulk geometry, and that it has the same expression in all cases, is a natural consequence of the ultra-local character of the map. Thus, in a sense, the dynamics occurring at the boundary are reproduced throughout the bulk geometry (in a slightly distorted manner). At first sight, this enstrophy construction in the bulk seems to indicate that one could define an interesting local quantity from the spacetime perspective in the bulk. However, since this quantity is defined on a distorted (and dynamical) surface, one should exercise caution, especially in highly distorted scenarios.
FIG. 9. Contour plots for Ψ1, (Ψ2) and Ψ3 (left to right) illustrating the vortex structure from the horizon to the boundary. Note that the "conical" structure is a result of the dependence upon the radial coordinate of the plotted quantities (see Eq. (B6)). (In contrast to Fig. 8, we are not rescaling here by any powers of r.) FIG. 10. Radial profile of the "enstrophy" at the horizon, for a vortex solution, as calculated by Eq. (A24), the square of the Pontryagin density [12], and [ (Ψ2)] 2 [49] (labelled Z1, Z2 and Z3, respectively). Each curve has been rescaled by a trivial constant factor for easier comparison.
illustrates the radial profile of the three quantities, showing good agreement.