Holographic Turbulence in a Large Number of Dimensions

We consider relativistic hydrodynamics in the limit where the number of spatial dimensions is very large. We show that under certain restrictions, the resulting equations of motion simplify significantly. Holographic theories in a large number of dimensions satisfy the aforementioned restrictions and their dynamics are captured by hydrodynamics with a naturally truncated derivative expansion. Using analytic and numerical techniques we analyze two and three-dimensional turbulent flow of such fluids in various regimes and its relation to geometric data.


Introduction and Summary
Turbulent fluid flow is prevalent in every-day phenomena. Yet, while common, it is difficult to characterize it in a quantitative manner. Recently, the gauge-gravity duality has provided a novel geometric means to study relativistic fluid flow of strongly coupled, large N gauge theories, [1]. Indeed, in a suitable parameter range, and given appropriate initial conditions, such flows have been shown to exhibit (decaying) turbulent behavior [2,3]. This development hints at a fascinating set of connections between turbulence and the dynamics of black hole horizons in general relativity.
A connection between turbulence and event horizons of asymptotically AdS black holes opens the possibility of using geometric tools to study turbulence. Indeed, various pioneering works have suggested a connection between the dynamics of the event horizon and a Kolmogorov-type cascade [2][3][4]. Unfortunately, the structure of time-dependent black hole horizons is not fully understood, which is one of the reasons this research program has not yet come to fruition. In what follows we continue exploring these connections using newly developed tools for studying gravity in a large number of dimensions. More specifically, we consider the turbulent behavior of event horizons of asymptotically AdS spaces, in the limit where the number of dimensions, d, is very large.
There is an extensive, ongoing, program to understand general relativity in the limit where the number of dimensions is large [5][6][7][8][9][10][11][12][13][14][15][16][17][18][19]. The works of [8,13,15,19,20] in particular address the large d limit of asymptotically AdS spaces. To obtain a large d limit, one considers a metric ansatz in d = n + p + 1 spacetime dimensions, where the dynamics depend on p + 1 dimensions. The remaining number of dimensions, n, is taken to be infinite. This allows one to relate the large d limit of gravity to gravity in p + 1 dimensions. Since our goal is to relate gravity at large d to hydrodynamic flow, we obtain the d → ∞ limit through a novel but circuitous path which will eventually allow a direct relation between horizon dynamics and fluid dynamics.
We start our analysis by considering the equations of motion of relativistic (conformal) hydrodynamics in the limit where the number of dimensions is large. Generically, relativistic hydrodynamics does not simplify in such a limit. Nevertheless, as we show in section 2, if we scale the time and space coordinates appropriately, tune the transport coefficients of the theory to lie in an appropriate subset of parameter space, and work in a suitable fluid frame, the hydrodynamic equations significantly simplify. In appropriate variables they are given by (2.18).
Given the simplicity of the hydrodynamic equations in an appropriately scaled coordinate system, we turn our attention to asymptotically AdS gravity in the same limit. We find that our scaling laws for the coordinates are precisely those introduced in [11] to study the Einstein equations in the large d limit. Additionally, the holographic transport coefficients satisfy precisely the relations required for the aforementioned simplification of the large d hydrodynamic equations. Indeed, as we show in section 3 the large d Einstein equations of [11] are precisely the large d hydrodynamic equations of motion, albeit in an unusual, but well motivated, fluid frame.
Since we were compelled to scale our coordinate system, the large d limit of the relativistic hydrodynamic equations of motion bear a striking resemblance to the (non-relativistic) Navier-Stokes equations. More precisely, they are a variant of the compressible version of the Navier-Stokes equations. As such, they may be studied using traditional tools of non relativistic fluid dynamics. As we show in section 4 when the flow is subsonic the Kolmogorov or Kraichnan scaling laws are expected to emerge and an inverse cascade should be exhibited when the system mimics two dimensional (sustained) turbulent fluid flow. In fact, even when the Mach number for the flow is not small we find evidence for an inverse cascade-like behavior.
Following the discussion in section 4 we turn to a numerical analysis of the equations of motion in section 5. By using an appropriate initial condition we generate flow which exhibits decaying turbulence. We then compare various stages of the flow to expectations made in section 4 for two and three spatial dimensions and to an analysis of holographic (decaying) turbulence in AdS 4 [2,3].
Given our improved understanding and analytic control over the relation of the fluid flow to the dynamics of the event horizon (for subsonic flows in particular), we study, in section 6, a proposed relation between the horizon curvature power spectrum and the hydrodynamic energy power spectrum [2]. While these quantities are linearly related in the regime of small Mach number, the relation between these two quantities becomes obscure as the Mach number increases.
Finally, we present a summary of our findings in section 7 where we also discuss possible extensions and puzzles associated with out results.

