Holographic turbulence in Einstein-Gauss-Bonnet gravity at large D

We study the holographic hydrodynamics in the Einstein-Gauss-Bonnet (EGB) gravity in the framework of the large D expansion. We find that the large D EGB equations can be interpreted as the hydrodynamic equations describing the conformal fluid. These fluid equations are truncated at the second order of the derivative expansion, similar to the Einstein gravity at large D. From the analysis of the fluid flows, we find that the fluid equations can be taken as a variant of the compressible version of the non-relativistic Navier-Stokes equations. Particularly, in the limit of small Mach number, these equations could be cast into the form of the incompressible Navier-Stokes equations with redefined Reynolds number and Mach number. By using numerical simulation, we find that the EGB holographic turbulence shares similar qualitative feature as the turbulence from the Einstein gravity, despite the presence of two extra terms in the equations of motion. We analyze the effect of the GB term on the holographic turbulence in detail.


Introduction
In recent years, it has been found that black hole physics simplify significantly in the limit that the number of spacetime dimensions D is very large [1][2][3][4]. The key feature in the large D limit is that the gravitational field of a black hole is strongly localized near its horizon due to the very large radial gradient of the gravitational potential. By performing the 1/D expansion in the near region of the black hole, the Einstein equations can be reduced to several effective equations which capture the dynamics of the black hole. As a consequence one can solve these equations to construct various black hole solutions, including black string and black ring with different asymptotic behaviors and topologies. Furthermore one can study the linear stability of these solutions perturbatively to obtain the quasinormal modes, and even study the non-linear evolutions of the solutions numerically [5][6][7][8][9][10][11][12][13][14].
The fact that the dynamics of the black holes are captured by the effective equations of the horizons is reminiscent of the membrane paradigm [15][16][17][18][19][20], initialized in [21]. 1 One essential feature of the membrane paradigm is that the evolution of the horizon is like a viscous fluid. In [23] it was found that the large D effective equations for (AdS or asymptotically flat) black branes can be interpreted as the equations for dynamical fluid. For AdS black branes the transport coefficients were found to match well with the result obtained from the holographic study. More importantly, the study in [23] indicates that JHEP01(2019)156 the hydrodynamical gradient expansion is truncated at a finite order and the higher order transport coefficients are vanishing in the D → ∞ limit.
On the other hand, according to the AdS/CFT correspondence, the geometry of dynamical black holes in asymptotically anti-de Sitter spaces (AAdS) can be described in terms of a conformal fluid living on the conformal boundary of the spaces [24]. This fluid/gravity duality not only provides a geometric way to study fluid dynamics but also help to reveal new phenomena that never occurs in gravity. For instance, it is well-known that the turbulence is a ubiquitous property of the fluid if the Reynolds number is sufficient large. The presence of the turbulence in hydrodynamics indicates a similar turbulent behavior should appear for the AdS black holes. As expected, it was found that by numerical simulation a perturbed black brane in AAdS can exhibit turbulent behavior when the Reynolds number of the fluid counterparts is sufficiently large [25][26][27]. Within the fluid/gravity duality it is possible to use the geometric tool to study the turbulence, such as the geometric interpretation of turbulence [25]. Therefore, the natural questions are what is the relation between the large D effective equations for the AdS black branes and the large D limit of the conformal fluid living on the AdS boundary, and can the fluid equations exhibit turbulent behavior in a certain regime? Indeed, as demonstrated by a recent work [28], the large D effective equations for the AdS black branes are equivalent to the large D limit of relativistic hydrodynamics which is naturally truncated at the second order in derivatives. Actually, in deriving the large D limit of the relativistic hydrodynamics, one need to decompose D = n + q + 1 and consider the dynamical variables depending only on q + 1 coordinates. Moreover one has to rescale the time and space coordinates appropriately such that the hydrodynamics on q + 1 dimensional system is simplified significantly in the large n limit. Such scaling laws for the coordinates are precisely the ones used in studying the large D limit of the gravity. More interestingly the effective equations for the black brane in the large D limit are exactly the same as the hydrodynamic equations of motions in a fluid frame. These equations are a variant of the compressible version of the non-relativistic Navier-Stokes equations. Then by using analytic and numerical techniques the authors in [28] analyzed two and three-dimensional turbulent flow of the fluid in various regime and the relation with the geometry of the black branes. In this paper we would like to extend the study to the holographic hydrodynamics in the framework of Einstein-Gauss-Bonnet (EGB) gravity at large D.
The Einstein-Gauss-Bonnet gravity provides a nice platform in studying the holographic hydrodynamics. According to the holographic dictionary, the higher-derivative terms in gravity may come from the stringy correction or the string interaction. The Gauss-Bonnet term, including the quadratic terms of the curvature tensors, appears as the leading order correction in the low energy effective action of the heterotic string theory [29,30]. Another appealing feature of the EGB gravity is that its equations of motion remain second order such that the fluctuations around the vacuum do not have ghostlike mode. The EGB gravity has been well studied in the holographic hydrodynamics, in particular on the Kovtun-Son-Starinets(KSS) viscosity bound [31] η s ≥ 1 4π .

