Holographic Turbulence in a Large Number of Dimensions
Moshe Rozali Evyatar Sabag Amos Yarom
hep-th gr-qc nlin.CD physics.flu-dyn


Unconditional violation of the shot noise limit in photonic quantum metrology
    Geoff J. Pryde
     
=============================================================================





§ 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, <cit.>. Indeed,  in a suitable parameter range, and given appropriate initial conditions, such flows have been shown to exhibit (decaying) turbulent behavior <cit.>.  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 <cit.>. 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 <cit.>. The works of <cit.> 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 <ref>, 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 (<ref>).

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 <cit.> 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 <ref> the large d Einstein equations of <cit.> 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 <ref> 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 <ref> we turn to a numerical analysis of the equations of motion in section <ref>. 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 <ref> for two and three spatial dimensions and to an analysis of holographic (decaying) turbulence in AdS_4 <cit.>.

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 <ref>, a proposed relation between the horizon curvature power spectrum and the hydrodynamic energy power spectrum <cit.>. 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 <ref> 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

    T^μν=T^μν_(0)+T^μν_(1)+T^μν_(2)+𝒪(∂^3) ,

where

    T^μν_(0) =   ϵ(T) u^μu^ν + ϵ(T)/d-1 P^μν ,
                
            T^μν_(1)  =  -2η(T) σ^μν ,
    
            T_(2)^μν =   λ_1(T) u^λ𝒟_λσ^μν+λ_2(T) (σ^λμσ_λ^  ν-σ^αβσ_αβ/d-1P^μν)   
        +λ_3(T) (ω^μλσ_λ^  ν+ω^νλσ_λ^  μ)+λ_4(T)(ω^μλω^ν_λ  -ω^αβω_αβ/d-1P^μν)  ,

and

    P^μν    =  g^μν+u^μu^ν,
    σ^μν    =  1/2P^μαP^νβ(∂_αu_β+∂_βu_α)-1/d-1P^μν∂_αu^α,
    ω^μν    =  1/2P^μαP^νβ(∂_αu_β-∂_βu_α), 
    𝒟_λσ^μν    =  P^μαP^νβ∂_λσ_αβ+∂_αu^α/d-1σ^μν ,

and we have set the speed of light to unity.

See <cit.> for details.

Following the standard reasoning behind the construction of effective field theories, the constitutive relations (<ref>) 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 (<ref>) 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 (<ref>). 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 

    ∂_μT^μν=0 ,
  
which need to be solved perturbatively in a derivative expansion. For instance, at leading order in derivatives one has

    u^μ∂_μϵ + d ϵ∂_μu^μ    = 0 
            
    d/d-1ϵu^μ∂_μu^ν  + P^μν∂_μϵ    = 0 .

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 (<ref>) for T_(1)^μν and T_(2)^μν  by an appropriate field redefinition of T and u^μ. The constitutive relations (<ref>) were written in the Landau frame, defined via

    u_μT^μν = u_μT_(0)^μν  = - ϵ(T) u^ν .

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

    ϵ_new(T_new)    = ϵ(T) - k_1 u^μ∂_μϵ(T)   + 𝒪(∂^2)  
    ϵ(T) u^μ_new    = ϵ(T) u^μ - c_1 ∂^μϵ(T) +  𝒪(∂^2)  ,

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 (<ref>) into the conservation equation (<ref>), are lengthy and cumbersome. While the zeroth order equations of motion (<ref>) 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 𝒪(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)β_conventional^i.


With this parameterization in mind, let us now go back to (<ref>). 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

    ds^2 = - dt^2 + δ_ab/ndζ^a dζ^b + dχ⃗_⊥^2/n

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

    u^μ = (1,β^a(t,ζ), 0⃗)/√(1-β^b(t,ζ)β^c(t,ζ)δ_bc/1+p+n) .

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 (<ref>) evaluate to

    σ^μν    =  n/2δ^μ aδ^ν b(∂_aβ_b + ∂_bβ_a) + 𝒪(n^0),
    ω^μν    =  n/2δ^μ aδ^ν b(∂_aβ_b - ∂_bβ_a) + 𝒪(n^0), 
    
            u^λ𝒟_λσ^μν    =  n/2δ^μ aδ^ν b(∂_t+β^c∂_c)(∂_aβ_b+∂_bβ_a)+𝒪(n^0) .