The Large d Limit of Hydrodynamics
Hydrodynamics is a universal low energy effective description of many-body systems. Following a Wilsonian type of coarse graining, the equations of motion of relativistic hydrodynamics may be characterized by a handful of fields: a unit normalized velocity field u µ (x), a temperature field T (x) and, in the presence of charged matter, a chemical potential µ(x). When working in the limit where gradients of the hydrodynamic fields are small compared to the inverse mean free path, all observables may be expressed as local functions of the hydrodynamic fields. These expressions are referred to as constitutive relations.
In a d spacetime dimensional conformal theory in flat space (and in the absence of charge, anomalies or parity breaking terms) the constitutive relations for the stress tensor take the form and and we have set the speed of light to unity. See [21,22] for details. Following the standard reasoning behind the construction of effective field theories, the constitutive relations (2.2) are the most general ones allowed, compatible with the conformal symmetry of the underlying theory, with a local version of the second law, and taking into account useful features of the derivative expansion. Indeed, in writing down (2.2) we have used the equations of motion at order n − 1 in derivatives to simplify the constitutive relations at order n. In addition, we have used a particular definition of the fluid velocity and energy density, usually referred to as the Landau frame, in order to further simplify the relations (2.2). Since these features of the derivative expansion will be essential in our study of the large d limit of hydrodynamics, we describe them in some detail in what follows.
In the absence of charge, the equations of motion for the hydrodynamic fields are nothing but the conservation equations for the stress tensor which need to be solved perturbatively in a derivative expansion. For instance, at leading order in derivatives one has Therefore, when constructing T µν (1) it is possible to replace u µ ∂ µ and u µ ∂ µ u ν with ∂ µ u µ and P µν ∂ µ . The same reasoning applies to the construction of T µν (2) ; at n'th order in the derivative expansion, the available tensor structures from which the stress tensor may be constructed is reduced due to the equations of motion at lower orders. It is possible to modify the relations (2.2) for T µν (1) and T µν (2) by an appropriate field redefinition of T and u µ . The constitutive relations (2.2) were written in the Landau frame, defined via From a physical standpoint the Landau frame is defined such that u µ points in the direction of the energy flux, and the relation between (T ) and T does not get modified by derivative corrections.
When the system is out of equilibrium, the Landau frame is only one of many choices for possible definitions for T and u µ . Indeed, working perturbatively in the derivative expansion one may always redefine the temperature and velocity field via, e.g., where the c's and k's are real numbers. The equations of motion for the velocity field and temperature, which follow from inserting the constitutive relations (2.2) into the conservation equation (2.4), are lengthy and cumbersome. While the zeroth order equations of motion (2.5) are relatively short and simple, writing the first order equations in terms of u µ is a formidable task, let alone the second order equations. In the remainder of this section we will study the equations of motion in the limit where the number of dimensions d becomes very large. We will see that with some restrictions, to be elaborated on below, the large d limit of these hydrodynamic equations become easier to handle.
Before working out the hydrodynamic equations, let us first consider thermodynamic equilibrium in the large d limit. In equilibrium, the temperature and velocity field are constant so that the stress tensor is given exactly by T µν (0) . Since the velocity is constant, we may conveniently choose a coordinate system where u µ ∂ µ = ∂ t . If we assume, without loss of generality, that the energy density is O(d 0 ) and insist that the pressure term, proportional to P µν , not be suprressed relative to the energy density, then we must work in a coordinate systems scaled so that the spatial coordinates carry an extra factor of 1/ √ d relative to the time coordinate: ds 2 = −dt 2 + i (dx i ) 2 /d. For the same reasoning, we parameterize the velocity field of a boosted thermally equilibrated system by u µ = (1, β i )/ 1 − β i δ ij β j /d. Note that our definition of β i differs from the conventional definition by a factor of √ d, β i = √ dβ i conventional . With this parameterization in mind, let us now go back to (2.4). Conservation of energy and momentum provide d equations for the d hydrodynamical variables T and u µ . In order for the large d limit to provide a useful proxy to a finite p + 1 dimensional system we must take the large d limit in such a way that the number of dynamical equations remains p + 1. To do so, we use an ansatz where the dynamical variables do not depend on n out of the d = 1 + p + n coordinates. Thus, we use where χ ⊥ denotes the n 1 coordinates on which T and u µ do not depend, and Latin indices run over the remaining p spatial components. The velocity field in this coordinate system is given by (2.9) -4 -In the remainder of this work we will raise and lower Latin indices with δ ab and δ ab respectively. Keeping with this notation, the tensor structures (2.3) evaluate to (2.10) Inserting (2.10) into (2.2), we find that in order for the shear viscosity and second order transport coefficients to contribute to the dynamics they must both scale as 1/d relative to . Recall that dimensional analysis and conformal invariance imply that for any finite d where the proportionality constants are dimensionless numbers. Thus, we may write where ξ is an emergent length scale which is determined by requiring that 1/d ξ remain finite in the large d limit.
Thus, in what follows we will use the parameterization where h 0 and i are dimension independent numbers. In the majority of this work we will set ξ = 1. It will be reintroduced where necessary. The large d stress-energy tensor is now given by with, The equations of motion (2.4) which result from the constitutive relations (2.15) are still somewhat unwieldy. To simplify them further we will use our freedom to perform field redefinitions of the velocity and temperature fields in order to remove all third order derivative terms from the equations of motion. The most general frame transformation we can carry out at second order in derivatives and in the large d limit is given by Inserting (2.16) into (2.14) and utilizing the derivative expansion, c.f., (2.5), we find that if then the second order terms in the equations of motion vanish (i.e., no contributions from second order terms in the stress tensor), provided we choose a frame where In what follows we will refer to (2.17b) as the large d frame. To iterate our finding in the previous paragraph, when (2.17) are satisfied, the equations of motion take the form and c 1 is given by the quadratic equation in (2.17b). The stress tensor whose divergence gives these equations of motion is given by (2.20) We emphasize that even though ∂ µ T µν (2) = 0, the equations of motion are only second order in derivatives due to cancellations coming from the perturbative expansion.