JHEP01(2019)156
In [32,33], it was shown that the KSS bound was violated in the EGB gravity where λ GB is the Gauss-Bonnet coefficient in five dimensions. The causality constraints require that in five dimensions [34][35][36] In general dimensions D, the viscosity satisfies (1.5) The above constraints comes either from the requirement that the energy flux in the scattering experiment is everywhere positive, or equivalently from the causal propagation of fluctuations.
The large D study of the EGB gravity [39][40][41] shares many similar features with that of the Einstein gravity. For example, there are also two classes of quasinormal modes for the EGB black holes: one consists of the decoupled modes, which characterize the dynamics of the black hole, and the other one consists of featureless non-decoupling modes [39]. For the EGB black strings (branes), the final states of the non-linear evolution are nonuniform black strings if the strings are thin enough. This is qualitatively the same as the behavior of black strings in the Einstein gravity. Besides, the large D effective equations for the (asymptotically flat) EGB black branes can be interpreted as the fluid equations as well [41]. So for the AdS black branes in the EGB gravity one may expect that the large D effective equations are also of the equations for a dynamical fluid, which allows us to study the holographic turbulence in the EGB gravity.
In the paper, based the work of [28] we will study the holographic turbulence in the EGB gravity by using the large D expansion method. Unlike the case in the Einstein gravity, up to second order the hydrodynamics dual to the EGB gravity is unknown. In this work, we focus on the holographic EGB hydrodynamics. As we show in section 2, if we take the same scalings for the coordinates as those used in constructing the EGB black strings [41] at large D, then we obtain the effective equations which have the form of the hydrodynamical equations naturally truncated at second order in derivatives. Moreover, the transport coefficients are the same as those obtained from AdS/CFT [37].
In section 3 we give an analytical discussion on the large D EGB fluid flows. We show that in the small Mach number limit, the EGB equations could be reduced to the incompressible Navier-Stokes equations with modified Reynolds and Mach number. However, if the Mach number is not small, the extra terms coming from the GB term may play an important role, especially in 2D flow. Then in section 4 we numerically solve the equations JHEP01(2019)156 of motion and analyze the effect of the GB term on the turbulence. We surprisingly find that the extra terms play negligible role in the evolution of the turbulence. We find that the 2D EGB turbulence has qualitatively similar behavior as the one from the Einstein gravity. For example, we notice an inverse and direct cascade in the turbulence, with the energy spectrum of the direct cascade obeying a k −4 power law. In section 5 we verify that the relation between the horizon curvature power spectrum and the hydrodynamic energy power spectrum proposed by [25] holds for the EGB gravity. We end with a summary and some discussions in section 6.
2 Equations of motion of relativistic hydrodynamics

1/n expansion of the EGB equations
The action of the D dimensional EGB gravity with a negative cosmological constant Λ = −(D − 1)(D − 2)/2 (we have set the radius of AdS space to unity) is given by where α is the GB coefficient. Here we follow the conventions in [39][40][41]. The GB coefficient α is related to the one used in [37] by As the GB term appears in the low energy effective action in the heterotic string theory [30], the GB coefficient is usually taken to be positive. However, just taking the EGB theory as a gravity theory, there is no reason to restrict the GB coefficient to be positive. Actually from the causal constraints [37] it can be negative as well.
There are more stringent constraints on the value of the GB coefficient than (1.5) by studying the scattering of a probe graviton from the shock wave. As shown in [38], in a theory with higher curvature terms, when the impact parameter is small enough, there could be time-advance which overwhelms the Shapiro time-delay and leads to the violation of causality. The higher derivative terms should be small in order not to violate the causality. For the GB term, in the flat spacetime, the GB coefficient should be at the same order of l 2 p . In the AdS spacetime, there is a AdS radius and the GB coefficient in units of AdS radius should be small as well. For example, in the well-known AdS 5 /SYM 4 correspondence, the GB coefficient should be smaller than (λ) −1/2 up to a undetermined factor, where λ is the t' Hooft coupling. The picture is that the pure gravity with higher derivative terms could lead to causality violation, but at a very short scale the other string states may help evade the violation.
According to the study in [38], it is impossible to have a consistent holographic theory of quantum gravity with a finite α unless an infinite tower of higher spin fields are included.

