Linearized fluid/gravity correspondence: from shear viscosity to all order hydrodynamics

In ref. \cite{1406.7222}, we reported a construction of all order linearized fluid dynamics with strongly coupled $\mathcal{N}=4$ super-Yang-Mills theory as underlying microscopic description. The linearized fluid/gravity correspondence makes it possible to resum all order derivative terms in the fluid stress tensor. Dissipative effects are fully encoded by the shear term and a new one, emerging starting from third order in hydrodynamic derivative expansion. In this work, we provide all computational details omitted in \cite{1406.7222} and present additional results. We derive closed-form linear holographic RG flow-type equations for momenta-dependent transport coefficient functions. Generalized Navier-Stokes equations are shown to emerge from the constraint components of the bulk Einstein equations. We perturbatively solve the RG equations for the viscosity functions, up to third order in derivative expansion, and up to this order compute spectrum of small fluctuations. Finally, we solve the RG equations numerically, thus accounting for all order derivative terms in the boundary stress tensor.


Introduction
In heavy ion collisions at RHIC and LHC, a novel state of QCD matter, quark-gluon plasma, is created. The quark-gluon plasma produced at RHIC was discovered to behave like a nearly perfect fluid reflecting strongly coupled regime of QCD [2,3], with relativistic hydrodynamics found to provide an accurate description of the plasma fireball expansion. The hydrodynamic evolution of the quark-gluon plasma is characterized by a set of transport coefficients, which have to be computed from the microscopic QCD. However, strongly coupled nature of this system prevents from a first principe analytic calculation of these coefficients. While lattice methods are quite reliable in studying QCD thermodynamics [4][5][6], they usually fail in extracting transport coefficients due to limited applicability to realtime dynamics. Therefore, various microscopic models are indispensable to understand transport properties of this QCD matter.
Fluid dynamics [7,8] is an effective description of most interacting quantum field theories at long wavelengths. The description is of statistical nature: it captures collective dynamics of a large number of microscopic degrees of freedom. The collective variables suitable for such a description are local densities of conserved charges, local fluid velocity and temperature. Hydrodynamic equations are local conservation laws for corresponding currents, which are specified via constitutive relations. As a low energy effective field theory, fluid dynamics describes near-thermal equilibrium systems and is naturally defined in terms of derivative expansion of the local fluid mechanical variables. Up to a finite number of transport coefficients, the derivative expansion at any given order is completely fixed by thermodynamic and symmetry considerations. Transport coefficients, such as viscosities and conductivities, must be determined from either experimental measurements or theoretical computations in the underlying microscopic theory.
The stress tensor T µν of a relativistic fluid is usually presented as a sum of two terms where T Ideal µν corresponds to ideal fluid dynamics and for conformal fluids has the form 1 where η µν is the four dimensional Minkowski metric, the fluid's four velocity is u µ and, in our notations, the temperature is T Diss µν accounts for dissipative effects and is chosen to be nonzero only for spatial components.
Up to the first velocity gradient it is given by the Navier-Stokes term 2 where η 0 is the shear viscosity coefficient and the tensor σ ij has the form 3 (1.5) Theoretical foundations of relativistic dissipative fluid dynamics are not yet fully established. The relativistic Navier-Stokes equations are a-causal and unstable [9][10][11][12]: the irreversible currents are linearly proportional to the thermodynamic forces, which have instantaneous influence on the currents, obviously violating causality. These problems can be solved by introducing retardation into the definitions of the irreversible currents [13,14], leading to equations of motion for these currents, which thus become independent dynamical variables. Theories of this type are known as causal relativistic dissipative fluid dynamics. Causality usually also implies stability [15]. To obtain a causal formulation, one needs to include higher order terms in the gradient expansion of the currents, in which case additional transport coefficients arise. Truncation at any fixed order would presumably lead to violation of causality that can be fully restored only at infinite order, which we refer to as an all-order gradient resummation. Non-trivial physical consequences imposed by causality of relativistic fluids were investigated in [9][10][11][12][15][16][17] and lead to certain constraints on possible values of higher order transport coefficients.
AdS/CFT correspondence [18][19][20][21] emerged over a decade ago as a standard tool and a model playground for addressing strongly coupled dynamics of gauge theories. The 1 The central charge of boundary CFT was factored out by a suitable normalization of the Boltzmann constant kB. 2 The bulk viscosity vanishes for the conformal fluids to be discussed below. 3 The linearization (2.7) was assumed in writing down the shear tensor σij.
AdS/CFT correspondence relates holographically a large N strongly interacting quantum field theory with dynamics of classical gravity in (asymptotically) AdS spacetime. Fluid/gravity correspondence [22][23][24][25] is a long wavelength limit of the AdS/CFT correspondence: it gives a map between black holes in asymptotically AdS spacetime and fluid dynamics of a strongly coupled boundary field theory. The most celebrated prediction of the fluid/gravity correspondence is a ratio of the shear viscosity η 0 to the entropy density s [22][23][24], η 0 s = 1 4π . (1.6) The ratio (1.6) is universal for a large class of strongly coupled gauge theory plasmas for which holographic duals are governed by Einstein gravities in asymptotically AdS spacetimes. Remarkably, this small ratio is quite close to the values extracted for QCD plasma from the RHIC experiments. Universality of this ratio was further proved in [26,27] and in [28] using black hole membrane paradigm. It was later found that university of the ratio (1.6) gets violated by either modifying Einstein gravity [29][30][31][32][33][34][35][36][37][38] or breaking isotropic invariance among spatial directions [39][40][41][42]. We refer the reader to [43] for a comprehensive review on early literature.
A specific realization of the fluid/gravity correspondence, and the one we will closely follow below, was established in [44]: it provides a systematic framework to construct a universal nonlinear fluid dynamics, order by order in the boundary derivative expansion. The stress tensor for the dual conformal fluid was explicitly constructed in [44] up to second order in derivative expansion, which is in perfect agreement with a general form of second order conformal hydrodynamics as analyzed in [45]. Computations of [44] were later generalized to conformal fluids in flat [46] and also weakly curved [47] background manifolds of various dimensions. Forced fluids in a weakly curved manifold were examined in [48] by studying long wavelength solutions of Einstein-dilaton gravity with negative cosmological constant. Further developments can be found in reviews [49,50] and references therein.
All order or resummed hydrodynamics is found to accommodate certain contributions, which are not present in a strict low momentum approximation. To avoid any confusion we thus clarify the terminology: by hydrodynamics we mean an effective theory given by a constitutive relation for the stress tensor in terms of temperature, fluid velocity and their gradients only. Navier-Stokes or any "unresummed" hydrodynamics involves only a finite number of gradients, while all order or "resummed" hydrodynamics means infinite number of gradients, but no other degrees of freedom.
The higher order derivative expansion generally includes two types of terms: nonlinear in the fluid velocity, like (∇u) 2 , and linear ones like ∇ · · · ∇u. The nonlinearities are significant when amplitudes of local fluid mechanical variables are large. However, even for fluid perturbations with small amplitudes, one can get large contributions from the linear terms when the momenta associated with the fluid perturbations are large. Given that these two types of terms are controlled by different parameters, it is possible to separate these contributions and to have the linear terms under theoretical control. We will focus on those linear terms in the rest of this paper.
To collect all order linear terms in a self-consistent manner, the shear viscosity η 0 was generalized into momenta-dependent function η(ω, q 2 ). This viscosity function is expressed in momentum space which follows from the replacement ∂ µ → (−iω, i q) in the linear gradient expansion of T µν . By postulating a constitutive relation for T µν in terms of the shear viscosity function η(ω, q 2 ), the authors of [52] attempted to extract η(ω, q 2 ) from thermal correlators of the stress tensor computed on the gravity side [56,57]. While certain progress was achieved in [52], its prime goal of complete determination of the viscosity function was not reached. It was realized that even in the case of linearized hydrodynamics, knowledge of retarded correlators contains insufficient information about all transport properties of the system. In ref. [1], we succeeded to completely solve this problem and below we provide all the details related to our computations.
In ref. [1], we reported progress, achieved via linearizing fluid/gravity correspondence, in consistently generalizing relativistic hydrodynamics to all orders. Upon linearization, perturbative computations of [44] can be straightforwardly extended to arbitrary order in the boundary derivative expansion. Our procedure is, however, slightly different from that of [44]. Particularly, the dissipative part in the stress tensor is collected in a unified way, rather than being determined order by order in derivative expansion. The Einstein equations in the bulk are split into two sets: dynamical equations and constraints. It turns out that in order to derive transport coefficient functions, it is sufficient to solve dynamical components of the bulk Einstein equations only. By solving only those we construct an "off-shell" fluid stress tensor. The remaining constraint components of the Einstein equations are equivalent to conservation laws of thus constructed fluid stress tensor and lead to generalized Navier-Stokes equations. It is worth emphasizing that the bulk dynamics is not absorptive, rather the bulk acts as a non-linear dispersive medium. Dissipative effects emerge via absorptive boundary conditions at the black hole horizon.
We find that the dissipative part of the stress tensor has the following form where π ij is a third order tensor structure and ζ(ω, q 2 ) is a new viscosity function, which in [52], apparently incorrectly, was argued to be zero. We here express the viscosity functions in momentum space but with tensors σ ij and π ij formulated as explicit derivatives of the fluid velocity. Later, the notations η(∂ v , ∂ 2 ) and ζ(∂ v , ∂ 2 ) will also be used to denote the viscosity functions. All these expressions are interchangeable under the above rule of replacement. We will be working using dimensionless units for all the momenta, choosing units such that the temperature is normalized to πT = 1. So all the physical momenta should be understood as in units of π T : ωπT and q i πT .
In the hydrodynamic limit ω, q i 1, η(ω, q 2 ) and ζ(ω, q 2 ) are expandable in power series and a perturbative analysis to be presented below reveals a few first terms in the expansions, where in our units, the first term in η corresponds to the ratio (1.6), while the second term accounts for the relaxation time [44,45,[58][59][60]. The remaining two terms in equation (1.9) are new third order transport coefficients. In section 4, η(ω, q 2 ) and ζ(ω, q 2 ) are computed numerically to all orders, and we will also present expansions of these viscosity functions up to fifth order. The remaining part of this paper is organized as follows. In section 2, we outline the derivation of the fluid dynamics from the bulk gravity. The main results are formulated as closed linear holographic RG flow-type equations for the generalized transport coefficient functions and generalized Navier-Stokes equations. We then compute dispersion relations for sound and shear modes. In section 3, to make comparison with previous studies in the literature, we perturbatively solve these holographic RG flow-type equations and obtain the fluid stress tensor up to third order in the derivative expansion. In section 4, we numerically solve the RG flow-type equations and obtain generalized transport coefficient functions, extending the perturbative analysis of section 3 to very large momenta. We also extract the viscosity functions via an approximate matching scheme, and find perfect agreement with the numerical results. Summary and discussion can be found in section 5.
2 Fluid dynamics from the bulk gravity

Linearized fluid/gravity correspondence
The fluid/gravity correspondence makes it possible to construct the fluid stress tensor and prove its conservation law (Navier-Stokes equations) by solving the bulk Einstein equations in asymptotically AdS spacetime, in the long wavelength limit. We start by considering a universal sector of the AdS/CFT correspondence: classical Einstein gravity with a negative cosmological constant in five dimensional spacetime, where the AdS radius is set to one for convenience. The Einstein equations which follow from the action (2.1) are We use upper case Latin indices {M, N, · · · } and lower case Greek indices {µ, ν, · · · } to denote bulk and boundary directions, respectively. Lower case Latin indices {i, j, · · · } will be used to specify spatial directions along the boundary.
Besides pure AdS 5 spacetime which is dual to the vacuum state of the boundary CFT, the action (2.1) also admits a 4-parameter family of Black Hole solutions, and f (r) = 1 − 1/r 4 . Hawking temperature of this black hole is which is identified as the temperature of the dual CFT. The operator P µν = η µν + u µ u ν acts as a projector onto spatial directions. Notice that so far the parameters β i and b are held constant, so that the line element (2.3) does form a class of solutions to the Einstein equations (2.2). As pointed out in [44], the line element (2.3) defines a metric of a uniform black brane written in the ingoing Eddington-Finkelstein coordinate, moving at velocity β i along spatial direction x i . Discussing fluid dynamics we closely follow [44] and promote the constant parameters β i and b into arbitrarily slowly varying functions of boundary coordinates x α , In general, (2.6) no longer solves the Einstein equations (2.2). The method developed in [44] is to add suitable corrections to (2.6), so that (2.2) is satisfied by the new line element. The corrected metric is not easily found for a general configuration. The authors of [44] introduced a systematic way to construct the corrected metric, expanding in the long wavelength limit. A scale associated with this expansion should be much larger than a characteristic scale of the system, such as inverse of the temperature 1/T . That is the velocity and temperature fields are assumed to vary slowly on this scale, admitting a gradient expansion around some arbitrarily chosen point, such as spacetime origin x α = 0. The Einstein equations (2.2) for the metric are then solved order by order in the boundary derivative expansion. As has been pointed out in the Introduction, the dual metric was constructed up to second order in the velocity gradient, including nonlinear terms quadratic in the velocity gradient.
Our main goal is to perform a summation over all higher order derivative terms in the boundary stress tensor, but as stressed employing linear approximation. Rather than resorting to order by order derivative expansion of [44], we instead linearize the problem in perturbations of the fluid mechanical variables u µ (x α ) and b(x α ). More specifically, the fluid velocity and temperature parameters are expanded as where, as in [44], we multiply β i (x α ) and b 1 (x α ) by a small number , which is an order counting parameter. Below we are going to systematically trace only the terms linear in and set to one in the final expression of the fluid stress tensor. b 0 corresponds to the equilibrium temperature of the dual fluid system while b 1 (x α ) accounts for the dissipative corrections. In what follows, we use conformal symmetry to set b 0 to one.
Substituting (2.7) into (2.3), the "seed" metric, i.e., a linearized version of (2.6) is where the first line is the line element of AdS 5 black hole written in the ingoing Eddington-Finkelstein coordinate. As has been explained above, the metric (2.8) does not solve the Einstein equations (2.2) at order . The term linear in in (2.8) is only a part of the metric we are after, and additional corrections at this order have to be introduced. The full metric is formally written as where g (0) is the first line of (2.8) and g (1) in corresponds to the term linear in of (2.8). The term g (1) corr is the added correction, whose form has to be determined via solving the bulk Einstein equations (2.2).
In order to proceed with the computations, we fix gauge. Following [44], we work in the "background field" gauge, (2.10) We pause to explain the results of the above gauge condition. The most general form of the undetermined metric correction g (1) corr could be parameterized as, The condition g rr = 0 implies g rr = 0. The metric components g rµ can be read off from eqs. (2.8) and (2.11), g rµ = (1, − β i ) + (g rv , g ri ) .
(2.12) Therefore, the second condition g rµ ∝ u µ amounts to requiring that In other words, up to O( ), the vector component g ri will be set to zero. The last condition in (2.10) gives a constraint, Under the gauge condition (2.10), the line element (2.11) for g (1) corr can be rewritten in the similar fashion as that of [44], where the trace part of g ij is explicitly denoted as scalar function h. Therefore, α ij is a symmetric traceless tensor of rank two. Notice that all the metric components {h, k, j i , α ij } are functions of the bulk coordinates {x α , r}. We shall find that these functions are functionals of the fluid velocity β i , which we leave as an undetermined parameter. Their precise forms have to be determined by solving the bulk Einstein equations (2.2), supplemented with appropriate boundary conditions to be discussed in subsection 2.2.
Once the dual metric is found, the fluid stress tensor of the boundary CFT can be computed via holographic dictionary [19,20]. The dual fluid system is defined on the r = ∞ hypersurface. However, in an asymptotically AdS spacetime, holographic renormalization is needed to remove divergences near conformal boundary r = ∞. To proceed, we consider a hypersurface Σ at constant r. The outgoing normal vector n M to Σ is The induced metrics γ M N and γ M N on the hypersurface Σ are constructed as The extrinsic curvature tensor K M N of the hypersurface Σ is Using the formula of [61,62], the stress tensor for the dual fluid is where G µ ν is the Einstein tensor constructed from the induced metric γ µν and K ≡ γ µν K µν . The last two terms in eq. (2.19) are the counter-terms needed to remove divergences near the conformal boundary r = ∞.
Applied to the metric (2.15), the fluid stress tensor can be expressed in terms of the functions {h, k, j i , α ij }. We here record all the components ofT µ ν , where the notation ∂ (i β j) etc stands for symmetrization over the indices i and j. It is important to notice that in the above expression forT µ ν , we have already dropped terms, which explicitly vanish as r → ∞.

Dynamics: the bulk gravity and the boundary hydrodynamics
In this subsection we write down dynamical equations which determine the functions {h, k, j i , α ij } in the bulk. Having these functions at hand, we extract the fluid stress tensor via eq. (2.20). To proceed, we have to specify proper boundary conditions for these metric functions. The first one is a regularity requirement for the metric over the whole range of r, in particular at unperturbed horizon r = 1. This follows from our choice of the ingoing Eddington-Finkelstein coordinate in which the metric is free of any coordinate singularity. The second boundary condition comes from asymptotic considerations. In the present paper we restrict our analysis to boundary fluid dynamics in flat spacetime with the Minkowski metric η µν . Therefore, we require that the metric correction does not modify the asymptotic features of the metric (2.6). This condition tightly constrains the large r behavior for metric functions {h, k, j i , α ij }. Specifically, as r → ∞, their behaviors should be restricted as Yet some integration constants remain unfixed due to a freedom of defining fluid velocity. We follow [44] and choose a frame for the dual fluid system. We will work in the "Landau frame" defined by u µ T µν We are now in the position to study dynamics of the bulk gravity. As pointed out in section 3.2 of [44], there is one redundancy among the total 15 components of the Einstein equations (2.2). Similarly to [44], we classify the remaining 14 equations into constraint equations and dynamical ones. In order to construct the fluid stress tensor, we only solve the dynamical equations. This will lead us to an "off-shell" fluid stress tensor with undetermined fluid velocity but with the transport coefficient functions fixed. The constraint equations will be shown to be equivalent to the conservation law of thus constructed fluid stress tensor.
The first equation we are to consider is E rr = 0, where s 0 and s 1 are arbitrary functions of the boundary coordinates x α . A nonzero function s 0 would violate the asymptotic requirement for h as specified in eq. (2.21). In addition, s 1 = 0 will result in T 00 Diss = 0 as can be seen from eq. (2.20). Therefore, the constraint on the asymptotic behavior at the infinity and the "Landau frame" convention enforce h = 0.
We proceed by considering E rv = 0, Clearly, the function k cannot be found until the vector j i and tensor α ij are computed and we will postpone integration of eq. (2.25) until first solving for j i and α ij . Fortunately, the dynamical equations for j i and α ij can be disentangled from k. The dynamical equation for j i follows from E ri = 0, The diagonal and off-diagonal components of α ij should be treated separately. We first consider the off-diagonal components emerging from E ij = 0 with i = j, The diagonal components of α ij are coupled with the function k. We present the dynamical equation for α 11 , where the term ∂ 2 r k will be eliminated using eqs. (2.25) and (2.26), The equation for α 11 can be put into a new form, (2.30) Similar equations hold for α 22 and α 33 . Remarkably, the equations for the off-and diagonal components of α ij can be combined into a unified form, We are to solve these second order partial differential equations (2.26) and (2.31). As has been outlined above, j i and α ij are linear functionals of β i . They can be uniquely decomposed as where σ ij and π ij are defined in (1.5) and (1.8). The above decomposition is obviously inspired by the special structure of the source terms in (2.26) (2.33) The dynamical equation (2.25) is rendered into the form, In section 4 we will study these asymptotic behaviors by exploring the equations (2.33) near the boundary.
Here we summarize these studies, where the momenta-dependent functions η(ω, q 2 ) and ζ(ω, q 2 ) are formally introduced as coefficients in the asymptotic expansion. The boundary conditions (2.21) and (2.22) have been already applied while deriving eq. (2.35). Asymptotic behaviors as r → ∞ by themselves cannot fix the coefficients η(ω, q 2 ) and ζ(ω, q 2 ). To find them, we have to solve the RG equations (2.33) in full starting from the horizon and integrating up to the boundary. The regularity requirement at r = 1 is found to be sufficient to fix η and ζ uniquely. We postpone completing this computation until sections 3 (analytic) and 4 (numeric). Meanwhile, based on (2.35), as r → ∞ the components j i and α ij behave as (2.36) The large r behavior of the function k follows from first integrating (2.34) over r and then making use of eq. (2.35),  7). This establishes, fully and uniquely, the relation (1.7) between the dissipative part of the stress tensor and the large r behaviors of the functions c and d, as encoded in the viscosity functions η(ω, q 2 ) and ζ(ω, q 2 ).