Inserting (<ref>) into (<ref>), 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

    ϵ∝ T^d ,
                η∝ T^d-1 ,
                λ_i ∝ T^d-2 ,

where the proportionality constants are dimensionless numbers. Thus, we may write

    η∝1/dϵ^1-1/d
    
    	        = ϵξ/d(ϵξ^d)^-1/d
    ϵξ/d

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

    η = h_0/nξϵ + 𝒪(n^-2)
                λ_i = 2 ℓ_i/nξ^2 ϵ + 𝒪(n^-2) ,

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

    T^μν = ϵ([                                      1                                    β^b;                                    β^a β^aβ^b+δ^ab+T̃^ab_(1)+T̃^ab_(2)+𝒪(∂^3) ])+𝒪(n^-1)

with,

    T̃^ab_(1) =     -h_0(∂^aβ^b+∂^bβ^a),
    T̃^ab_(2)  =    ℓ_1(∂_t+β^c∂_c)(∂^aβ^b+∂^bβ^a)
                    +ℓ_2/2(∂^aβ^c+∂^cβ^a)(∂^bβ_c+∂_cβ^b)  
        +ℓ_3/2( (∂^ aβ^c-∂^cβ^ a)(∂^bβ_c+∂_cβ^b) + (a↔b))
                    +ℓ_4/2(∂^aβ^c-∂^cβ^a)(∂^bβ_c-∂_cβ^b) .



The equations of motion (<ref>) which result from the constitutive relations (<ref>) 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

    ϵ→    ϵ + k_1ϵ∂_bβ^b 
                    +  k_2∂_b∂^bϵ + k_3ϵ(∂_bβ^b)^2 + k_4ϵ∂_bβ^c∂_cβ^b + k_5∂_bϵ∂^bϵ/ϵ + k_6ϵ∂_bβ_c∂^bβ^c  + 𝒪(∂^3)
    β^a →    β^a + c_1∂^aϵ/ϵ 
                     + c_2∂_bβ^b∂^aϵ/ϵ + c_3∂^aβ^b∂_bϵ/ϵ + c_4∂^a∂_bβ^b + c_5∂_bϵ/ϵ∂^bβ^a + c_6∂_b∂^bβ^a + 𝒪(∂^3)  .


Inserting (<ref>) into (<ref>) and utilizing the derivative expansion, c.f., (<ref>), we find that if



    ℓ_2 + 2 ℓ_3 + ℓ_4 = 0 ,
                and    ℓ_1 - ℓ_2 - ℓ_3 = 0 ,

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

    k_2 = ℓ_1 + ℓ_4
                
            -2 c_1 h_0 = c_1^2 + 2 ℓ_1
                
            c_5 = 2 (ℓ_1 + ℓ_4) 
                
            c_6 = ℓ_1 + ℓ_4 ,


and the remaining c_i's and k_i's vanish.

In what follows we will refer to (<ref>) as the large d frame. To iterate our finding in the previous paragraph, when (<ref>) are satisfied, the equations of motion take the form 
    ∂_t ϵ + c_1 ∂_a ∂^a ϵ    = - ∂_a j^a 
    ∂_t j^a - h_0  ∂_b ∂^b j^a     = -∂^a ϵ - ∂_b ( j^a j^b/ϵ) - 
                     (c_1/2 - ℓ_1/c_1) ∂_b(j^a/ϵ∂^bϵ - ϵ∂^a(j^b/ϵ))

where we have defined

    j^a = β^a ϵ ,

and c_1 is given by the quadratic equation in (<ref>).
The stress tensor whose divergence gives these equations of motion is given by

    T_(0)^μν   = [                 ϵ               j^a;               j^b δ^abϵ + j^a j^b/ϵ ]    
             T_(1)^μν = [                                                                     0                                                             c_1 ∂^a ϵ;                                                             c_1 ∂^b ϵ (c_1 + h_0) (  j^a/ϵ∂^b ϵ + j^b/ϵ∂^a ϵ)  -h_0 (  ∂^a j^b  + ∂^b j^a ) ]
    
            T_(2)^μν    = [               0               0;               0 c_1^2 ∂^a ∂^b ϵ ] + (ℓ_1 + ℓ_4) ∂_c ∂^c T_(0)^μν .

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 (<ref>), 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 (<ref>) and (<ref>). Both equations impose restrictions on the transport coefficients when the number of dimensions is large. Equation (<ref>) implies that the four transport coefficients of the theory are restricted to lie on a plane in parameter space while equation (<ref>) 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 (<ref>) and (<ref>) 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 (<ref>) are valid in a holographic setup. There, we find that

    λ_1-λ_2 - λ_3 = 0
                λ_2 + 2λ_3 -λ_4 = 𝒪(d^-3)