JHEP01(2019)156
More precisely, the GB coefficient is closely related to the scale that the higher spin fields appear at, and so should be very small, depending on the details of quantum gravity. In the large D expansion, the real parameter appearing in the discussion is λ D GB = α(D − 3)(D − 4). In the following discussion, we consider the GB coefficient λ D GB satisfying the constraints (1.5). It is finite, while the coefficient α is actually of order 1/D 2 , which could be very small. We expect that the value of α is not conflict with the study in [38]. 2 From the action, we obtain the equations of motion for the metric where The black brane solution in EGB gravity can be written as [42]

8)
r + is the horizon radius, andα is exactly the Gauss-Bonnet coefficient used in [37]. Note that in the above expression L c is introduced as the effective radius of the AdS space in EGB gravity, since as r → ∞, h(r) → 1/L 2 c , with (2.9) From above we can find that there exists an upper bound for the GB coefficientα, i.e.
beyond which the theory is not well-defined. On the other hand, the causality constraint (1.5) requires that in the large D limit, It turns out that the upper bound is the same as (2.10).
Using the Eddington-Finkelstin coordinates the solution (2.7) can be written as

JHEP01(2019)156
L c can be absorbed into the redefinition of x i , i.e. x i → x i /L c , such that L c does not appear in the metric explicitly. This form of the metric leads to the following metric ansatz where a, b, c = 1, · · · q, x is an n dimensional vector and D = n + q + 2. Note that as in [8,28] we keep q finite and let n → ∞. In this case most of the spatial directions remain translational invariance, so that we only study the deformations of the black brane that depend only on a finite number of coordinates z a . Moreover, in the following we use 1/n rather than 1/D as the expansion parameter. In order to perform the 1/n expansion properly we need to specify the large D scalings of the metric functions. According to the properties of the boundary hydrodynamics [37], the speed of sound of long-wavelength perturbations is of O(1/ √ n), this indicates that to capture the physics of the long-wavelength perturbations, we should rescale In addition, we consider small velocities O(1/ √ n) along the brane directions. Therefore, the large n scalings of the metric functions are respectively (2.14) These scalings are exactly the same as the one in discussing the black string in the large D EGB theory [41]. We can rewrite the metric ansatz (2.13) as where such that the derivative with respect to R is finite in the large D limit, where r 0 is a horizon length scale which can be set to be unity r 0 = 1. To solve the EGB equations we need to specify the boundary conditions at large R, which are On the other hand, the solutions have to be regular at the horizon. At the leading order of the 1/n expansion, the EGB equations only contain R-derivatives so they can be solved by performing R-integrations. After imposing the boundary conditions the leading order solutions are obtained as 20) and to simplify the expressions we have introduced the quantities b and β As shown in [23] and [41], the 1/n terms inG ab are obtained at the next-to-leading order in the 1/n expansion of the EGB equations. It must be included since it also appears in the EGB equations at the leading order of the 1/n expansion. m(v, z a ) and p a (v, z a ) are the integration functions of R-integrations of the EGB equations. At the next-to-leading order of the 1/n expansion, m and p a must satisfy the effective equations In the limitα → 0, i.e. β → 1, the above equations reproduce the ones in the Einstein gravity [23].

Linear analysis
Considering a small perturbation around the static uniform black brane solution, i.e. m = m 0 = const., p a = 0, and plugging these into the effective equations (2.23) and (2.24), we can read the quasinormal mode frequencies of the shear mode and the sound mode.
Shear mode. The frequency is As we will see in the next subsection, this dispersion relation is related to the transport coefficient of the viscous hydrodynamics.

JHEP01(2019)156
Sound modes. The frequencies are It is easy to see that the perturbation is always stable. Up to the leading order of k, one finds ω ± = ± 2 β+1 k. As we will see in the next subsection, this coefficient corresponds to the speed of sound of the long-wavelength perturbation of the fluid.