Holography in a Large Number of Dimensions
Our final result, equation (2.18), is, evidently, a simplified version of the full equations of motion. In order to obtain it we had to restrict our transport coefficients to satisfy (2.13) and (2.17a). Both equations impose restrictions on the transport coefficients when the number of dimensions is large. Equation (2.17a) implies that the four transport coefficients of the theory are restricted to lie on a plane in parameter space while equation (2.13) implies a particular scaling of the transport coefficients with d when d is large. Since the large d limit is used as a proxy for finite dimensional system one may, naively, enforce by hand a 1/d scaling of transport coefficients for large d so that both (2.13) and (2.17a) will be satisfied. Notwithstanding, it seems likely that such a behavior will affect the radius of convergence of the large d expansion. It is then somewhat surprising that the relations (2.17a) are valid in a holographic setup. There, we find that see [21,22]. Note that the first of the two relations is satisfied in any number of dimensions [21] even in the presence of matter [23] and also in other instances [24]. As we will now demonstrate the constraints (2.13) are also satisfied in a holographic context.
The constraints (2.13) and (2.17a) guarantee that the equations of motion will take the simplified form (2.18). Therefore, if they are satsified, and we work perturbatively in the derivative expansion, equations (2.18) must emerge from a holographic analysis [1]. We find that not only do (2.18) emerge, they are also valid to all orders in the derivative expansion. Put differently, the Einstein equations at large d are equivalent to large d hydrodynamics naturally truncated at second order in derivatives. In what follows we will show this property of the Einstein equations explicitly.
Consider the D = d + 1 space-time dimensional Einstein-Hilbert action where we have set the radius of AdS space to unity. Following [11], we wish to take the large D limit of this action while retaining translation invariance in most of the spatial directions. To this end, we use the ansatz where a, b = 1, . . . , p, D = p + n + 2 and x ⊥ is n dimensional. To obtain a boundary metric as in (2.8) we make the replcaements and choose boundary conditions such that (3.6) In addition, we rescale the radial coordinate such that is finite in the large d limit. Note that r 0 serves as a reference length scale. In what follows we will set r 0 = 1 for clarity. In the new coordinate system the line element takes the form (3.8) The scaling used to obtain (3.8) is precisely that of [11]. See also [15]. Inserting (3.8) into the Einstein equations and solving order by order in n we find that where the functions a and f a must satisfy the constraint equations As expected, equations (3.11) and (2.18) match as long as we set and make the identifications Thus, in the large D limit the Einstein equations precisely reproduce the large d limit of the equations of motion for relativistic hydrodynamics. The associated stress tensor is given in (2.20).
Before proceeding let us pause to consider the relation between the scale r 0 associated with the horizon, (3.7), and the scale ξ generated by taking the large d limit of hydrodynamics (2.12). Recall that for any finite d the (local) Hawking temperature associated with the event horizon is , (3.14) and that in the absence of charge the entropy density is given by where P is the pressure and = a the energy density. Using the holographic result for the shear viscosity to entropy density ratio [25] we find Comparing (3.16) to (2.12) and using (3.12) we obtain Given (3.17) the expression for the entropy density (3.15) takes the form in the large d limit and in our current conventions. 1 Since we have a good handle over the dynamics of the area element of the event horizon we can test various proposals regarding its behavior in the presence of turbulence, and analyze the geometric structure associated with turbulent flows.