see <cit.>. Note that the first of the two relations is satisfied in any number of dimensions <cit.> even in the presence of matter <cit.> and also in other instances <cit.>. As we will now demonstrate the constraints (<ref>) are also satisfied in a holographic context. 


The constraints (<ref>) and (<ref>) guarantee that the equations of motion will take the simplified form (<ref>). Therefore, if they are satsified, and we work perturbatively in the derivative expansion, equations (<ref>) must emerge from a holographic analysis <cit.>. We find that not only do (<ref>) 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

    S = ∫√(g)(R + (D-1)(D-2) ) d^Dx ,

where we have set the radius of AdS space to unity.
Following <cit.>, 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

    ds^2=2dt(-A(t,r,z^i)dt+dr-F_a(t,r,z^i)dz^a)+G_ab(t,r,z^i)dz^adz^b+G_⊥(t,r,z^i)dx⃗_⊥^2

where a,b=1,…,p, D=p+n+2 and x⃗_⊥ is n dimensional. To obtain a boundary metric as in (<ref>) we make the replcaements

    z^a=ζ^a/√(n)    x⃗_⊥=χ⃗_⊥/√(n)

as well as

    F̃_a = F_a/√(n)    G̃_ab = G_ab/n    G̃_⊥ = G_⊥/n ,

and choose boundary conditions such that 

    A=r^2(1/2+𝒪(r^-n))
              F̃_a=r^2/n(𝒪(r^-n))
              G̃_ab=r^2/n(δ_ab+𝒪(r^-n))
              G̃_⊥=r^2/n(1+𝒪(r^-n)) .

In addition, we rescale the radial coordinate such that

    R=(r/r_0)^n

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

    ds^2=2dt(-A(t,R,ζ^c)dt+R^1/n-1/ndR-F̃_a(t,R,ζ^c)dζ^a)+G̃_ab(t,R,ζ^c)dζ^adζ^b+G̃_⊥(t,R,ζ^c)dχ⃗_⊥^2 .


The scaling used to obtain (<ref>) is precisely that of <cit.>. See also <cit.>. Inserting (<ref>) into the Einstein equations and solving order by order in n we find that

    A =1/2-a(ζ^μ)/2R + 𝒪(n^-1)
                F̃_a  =f_a(ζ^μ)/nR+ 𝒪(n^-2)
       

    G̃_ab=δ_ab/n + 𝒪(n^-2)
                G̃_⊥= 1/n + 𝒪(n^-2)

where the functions a and f_a must satisfy the constraint equations

    ∂_ta-∂_b∂^ba     =     -∂_bf^b
    ∂_tf_a-∂_b∂^bf_a    =     -∂_aa-∂_b(f_af^b/a) .


As expected, equations (<ref>) and (<ref>) match as long as we set

    h_0=1 ,
                
            c_1=-1 ,
                ℓ_1= 1/2

and make the identifications

    a = ϵ    
            f_a = j_a .

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 (<ref>).

Before proceeding let us pause to consider the relation between the scale r_0 associated with the horizon, (<ref>), and the scale ξ generated by taking the large d limit of hydrodynamics (<ref>). Recall that for any finite d the (local) Hawking temperature associated with the event horizon is

    1/a^1/d r_0 = d/4π T ,

and that in the absence of charge the entropy density is given by

    s = ∂ P/∂ T = a/T

where P is the pressure and ϵ=a the energy density. Using the holographic result for the shear viscosity to entropy density ratio <cit.> we find

    η = a^1-1/d/d r_0 .

Comparing (<ref>) to (<ref>) and using (<ref>) we obtain

    ξ = 1/r_0 .




Given (<ref>) the expression for the entropy density (<ref>) takes the form

    s = 4 πξ/d a