Dynamical fluid
In [23], it was shown that the effective equations for the black branes can be interpreted as the equations for a dynamical fluid. This turns out to be true for the asymptotically flat EGB black branes as well [41]. Here we show that the effective equations (2.23) and (2.24) can be actually transformed into the fluid equations. Firstly, we introduce v a by with m and v a being taken as the energy density and the velocity of the fluid flows. Furthermore, we notice that (2.24) can be written in terms of v a as This equation is for the momentum conservation, with τ ab being the stress tensor. Therefore, we may interpret the effective equations as the equations of motion for non-relativistic, 3 compressible fluid naturally truncated at second order in derivatives. As shown in [28], the Einstein equations at large D are equivalent to the large D hydrodynamics within the fluid/gravity duality. This statement holds for the EGB gravity as well: the EGB equations at large D are equivalent to the large D hydrodynamics describing the conformal fluid living on the conformal boundary of the AdS space. Although the relativistic hydrodynamics that dual to AdS black branes in the EGB gravity up to second order is unknown at present, we can still find some hints by comparing the properties of the fluid flows (2.30) and (2.31) with the results obtained from AdS/CFT [37]. From (2.32) we find the pressure In fact, since the physical velocity for the fluid flow is v i = v i / √ n, the fluid flow is non-relativistic in the large n limit. In [28], due to the extra factor 1/ √ n carried by the velocity filed, the large D hydrodynamics is shown to be non-relativistic.

JHEP01(2019)156
which gives the equation of state of the black string. The speed of sound of the longwavelength perturbations is then Due to the scaling relation we used in the metric ansatz, the physical speed of sound is c s = c s / √ n, which is small in the large n limit. Moreover, the shear viscosity is then the ratio of the shear viscosity to the entropy density is given by with the entropy density s = 4πm. Both this result and the speed of sound are in accord with the results found previously in [37] by taking the large D limit and using our convention. 4 3 General analysis of the large D EGB fluid flows In this section we give a general discussion on the holographic EGB fluid flows. We will relate the equations of motion (2.30) and (2.31) to the compressible Navier-Stokes equations.
In terms of In fact, if we use the convention in [37] and for simplicity only consider one momentum along the brane direction, then we find the effective equations In terms of pz = 1 β mvz + 1+β 2β ∂zm, the effective equations have the similar form as (2.30) and (2.31). In this case the ratio of the shear viscosity to the entropy density is found to be η s = 1 4π β 2 +1 2 which is exactly the same as the large D limit of (1.4). 5 Note that ∇ f · ∇m is different from ∇m · ∇ f : the former one is written as ∂af b ∂ b m but the latter one

JHEP01(2019)156
It is convenient to introduce the following dimensionless variables, where L 0 is a characteristic length scale of the system, U is a characteristic velocity and E is a characteristic energy density. In terms of these dimensionless variables, (3.2) and (3.3) become of the forms

5)
are the Reynolds number and the Mach number. Although the above equations are slightly more complicated than the ones appearing in the Einstein gravity, they can actually be taken as a variant of the compressible version of the non-relativistic Navier-Stokes equations. This can be seen by taking the incompressible limit and redefining the Reynolds number and the Mach number. Let us first consider the simple case when the Mach number is small, and the fluid flow is nearly incompressible. We can expand u and in series of M , i.e.
Then from the above equations we find that (0) and (1) are constant, and u (0) and (2) satisfy the following equations, with (0) being the density and (2) being the pressure. If we redefine the Reynolds number and the Mach number as then instead of (3.9), u (0) and (2) satisfy where η and c s are the shear viscosity and the sound velocity. From the relations (2.35) and (2.34), it is easy to see that they are exactly the same as (3.10).
Up to now, we have found that in the small Mach number limit, the equations of motion for the holographic fluid flow in the EGB theory take exactly the form of the Navier-Stokes equations for an incompressible fluid. The impact of the GB term on the fluid flow is totally encoded in the redefinition of the Reynolds number and the Mach number. The incompressible flow has been well studied (see for example in the textbook [43]), and here we just give a brief review on the main results for the freely decaying fluid turbulence.
One may define the kinetic energy and the enstrophy by where is the vorticity 2-form. It is illuminating to write the kinetic energy and the enstrophy in terms of the energy spectrum E(k), The energy spectrum E(k) provides a measure of how the energy is distributed across the various eddy size ∼ 1/k. From the equations (3.11), one can get the evolution equations for E I and Ω I where σ (0) and P I are the non-relativistic shear tensor and the palinstrophy, defined by As Ω I ≥ 0, the energy E I is always decreasing. For the enstrophy equation (3.17), the two terms on the right-hand side correspond to the generation of the enstrophy via the vortex line stretching and the destruction of the enstrophy by the viscous forces. In the large Reynolds number limit, there is an approximate balance between the production and the