Analysis of Large d Fluid Flows
Thus far we have described hydrodynamics using the variables and β a . Since we are working in the large d frame (2.17b) these variables are somewhat different from the energy density and velocity field as described in the Landau frame. Indeed, and β a coincide with the energy density and fluid velocity of the Landau frame only in equilibrium. Otherwise, they are related to the Landau frame variables 1 Note that using the Bekenstein-Hawking formula S = A/4G N we obtain s = through (2.16). Nevertheless, if we rewrite the equations of motion (2.18) in terms of and β a = j a / , we find In order to understand the dynamics associated with equations (4.1) and their relation to the compressible Navier-Stokes equations, it is useful to switch to dimensionless variables, where L 0 is a characteristic length scale of the system, √ nU a characteristic velocity , and E a characteristic energy density. The factor of √ n in the definition of the characteristic velocity arises from the scaling we used for β a in (2.9). In terms of these dimensionless variables, equations (4.1) take the form where, in analogy with the compressible Navier-Stokes equations, we have defined a Reynolds number Re = √ nL 0 U and a Mach number M = √ nU . If we reinsert the dimensionful parameters c and ξ we find where c s is the speed of sound. When the Mach number is small we may expand u and p in powers of the Mach number, On a manifold with no boundary p (0) and p (1) are constant and u (0) and p (2) satisfy the incompressible Navier-Stokes equation, with p (0) the density and p (2) the pressure. Equations (4.6) are precisely the Navier-Stokes equations for an incompressible fluid. Thus, in the limit of small Mach number we may use the exhaustive machinery developed for incompressible flow to study solutions of (4.6). While this is textbook material (see for instance [26]) let us remind the reader of the salient features of such flows. We define the total energy and enstrophy via with ω (0)ij = ∂ i u (0) j − ∂ j u (0) i the vorticity two-form. The evolution equations for E I and Ω I are given by where we have introduced the non relativistic shear tensor, σ ij (0) and the Palinstrophy, P I , We once again emphasize that the energy density is given by T 00 via (2.14) and that E I is the would be energy of the analog Navier-Stokes equation in the incompressible limit. Nevertheless, we shall, with some abuse of language, refer to E I as the energy and to Ω I as the enstrophy. Note that Ω I ≥ 0. Thus, the energy E I can only decrease. On the other hand, Ω I itself can increase or decrease depending on the sign of the first term on the right hand side of (4.8b). The dynamics associated with this term is often referred to as vortex stretching. There is ample phenomenological evidence and various imperfect arguments that vortex stretching increases the enstrophy for turbulent incompressible flow. Indeed, one usually posits that, for sustained turbulent flow, lim Re→∞ Ω I /Re = e 0 > 0. The appearance of an emergent scale e 0 at large Reynolds number implies an energy cascade. To see this consider the two point function and its Fourier transform We define the energy density E(k) via where dθ k is a solid angle in momentum space, viz. d p k = k p−1 dθ k dk. With this definition, and also If E I is to decrease while Ω I is to remain constant then E(k) must distribute itself in such a way that energy will flow from lower momentum modes to higher ones. This process is referred to as the Kolmogorov cascade. If we constantly supply energy into the system and assume that E(k) depends on e 0 and k we find the celebrated −5/3 law, Various numerical and experimental verifications of (4.15) can be found in [27,28].
An exception to (4.15) arises when p = 2 [29][30][31]. In two spatial dimensions we may treat ω ij as a volume form. It is then straightforward to show that the vortex stretching term vanishes for incompressible flow. In this case the enstrophy is bound from above by its value at t = 0 so that energy is conserved when the Reynolds number becomes large. In this case the palinstrophy plays the same role that enstrophy played in three dimensional flow. One then expects that the energy distribute itself towards lower wavenumber. If energy is continuously pumped into large scales via a driving force then in order for the system to remain in a steady state some type of large scale friction needs to be introduced into the system. This friction introduces an energy scale and yields a power law behavior as in (4.15) referred to as the inverse cascade. In addition to the inverse cacade there is a direct cascade associated with the enstrophy production term. An analysis similar to the one that led to (4.15) implies See [32] for various discussions.
In the current work we will study decaying turbulence where the typical velocity scale U and length scale L 0 vary with time. In the context of the Navier-Stokes equation, a power law behavior associated with a direct cascade, c.f., (4.15), was observed in numerical simulations of decaying turbulence in p = 3 dimensions [33][34][35]. In p = 2 dimensions numerical simulations of the Navier-Stokes equations usually lead to either an inverse or direct cascade. Numerical simulations of decaying turbulence with no slip boundary conditions may exhibit both a direct cascade and an inverse cascade, as in (4.15) and (4.16) [36,37]. In periodic domains a direct cascade as in (4.16) is often observed [38,39]. Recent simulations exhibit both a direct and an inverse cascade in such scenarios [40].
Let us return our attention to (4.3). Using intuition gained from the analysis at small Mach number we focus on the evolution of the energy and enstrophy. We find that our equations of motion where we now define and Note that we have used the letter p for both the number of spatial dimensions and the dimensionless energy density. In the equations above and in the remainder of this section the number of dimensions p will appear only in the measure d p x. Similar to the case of small Reynolds number, we find that E C is approximately conserved for large M and large Re. Furthermore, similar to (4.8b), we find the following equation for the dynamics of the enstrophy with We will refer to the first term on the right hand side of (4.20) as a vortex stretching term. This term vanishes for p = 2 (as does the last term on the right hand side of (4.20)). Thus Ω C is conserved at high Reynolds number and we may expect an inverse cascade in such a configuration. When Re is finite, it is difficult to determine the sign of the rate of change of Ω C and therefore difficult to assess a priori whether a scaling regime exists or not. If it does exist, we may expect the same power law behavior as in (4.16). Since (4.20) is less manageable than (4.8b), to make headway we must resort to numerical methods. We carry out such an analysis in the next section.