in the large d limit and in our current conventions.[Note that using the Bekenstein-Hawking formula S=A/4 G_N we obtain s = r_0^p+n/4 G_N n^n/2 a. Equation (<ref>) agrees with this result up to an overall normalization which was chosen when we set ϵ = a in (<ref>).]  
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 (<ref>) 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 through (<ref>). Nevertheless, if we rewrite the equations of motion (<ref>) in terms of ϵ and β^a = j^a/ϵ, we find

    ∂_tϵ + (β⃗·∇⃗) ϵ + ϵ∇⃗·β⃗    = ∇^2 ϵ
    ∂_tβ⃗+(β⃗·∇)β⃗ + ∇⃗ϵ/ϵ    =  ∇^2 β⃗ +  2(∇⃗ϵ/ϵ·∇⃗)β⃗ .



In order to understand the dynamics associated with equations (<ref>) and their relation to the compressible Navier-Stokes equations, it is useful to switch to dimensionless variables,

    u⃗ = β⃗/√(n) U ,
              
            p = ϵ/E ,
              ∂/∂ t→L_0/√(n) U∂/∂ t ,
              ∇⃗→ L_0 ∇⃗,

where L_0 is a characteristic length scale of the system, √(n)U 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 (<ref>). 
In terms of these dimensionless variables, equations (<ref>) take the form

    ∂_tp +(u⃗·∇⃗)p + p∇⃗·u⃗   =  1/Re∇^2 p  
    ∂_tu⃗+(u⃗·∇⃗)u⃗ + ∇⃗ p/M^2 p    = 1/Re∇^2u⃗ + 2/Re(∇⃗ p/p·∇⃗)u⃗ ,

where, in analogy with the compressible Navier-Stokes equations, we have defined a Reynolds number Re=√(n)L_0U and a Mach number M=√(n)U.
If we reinsert the dimensionful parameters c and ξ we find

    Re=√(n)L_0U/ξ c=L_0U/ξ c_s    
            M = √(n)U/c=U/c_s .

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,


    u⃗ = ∑_n=0 M^n u⃗_(n)    
            p = ∑_n=2 M^n p_(n) .

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,

    ∇⃗·u⃗_(0) = 0
                ∂_tu⃗_(0)+(u⃗_(0)·∇⃗)u⃗_(0) + ∇⃗ p_(2)/ϵ_(0)  = 1/Re∇^2 u⃗_(0) .

with p_(0) the density and  p_(2) the pressure.
Equations (<ref>) 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 (<ref>). While this is textbook material (see for instance <cit.>) let us remind the reader of the salient features of such flows. 

We define the total energy and enstrophy via

    E_I =1/2∫  |u⃗_(0)|^2 d^px ,
                Ω_I = 1/2∫ω_(0)ijω_(0)^ij d^px

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



    ∂/∂ t E_I    = -1/ReΩ_I 
    ∂/∂ tΩ_I    = ∫ω_(0)^ijω_(0)jkσ_(0)^ki d^px - 1/Re P_I


where we have introduced the non relativistic shear tensor, σ_(0)^ij and the Palinstrophy, P_I,

    σ_(0)ij = ∂_i u_(0) j + ∂_j u_(0) i ,
                
            P_I = 1/2∫∂_k ω_(0)ij∂^k ω_(0)^ij d^dp .

We once again emphasize that the energy density is given by T^00 via (<ref>) 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 (<ref>). 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

    Q(r⃗) = ∫u⃗_(0)(x⃗) · u_(0)(x⃗ + r⃗) d^p x/∫ d^p x

and its Fourier transform

    Q(k⃗) = 1/(2π)^p∫ Q(r⃗) e^-i k⃗·r⃗ d^pr .

We define the energy density E(k) via

    E(k) = 1/2∫Q(k⃗) k^p-1 dθ_k

where dθ_k is a solid angle in momentum space, viz. d^pk = k^p-1 dθ_k dk. With this definition,

    E_I  = ∫_0^∞ E(k) dk ,

and also

    Ω_I = ∫_0^∞ k^2 E(k) dk .

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,

    E(k) ∼ e_0 k^-5/3 .

Various numerical and experimental verifications of (<ref>) can be found in <cit.>.

An exception to (<ref>) arises when p=2 <cit.>. 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 (<ref>) 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 (<ref>) implies

    E(k) ∼ w_0 k^-3 .

See <cit.> 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., (<ref>), was observed in numerical simulations of decaying turbulence in p=3 dimensions <cit.>. 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 (<ref>) and (<ref>) <cit.>. In periodic domains a direct cascade as in (<ref>) is often observed <cit.>. Recent simulations exhibit both a direct and an inverse cascade in such scenarios <cit.>.