JHEP01(2019)156
viscous dissipation of the enstrophy. Thus, to keep Ω I a constant, the energy spectrum has to be distributed in such a way that the energy will flow from the lower momentum modes to the higher momentum modes. Under Kolmogorov's Similarity Hypothesis, the energy spectrum behaves as which is known as the Kolmogorov's five-thirds law, where ε denotes the rate of energy dissipation at small scales. The turbulence in two spatial dimensions is a little bit special. In this case, the vortex stretching term vanishes for an incompressible flow, and the enstrophy declines monotonically as the flows evolves so that it is bounded from above by its initial value. In the large Reynolds number limit, the energy is conserved. It turns out that there could be another kind of energy cascade, the so-called inverse energy cascade, in which the energy flow to the lower momentum modes. Besides, before the turbulence adjusts to its fully developed state, there is a direct cascade of the enstrophy from the large scales down to the small scales due to the continual filamentation of the vorticity. In the inertial range, according to Batchelor [44], the energy spectrum for the direct cascade could be where δ is the rate of dissipation of the enstrophy. This is the two-dimensional analogue of the Kolmogorov's five-third law. However, unlike Kolmogorov's five-thirds law, the direct cascade law (3.20) is far from being well accepted. In fact, the numerical simulations suggest that E ∼ k −n , where n is typically a little larger than 3. For more details, one can see [43] and the references therein. On the other hand, according to Kraichnan [45], the possibility of an inverse cascade depends on the relative strengths of the non-linear interactions between scales. If the input of the energy is not sustained for sufficiently long time, the inverse cascade cannot develop and there appears only the direct cascade in the flow. Thus, unlike the case of the forced two-dimensional flows, the inverse cascade with a Kolomogrov-like energy scaling ∼ k −5/3 is puzzling in the decaying case [46]. It is illuminating to study the equations (3.5) and (3.6) beyond the small Mach number limit. We may define the energy and enstrophy via Here and in the following we need the vorticity and the shear tensor From the equations (3.5) and (3.6), we have ∂ v Ω C = 1 ω ij ω jk σ i k − δ i k ∇ · u d q x + (terms proportional to 1/R e ) . (3.24)

JHEP01(2019)156
There are two remarkable points on these two equations. The first is that the effective Reynolds number R e and Mach number M appear in the equations. This suggest that the effect of the GB term could be largely encoded by these two parameters. The second is that the first term on the right hand side in the equation (3.24) could be taken as a vortex stretching term, by which the energy is passed down the cascade to the small scales. The other terms are proportional to 1/R e , and play the similar role as the palinstrophy that characterizes the viscous dissipation at small scales. If the Mach number is not small, the equations (3.5) and (3.6) are more complicated than the ones appearing in the Einstein gravity. Actually there are two more terms proportional to (β − 1) in (3.6), whose physical implication is not clear. However, both terms are inversely proportional to 1/R e , and consequently when the Reynolds number is very large, they may not play significant role on the dynamics of the flow except the effect on the viscous dissipation. In three spatial dimensions, if the Reynolds number R e is large, the Kolmogorov cascade should appear in the turbulence as the vortex stretching term dominates. In two spatial dimensions, the vortex stretching term is vanishing. Therefore one may expect the inverse energy cascade in this case. However, due to the presence of two extra terms induced by the GB term, it is not clear if the inverse energy cascade do appear or present different feature. In order to understand the effects of the GB term, we need to do numerical simulation.

Numerical study of turbulent flows
In this section by numerically solving the equations of motion (2.23) and (2.24) we will give a detailed analysis for the turbulent flows in the EGB theory. In particular, we focus on the distinctions between the Einstein gravity and the EGB theory.
As in [28], we consider the flows in a toroidal domain of size L. To solve the equations of motion (2.23) and (2.24) we use the Fourier spectral method in the spatial directions and the fourth order Runge-Kutta method for time evolution. The initial conditions are taken to be shear flows with small random perturbations, since in the inviscid case, i.e. η → 0, they are stationary solutions of the fluid equations. Explicitly, the initial conditions are for q = 2, and