Analysis of Turbulent Flows
We now turn to a numerical analysis of turbulent flows. While the analysis of section 4 was valid for unbound domains, our simulations focus on flows in a bound, toroidal domains of length Lξ in each direction. We comment on the influence of these boundary conditions on our results when relevant. We solved the equations of motion (3.11) using a variety of methods including Fourier spectral methods and finite differencing in the spatial directions. We evolved the variables a and f a = u a forward in time using third order Adams-Bashforth or explicit Runge-Kutta. The initial conditions we used for our simulations were "perturbed shear flows", i.e. constant density flows where the velocity field is perpendicular to its gradient. Such flows solve (4.3) in the limit of infinite Reynolds number.
More specifically, for p = 2 we used whereas for p = 3 we used Here E is an energy scale which drops out of the equations of motion. The parameter a 0 in both (5.1a) and (5.1b) is constant in space and the δf i denote perturbations of the form where m, A i, m and ∆φ i, m are chosen from a random sample. The sum in (5.2) included 10 to 100 random modes. Typical simulations of two dimensional flow included 20 modes and typical three dimensional flow included 40 or 80 modes. The components of the wavenumber m was chosen from a uniform distribution of integers running from 1 to 64 for a typical two dimensional flow and 1 to 16 or 1 to 32 for a typical three dimensional flow. The amplitude A i, m was chosen from a uniform distribution ranging from 0 to A with A a pre-determined parameter. The phase ∆φ i, m was chosen from a uniform distribution ranging from 0 to 2π. Since our setup does not involve a driving force, the flow we generate does not reach steady state turbulence. Nevertheless, we may associate a Reynolds number and a Mach number to the initial conditions of the flow. Choosing U = max(u) − u and L 0 = Lξ/n, we find that If the initial Reynolds number is sufficiently large we expect a turbulent instability to emerge. From a numerical standpoint these instabilities are trigered by the modes characterized by A i, m in (5.1). For sufficiently large Re we found that numerical roundoff error was sufficient to trigger the instability even for A i, m = 0. In presenting our results we find it convenient to use scaled coordinates for time and energy power spectrum respectively. Typical results of two dimensional and three dimensional fluid flow with relatively high Reynolds number can be found in figures 1 and 2.
In both two and three dimensions we found that the dynamical behavior of the fluid can be characterized by three stages. An initial stage where an instability drives the shear flow into a chaotic configuration. Typical flow in this stage can be found in the first row of figures 1 and 2. Once the instability sets in, we observe a turbulent regime where the energy power spectrum exhibits power law behavior. Typical turbulent behavior can be seen in the second row of figures 1 and 2. The distinction between two dimensional flow and three dimensional flow is most apparent (visually) in the turbulent phase: in two dimensional flow large scale structure is formed while in three dimensional flow, vortices break up into smaller ones. In the final stage the flow decays into an equilibrated configuration. This behavior can be observed in the last row of figures 1 and 2. In the remainder of this section we will discuss each of these regimes in detail.

Initial Phase: Onset of Instability
The initial phase is characterized by the development of instabilities. To quantify these instabilities, and later also turbulence, we analyze the energy spectrum E C , defined similar to E I in (4.13), x) e −i k· x . When the Reynolds number is small we find that the shear flow quickly decays to an equilibrated configuration. When the Reynolds number is sufficiently large we find that instabilities dominate the spectrum and lead to turbulent berhavior. Typical behavior of turbulent instabilities in two dimensional fluid flow can be observed in the top right panel of figure 3. This particular run was carried out with Re = 1562.5, M = 0.5, n = 32 and an amplitude A = 10 −5 for the noise. We have carried out similar runs with Mach numbers ranging from M = 0.005 to M = 2 and Reynolds number of up to Re ∼ 1562. As the Mach number increased the amplitude of the random noise required to generate turbulence became higher; for low Mach number numerical noise was sufficient to generate the instability while for M = 2 we had to set A = 0.01 to generate a turbulent instability. It is possible that at very high Mach numbers and the initial conditions (5.1) a turbulent instability will not form.
Going back to figure 3, around t = 10τ the initial disturbance with k = 32 × 2π/L has decayed and an unstable mode with k = 22 × 2π/L can be observed. A more detailed analysis of the spectrum shows that the unstable mode is a shear mode orthogonal to the first. Indeed, in figure 4 we plot the Fourier decomposition of |u p | 2 at different times. At t = 6τ , the instability has not set in and |u p | receives support from the initial shear mode which we have set as input into the system. At t = 12τ some remnants of the initial shear mode are observed (in red), but most of the support for |u p | comes The instability becomes apparent around t ≈ 10τ , followed by a chaotic behavior generating small eddies shown at t = 20τ . The inverse cascade is apparent in the "merging" of co-rotating eddies to form larger ones starting around t ≈ 20τ and ending at t ≈ 960τ . The termination of the inverse cascade is caused by the finite length of the box. The final stage shown at t = 960τ is of slow decay of two counter rotating eddies. from a transverse shear mode. The fact that instabilities with lower wavenumber become excited is suggestive of the inverse cascade present in two dimensional flow. In three dimensions, the instabilities are associated with higher wave numbers. Typical behavior for instabilities in three dimensions can be found in the bottom right corner of figure 3. The data for the bottom right plot in figure 3 was obtained for a flow with Re = 162.5, M = 1 and an initial perturbation associated with the fourth harmonic, n = 4. We have carried out similar runs with Mach numbers of up to 10. Here, as opposed to the two dimensional case, a turbulent instability seems to emerge even for high Mach numbers. For the initial conditions associated with figure 3, around t = 0.1τ we observe a growing unstable mode which is roughly double the wavelength of the original.
A typical Fourier decomposition of the velocity field for three dimensional flow can be found in figure 5 where we show constant k z slices of the spectrum of | u| at different times. The spectrum of the instabilities suggests that the dominant unstable modes take the form u i ∼ cos 2πn L (i ± j) where i, j ∈ x, y, z, j is the initial mode direction of f i in (5.1). (For example: if the initial shear flow was of . At t = 6τ (left) the initial disturbance at kx = 32 × 2π/L, ky = 0 dominates the spectrum. At t = 12τ an instability with kx ∼ 0 and ky = 22 × 2π/L, orthogonal to the initial disturbance but with approximately half the initial wavenumber has set in and dominates the spectrum. Additional weaker instabilities with kx ∼ 20 × 2π/L and ky ∼ 0 can also be observed.
the form f x = Ec s cos 2πn L y , then the associated unstable mode is u x ∼ cos 2πn L (x ± y) . ) The authors of [3] define a critical Reynolds number as the minimal Reynolds number required for a turbulent instability to exist. Put differently, having the Reynolds number above the critical one is a necessary condition for the appearance of turbulence. An analysis similar to that of [3] implies that the critical Reynolds number associated with our initial conditions is of the order of 15, roughly two orders of magnitude smaller than the initial Reynolds number needed for turbulence to develop. We provide more details regarding the computation of the critical Reynolds number in appendix A.