Let us return our attention to (<ref>). 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 (<ref>) imply

    ∂/∂ t E_C  = 1/M^2∫ p    (∇⃗·u⃗)   d^px  -1/4 Re∫ p  ( ω_ijω^ij + σ_ijσ^ij)d^px

where we now define

    E_C = ∫ p   |u⃗|^2 d^px

and 

    ω_ij = ∂_i u_j - ∂_j u_i
                σ_ij = ∂_i u_j + ∂_j u_i  .

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

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 (<ref>), we find the following equation for the dynamics of the enstrophy

    ∂/∂ tΩ_C = ∫1/pω_ijω^jk( σ_k^i - δ_k^i∇⃗·u⃗) d^px 
            + 1/Re∫|∇⃗ p |^2/p^3ω^2 
            - p ∂_k ω_ij/p∂_k ω_ij/p
    
            - 2/p^2σ_klω^ lj( ∂_j ∂_k p - ∂_j p ∂_k p/p)
            - 2/p^2ω_klω^ lj( ∂_j ∂_k p - ∂_j p ∂_k p/p - δ_jk( ∇^ 2p - |∇⃗ p |^2 )) d^px ,

with

    Ω_C = ∫ω_ijω^ij/p d^px .

We will refer to the first term on the right hand side of (<ref>) as a vortex stretching term. This term vanishes for p=2 (as does the last term on the right hand side of (<ref>)). 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 (<ref>). Since (<ref>) is less manageable than (<ref>), 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 <ref> 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 (<ref>) 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 (<ref>) in the limit of infinite Reynolds number.

More specifically, for p=2 we used



    f_x =E c_s δ f_x(x⃗),     
    	f_y =E c_s cos(2π n/L x)+δ f_y(x⃗),     
    	a = a_0 E

whereas for p=3 we used

    f_x    =     E c_s cos(2π n/L y) + δ f_x(x⃗),     
    	f_y  =E c_s cos(2π n/L z) + δ f_y(x⃗)
    
    	f_z    =     E c_s cos(2π n/L x) + δ f_z(x⃗),     
    	a  = a_0 E .


Here E is an energy scale which drops out of the equations of motion. The parameter a_0 in both (<ref>) and (<ref>) is constant in space and the δ f_i denote perturbations of the form

    δ f_i(x⃗)=∑_m⃗ A_i,m⃗Ec_scos(Δϕ_i,m⃗+2π(m⃗·x⃗)/L)

where m⃗, A_i,m⃗ and Δϕ_i,m⃗ are chosen from a random sample. The sum in (<ref>) 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

    Re = L/a_0 n ,
    	    
    	M=1/a_0 .

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 (<ref>). 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

    τ = L_0/U = Re/M^2ξ/c_s  
    	e=a U^2/L_0^p-1= M^p/Re^p-1E/ξ^p-1 ,

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





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

(<ref>),

    E_C(t,k) = ∂/∂ k∫_|k⃗'|≤ kd^pk'/(2π)^p|u⃗_p(t,k⃗')|^2

with u⃗_p(t,k⃗) = ∫ d^px√(p(t,x⃗))u⃗(t,x⃗)e^-ik⃗·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 <ref>. 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 (<ref>) a turbulent instability will not form.

Going back to figure <ref>, 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 <ref> 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 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 <ref>.  The data for the bottom right plot in figure <ref> 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 <ref>, 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 <ref> 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 (<ref>). (For example: if the initial shear flow was of the form f_x=E c_s cos(2π n/L y), then the associated unstable mode is u_x  ∼cos( 2π n/L(x ± y)). )


The authors of <cit.> 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 <cit.> 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 <ref>.






 §.§ 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 (<ref>). We present a detailed analysis of the power spectrum during the turbulent phase in appendix <ref>. 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 <ref> and <ref> 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 <ref>. 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 <cit.>. 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. 






 §.§ 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 ℝ^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 <cit.>. 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 times

    ϵ(x,y)     =    ℰ_0 
    
    	u_x(x,y)     =     U_x 0e^- (2π/L)^2tcos((2π/L)y+ϕ_1)
     
    	u_y(x,y)     =     U_y 0e^- (2π/L)^2tcos((2π/L)x+ϕ_2) ,

where ℰ_0, U_x 0 and U_y 0 are constants.
One can check that (<ref>) solves (<ref>) for low Mach number up to exponentially suppressed corrections. 
Typical flow of the form (<ref>) is depicted in the lower right corner of figure <ref>. A spectral analysis of the flow can be found in figure <ref>.

To get a quantitative handle on (<ref>) 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)^2t. 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.