JHEP01(2019)156
In the above expressions, A i, m denotes the amplitude of the perturbation and is chosen from a uniform distribution ranging from 0 to a small value. The phase δφ i, m is chosen from a uniform distribution ranging from 0 to 2π. m is the wavenumber in unit of 2π/L. For more details about the initial conditions, one can see [26][27][28].
Due to the decay caused by the shear viscosity, the shear flow we study is not steady. The Reynolds number, as a useful dimensionless quantity to predict the stability of the steady flow, may still be useful in our case. Thus the stability of the flow depends only on the instantaneous value of the Reynolds number [27].
For our flow, we follow the definition in [28] where the characteristic velocity appearing in (3.10) is given by 6 U characterizes the velocity fluctuation. For the initial conditions of the flow, it is easy to obtain the Reynolds number and the Mach number where we have used L 0 = L n I . From the above expressions we can see the effect of the GB term on the initial conditions. On the one hand, the presence of the GB term always lowers the initial Reynolds number. On the other hand, for a positive GB coefficientα (i.e. 0 ≤ β < 1), the presence of the GB term always lowers the Mach number, but for a negativeα (i.e. 1 < β ≤ 2) it would always lead to a larger Mach number. Similar to the steady flow case, we expect that if the initial Reynolds number is sufficiently small the flow is stable against small perturbations, and the final state is laminar. While for a sufficiently large initial Reynolds number, the flow is unstable and we should expect the emergence of the turbulent flow. Between these two extremes, for an intermediate initial Reynolds number, there exists a critical Reynolds number [27] beyond which the transition from the laminar flow to the turbulent flow may occur, as we will see in the following.
As a preview, we show the typical simulations of two dimensional fluid flow in figures 1, figure 2 and figure 3 and the simulations of the three dimensional fluid flow in figure 4 with large Reynolds numbers for different values of β. Qualitatively, the behavior of the fluid flows in the EGB theory is very similar to the one in the Einstein gravity. It can be described by three phases: the initial growth of the instabilities, the turbulent regime in which the inverse/direct energy cascade is observed, and the late time decay into an equilibrium. An obvious distinction between these two theories is that, for a positive GB coefficient the time for the turbulence to emerge is later, but for a negative GB coefficient the time is earlier, under the same initial conditions (L, m 0 and n I ). This is largely because the presence of the positive GB coefficient lowers both the initial values of the Reynolds number and the Mach number as can be see from (4.5), thus it costs more time for the flow to reach the turbulent phase. In contrast, the presence of the negative GB coefficient lowers the initial value of the Reynolds number but increase the initial value of the Mach number. As a consequence, it costs less time for the flow to reach the turbulence phase. Figure 1. The vorticity field ω = ∂ x p y − ∂ y p x for a two dimensional fluid flow in the Einstein gravity (β = 1), with initial Reynolds number R e0 = 937.5, m 0 = 2, initial mode n I = 8, and A = 10 −5 at various times. We can see the initial growth of the instability, and later on the formation of filaments which indicates the direct cascade, and then the turbulent regime at which a large scale structure is formed and the inverse energy cascade is apparent. Finally, we can see the slow decay of two counter-rotating vortices. Note that here and in the following we employ τ as the units of time with τ = L 0 m 0 .    the Einstein gravity, a positive GB coefficient has a larger evolution rate while a negative GB coefficient has a smaller evolution rate. This may be ascribed to the distinction of the dynamical equation. After we adopt the definition (3.10), the dynamical equation has two extra terms proportional to the inverse of the Reynolds number. These two terms induce viscous forces and affect the evolution of the fluid flow. As we will see in the next subsection, from the linear analysis, the positive GB coefficient has the smaller viscous effect and the larger evolution rate, while the negative GB coefficient has the larger viscous effect and the smaller evolution rate.

JHEP01(2019)156
In the following we will pay our attention to the effects of the GB term on the fluid flows in these different stages.

Stability of shear flow
In this subsection we analyze the early stage of the flow. Due to the restriction of our computational resource, in the following we mainly show our results of the two dimensional flows. As we mentioned before, depending on the initial Reynolds number, the small random perturbations of the shear flow grow or decay with time. The growth of the perturbations can be quantified by the growth of p x field [27] or as in [28] by the energy spectrum E C whose definition is 7 In what follows we will take either one of the two quantities when necessary. 7 The definition of the energy spectrum is a little different from the one in [28]: here we use p rather than u in EC . We find that this would be easier for numerical analysis. It works pretty well.  At low values of R e0 the flow is laminar, as shown in figure 6. We can see that the L 2 norms of the components of p i ,