Generalized Navier-Stokes equations and spectrum of small fluctuations
Generalized, all order Navier-Stokes equations follow from the bulk constraints E vv = E vi = 0. We find it more convenient to study suitable linear combinations of these constraints and dynamical equations. More specifically, the combination r 2 f (r)E vr + E vv = 0 states (2.38) The second constraint r 2 f (r)E ri + E vi = 0 yields (2.39) Taking the large r limit of (2.38) and (2.39) results in Eqs. (2.40) are equivalent to the boundary stress tensor conservation law ∂ µ T µν = 0, which determines the temperature and velocity profiles as functions of time, provided initial configuration is specified. We close this section by studying the spectrum of small fluctuations of the local fluid mechanical variables. We consider a plane wave ansatz for the velocity β i (x α ) and temperature b(x α ), Similarly, we obtain dispersion relation for sound mode by taking q β, Notice that the second viscosity function ζ(ω, q 2 ) does not contribute to the shear mode. The dispersion relations (2.42,2.43), are exact, namely valid for any value of ω and q. As has been already discussed in [52], beyond small ω, q limit, these equations generate infinitely many solutions, which are the quasinormal modes of the dual theory.
In the next section, we will determine the viscosity functions η(ω, q 2 ) and ζ(ω, q 2 ) in the hydrodynamic limit ω, q i 1. The results have been quoted in eq. (1.9). With these expressions at hand, we consider corrections to dispersion relations due to the higher order derivative terms. Solving the dispersion equations (2.42) and (2.43) perturbatively, we obtain (2.44) These results can be compared with those obtained in [44,45]. The shear mode dispersion relation agrees with the quasi-normal mode computations of [45] up to q 4 terms. This further validates the argument of [44,45] that the third order derivative terms in the fluid stress tensor are crucial in correctly producing the shear wave mode at order q 4 . However, the sound wave dispersion relation presented here is different from the one obtained in [44,45]: it is being corrected by the third order derivative in the fluid stress tensor. The analytical expressions for the dispersion relations should agree at small ω, q i with numerical results of [63].