Turbulent Phase
Once the instability of the initial phase is strong enough, it drives the fluid to a turbulent regime which we characterize using the energy density power spectrum, E C , defined in (5.5). We present a detailed analysis of the power spectrum during the turbulent phase in appendix B. Here we confine ourselves to summary of the salient features of that analysis. Representative plots of E C for two and three dimensional fluid flow can be found in figures 6 and 7 respectively.
For two dimensional fluid flow we found a consistent k −4 power law for Mach numbers between 0.005 and 2 and Reynolds number between 781 and 1562. As the Reynolds number decreases the range for which power law behavior is observed becomes smaller and vanishes completely around Re ∼ 390. See appendix B. While it is clear that lower modes get populated as time progresses, indicative of an inverse cascade, the expected k −5/3 and (or) k −3 law is absent from the simulations we have studied. Early simulations of non-relativistic decaying turbulence displayed similar scaling [38,39]. Perhaps increasing the initial Reynolds number of our flow will ameliorate this problem, as it does in the non-relativistic case.
In contrast to two dimensional fluid flow, in three dimensions we found remarkable agreement with an E C ∼ k −5/3 power law. Our runs include Mach numbers between M = 0.1 and M = 10 and Reynolds numbers between Re = 81.25 and Re = 750, and initial data involving initial modes n = 1 and n = 4. Simulations with initial modes with n > 4 and reasonably high (initial) Reynolds number are expensive. Indeed, for the n = 4 run the initial Reynolds number was rather low, Re = 162.5, which apparently manifested itself as a visible λ = L/4 periodic behavior of the flow even in the "turbulent" regime where Kolmogorov scaling was observed. Figure 6. The energy spectrum for two dimensional flow with an initial Reynolds number of Re = 1562.5 initial Mach number given by M = 0.5 and initial mode n = 32. The right plot provides an overlay of the energy spectrum at various times. An emergent power law EC (t, k) ∼ k −4 is observed from k ∼ 10 × 2π/L to k ∼ 100 × 2π/L. By restricting oneself to lower wavenumber it is possible to fit the data to other power law behavior though the merit in doing so is unclear.