JHEP01(2019)156
decay exponentially both for the Einstein gravity and the EGB theory. From the linear analysis in section 2.2, we know that the presence of the positive GB coefficient lowers the quasinormal mode frequency (2.27) and so the viscosity (2.35), thus the L 2 norms of p i decay slower than the ones in the Einstein gravity. In contrast, the presence of the negative GB coefficient increases the decay rate so that we have steeper lines in the right panel of figure 6. Moreover, we find that for larger Mach numbers the results are essentially the same as the ones of small Mach number case. At a large enough R e0 we observe the turbulent instabilities, as shown in figure 7. At the first glance the energy spectrums in the Einstein gravity and the EGB theory are very similar. In both cases we find the decay of the initial disturbance at k = 8 × 2π L and then an unstable mode emerging at a lower wavenumber k ∼ 3 2 × 2π L , which may indicate the appearance of the inverse cascade. However, the direction cascade is not apparent, this may attribute the insufficient precision of our numerical simulation. As we mentioned before since the flow in the EGB theory has a lower Reynolds number, the initial JHEP01(2019)156 For the positive GB coefficient (β = 0.5), the decay rate of log U is smaller and the maximum of log L 2 (p x ) appears at later time. In contrast, for the negative GB coefficient (β = 2), the decay rate is larger and it takes less time for log L 2 (p x ) to reach the peak.
disturbance decays slower and the new unstable mode grows slower, comparing with the Einstein gravity.
For an intermediate value of the initial Reynolds number, as shown in figure 8, the initial small perturbations characterized by L 2 (p x ) grow exponentially to reach a maximum that is smaller than L 2 (p y ), and then decay exponentially. From figure 8 we can find the exponential decay of the characteristic velocity U defined in (4.4). From the definition (3.10) the Reynolds number varies with time. In this case we can define a critical Reynolds number R c as in [27]. Initially, we have R e0 > R c , thus the unstable mode L 2 (p x ) grows exponentially. As time progresses, R e decreases, then when R e < R c , L 2 (p x ) is not growing any more, instead it decays exponentially. Therefore, at the peak of L 2 (p x ) we may identify R e = R c .
We have identified the critical Reynolds number R ec for different initial Reynolds numbers, different initial Mach numbers and different β's. Firstly, for the Einstein gravity β = 1, we find that for a given Mach number the larger is initial Reynolds number, the larger is the critical Reynolds number, as shown in table 1. Secondly, for a given initial Reynolds number, we find that the larger is the Mach number, the smaller is the critical Reynolds number, as shown in table 2. Finally, we can see from table 3 that for the fixed initial Reynolds number and Mach number, the critical Reynolds numbers of different GB coefficients are pretty close.

Turbulent regime
If the initial Reynolds number is sufficiently large, the initial instability grows until it reach the same amplitude as that of the initial shear mode, then the fluid goes into the turbulent regime, as shown in figure 9. In this stage, the small eddies merge into vortices, which continue to merge into increasingly large vortices until two counter-rotating vortices are formed, as shown in figures 1 and 2. This process is the so-call inverse energy cascade. Besides, there is a direct cascade caused by the formation of the filaments of the vorticity. These two cascades can be characterized by the shape of the energy spectrum E C in the inertial range. The detailed analysis of the energy spectrum during the turbulent phase can be found in [28]. Here we only want to emphasize that for the two dimensional fluid flow in the EGB theory the energy spectrum of the direct cascade satisfies a k −4 power law for a small Mach number and a large Reynolds number, as shown in figure 10. Similar    Table 3. The critical Reynolds numbers of the fluid flow in the EGB gravity with different GB coefficients, here the initial Reynolds number is fixed R e0 = 160. to the Einstein gravity, the k −5/3 law for the inverse cascade is absent from the energy spectrum analysis in the case of the EGB theory. From the analysis in [46], it was found that the range at which the k −5/3 law applies is very short and the distinction from the k −3 law demands high numerical precision.

JHEP01(2019)156
For three dimensional fluid flow, the celebrated k −5/3 power law has been found in [28], even the Mach number is large. Due to the limited computational power, we are not able to carry out the numerical analysis on three dimensional fluid flow. Nevertheless, we expect the same energy cascade happens in the EGB turbulence as well.