Viscosity functions in the hydrodynamic regime: Perturbative analysis
In this section we solve eqs. (2.33) and (2.34) perturbatively in the hydrodynamic regime ω, q i 1. We then compute the fluid stress tensor using thus constructed perturbative metric correction, up to third order in the derivative expansion. Up to second order, we reproduce some well-known results from the literature, validating correctness of our formalism. We also compute new third order transport coefficients. Finally, we give a formal construction to any order in the derivative expansion, which is fully consistent with the numerical analysis of section 4.
To trace the order in the derivative expansion, we multiply ω and q i by a small parameter λ ω → λω, q i → λq i .  Correspondingly, the metric correction is counted by powers of λ, which follows from (3.2). Notice, however, that different orders in the coefficients a n , · · · , d n mix due to explicit derivatives in the decomposition (2.32). It is straightforward to write down equations for a n , · · · , d n at each order in λ. In what follows, we explicitly solve these equations imposing the boundary conditions discussed in section 2.

Perturbative solution to the metric correction
Metric correction at zeroth order. To the lowest order, we have equation for a 0 only, whose generic solution is a 0 (ω, q, r) = Cr 4 + C . (3.5) In the above solution, the constant C multiplies a non-normalizable mode, which deforms the metric of the boundary field theory. It has to vanish by the condition (2.21). The remaining constant C corresponds to a shift of the fluid velocity and is set to zero by the "Landau frame" convention (2.22). Therefore, there exists no nontrivial solution for a 0 , which is consistent with intuition that metric correction appears starting from first order in the derivative expansion.
Metric correction at first order. Up to the first order in the derivative expansion, it is sufficient to consider the system of differential equations, where a 0 = 0 was already used. We study (3.6) as an example of how to fix integration constants. We first rewrite equation for c 0 as (r 7 − r 3 )∂ 2 r c 0 + (5r 6 − r 2 )∂ r c 0 + 3r 4 = 0 =⇒ r 2 ∂ r (r 5 − r)∂ r c 0 + 3r 4 = 0. (3.7) First integration is done from 1 to r: where the integration constant c 0 is set to zero by the regularity condition at r = 1. Then, as r → ∞, the right-handed side of (3.8) falls off like ∼ 1/r 2 , so that it is valid to integrate the above equation from r to ∞, To fix the integration constant c 0 , we consider the large r behavior of F (r), Therefore, to keep the asymptotic requirement for α ij as given in (2.21), the integration constant c 0 = 0. So, c 0 (ω, q, r) = F (r).