See figure <ref>.


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 (<ref>). For the runs we have carried out in three dimensions with n=1 and n=4 (see (<ref>)) we find that the late time behavior takes the form

    u_x    =     U_x e^- (2π n/L)^2tcos(2π n/L y),     
    	u_y  = U_y e^- (2π n/L)^2tcos(2π n/L z)
    
    	u_z    =     U_z e^- (2π n/L)^2tcos(2π n/L x),    ϵ  = ℰ_0 .

where n is the 
mode number injected into the system in (<ref>).
Typical flow of the form (<ref>) is exhibited on the bottom row of figure <ref>. A k-space view of the velocity and energy density fields matching (<ref>) can be found in figure <ref>. A possible explanation for (<ref>) 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 <cit.> and there are numerous suggestions for identifying an appropriate geometric quantity which captures the turbulent behavior of the dual fluid <cit.>. In what follows we will focus on the work of <cit.> where the authors argue that the horizon power spectrum, 𝒜, to be defined shortly, is proportional to the energy power spectrum defined in (<ref>) and (<ref>). 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 

    ds^2=dt (-(1-a(ζ^μ)/R)dt+2/nRdR - 2f_a(ζ^μ)/nRdζ^a)+δ_ab/ndζ^adζ^b+1/ndχ⃗_⊥^2 .

The only null surface associated with the metric (<ref>) 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

    Θ_MN≡Π_ M^PΠ_ N^Q∇_Pn_Q ,     Π_ N^M≡δ_ N^M+ℓ^Mn_N

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 <cit.> we consider the rescaled traceless horizon curvature, defined by

    θ_ j^i≡√(γ/κ^2)Σ_ j^i ,     Σ_ j^i≡Θ_ j^i-1/d-2Θ_ n^nδ_ j^i,

where i,j run over the n+p spatial dimensions of the horizon, √(γ)≡√((g_ij)) is the area element on a spatial slice of the event horizon, and κ is defined by the geodesic equation n^M∇_Mn_Q=κ n_Q. 
The horizon curvature power spectrum is defined by

    𝒜( t,k ) ≡∂/∂ k∫_|𝐤'|≤ kd^pk'/(2π)^pθ̃_  j^*i( t,𝐤') θ̃_ i^j( t,𝐤') ,

where

    θ̃_ j^i( t,𝐤) ≡∫d^px θ_ j^i( t,𝐱) e^-i 𝐤· 𝐱

and p is the number of spatial dimensions where the dynamics take place. 

An explicit computation gives us

    n_M dx^M = dR - ∂_t a dt - ∂_b a dζ^b ,
    	    ℓ_M d x^M = - 1/nRdt  .

and 

    √(γ)=a(ζ^μ)n^-(n+p)/2 ,
    	    κ=a(ζ^μ)n^2/2

from which

    θ_ j^i = aδ^il/n^(p+n)/2(∂_j(f_l/a-∂_la/a)+∂_l(f_j/a-∂_ja/a))

follows. Substituting a=ϵ, f_a=ϵβ_a and taking the incompressible limit, we find

    θ_ j^i=ϵ/n^(p+n)/2(∂_jβ^i+∂^iβ_j) .


The Fourier transform θ̃ will involve a convolution of ϵ̃ and β̃ which are the Fourier transform of ϵ and β respectively, a 'la (<ref>). However, if we are working in the limit of small Mach number then ϵ is approximately constant and we find

    θ̃_ j^i ϵ/n^(p+n)/2(i k_jβ̃^i+i k^iβ̃_j) .

Inserting (<ref>) into (<ref>) and comparing to (<ref>) we find that at low Mach number,

    𝒜(k)/E(k) = 2 ϵ/n^(p+n) k^2

as predicted in <cit.>. If, on the other hand, the Mach number is not small a relation of the form 𝒜∼ k^2 E will not hold. 