Late time decay
The final phase of the turbulent flow is the formation of two counter-rotating vortices, as shown in figures 1, 2 and 3. We find that the late time decay takes the following form m(x, y) = m 0 , where m 0 , p x0 and p y0 are constant. The above behavior can be understood as follows.
Since in this stage the velocity field becomes small, thus the linear analysis for the fluid flows is applicable. The decay rate is just the damping part of the quasinormal mode of the black brane. On the other hand, due to the inverse cascade at late time the lowest mode dominates the flow, as shown in figure 11.

Geometric interpretation of turbulence
In this section, we would like to see if the horizon power spectrum defined in [25] is applicable to the turbulence of the holographic EGB fluid.

JHEP01(2019)156
Thus in the low Mach number case the emergent k 2 relation holds for EGB gravity as well, with the coefficient embodying the effect of the GB term. As in the case of the Einstein gravity [28], if the Mach number is not small the simple relation A ∼ k 2 E C is expected to be broken.

Summary
In this paper we studied the holographic hydrodynamics of the fluid flow in the Einstein-Gauss-Bonnet gravity by using the large D effective theory. After performing the 1/D expansion and integrating out the radial direction of the EGB equations, we obtained the effective equations for the EGB-AdS black branes. We found that the effective equations can be easily turned into the form of the equations for a dynamical fluid, the properties of which are consistent with the results obtained from the other studies in AdS/CFT. In the small Mach number limit, we found that the holographic hydrodynamic fluid equations in the EGB gravity could reduce to the Navier-Stokes equations for the incompressible fluid, with the modified Reynolds number and the modified Mach number Therefore, as in the case of the Einstein gravity, the Kolmogorov scaling laws are expected to emerge for three dimensional fluid flow, and an inverse/direct cascade should exhibit for two dimensional turbulent fluid flow. Compared to the hydrodynamical equations derived in the Einstein gravity, there are two extra terms in the equations of motion (3.6), being proportional to 1/R e . In three dimensions, the vortex stretching term has the dominant influence on the kinetic energy transfer between different scales, and it leads to the Kolmogorov scaling. In this case, the two extra terms play minor role. In two dimensions, as the vortex stretching term is vanishing, and the evolution of the enstrophy is determined by the palinstrophy-like terms, which include the contribution from the extra terms. One may expect that the 2D turbulence may present different feature due to the extra terms. Surprisingly, we found quite similar behaviors for the 2D holographic turbulence in the EGB theory as the one in the Einstein gravity. We generated the turbulent flow in a toroidal domain by numerically solving the equations of motion. The evolution of the flow can be divided into three stages: the initial growth of the instabilities, the turbulent regime and the late time decay. Qualitatively, in all these stages the behavior of the EGB fluid flow is similar to the one in the Einstein gravity: the initial instability grows accompanied with the formation of filaments, then the flow enters the turbulence regime in which the large scale structure is formed, and finally the formed counter-rotating vortices decay slowly. The formation of the large scale structure indicates an inverse cascade and the formation of filaments indicates the direct cascade. We found that the energy spectrum of the direct cascade obeys a k −4 power law but the k −5/3 law for the inverse cascade is absent. A possible reason for this could be the insufficient numerical precision [46].

JHEP01(2019)156
Despite the qualitative similarity, the EGB turbulence presents some different features. The effect of the GB term is twofold. On the one hand, the GB term affects the initial Reynolds number and Mach number. If the initial conditions are fixed, for a positive GB coefficient, both the Reynolds number and Mach number are smaller than those of the fluid in the Einstein gravity, such that the EGB fluid needs more time to reach the turbulent phase. However, for a negative GB coefficient, the Reynolds number is smaller but the Mach number is larger, so that the fluid needs less time to reach the turbulent phase. On the other hand, if we keep the initial Reynolds number and Mach number fixed, the extra terms in the dynamical equation affect the dissipation rate and so the evolution rate of the fluid flow. Moreover, we found that the critical Reynolds number is not sensitive to the value of the Gauss-Bonnet coefficient.
Due to the simplification at large D, we showed analytically that in the regime of small Mach number the proposed relation between the horizon curvature power spectrum and the hydrodynamic energy power spectrum A/E C ∝ k 2 [25] still holds in the EGB gravity.
In this work, we only did spectrum analysis for the 2D holographic EGB turbulence. It would be interesting to consider the 3D case. Even though a Kolmogorov cascade is expected in 3D case, it would be nice to show it explicitly.