(3.11)
A remark about the integration constant c 0 is worthy. In principle, the outer integral in eq. (3.9) might also be done from 1 to r, but with a new integration constant different from c 0 . This constant would have to be determined by the same asymptotic considerations. The final result for c 0 is still given by eq. (3.11). We proceed with a 1 , =⇒ a 1 (ω, q, r) = −iωr 3 + A(ω, q 2 )r 4 + A (ω, q 2 ). (3.12) The functions A(ω, q 2 ) and A (ω, q 2 ) should equal zero, following the same arguments as below eq. (3.5). We arrive at a first nontrivial expression for the function a, a 1 (ω, q, r) = −iωr 3 .

(3.13)
Up to the first order in the derivative of the fluid velocity, equation (2.34) simplifies to ∂ r k (1) = 2r 2 ∂β, (3.14) and the result for k (1) is where we use the convention (2.22) to set the integration constant in (3.15) to zero. We summarize the metric correction at the first order in the derivative of the fluid velocity, Metric correction at second order. The analysis above can be straightforwardly extended to the second order. The relevant equations are where the expression (3.13) for a 1 has been used to simplify the equation for c 1 . The function b 0 obeys the same equation as a 0 except for the source term, where the large r behavior for the source term is ∼ 1/r 3 . The generic solution for b 0 can be obtained by integration over r, The outer integral is taken from 1 to r to make it well-defined. The integration constants b 0 and b 0 will be again fixed by the asymptotic conditions and "Landau frame" convention as done for a 0 , The solution for b 0 now reads (3.20) The equation for c 1 can be solved similarly to its zeroth order counterpart c 0 . We present final results, dy 2y 3 ∂ y c 0 (y) − y + 3y 2 c 0 (y) , (3.21) which, as r → ∞, falls off as Since the source term in the equation for a 2 decays rather rapidly in the large r regime, solution for a 2 is where the integration constants are fixed in a similar manner as in the above cases of a 0 and a 1 . The second order correction k (2) is solved by, where the integration constant is again fixed by the condition (2.22).
We are led to the large r behavior for the metric correction at second order,