Final Phase
Since there is no driving force the fluid is expected to reach equilibrium at late times. In two dimensional non-relativistic and incompressible fluid flow on R 2 the late time behavior of the velocity field is given by the Oseen Vortex solution which is an attractor of the Navier-Stokes equation [41]. Since we are placing our fluid on a torus the late time behavior of the fluid is somewhat different. In particular, we find that due to the inverse cascade the lowest lying mode dominates the flow, such that at late where E 0 , U x 0 and U y 0 are constants. One can check that (5.6) solves (4.3) for low Mach number up to exponentially suppressed corrections. Typical flow of the form (5.6) is depicted in the lower right corner of figure 1. A spectral analysis of the flow can be found in figure 8. To get a quantitative handle on (5.6) we have fitted the late time behavior of u x and u y to an exponential decay law u i ∼ e −αi(L)t = e −νi( 2π L ) 2 t . Let us denote the L 2 norm of a quantity X by L 2 (X). We evaluate α i (L) by fitting the dependence of L 2 (u x ) and L 2 (u y ) on time to a power law fall off, obtaining ν x = 1.004 ± 0.066, ν y = 1.002 ± 0.053. (5.7) See figure 9.
Since simulating three dimensional fluid flow at late times is expensive, we have not carried out a full analysis of its late time behavior as we have done for two dimensional flow. There is some indication that the late time solution will be dominated by the modes associated with the initial condition (5.1). For the runs we have carried out in three dimensions with n = 1 and n = 4 (see (5.1)) we find that the late time behavior takes the form (5.8) Figure 9. The decay rate of the velocity fields ux and uy, denoted by αx (left) and αy (right) at late times as a function of (2π/L) 2 with L the length of the torus. Box lengths ranged from L = 5000 to L = 15000. The solid lines denote a linear fit where n is the mode number injected into the system in (5.1b). Typical flow of the form (5.8) is exhibited on the bottom row of figure 2. A k-space view of the velocity and energy density fields matching (5.8) can be found in figure 10. A possible explanation for (5.8) may be that the initial perturbation still holds most of the energy at late times and due to the direct cascade is slowest to decay. It is also possible that our three dimensional simulations do not exhibit full turbulent behavior which would be in line with the periodic behavior we observe for three dimensional flow at n = 4, as mentioned earlier.

Geometrizing Turbulence
As should be clear from our discussion so far, a full understanding of turbulence, even in the incompressible limit, is far from complete. It would be favorable if it were possible to utilize the geometric tools available from numerous studies of black hole dynamics to address turbulence. The AdS/CFT correspondence opens the possibility for such a procedure [1] and there are numerous suggestions for identifying an appropriate geometric quantity which captures the turbulent behavior of the dual fluid [2,4,42]. In what follows we will focus on the work of [2] where the authors argue that the horizon power spectrum, A, to be defined shortly, is proportional to the energy power spectrum defined in (4.7) and (4.13). By appealing to the large d limit we obtain an analytic handle over such a relation and can study its regime of validity.
Recall that the leading order contribution to the black brane metric is given by The only null surface associated with the metric (6.1) which agrees with a black brane topology is given by R = a (t, ζ a ) .
which must therefore be identified with the event horizon. The extrinsic curvature on the event horizon is given by where Latin indices run over all d = n + p + 2 dimensions, n M is the null normal to the horizon, and M is an auxiliary null vector, which satisfies M n M = −1. Following [2] we consider the rescaled traceless horizon curvature, defined by where i, j run over the n + p spatial dimensions of the horizon, √ γ ≡ det(g ij ) is the area element on a spatial slice of the event horizon, and κ is defined by the geodesic equation n M ∇ M n Q = κn Q . The horizon curvature power spectrum is defined by (6.6) and p is the number of spatial dimensions where the dynamics take place. An explicit computation gives us and follows. Substituting a = , f a = β a and taking the incompressible limit, we find The Fourier transformθ will involve a convolution of˜ andβ which are the Fourier transform of and β respectively, a 'la (6.6). However, if we are working in the limit of small Mach number then is approximately constant and we find θ i j n (p+n)/2 ik jβ i + ik iβ j . (6.11) Inserting (6.11) into (6.5) and comparing to (5.5) we find that at low Mach number, as predicted in [2]. If, on the other hand, the Mach number is not small a relation of the form A ∼ k 2 E will not hold. Indeed, in figure 11 we show typical results for the horizon power spectrum for high and low initial Mach number in two and three dimensions. As expected, A ∼ k 2 E holds only for very low Mach number. A more refined analysis can be found in appendix C.