Indeed, in figure <ref> we show typical results for the horizon power spectrum for high and low initial Mach number in two and three dimensions. As expected, 𝒜∼ k^2 E holds only for very low Mach number. A more refined analysis can be found in appendix <ref>.






§ 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 <cit.>. 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 <cit.>. 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 <cit.>, 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 <cit.>. 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.













§ ACKNOWLEDGEMENTS

The work of AY and ES is supported by an ISF grant and an ISF-UGC grant. The work of MR is supported by a Discovery grant from NSERC.








§ CRITICAL REYNOLDS NUMBER



Following <cit.> one may quantify the critical Reynolds number, above which the instability grows. 
Given the initial conditions (<ref>), we define a critical Reynolds number, Re^c, for these conditions, as the instantaeous Reynolds number at the moment where the amplitude of the unstable mode reaches a maximal value which is lower than the amplitude of the initial shear mode. That is, starting with (<ref>) we look for a maximum of the amplitude of the growing perturbed mode. We have carried out the analysis for several values of n (mode number), a_0 (Mach number) and L (box size). 

Our results for two dimensional flow can be found in figure <ref>. We find


    Re^c = 15.38 ± 0.81.

independent of n for L=10^6 and M=0.01 and

    Re^c = 15.33 ± 1.38.

for L=10^4 and M=1.
As described in the main text, once the Mach number is much larger than 1, it is difficult to generate a turbulent instability using the initial condition (<ref>).

Results for three dimensional flow can be found in figure <ref>. 
Here, we find that 

    Re^c = 14.7 ± 1.3.

for M=0.01 and L=10^4 and

    Re^c = 12 ± 1.3

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 (<ref>) 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.






§ DETAILED ANALYSIS OF THE POWER SPECTRUM



As discussed in section <ref>, 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 <ref> we have plotted the energy power spectrum for fixed initial Reynolds number and varying Mach number. In figure <ref> 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 <ref> and <ref>. In figure <ref> we have plotted E_C for a variety of initial Reynolds and Mach numbers and an initial disturbance (<ref>) with n=4. Figure <ref> 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 <cit.>.







§ ANALYSIS OF THE HORIZON AREA POWER SPECTRUM



In (<ref>) 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 (<ref>). Once the Mach number is small enough the flow becomes incompressible and we may approximate (<ref>) by (<ref>)  which yields 𝒜/E_C ∼ k^2. A numerical analysis of 𝒜/E_C for various initial Mach and Reynolds number is displayed in figure <ref>. As expected, the dependence of 𝒜/E_C on k seems to deviate from k^2 once the Mach number becomes large. 


To quantify the deviation of 𝒜/E_c from a k^2 behavior we go back to the compressible contributions to the traceless expansion (<ref>).
Let us denote the Fourier transform of a quantity X by ℱ(X). Then, in order for the incompressible terms in (<ref>) to dominate over the compressible ones, we need ∂/∂ k∫ dk^p|ℱ[∂_j(f_ℓ/a)]|^2 to be much larger than ∂/∂ k∫ dk^p|ℱ[∂_j(∂_ℓa/a)]|^2. In order quantify this relation we define the following matrix,

    M_jℓ = ∂/∂ k∫ dk^p|ℱ[∂_j(∂_ℓa/a)]|^2/∂/∂ k∫ dk^p|ℱ[∂_j(f_ℓ/a)]|^2,

which should become negligible whenever 𝒜/E_C ∝ k^2 is valid. In figure <ref> we have plotted M_ij for a two dimensional flow with Mach numbers ranging from 5× 10^-3 to 2. While M_j=ℓ remains small, M_j≠ l becomes non negligible at large Mach number.


A similar analysis for three dimensional systems can be found in figures <ref> and <ref> which display 𝒜/E_C for flows with initial conditions (<ref>) and n=4 or n=1 respectively. For both types of flows a k^3 scaling law is observed, even for the lowest Mach number which we could numerically access, M=0.1.

Extrapolating the data regarding the ratios of ∂/∂ k∫ dk^p|ℱ[∂_j(f_l/a)]|^2 to ∂/∂ k∫ dk^p|ℱ[∂_j(∂_la/a)]|^2 from figure <ref> we estimate that in order for the flow to become incompressible we need M ∼ 10^-4 with Re=750 and n=1 or Re=162 and n=4. Three dimensional flows with such a low Mach number are expensive and will not be covered in this work.







JHEP