(3.25)
Metric correction at third order. In order to extend previous perturbative analysis to O(∂ 3 ), we consider the following system of differential equations,

(3.26)
Since the equation for d 0 is of the same structure as that of c 0 , we can solve for it in the very same way as we did for c 0 : yc 0 (y) dy, (3.27) where the integration constants are fixed by the boundary conditions (2.21) and (2.22). In the large r limit, d 0 behaves as (3.28) It is straightforward to integrate the remaining equations in (3.26) over r and fix the integration constants in the same way as has been done for the lower order counterparts. For brevity, we only present the final results, (3.29) We also have a third order version of eq. (2.34), The large r behavior for metric correction at order O(∂ 3 ) is,

(3.31)
Metric correction at O(∂ n ) with n ≥ 4. We end with a formal argument towards arbitrarily higher order metric correction in the derivative expansion. For n ≥ 4, we have the following system of recursive differential equations, r c n−1 + (5r 6 − r 2 )∂ r c n−1 − 2iωr 5 ∂ r c n−2 − r∂ r a n−1 + a n−1 − 3iωr 4 c n−2 , 0 =r∂ 2 r a n − 3∂ r a n − q 2 r 3 ∂ r c n−2 , where d −1 should be understood as null. It can be shown by induction that large r asymptotic behaviors of the coefficients a n , b n , c n , d n (n ≥ 4) are of universal form, a n (ω, q, r) → S n a (ω, q) r (3.33) The n-th order counterpart of (2.34) is The functions S n a etc. are to be determined by solving the recursive equations (3.32), similarly as we did for the lower order metric corrections. Generically, they will take a form of fixed order polynomials in ω and q, S n = m=n m=0 ρ m ω m q n−m . Although we are not able to give exact analytical expressions for S n a etc., the formal analysis presented here is useful in obtaining the general structure of T µν up to arbitrary order in the derivative expansion. At any order O(∂ n ) with n ≥ 4, the components j i and α ij fall off as in the large r regime.