Summary and Outlook
In this paper we have discussed relativistic hydrodynamics in the limit where the number of dimensions is large, revealing the simplifications that occur in this limit. We have focused our attention particularly on holographic theories whose dynamics surprisingly follow these simplified equations. We have analyzed turbulent flows of these dynamical systems, and their relation to the geometry of black hole horizons. We have seen that three dimensional flows exhibit a Kolmogorov cascade with the expected power law behavior, k −5/3 , for a wide range of Mach numbers. While promising, we remind the reader that the initial conditions for our three dimensional simulations where carried out with a low wavenumber (n = 1 and n = 4). Higher wave numbers would require more computational resources. Moreover, in the n = 4 simulations we observed a periodicity with wavelength λ = L/4 throughout the flow. We expect that such behavior will disappear for sufficiently high initial Reynolds number which we have not yet reached.
In contrast to three dimensional flows, two dimensional flows did not exhibit either a k −5/3 power law or a k −3 power law, albeit displaying a tendency to form large scale structures, associated with the expected inverse cascade. The behavior we observe is probably due to an insufficiently high Reynolds number. It would be worthwhile to improve on this point; simulations of decaying turbulence which exhibited both a direct and inverse cascade usually require a grid much larger than the one we have used [40]. Preliminary runs do indicate that increasing the Reynolds number may result in a canonical power law.
Recall that the power law behavior of the energy power spectrum has been evaluated for sustained turbulence in which case a driving force continuously supplies energy into the system. It would be interesting to study sustained turbulence in our setup, which would allow for more robust scaling relations to be observed in the steady state. In order to reach that steady state, one would need to add stochastic force to supply energy at the appropriate range of wave numbers, and for the scenario of inverse cascade to include friction as a sink of energy at large scales (keeping in mind that in the direct cascade energy is dissipated at small scales.) In addition to the Kolmogorov scaling discussed in this work, steady state turbulence would allow us to investigate real space scaling relations of the sort recently described in [43]. In the present context of decaying turbulence, those relations are not stable enough to be clearly visible, but we expect that in sustained turbulence we would be able to investigate them more reliably. Similar comments apply to the holographic study of superfluid turbulence. While holographic superfluidity was studied in the large dimension limit [44], the reduction to horizon equations does not hold in the presence of a scalar field. We hope to return to this problem in the future.
While in principle these additional elements are possible in the large dimension limit, it seems we are no longer afforded the simplifications of that limit, with those additional elements. The dimensional reduction of the equations, allowing us to discuss the horizon fluid first and deduce the resulting spacetime subsequently, is no longer in effect when adding external dials, since those are imposed at infinity. Nevertheless, perhaps there is another simplifying limit that would allow for investigation of sustained turbulence in the present context.
Finally, using the simplicity of the large d black brane metric we were able to compare our analytic expression for the area power spectrum to the energy power spectrum. We have found that the power law behavior of these two quantities is related only at very low Mach number at which point there is very good agreement with [2]. It would be interesting to identify a geometric entity which captures the Kolmogorov power law behavior. With such a quantity at hand one may be able to use the powerful tools of general relativity to compute this quantity in a more precise manner.
for M = 100 and L = 1. While it is relatively simple to generate three dimensional flow with large Mach number, it is difficult to numerically simulate a flow with initial conditions (5.1) and large n. The main complication in simulating large n is related to the direct cascade where the initial disturbance is pushed to large wave numbers. The fluids tendency to populate modes with large wavenumber require a large number of grid points to simulate.

B Detailed Analysis of the Power Spectrum
As discussed in section 5.2, given a sufficiently high Reynolds number the flow eventually reaches a turbulent regime in which the energy power spectrum E C has power law behavior. In two dimensional decaying turbulence we have observed a k −4 power law behavior for a variety of initial Mach numbers and initial Reynolds numbers. In figure 14 we have plotted the energy power spectrum for fixed initial Reynolds number and varying Mach number. In figure 15 we have plotted the energy power spectrum for fixed Mach number and varying Reynolds number. As the initial Reynolds number decreases, the time t at which the power law behavior may be observed increases, and the size of the inertial range decreases.
A spectral analysis of E C at intermediate times, in three spatial dimensions, can be found in figures 16 and 17. In figure 16 we have plotted E C for a variety of initial Reynolds and Mach numbers and an initial disturbance (5.1) with n = 4. Figure 17 describes fluid flow for the same initial Reynolds and Mach numbers but an initial disturbance with n = 1. Somewhat surprisingly, the Kolmogorov power law E C ∼ k −5/3 , seems to be robust and holds also for large Mach number when the incompressible approximation is no longer valid [45].   and various initial Reynold number. We note that the size of the inertial range, the power law behavior and the times at which power law behavior is observed are dependent on the Reynolds number. This behavior should be contrasted with that exhibited in figure 14 .

C Analysis of the Horizon Area Power Spectrum
In (6.12) we have argued that the ratio of the horizon area power law spectrum to the energy spectrum grows like k 2 only for flows with low Mach number. Recall that the traceless part of the expansion, θ i j is given by equation (6.9). Once the Mach number is small enough the flow becomes incompressible and we may approximate (6.9) by (6.10) which yields A/E C ∼ k 2 . A numerical analysis of A/E C for various initial Mach and Reynolds number is displayed in figure 19. As expected, the dependence of A/E C on k seems to deviate from k 2 once the Mach number becomes large.
To quantify the deviation of A/E c from a k 2 behavior we go back to the compressible contributions to the traceless expansion (6.9). Let us denote the Fourier transform of a quantity X by F(X). Then, in order for the incompressible terms in (6.9) to dominate over the compressible ones, we need The results on the right are averaged in jumps of four k points, since the initial mode n = 4 sets the jump scale. One can see the Mach number has little to no effect on the spectrum. All the plots display a weak signature of the expected power law EC (t, k) ∼ k −5/3 . The lowering of the Reynolds number has a significant effect on the power law scaling range, shortening it from k ∈ (10, 30) when Re = 162.     figure 21. The right plots shows the typical relation for the diagonal terms M j= . The left plots depicts the typical relation for the off diagonal terms M j = . The variance between Mach numbers, rather high or low appears almost linear on the log scale for M j = , and is drastically different from the two dimensional flows on figure 18. This is responsible for the lack of accuracy of the relation (6.12). The aforementioned variance between Mach numbers leads to the conclusion that a much lower Mach number than M ∼ 0.1 is required in order to drop high order derivative terms in (6.9).