Fluid stress tensor up to third order and beyond
With the perturbative solutions to the metric correction at hand, we proceed by computing the fluid stress tensor (2.20). Up to the second order O(∂ 2 ), which is exactly the results obtained in [44] when linearized as in (2.7). Let us write the fluid stress tensor as a formal derivative expansion, where the zeroth order T 6∂ 2 + 6π − π 2 + 12 2 − 3 ln 2 + ln 2 2 ∂ 2 v σ ij where the tensor structure π ij appears for the first time. From the expression for the stress tensor we can immediately read off the viscosity functions η and ζ as quoted in (1.9). Substituting eq. (3.35) into eq. (2.20), we arrive at a general form of T µν at arbitrary order in the gradient expansion,

All order linearized hydrodynamics
To fully account for all order derivative terms in the fluid stress tensor, we resort to numerical techniques for solution of the RG equations (2.33), extending validity of the above discussed hydrodynamic regime to large momenta. It is convenient to rescale the functions a(r) and b(r) a(ω, q, r) = r 4ã (ω, q, r), b(ω, q, r) = r 4b (ω, q, r), (4.1) and also use u-coordinate instead of r In u coordinate, the horizon is located at u = 1 while the conformal boundary is at u = 0.
In what follows, we also use notationsc(u) = c(r) andd(u) = d(r) to stress that they are functions of u. Equations (2.33) become where prime denotes derivative with respect to u. The problem of resumming all order derivative terms in the boundary stress tensor is reduced to a boundary value problem of the system of ordinary differential equations (4.3).
In the rest of this chapter we will solve this problem by two methods. The first method will be fully numerical while the second one is an approximate analytic scheme. Both methods are demonstrated to converge to the same results.

Numerical results from the shooting method
We have to impose boundary conditions both at the horizon and asymptotic infinity. We apply a shooting method to solve the system (4.3). The main idea behind the shooting method is to reduce the boundary value problem to an initial value problem for a system such as (4.3). One starts from a trial solution (initial condition) at one boundary (horizon) and integrates the system until the other boundary. Then, thus obtained solution should be matched with boundary conditions at the end of the integration. That would not happen for an arbitrary trial initial condition: the trial solution has to be fine tuned in order for the boundary conditions at the end of the integration to be satisfied. This fine tuning problem can be turned into an optimization procedure.
We are now to discuss an implementation of this method for the system (4.3) given the boundary conditions presented in section 2. In order to fully find a solution for four second order differential equations, we have to specify overall eight boundary conditions. The regularity requirement at horizon, boundary conditions (2.21) and the Landau frame convention (2.22), indeed do provide precisely eight conditions: two at the horizon and six at the conformal boundary.
Series solution near the horizon. We start from the regularity requirement at the unperturbed horizon u = 1. To have a regular black hole solution near u = 1, the functions a,b,c,d have to be Taylor expandable, Series solution near the conformal boundary. We turn to discuss near u = 0 behavior for these functions. At u = 0, the characteristic indices for the system (4.3) are 0 and 4. Series solution then takes the form, where the subscript "b" marks the asymptotic infinity u = 0. The logarithmic branch, whose coefficients are labeled with the subscript "L", is necessary due to the fact that the difference between two characteristic indices is integer. Similarly to the near horizon expansion, by substituting (4.5) into (4.3), all the coefficients of (4.5) are related to the following eight The boundary condition (2.21) implies that while the "Landau frame" convention (2.22) constrains two more expansion coefficients leaving only two undetermined coefficients, C 4 b and D 4 b , which have to be determined through dynamical evolution from the horizon.
With the conditions (4.6) and (4.7) at hand, the logarithmic branch in (4.5) vanishes identically and the large r behavior for these functions is thus (4.8) In terms of functions a, b, c, d, we have (4.9) The coefficient functions C 4 b and D 4 b can be now identified with the viscosity functions η(ω, q 2 ) and ζ(ω, q 2 ) Our problem is now mapped into finding such six near horizon expansion coefficients b to vanish in accordance with the above boundary conditions. Once this is achieved, the coefficients C 4 b and D 4 b should be read off from the final solution. Therefore, the boundary value problem for the system (4.3) is reduced to the problem of root-finding or optimization in numerical analysis, The procedure is repeated for each value of ω and q.
Further numerical details and results. The near-horizon series solution (4.4) makes it possible to evaluate ã,ã ,b,b ,c,c ,d,d (ω, q, u) (4.12) at some point u = 1 − u , close to the horizon (u 1). That helps to avoid a numerically problematic region near horizon where the system of equations has a singular point. In our shooting routine, we integrate the system (4.3) from this near-horizon point u = 1 − u till some point u = u , close to asymptotic infinity u = 0.
We use Newton's method [64] for root-finding and find it works rather well. Efficient initial guesses for the trial values of A 0 h , A 1 h , B 0 h , B 1 h , C 0 h and D 0 h are provided by linearly extrapolating previously computed roots along the ω or q 2 -axis. The numerical procedure is started from the point ω = q 2 = 0, known exactly from section 3 where C is the Catalan constant with approximate value C ≈ 0.915966. We set u to 10 −5 and checked stability of the results with respect to variations of u . Our numerical results for the viscosity functions η(ω, q 2 ) and ζ(ω, q 2 ) are shown in Figure 1. The real parts of η and ζ decrease with momenta until reaching minima around points ω ≈ 3.0, q 2 = 0 and ω ≈ 1.9, q 2 = 0 , respectively. A sign of damped oscillations is observed in the results, while eventually, the real parts vanish at very large ω and/or q 2 . The imaginary parts of the viscosities first increase from zero up to some maxima around ω ≈ 1.7, q 2 = 0 for η and ω ≈ 1.0, q 2 = 0 for ζ. With further increase of the momenta, the imaginary parts decrease reaching zero at large momenta.
Vanishing of transport coefficients at very large momenta is well anticipated: there should be no response at very short times or distances. This point is critical for the generalized relativistic hydrodynamics to be causal. To further confirm our observation, we focus on large momenta behavior for the viscosity functions. In Figures 2 and 3, we show our results for very large ω or q 2 . The imaginary parts of η and ζ are identically zero at ω = 0.  This is obvious from the eqs. (4.3), which have no imaginary terms at ω = 0. Vanishing of the viscosity functions at large ω (and/or q 2 ) is an important factor for reliably addressing early time stage in heavy ion collisions [51,55,65]. The imaginary part of η never turns negative. This is indeed a necessary condition for stability of the fluid equations, as is seen from the shear mode dispersion relation (2.42). In contrast, the imaginary part of ζ does change sign at intermediate values of ω and q 2 as seen in Figure 1. However, negative ζ does not mean an instability: it only contributes as a practically negligible correction to the dispersion relation of sound mode (2.43), but never turns any of these modes into unstable. Indeed, the absolute value of ζ(ω, q 2 ) is always highly suppressed as compared to that of η(ω, q 2 ). Furthermore, the inequality q 2 ζ(ω, q 2 ) η(ω, q 2 ) is valid in all the kinematic range covered by our numerical results. Therefore, for any practical applications and hydrodynamic modelings, it is probably always safe to ignore the viscosity function ζ(ω, q 2 ).

Approximate results from the matching method
In this subsection we provide an alternative approach to solving equations (4.3) based on an analytic approximate scheme. This provides us with a possibility to check the numerical results of the previous subsection. The main idea is to adopt a matching method, which was introduced in [66] in order to provide analytical evidence for condensation phenomena in a holographic superconductor model [67]. The method is based on the series expansions (4.4) and (4.5), which not only exactly solve the system (4.3) but also should match over the whole regime of u ∈ [0, 1].
The approximation we are to employ is a truncation of the series (4.4) and (4.5): we will keep up to eleven terms in each expansion, i.e., order u 10 and (u − 1) 10 , respectively. While the truncated series would not any more solve the system (4.3) exactly, keeping enough terms guarantees accurate solutions near the horizon and conformal boundary. We then match the truncated series solutions at an intermediate point such as u = 1/2. Taking a(ω, q, u) as an example, the matching of its value and first order derivative at u = 1/2 results in, where the boundary conditions (4.6) and (4.7) should be imposed. This method casts the system (4.3) into algebraic equations for these expansion coefficients. On the one hand, having kept a large number (eleven) of terms in the expansions, we achieve stable and highly accurate results. The viscosity functions obtained from the matching scheme are displayed in Figure 4. On the other hand, the large number of terms kept in (4.14) lead to analytical but rather lengthy and not very illuminating expressions for the viscosity functions. They are of the type where P s and Qs are finite order polynomials in ω and q 2 . Below, we only quote a few terms in hydrodynamic expansion of these polynomials, Q ζ (ω, q 2 ) =1.00 − 3.13iω − 4.93ω 2 + 0.526q 2 − 1.44iωq 2 + 5.20iω 3 − 2.01ω 2 q 2 + 4.11ω 4 + 0.107q 4 + 1.89iω 3 q 2 − 2.59iω 5 + · · · . (4.16) The structure (4.15) implies that the viscosities have a number of poles (zeros of Qs) and this is quite consistent with the arguments made in [52] that exact η(ω, q 2 ) in fact has infinitely many poles.

Summary and discussion
In this paper we have provided all the details and expanded presentation of the results advertised in our short publication [1]. As a further development of the ideas put forward in [51], we have consistently determined the energy-momentum stress tensor of a weakly perturbed conformal fluid, whose underlying microscopic description is a strongly coupled N = 4 super-Yang-Mills theory at finite temperature. The results were derived by linearizing the fluid/gravity correspondence. We have included all order derivative terms in the boundary stress tensor, limiting the study to small amplitude perturbations only. We have found that all order dissipative terms in the fluid stress tensor are fully accounted for by two (generalized) momenta-dependent viscosity functions η(ω, q 2 ) and ζ(ω, q 2 ). η(ω, q 2 ) is a transport coefficient in front of the shear tensor σ ij while ζ(ω, q 2 ) is a coefficient of a new tensor π ij , which is given in terms of third order derivative of the fluid velocity. As one of our main results, we have derived a closed-form linear holographic RG flowtype equations (2.33) for the viscosity functions. Intriguingly, an analogous holographic RG flow equation for electrical conductivity obtained in [28] is a nonlinear one. The constraint components of the bulk Einstein equations (2.2) have been shown to generalize the Navier-Stokes equations, consistently with the conservation laws of the fluid stress tensor. We have analytically computed the viscosity functions, up to third order in the hydrodynamic gradient expansion, and the dispersion relations for the shear and sound waves. These third order corrections are needed in order to correctly reproduce the dispersion relation for the shear wave up to order O(q 4 ), as emphasized in [44,45].
To include all order dissipative effects in the fluid stress tensor, we have solved numerically the RG flow-type equations (2.33). The numerical results for the viscosity functions are displayed in Figure 1. Based on our numerical calculations, we have been able to significantly extend knowledge about the viscosity functions in the hydrodynamic limit, providing in (4.17) an expansion up to fifth order. Consistently with physical intuition, the viscosities vanish at very large momenta as seen in Figures 2 and 3. We have verified our results by solving the equations (2.33) using an alternative method.
Importantly, the hydrodynamic theory constructed in this work is causal and should be free of any instabilities if implemented as a dynamical model for plasma evolution. The viscosity function encodes an infinite set of quasi-normal modes [63] including corresponding residues.
Obviously, this is not QCD, but we hope that some generic features about momentadependent transport coefficients and high gradient structures have been captured by our results. They can help in building new models of causal relativistic hydrodynamics, beyond Navier-Stokes or Israel-Stewart.
The stress tensor computed in this work resums all order gradients linearized in the amplitude of fluid velocity and is not readily applicable to phenomena where nonlinearities are important, such as Bjorken flow. Nevertheless, the results reported here might be useful both in estimating phenomenological roles and sizes of higher gradient effects as suggested in [51,52] and in direct studies of experimentally observed phenomena where linear dissipative terms play an important role. Sonic booms created by jets or heavy quarks, fluctuations in the flows and correlations are examples of applications which we have in mind.
A question of convergence of the gradient expansion has been raised in [55], and it has been argued there that radius of convergence is in fact zero presumably due to a factorial growth in the number of terms at high orders. We have not observed any convergence issues, thus indirectly confirming the conclusions of [55] about a nonlinear origin of the convergence problem.
Throughout the paper, we have been referring to equations (2.33) as holographic RG flow-type equations for the viscosity functions. Indeed, the radial coordinate r is frequently associated with a scale of RG flow of the boundary CFT. While we do have an evolution in r, we identify physical quantities (viscosities) only at the infinite boundary. In the spirit of holographic Wilsonian RG [68][69][70][71][72][73], it would be proper to introduce a finite cutoff along radial direction and define associated physical quantities, such as η(r, ω, q 2 ) and ζ(r, ω, q 2 ), at the cutoff surface. That would result in an RG evolution of these momenta-dependent coefficients, thus extending previous results on RG flows of the shear viscosity coefficient η 0 . So far, the university of the ratio (1.6) was clarified as a consequence of no RG flow from the horizon to the boundary of the Navier-Stokes hydrodynamics [28].
As further development of this project, we plan to extend our present study to conformal fluids in a weakly curved background manifold, with all order derivative terms resummed. Metric perturbations at the boundary may be taken into account following [47,48]. Additional transport coefficient functions associated with the boundary curvature are expected to emerge [52].