Second-order transport, quasinormal modes and zero-viscosity limit in the Gauss-Bonnet holographic fluid

Gauss-Bonnet holographic fluid is a useful theoretical laboratory to study the effects of curvature-squared terms in the dual gravity action on transport coefficients, quasinormal spectra and the analytic structure of thermal correlators at strong coupling. To understand the behavior and possible pathologies of the Gauss-Bonnet fluid in $3+1$ dimensions, we compute (analytically and non-perturbatively in the Gauss-Bonnet coupling) its second-order transport coefficients, the retarded two- and three-point correlation functions of the energy-momentum tensor in the hydrodynamic regime as well as the relevant quasinormal spectrum. The Haack-Yarom universal relation among the second-order transport coefficients is violated at second order in the Gauss-Bonnet coupling. In the zero-viscosity limit, the holographic fluid still produces entropy, while the momentum diffusion and the sound attenuation are suppressed at all orders in the hydrodynamic expansion. By adding higher-derivative electromagnetic field terms to the action, we also compute corrections to charge diffusion and identify the non-perturbative parameter regime in which the charge diffusion constant vanishes.

3 Gauss-Bonnet transport coefficients from fluid-gravity correspondence

Introduction
Gauge-string duality has been applied successfully to explore qualitative, quantitative and conceptual issues in fluid dynamics [1][2][3][4][5][6][7][8][9][10]. Although the number of quantum field theories with known dual string (gravity) descriptions is very limited, their transport and spectral function properties at strong coupling can in principle be fully determined, thus giving valuable insights into the behavior of strongly interacting quantum many-body systems. Moreover, dual gravity methods can be used to determine coupling constant dependence of a variety of physical quantities with an ultimate goal of interpolating between weak and strong coupling results and describing, at least qualitatively, the intermediate coupling behavior in theories of phenomenological interest [11][12][13][14][15][16].
For generic neutral fluids, there are two independent first-order transport coefficients (shear viscosity η and bulk viscosity ζ), and fifteen second-order coefficients 1 (see e.g. [21]). For Weyl-invariant or "conformal" fluids, the additional symmetry constraints reduce the number of transport coefficients to one at first order (shear viscosity η) and five at second order 2 (usually denoted τ Π , κ, λ 1 , λ 2 , λ 3 ). The coefficients η, τ Π , λ 1 , λ 2 are "dynamical", whereas κ and λ 3 are "thermodynamical" in the classification 3 introduced in ref. [19]. In the parameter regime where the dual Einstein gravity description of conformal fluids is applicable (e.g. at infinite 't Hooft coupling λ = g 2 Y M N c and infinite N c in theories such as N = 4 SU (N c ) supersymmetric Yang-Mills (SYM) theory in d = 3 + 1 dimensions), the six transport coefficients (in d space-time dimensions, d > 2) are given by [22] η = s/4π , (1.1) , (1.5) where s is the entropy density, ψ(z) is the logarithmic derivative of the gamma function, 1 The existence of a local entropy current with non-negative divergence implies η ≥ 0, ζ ≥ 0 [17] and constrains the number of independent coefficients at second order to ten [18]. Alternatively, independent "thermodynamical" [19] terms in the hydrodynamic expansion can be derived from the generating functional without resorting to the entropy current analysis [4,5]. A computerized algorithm determining all tensor structures appearing at a given order of the hydrodynamic derivative expansion has been recently proposed in ref. [20]. Modulo constraints potentially arising from the entropy current analysis (not attempted in ref. [20]), it identifies 68 new coefficients for non-conformal neutral fluids and 20 coefficients for conformal ones at third order of the derivative expansion. 2 There are no further constraints in addition to η ≥ 0 coming from the non-negativity of the divergence of the entropy current in the conformal case [18]. 3 Essentially, the coefficient is called "dynamical" if the corresponding term in the derivative expansion vanishes in equilibrium, and "thermodynamical" otherwise.
In the regime of strong coupling, theories with gravity dual description appear to exhibit robust properties of transport coefficients and relations among them. One of such properties is the universality of shear viscosity to entropy density ratio η/s = 1/4π in the limit described by a dual gravity with two-derivative action [41], [42], [43][44][45]. Another one seems to be the Haack-Yarom relation: following the observation in ref. [46], the linear combination of the second-order transport coefficients 4 was proven to vanish in all conformal theories dual to two-derivative gravity 5 [47]. Eqs.
(1.2), (1.4), (1.5) show this explicitly. Somewhat surprisingly, the Haack-Yarom relation continues to hold to next to leading order in the strong coupling expansion, at least in N = 4 SU (N c ) supersymmetric Yang-Mills theory in d = 3 + 1 dimensions in the limit of infinite N c [30], in theories dual to curvature-squared gravity [30], in particular, in the Gauss-Bonnet holographic liquid 6 (perturbatively in the Gauss-Bonnet coupling) [49]. It was shown recently that the result H = 0 continues to hold for non-conformal liquids along the dual gravity RG flow [50]. 7 It remains to be seen whether such robustness extends to higher-order transport coefficients and/or other properties of strongly coupled finite temperature theories and whether it is related to the presence of event horizons in dual gravity. 8 Monotonicity and other properties of transport coefficients are of interest for studies of near-equilibrium behavior at strong coupling, in particular, thermalization, and for attempts to uncover a universality similar to the one exhibited by the ratio of shear viscosity to entropy density. Monotonicity of transport coefficients or their dimensionless combinations may seem more exotic than the monotonicity of central charges [51,52] or the free energy [53,54], yet it is often an observed property, at least in a given state of aggregation [41,55].
In N = 4 SYM at infinite 9 N c , the shear viscosity to entropy density ratio appears to be a monotonic function of the coupling [41], with the correction to the universal infinite coupling result being positive [23,25], η s = 1 4π 1 + 15ζ(3)λ −3/2 + . . . . (1.8) Subsequent calculations revealed that the corrections coming from higher derivative terms in the gravitational action can have either sign [56,57]. For the action with generic curvature squared higher derivative terms where the cosmological constant Λ = −6/L 2 , the shear viscosity -entropy density ratio is 10 (1.10) The sign of the coefficient α 3 affects not only viscosity but also the analytic structure of correlators in the dual thermal field theory [12].
Corrections to Einstein gravity results computed from generic higher-order derivative terms in the dual gravitational action can be trusted so long as they remain (infinitesimally) small relative to the leading order result, as they are obtained by treating the higherderivative terms in the equations of motion perturbatively. This limitation arises due to Ostrogradsky instability and other related pathologies such as ghosts associated with higherderivative actions [58], [59], [60], [61] (see also refs. [62], [63] for a modern discussion of Ostrogradsky's theorem, and ref. [64] for an interesting historical account of Ostrogradsky's life and work). One may be tempted to lift the constraints imposed by Ostrogradsky's theorem by considering actions in which coefficients in front of higher derivative terms conspire to give equations of motion no higher than second-order in derivatives as happens e.g. in Gauss-Bonnet gravity in dimension D > 4 or, more generally, Lovelock gravity [65]. Gauss-Bonnet (and Lovelock) gravity has been used as a laboratory for non-perturbative studies of higher derivative curvature effects on transport coefficients of conformal fluids with holographic duals [56,57,[66][67][68][69][70][71][72], [12]. In particular, the celebrated result for the shear viscosity-entropy ratio in a (hypothetical) conformal fluid dual to D = 5 Gauss-Bonnet gravity [57], has been obtained non-perturbatively in the Gauss-Bonnet coupling λ GB . The result would imply that there exist CFTs whose viscosity can be tuned all the way to zero in the regime described by a dual classical (albeit non-Einsteinian) gravity. It was found, however, that 9 At large but finite Nc, and large λ, the result for η/s is also corrected by the term proportional to 39,40]. 10 All second-order transport coefficients for theories dual to the background (1.9) have been computed in ref. [30].
for λ GB outside of the interval the dual field theory exhibits pathologies associated with superluminal propagation of modes at high momenta or negativity of the energy flux in a dual CFT [57, 66-68, 73, 74]. For Gauss-Bonnet gravity in D dimensions (D ≥ 5), the result (1.11) generalizes to [68,75] and the inequalities corresponding to eq. (1.12) become 11 [68,75] − (1.14) Given the constraints (1.14) and monotonicity of η/s in (1.13), one may conjecture a GB gravity bound on η/s [68,75], instead of the Einstein's gravity bound η/s ≥ 1/4π. For 3 + 1-dimensional CFTs, the GB bound would imply η/s ≥ (0.640)/4π [66]. Recently, the constraints (1.12) were confirmed and generalized to Gauss-Bonnet black holes with spherical (rather than planar) horizons by considering boundary causality and bulk hyperbolicity violations in Einstein-Gauss-Bonnet gravity [76]. Since these causality problems arise in the ultraviolet, one may hope that treating Gauss-Bonnet gravity as a low energy theory with unspecified ultraviolet completion would allow one to consider its hydrodynamic (infrared) limit without worrying about causality violating ultraviolet modes, i.e. that it is in principle possible to cure the problems in the ultraviolet without affecting the hydrodynamic (infrared) regime (one may also try to construct a theory with a low temperature phase transition breaking the link between the hydrodynamic IR and causality breaking UV modes [77]). However, a reflection on the recent analysis by Camanho et al. [78] of the bulk causality violation in higher derivative gravity seems to imply that, provided the relevant conclusions of ref. [78] are correct 12 , a reliable treatment of Gauss-Bonnet terms beyond perturbation theory for the purposes of fluid dynamics is not possible. The Einstein-Gauss-Bonnet action in D = 5 is given by where the cosmological constant Λ = −6/L 2 , and l GB is the scale of the Gauss-Bonnet term which a priori is not necessarily related to the cosmological constant scale set by L. As 11 Curiously, in the D → ∞ limit, the range (1.14) is −3/4 ≤ λGB ≤ 1/4. Note that the black brane metric is well defined for λGB ∈ (−∞, 1/4] for any D. We shall only consider D = 5 in the rest of the paper. 12 See refs. [76,[79][80][81][82][83][84] for recent discussions.
argued in ref. [78], the generic bulk causality violations in Gauss-Bonnet classical gravity can only be cured by including an infinite set of higher spin fields with masses squared m 2 s ∝ 1/λ GB l 2 GB . Integrating out these fields to obtain a low energy effective theory would lead to an infinite series of additional higher derivative terms in the gravitation action. Schematically, the modified action would have the form Considering a specific solution (e.g. a black brane whose scale is set by the cosmological constant) and rescaling the coordinates x →x = x/L leads to To suppress contributions (e.g. to transport coefficients) coming from the (unknown) terms with k > 1, one has to assumeλ GB 1. This is similar to the condition l s /L 1 in the usual top-down holography. Thus, generically one may expect results such as (1.11) to be potentially corrected by terms O(λ 2 GB ) and/or higher, and therefore be reliable only for λ GB 1. It seems, therefore, that one essentially cannot escape the Ostrogradsky problem (at least not in classical gravity) by engineering a specific higher-derivative Lagrangian with second-order equations of motion. An alternative view of the aspects of the analysis in ref. [78] has been advocated in refs. [79,81] (see also [80,82] and [76]). Our approach to these problems will be purely pragmatic: 14 we shall a priori ignore any existing or debated constraints on the Gauss-Bonnet coupling and explore the influence of curvature-squared terms on quasinormal spectra and transport coefficients for all range of the coupling allowing a black brane solution, i.e. for λ GB ∈ (−∞, 1/4] (see section 2). In particular, we are interested in revealing any generic features the presence of higher-curvature terms in the action may have (as pointed out in ref. [12], the spectra of R 2 and R 4 backgrounds exhibit qualitatively similar novel features not present in Einstein's gravity). We use the action (1.16) (with l GB = L) to compute transport coefficients, quasinormal spectrum and thermal correlators analytically and non-perturbatively in Gauss-Bonnet coupling, fully exploiting the advantage of having to deal with second-order equations of motion in the bulk. Different techniques will be used to compute Gauss-Bonnet transport: fluid-gravity duality, Kubo formulae applied to two-and three-point correlators, and quasinormal modes. We find that only the three-point functions method allows to determine all the coefficients analytically: other approaches face technical difficulties we were not able to resolve. In a hypothetical dual CFT, constraints on Gauss-Bonnet coupling considered e.g. in ref. [68] manifest themselves in the superluminal propagation of highmomentum modes for λ GB outside of the interval (1.14). In the far more stringent scenario of ref. [78], one may expect to detect anomalous behavior in the regime of small frequencies and momenta and in transport coefficients. Accordingly, we shall look for pathologies in the hydrodynamic behavior of the model at finite values of λ GB indicating the lack of ultraviolet completion and the potential need for corrections coming from the unknown terms in (1.18).
The full non-perturbative set of first-and second-order Gauss-Bonnet transport coefficients can be determined analytically and is given by 15 (1.21) where we have defined An alternative way of writing the Gauss-Bonnet second-order coefficients is given by Eqs.
(4.19) -(4.23). In the limit of λ GB → 0 (γ GB → 1), which corresponds to Einstein's gravity, one recovers the standard results for infinitely strongly coupled conformal fluids in 3 + 1 dimensions given by eqs. (A.1) and (A.2) in appendix A. The result for η was obtained in ref. [57] and the relaxation time τ Π was first found numerically in ref. [67]. Coefficients τ Π and κ were previously computed analytically in ref. [85], and we have reported λ 1 , λ 2 , λ 3 in ref. [86]. To linear order in λ GB , the results coincide 16 with those found in ref. [49].  Computing the energy-momentum tensor correlation functions in holographic models with higher-derivative dual gravity terms, one finds a new pole on the imaginary frequency 15 The shear viscosity η as a function of temperature and γGB is given in Eq. (2.46). 16 The notations used in ref. [49] are related to the ones in this paper by λ0 = ητΠ, δ = 4λGB and axis. This pole, first found in the quasinormal spectrum analysis of ref. [12], is moving from the complex infinity closer and closer to the origin as the parameter in front of the higherderivative term in the action (such as λ GB in eq. (1.16)) increases, and can be approximated analytically in the small-frequency expansion. The poles of this type appear to be generic in higher-derivative gravity: they are present in R 2 and R 4 gravity, and their behavior is qualitatively similar [12]. Another interesting feature of Gauss-Bonnet holographic liquid is the zero-viscosity limit. In ref. [87], Bhattacharya et al. suggested the existence of a non-trivial second-order non-dissipative hydrodynamics, i.e. a theory whose fluid dynamics derivative expansion has no contribution to entropy production while still having some of the transport coefficients non-vanishing. 17 For conformal fluids, the classification of [87] implies the existence of a four-parameter family of non-trivial non-dissipative fluids with η = 0 and non-vanishing coefficients τ Π , κ, λ 1 = κ/2, λ 2 and λ 3 . Given the result (1.11), the hypothetical theory dual to Gauss-Bonnet gravity in the limit of λ GB → 1/4 is a natural candidate for a dissipationless fluid (ignoring for a moment any potential corrections coming from (1.18)). In the limit of λ GB → 1/4 (γ GB → 0) we find [86] ητ Π = 0, (1.27) At first glance, this result realizes the dissipationless liquid scenario outlined in ref. [87]: the shear and bulk viscosities are zero while some of the second-order coefficients are not. However, the relationship κ = 2λ 1 , which is required for ensuring zero entropy production, does not hold among the coefficients in (1.27). We therefore conclude that the holographic Gauss-Bonnet liquid does not fall into the class of non-dissipative liquids discussed in ref. [87]. This may be a hint that the corrections from (1.18) must indeed be included. The paper is organised as follows. In section 2 we analyze the finite-temperature twopoint correlation functions of energy-momentum tensor in the theory dual to Gauss-Bonnet gravity as well as the relevant quasinormal modes in the scalar, shear and sound channels of metric perturbations, including the new pole on the imaginary axis at finite coupling λ GB . Kubo formulas determine the coefficients η, κ and τ Π . The shear channel quasinormal frequency is used to confirm the results for η and τ Π , and to find the third-order transport coefficient θ 1 . We discuss the limit λ GB → 1/4, where the full quasinormal spectrum can be found analytically, and the limit λ GB → −∞. In section 3, we apply the fluid-gravity duality technique to compute the Gauss-Bonnet transport coefficients. All coefficients except κ can be determined in this approach. However, due to technical difficulties, all of them with the exception of η can be found only perturbatively as series in λ GB . A more efficient method of three-point functions is considered in section 4, where all the coefficients are computed analytically and non-perturbatively, and we also discuss the monotonicity properties of the 17 The authors of [87] considered an effective field theory approach [88,89] to non-dissipative uncharged second-order hydrodynamics. The approach relies on a classical effective action and standard variational techniques to derive the energy-momentum tensor. It is thus unable to incorporate dissipation. The inclusion of dissipation into the description of hydrodynamics, using the same effective description, was analysed in [90,91]. coefficients and the zero-viscosity limit. Finally, in section 5 we discuss the influence of higher derivative terms on charge diffusion in the most general four derivative Einstein-Maxwell theory. Section 6 with conclusions is followed by several appendices: in Appendix A, a brief summary of second-order transport coefficients in N = 4 SYM at weak and strong coupling is given. A comparison of notations and conventions used in the literature on second-order hydrodynamics and, specifically, in the discussion of Haack-Yarom relation is given in Appendix B. In Appendix C we outline the procedure of setting the boundary conditions at the horizon in hydrodynamic approximation. Appendices D and E contain some technical results.
2 Energy-momentum tensor correlators and quasinormal modes of Gauss-Bonnet holographic fluid The coefficients of the four-derivative terms in the Gauss-Bonnet action (1.16) ensure that the corresponding equations of motion contain only second derivatives of the metric. The equations are given by The equations (2.1) admit a black brane solution 18 (2. 3) The arbitrary constant N GB will be set to normalize the speed of light at the boundary (i.e. in the dual CFT) to unity, and we henceforth use this value. The solution with r + = 0 corresponds to the AdS vacuum metric in Poincaré coordinates with the AdS curvature scale The parameter γ GB is defined in Eq. (1.25). We shall use λ GB and γ GB interchangeably, and set L = 1 in the rest of the paper unless stated otherwise. The Hawking temperature, the 18 Exact solutions and thermodynamics of black branes and black holes in Gauss-Bonnet gravity were considered in [92] (see also refs. [93][94][95][96][97]). entropy density and the energy density associated with the black brane background (2.2) are given, correspondingly, by The metric (2.2) is well defined for λ GB ∈ (−∞, 1/4] (or γ GB ∈ [0, ∞), with the interval of positive λ GB corresponding to the interval γ GB ∈ [0, 1)). We note that s/T 3 is a monotonically decreasing function of γ GB in the interval γ GB ∈ [0, ∞). The holographic dictionary relating the coupling λ GB of Gauss-Bonnet gravity in D dimensions to the parameters of the dual CFT has been thoroughly discussed in ref. [68] (see also the comprehensive discussion of the D = 5 case in ref. [98]). For a class of fourdimensional CFTs (usually characterized by the central charges c and a), there exists a parameter regime (e.g. λ N 2/3 c 1 [56,68]) in which the dual description is given by Einstein gravity with a negative cosmological constant plus curvature squared terms treated as small perturbations, so that e.g. the coefficient α 3 in the action (1.9) is α 3 ∼ (c − a)/c ∼ 1/N c 1, as in the discussion of the superconformal N = 2 Sp(N c ) gauge theory with four fundamental and one antisymmetric traceless hypermultiplets by Kats and Petrov 19 [56]. For finite λ GB , if a dual CFT exists at all, one may relate the Gauss-Bonnet coupling to the parameters characterizing two-and three-point functions of the energy-momentum tensor in the CFT [68]. In particular, the holographic calculation [68] gives the central charge c Note that the central charge is a monotonically increasing non-negative function of γ GB in the interval γ GB ∈ [0, ∞), with c = 0 at γ GB = 0 (i.e. at λ GB = 1/4). Generically, we may expect λ GB to be a function of both λ and N c at large but finite values of these parameters. We compute the retarded two-point functions G R µν,ρσ of the energy-momentum tensor in a hypothetical finite-temperature 4d CFT dual to the Gauss-Bonnet background (2.2) following the standard holographic recipe [99][100][101][102]. Gravitational quasinormal modes of the background corresponding to the poles of the correlators G R µν,ρσ [99,102] have been computed and analyzed in detail as a function of the Gauss-Bonnet parameter λ GB in ref. [12]. The quasinormal spectrum at λ GB = 1/4 is computed analytically in section 2.4 of the present paper.
The full gravitational action needed to compute the correlators contains the Gibbons-Hawking term and the counter-term required by the holographic renormalisation, (2.10) 19 Other examples, as well as the string theory origins of the curvature-squared terms in the effective action are discussed in ref. [98].
where S GB is the Gauss-Bonnet action (1.16), the modified Gibbons-Hawking term is given by 11) and the counter-term action is (see e.g. [103]) where (2.13) Here γ µν = g µν − n µ n ν is the induced metric on the boundary, n µ is the vector normal to the boundary, i.e. n µ = δ µr / √ g rr , R γ is the induced Ricci scalar and G µν γ is the induced Einstein tensor on the boundary. The extrinsic curvature tensor is K is its trace and the tensor J µν is defined as Similarly, J denotes the trace of J µν . Due to rotational invariance, we may choose the fluctuations h µν of the background metric to have the momentum along the z axis, i.e. we can set h µν = h µν (r)e −itω+iqz , which enables us to introduce the three independent gauge-invariant combinations of the metric components [102]-scalar (Z 1 ), shear (Z 2 ) and sound (Z 3 ): Throughout the calculation, we use the radial gauge h rµ = 0 and the standard dimensionless expressions for the frequency and the spatial momentum .
By symmetry, the equations of motion obeyed by the three functions Z 1 , Z 2 , Z 3 decouple [102]. Introducing the new variable u = r 2 0 /r 2 , the equation of motion in each of the three channels can be written in the form of a linear second-order differential equation where i = 1, 2, 3 and the coefficients A i and B i are given in Appendix D. For some applications, especially in fluid-gravity duality, it will be convenient to use yet another radial variable, v, defined by [57] so that the horizon is at v = 0 and the boundary at v = 1 − γ GB . The new coordinate is singular at zero Gauss-Bonnet coupling, λ GB = 0 (γ GB = 1), thus the results for λ GB = 0, which are identical to those of N = 4 SYM theory at infinite 't Hooft coupling and infinite N c , have to be obtained independently. On shell, the action (2.10) reduces to the surface terms, where the contribution from the horizon should be discarded [99], [104]. In terms of the gauge-invariant variables (2.16), (2.17) and (2.18), the part of the action involving derivatives of the fields can be written as where Z is the derivative of Z(u, ω, q) with respect to the radial coordinate. The functions A i include the boundary contributions from the parts S GB and S GH of the action (2.10), but not from S c.t. . The ellipsis in Eq. (2.23) stands for the boundary terms proportional to the products h µν ( , −ω, −q)h ρσ ( , ω, q) arising from all the three parts of the action (2.10).
In the following, we shall only need those terms in our discussion of the scalar sector 20 . The explicit expressions for A i are given by and Z i (u, ω, q) are the solutions to Eq. (2.20) obeying the incoming wave boundary condition at the horizon and normalized to Z where Z i (u, ω, q) are the incoming wave solutions to Eq. (2.20).

The scalar channel
In this section, we extend the analysis of the scalar sector of metric perturbations performed in ref. [57] to second order in the hydrodynamic expansion. To that order, the retarded two-point function of the appropriate components of the energy-momentum tensor obtained by considering a linear response to metric perturbation has the form [2] G R,lin.resp. xy,xy (2.28) Using dual gravity, we compute the retarded Green's function G R xy,xy analytically for w 1 and q 1, and read off the transport coefficients τ Π and κ by comparing the result with Eq. (2.28). A novel feature at finite γ GB is the appearance of a new pole of the function G R xy,xy (ω, q) in the complex frequency plane [12]. The pole is moving up the imaginary axis with γ GB increasing. It is entering the region w 1 at intermediate values of γ GB and thus is visible in the analytic approximation.
To compute the two-point function in the regime of small frequency, we need a solution of the scalar channel differential equation (2.20) for w 1 and q 1. Using the variable v defined by the relation (2.21) and imposing the in-falling boundary condition [99] by isolating the leading singularity at the horizon via where G is a function of w and q of the form in Eq. (2.29) is the normalization constant. To find a perturbative solution g(v) for w 1, q 1, we introduce a book-keeping expansion parameter µ [102] and write where the functions g n satisfy the equations The functions H n are determined recursively from G and g m with m < n by where n ≥ 1. At first order, g 0 = 1 and g −1 = 0 which gives H 1 = −iw. A solution to Eq. (2.36) can be written in the form where C n and D n are the integration constants. In particular, for n = 1 we have Factorization (2.29) implies that the functions g n must be regular at the horizon (at v = 0). In the case of g 1 , the regularity condition leads to C 1 = −iw/2. Furthermore, all g n with n > 1 must vanish at the horizon (see Appendix C). For n = 1, this amounts to setting D 1 = 0. Hence, to linear order in w and q we have Repeating the procedure, we find the function g 2 (v): The functions g (w) 2 and g (q) 2 appearing in Eq. (2.40) are given by lengthy but closed-form expressions. Even though we do not have a closed-form expression for the remaining integral in Eq. (2.40), this is irrelevant for the purposes of computing the two-point function in the hydrodynamic limit, since the existing expression for g 2 is sufficient for fixing both the boundary conditions on g 2 itself and for determining the near-boundary expansion of Z 1 . More precisely, the integral in Eq. (2.40) comes from the outer integration in (2.37) and does not affect the regularity at the horizon thus allowing to fix the integration constant C 2 . The integral in (2.40) can be evaluated order-by-order in the near-boundary expansion of the integrand and the constant D 2 can be re-absorbed into the integration constant.
The full on-shell action (2.22) including the contact terms is given by where we used the near-boundary regulator u = → 0. Here, the first term is minus the four-volume V 4 times the free energy density (i.e. the pressure P), where which is consistent with Eqs. (2.8) and (2.7). The ellipsis denotes higher-order terms in w and q and terms vanishing in the → 0 limit. The retarded two-point function G R xy,xy (ω, q) can then be computed by evaluating the boundary action (2.41). Using the solution (2.29) to first order in w and q (i.e. including only the function g 1 in the expansion (2.34)) we find The Green's function has a pole on the imaginary axis at The approximation in Eq. (2.44) assumes γ GB 1. The pole is absent from the spectrum at λ GB = 0 (γ GB = 1) or, rather, it is located at complex infinity. At non-vanishing λ GB of either sign, the pole moves up the imaginary axis with |λ GB | increasing. For positive λ GB , it reaches the quasinormal frequency value at λ GB = 1/4 in that limit, determined analytically in section 2.4. For negative λ GB , the pole moves up to the origin. Its location is correctly captured by the small frequency perturbative expansion of the solution g(v) only for sufficiently large γ GB (see Fig. 1 and ref. [12] for details).
A small frequency expansion of Eq. (2.43) is A comparison with Eq. (2.28) gives the familiar expression for pressure (2.42) and the shear viscosity [57]   The full quasinormal spectrum of metric fluctuations in the scalar channel as a function of γ GB has been analyzed in detail in ref. [12]. The spectrum qualitatively differs from the one at λ GB = 0 in a number of ways, depending on the sign of λ GB . For λ GB > 0, there is an inflow of new quasinormal frequencies (poles of G R xy,xy (ω, q) in the complex frequency plane), rising up from complex infinity along the imaginary axis. At the same time, the poles of the two symmetric branches recede from the finite complex plane as λ GB is increased from 0 to 1/4, and disappear altogether in the limit λ GB → 1/4. The spectrum in this limit coincides with the one obtained analytically at λ GB = 1/4 in section 2.4 of the present paper. For λ GB < 0, on the contrary, the poles in the symmetric branches become more dense with the magnitude of λ GB increasing, and the two branches gradually lift up towards the real axis. They appear to form branch cuts (−∞, −q]∪[q, ∞) in the limit γ GB → ∞. For small q and very large γ GB , this would imply accumulation of poles of the Green's function in the region |w| 1. We have not investigated this limit in detail. Also, as noted above, there is at least one new pole (seen in Fig. 1) rising up the imaginary axis. The residue and the position of the pole w g contribute to the shear viscosity and to the position of the corresponding transport peak of the spectral function. A qualitatively similar phenomenon has been observed in the case of N = 4 SYM at large but finite 't Hooft coupling [12].

The shear channel
The energy-momentum tensor two-point functions G zx,zx , G tx,tx , G tx,zx in the shear channel can be expressed through the single scalar function G 2 as explained in ref. [102]. For example, 21 where the ellipsis represents the contact terms. In holography, the function G 2 is determined by the solution Z 2 (u, ω, q) (2.27) of the equation (2.20) obeying the appropriate boundary conditions, and by the relevant part of the on-shell boundary action (2.23). The retarded correlators in the shear channel are characterized by the presence of the hydrodynamic diffusive mode whose dispersion relation is given by where θ 1 is the transport coefficient of the third-order hydrodynamics introduced in ref. [20].
Higher terms in the momentum expansion of the shear mode depend on the (unclassified) fourth-and higher-order transport coefficients. Since the Gauss-Bonnet fluid is Weylinvariant ("conformal"), we have ε = 3P and thus η/(ε+P ) = (1 − 4λ GB ) /4πT = γ 2 GB /4πT . In holography, the quasinormal mode (2.48) can be found analytically by solving the equation (2.20) perturbatively for w 1, q 1: The coefficient in front of the term quadratic in momentum coincides with the one predicted by hydrodynamics of the holographic Gauss-Bonnet fluid with known shear viscosity. Since the coefficient τ Π is also known (e.g. from Eq. (2.28)), the quartic term in (2.49) allows one to read off the coefficient θ 1 : (2.50) In the dissipationless limit γ GB → 0 we have θ 1 ∼ γ 3 GB → 0. In fact, it can be seen numerically [12] that the full shear mode (2.48) approaches zero in the limit γ GB → 0. At γ GB = 0 (λ GB = 1/4), this mode disappears from the spectrum altogether due to the vanishing residue which is consistent with our analytic results for the spectrum at λ GB = 1/4 in section 2.4.
The full quasinormal spectrum was investigated numerically and partially analytically in ref. [12]. Its behavior as a function of λ GB is qualitatively similar to the one in the scalar channel, with the exception of one curious phenomenon: at fixed q, the new pole rising up the imaginary axis with (negative) λ GB increasing in magnitude, collides with the hydrodynamic pole (2.48) at some λ GB = λ c GB (q), and the two poles move off the imaginary axis. This is interpreted as breakdown of the hydrodynamic regime at a given q = q c (λ GB ). Curiously, the range of applicability of the hydrodynamic regime (i.e. the range q ∈ [0, q c ]) increases with the field theory "coupling" (understood as the inverse of |λ GB |) increasing [12].
The retarded correlation functions of the energy-momentum tensor in the shear channel can be computed from the boundary action (2.23). For the function G 2 in Eq. (2.47) we find 22 In the hydrodynamic approximation, to first non-trivial order in w, q, with both w ∼ µ 1 and q ∼ µ 1 scaling the same way, the shear channel solution to Eq. (2.20) obeying the incoming wave boundary condition is 2 is the normalization constant. We note that in order to obtain the hydrodynamic dispersion relation (2.49) that includes information about the second and the third order transport coefficients, we need to find Z 2 to one order higher, but using the scaling ω ∼ µ 2 and q ∼ µ is sufficient to extract the diffusive pole.
For the correlation function G 2 in the regime w 1, q 1 we thus find the following expression where ω g = 2πT w g (see Eq. (2.44)). At vanishing Gauss-Bonnet coupling λ GB = 0 (γ GB = 1) one has |w g | → ∞ and we formally recover 23 the standard result for N = 4 SYM at infinitely strong 't Hooft coupling and infinite N c [100,102] but it should be noted that the formula (2.53) is accurate only for |w g | 1, i.e. for sufficiently large γ GB . The correlator (2.53) has two poles with the following dispersion relations, expanded to q 2 : The first is the usual diffusive pole, corresponding to quadratic part of the dispersion relation (2.49), while the second pole is a new non-hydrodynamic pole coming from complex infinity at non-zero λ GB . This pole moves up the imaginary axis with γ GB increasing and is responsible for the breakdown of hydrodynamics in the large γ GB limit for any fixed non-zero value of q (see ref. [12] for details). The above expression for the Green's function and the dispersion relations are only valid in a (double expansion) regime in which not only w ∼ q 1 but also γ GB 1. The latter condition is required for the gapped mode on the imaginary axis to satisfy |w| 1. 22 As in refs. [100,102], we ignore possible contact terms coming from Sc.t.. See remarks in Appendix A of ref. [102]. 23 Upon the identification N 2 c = 4π 2 /κ 2 5 .
Note also that the form of the dissipative corrections implies that γ GB q 1. Obviously, these restrictions are only necessary if we are interested in analytic expressions.
The location of the momentum density diffusion pole confirms the result (1.11) for the shear viscosity of Gauss-Bonnet holographic fluid. We note that in the limit λ GB → 1/4 (γ GB → 0) the residue of the diffusion pole vanishes. The full Green's function can be determined numerically. The corresponding spectral function in the shear channel for various values of γ GB has been computed numerically in ref. [12].

The sound channel
The correlation functions in the sound channel can be expressed through the single scalar function 24 G 3 [102]. For example, for the energy density two-point function in the conformal case we have The hydrodynamic modes in the sound channel are the pair of sound waves whose dispersion relation is predicted by relativistic hydrodynamics up to a quartic term in spatial momentum: where c s = 1/ √ 3 is the speed of sound, Γ = 2η/3(ε + P ), ε + P = sT in the absence of chemical potential, and τ Π , θ 1 , θ 2 are transport coefficients of the second-and third-order (conformal) hydrodynamics in four space-time dimensions.
Solving the equation (2.20) for Z 3 perturbatively for w 1, q 1, imposing the incoming wave boundary condition at the horizon and the Dirichlet condition at the boundary, we find the hydrodynamic quasinormal mode 25 24 See footnote 21. 25 Here it is tacitly assumed that γGB is small enough. For moderate and large γGB, in addition to the mode (2.58), there exists another mode moving up the imaginary axis with γGB increasing. This mode enters the hydrodynamic domain w 1, q 1 for γGB ∼ 2 − 4 and can be seen analytically, as discussed in ref. [12].
Comparing the expansion (2.58) to the prediction (2.57) of conformal hydrodynamics one finds the same expressions for the shear viscosity -entropy density ratio and the secondorder transport coefficient τ Π as the ones reported in Eqs. (1.19) and (1.20). This agreement is gratifying but more analytic work is needed to extend the expansion (2.58) to quartic order and determine the coefficient θ 2 of the third-order hydrodynamics. Other features of the quasinormal spectrum are qualitatively similar to the scalar case and are discussed in full detail in ref. [12].
The coefficients in front of the quadratic, qubic and possibly 26 quartic terms in the dispersion relation (2.57) vanish in the limit γ GB → 0. This limit is hard to study numerically but it is conceivable that the higher terms vanish as well leaving the linear propagating mode w = ±q/ √ 3. Such a mode, however, is absent in the exact spectrum at γ GB = 0 (see section 2.4).
To first order in the hydrodynamic expansion, the gauge-invariant mode is given by and we have used U 2 = u 2 +γ 2 GB −u 2 γ 2 GB . The correlation function G 3 can then be computed from giving As required by rotational invariance, G 1 (ω, 0) = G 2 (ω, 0) = G 3 (ω, 0) [102]. The contact term in the on-shell action (2.23) relevant for the computation of G tt,tt (ω, q) is (2.64) The full retarded energy density two-point function is then (2.65) 26 Possibly, because the expression for θ2 remains unknown.
The thermodynamic (equilibrium) contribution has been omitted from this expression. To this order in the hydrodynamic expansion, the spectrum contains three modes, (2.67) The first two are the attenuated sound modes (2.58) and the third mode is the gapped mode similar to those in the scalar and shear channels. As in the shear channel, these results require the following scalings to be respected: w ∼ q 1, γ GB 1 and hence, γ GB q 1. Second-order corrections to the two hydrodynamic sound modes were given by Eq. (2.58).
To study the spectrum beyond second-order hydrodynamics and investigate higher-frequency spectrum, we must again resort to numerics. We note that for better control over the numerics, it is useful to follow [106] and write we obtain the following expression convenient for the computation of quasinormal modes: The coefficient h can be found analytically, h = − (1 + γ GB ) 4 w 2 − q 2 2 /32. For a detailed discussion of the quasinormal spectrum, see ref. [12]. A comprehensive analysis of the large spatial momentum asymptotics similar to the one accomplished for the strongly coupled N = 4 SYM in refs. [107,108] would be of interest but has not been attempted neither in ref. [12] nor in the present paper.

Exact quasinormal spectrum at λ GB = 1/4
At λ GB = 1/4, the equations of motion (2.20) for all channels simplify drastically. They reduce to the following system Scalar channel: Shear channel: Sound channel: and not α 1,2 = {0, 2}, which are their values for any λ GB < 1/4 (and in fact for all fivedimensional bulk fluctuations dual to operators of conformal dimension ∆ = 4 of a 3 + 1dimensional boundary theory). The standard holographic dictionary then implies that at λ GB = 1/4 the dual theory operators scale as the energy-momentum tensor in six rather than four dimensions. Technically, the reason for this "dimensional transmutation" is related to the fact that the "standard" terms in the wave equations (2.20) most singular in the limit u → 0 are multiplied by the coefficients proportional to (1−4λ GB ) and thus vanish at λ GB = 1/4. At this value of the Gauss-Bonnet coupling, the theory becomes "topological gravity" [109] with a number of curious properties. 27 In particular, thermodynamic properties of the black brane solution at λ GB = 1/4 are different from the ones at λ GB < 1/4 [110]. The underlying physical reasons and significance of this limit are not entirely clear to us, although they might be related to the issues discussed in ref. [111] and refs. [112][113][114].
We note that the Gauss-Bonnet black brane metric is regular at λ GB = 1/4: Rescaling the coordinates t, x, y, z and the parameter L, it can be brought into the form Shear: 78) 27 We thank the referee for bringing refs. [109] and [110] to our attention.
where Ω ≡ −1 − iw 2 . Given the three solutions, the quasinormal spectrum is determined analytically by imposing the Dirichlet condition Z i (0) = 0 at the boundary. We find Scalar: Sound: where n 1 and n 2 are independent non-negative integers. The numerical study of the Gauss-Bonnet quasinormal spectrum in ref. [12] shows that in the limit λ GB → 1/4 the quasinormal frequencies approach the ones found above. The spectrum in the shear channel is q-independent. In the scalar and sound channels, for sufficintly large q the modes cross into the upper half plane of frequency thus signaling an instability. This is perhaps not surprising given the causality problems in the boundary theory observed for sufficiently large spatial momentum in ref. [57] and other publications. Finally, let us address the questions of what happens to the hydrodynamic poles in the limit of λ GB → 1/4 (γ GB → 0). By examining the limit of the sound correlator G tt,tt (ω, q) given by Eq. (2.65) computed for any generic value of γ GB (or the limit of G 3 given by Eq. (2.63)), we find a non-vanishing Green's function with an unattenuated sound mode, ω = ±q/ √ 3. On the other hand, the sound spectrum computed analytically at γ GB = 0 (cf. Eq. (2.81)) contains no such mode. This situation can be contrasted with the shear channel: there, the correlator G 2 (cf. (2.53)) vanishes in the same limit and there is no remaining diffusive mode in the spectrum. Consistently, the exact quasinormal spectrum at γ GB = 0 (cf. (2.80)) contains no mode at w = 0, either.
We do not have a full understanding of this phenomenon but can offer the following comments. Examine more closely the limit γ GB → 0 of the sound correlator G 3 (2.62). First, we notice that its Z 3 -independent prefactor gives different expressions depending on which of the two limits, γ GB → 0 or → 0, is taken first. Namely, where the ellipses denote terms subleading in the expansions of and γ GB around zero. Now, in the expansion around γ GB = 0, the Fröbenius series (2.68) becomes (2.84) By first taking and then γ GB to zero (the order of limits we took to find G 3 in Eq. (2.63)), one again recovers the leading order hydrodynamic expression (2.85) With the opposite order of limits, the prefactor (2.83) and the solution (2.84) yields where A and B depend on γ GB . What this expression reveals is that it is possible for the unattenuated sound mode to be a pole of the Green's function, having entered into the expression from the prefactor, not the ratio of B/A. Thus, such a pole would not appear as a part of the quasinormal spectrum.

The limit λ GB → −∞
It is tempting to investigate the limit λ GB → −∞ analytically to confirm the observations based on numerical simulations. However, taking this limit is problematic for two reasons. First, on a technical level, the equations of motion for fluctuations contain products of the type λ GB (r − r + ) which remain finite for r sufficiently close to the horizon r + , even at large |λ GB |. This can possibly be dealt with by a variable redefinition but the second problem is more serious. The Kretschmann curvature invariant evaluated on the black brane solution For λ GB ∈ [0, 1/4], the curvature singularity in Eq. (2.87) is at r = 0. However, for λ GB < 0 the curvature singularity is located at Thus, as λ GB is tuned from 0 to −∞, the curvature singularity moves continuously from r = 0 to the horizon r = r + and becomes a naked singularity 28 in the limit λ GB → −∞.
Because the classical background geometry is singular at the horizon, considering classical metric fluctuations in the limit λ GB → −∞ would be meaningless. In some sense, the need for an ultraviolet completion of gravity in this limit is in accord with the observations made in ref. [12] and in the present paper that the regime of large negative λ GB qualitatively corresponds to the regime of weak coupling in the field theory which generically requires the full dual stringy rather than dual gravity description. As a curious observation, we note the following. In the large (negative) λ GB expansion, the Ricci scalar evaluated on the solution (2.2) to leading order becomes 89) 28 The appearance of naked singularities in the solutions of Lovelock gravity has been investigated in ref. [115]. In fact, all three curvature scalars that appear in the Gauss-Bonnet term, R µνρσ R µνρσ , R µν R µν and R 2 , are singular at r = r + and scale as 1/λ GB , while their combination that appears in the action remains finite and independent of r: (2.91) As as result of these scalings, the Einstein-Gauss-Bonnet action (1.16) to leading order in λ GB reduces to the Gauss-Bonnet term and the cosmological constant Λ = −6/L 2 : (2.92) This theory has a black brane solution that coincides with the λ GB → −∞ limit of the solution (2.2), where we have introduced a rescaled radial coordinate r = (−λ GB ) 1/4r .

Gauss-Bonnet transport coefficients from fluid-gravity correspondence
From the analysis of quasinormal spectra and retarded two-point functions in section 2, we were able to determine non-perturbative expressions for the Gauss-Bonnet transport coefficients η, τ Π , κ (and also θ 1 of the third-order hydrodynamics). To find the remaining transport coefficients, one can use either the fluid-gravity correspondence or the Kubo formulae applied to three-point functions. In this section, we shall use the fluid-gravity methods [3,116]. Previously, fluid-gravity approach has been used to determine the shear viscosity [117] and second-order hydrodynamic coefficients [49] of Gauss-Bonnet holographic liquid perturbatively in λ GB . Fluid-gravity correspondence uses the fact that the bulk metric perturbations h µν source the energy-momentum tensor T µν in the generating functional of the boundary quantum field theory [118,119]. Gravitational bulk action should thus be able to capture all of the energy-momentum properties of the dual theory. The procedure for computing the holographic energy-momentum tensor, inspired by the old prescription of Brown and York [120], was proposed in ref. [121]. One expects then that in the appropriate variables a gradient expansion of the bulk metric should capture the hydrodynamic gradient expansion of the dual field theory's energy-momentum tensor. Following ref. [116], we write the Gauss-Bonnet black brane background solution (2.2) in the Eddington-Finkelstein coordinates, where N GB is given by Eq. (2.4). We set L = 1 for convenience and defined b ≡ 1/r + to be consistent with the notations used in ref. [116]. The function f (br) is The energy-momentum tensor is given by the expression where all the ingredients are defined just below Eq. (2.13). The next step is to boost the brane solution (3.1) along a space-time dependent velocity four-vector u a (x), where with i = 1, 2, 3 corresponding to the spatial boundary coordinates. Note that x a = (v, x, y, z) in Eddington -Finkelstein coordinates. The boosted black brane metric, which we denote by g (0) µν , becomes Generically, the metric (3.5) is no longer a solution of the Einstein-Gauss-Bonnet equations of motion (2.1). In fluid-gravity correspondence, assuming a slow-varying dependence of the coefficients on the coordinates x a and making a gradient expansion, one imposes the equations of motion (2.1) as the condition each term in the expansion must satisfy. We make a gradient expansions in the derivatives of the fields β i (x a ) and b (x a ) to second order, in agreement with the boundary theory's standard second-order hydrodynamic gradient expansion in velocity and temperature fields (see e.g. Appendix B). To second order, the metric will have the form where g (0) µν and g (1) µν are expanded up to terms involving two derivatives of b and β i inclusive. We shall use as a book-keeping parameter in the derivative expansion.
The procedure of solving equations (2.1) order by order is greatly simplified, if one notices that it is sufficient to solve the equations of motion locally around some point x a = X a . The global metric can be obtained from these data alone [116]. The local expansions of the fields b and β i are given by We choose to work in a local frame at the origin, X a = 0, where b 0 = 1 and Furthermore, it is consistent to choose a gauge with β i (1) = 0 at x a = X a [116].

First-order solution
The most general expression for the first-order metric g (1) µν can be conveniently written in a scalar-vector-tensor form where x i = (x, y, z), k 1 and h 1 are scalars, j i 1 is a three-vector and A ab is a tensor. As discussed above, we proceed by using the expanded forms of b and β i given in (3.7) and (3.8) to write the order-metric as g µν = g First, we solve the Dynamical equation 1 in (3.13) for h 1 (r). We then use Constraint 2 in (3.12) which relates k 1 (r) to h 1 (r) to solve for k 1 (r). Constraints 1 and 3 in (3.11) and (3.14) give Finally, we can solve the remaining Dynamical equations 2 and 3 in (3.15) and (3.16) to find j 1 (r) and the tensor A ab which contains information about shear viscosity.
The global first-order metric, g µν = g (0) µν + g (1) µν , can be written as a sum [116], where the six line elements A n are defined as

20)
The last step is to find the function F 0 (r) entering Eq. (3.20). This function is part of the tensor A ab satisfying Eq. (3.16). Explicitly, the second-order differential equation for F 0 is where F 1 (a, b, b ; c; w, z) is the Appell hypergeometric function of two variables and where 2 F 1 (a, b; c; z) is the Gauss hypergeometric function. The expansion of the Appell function at r → ∞ (explicitly written here for 0 < γ GB < 1) can be found by using the theorems in ref. [122]: (3.24) The result (3.24) allows us to find the expansion of F 0 (r) near the boundary, valid to order O(r −4 ) which is sufficient for the purposes of computing the boundary energy-momentum tensor. Substituting F 0 (r) into the first-order metric g µν (1) and computing the energy-momentum tensor (3.3) with the full first-order solution we recover the non-perturbative result for the shear viscosity η presented in (2.46).

Second-order solution
The second-order correction g (2) µν is computed in a similar way: first, we perturb g to second order in derivative expansion and then find g (2) µν requiring that the Einstein-Gauss-Bonnet equations of motion (2.1) are satisfied.
To find the second-order transport coefficients non-perturbatively, we would need to solve differential equations with the differential operator given by the left-hand side of Eq. (3.22) and the right-hand sides involving integrals over the Appell function (3.24). This program faces a certain technical challenge, and we were not able to find closed-form expressions for the transport coefficients in this way. It is possible, however, to obtain terms of the perturbative expansion of transport coefficients in γ GB and thus check the fully nonperturbative results (1.20), (1.22), (1.23) found by using the method of three-point functions (see ref. [86] and section 4), as well as perturbative results by Shaverin [49].
A convenient ansatz for the line element of the second-order metric g (2) µν is suggested by the tensor structure of second-order hydrodynamics (see e.g. Appendix B): 26) where x i = (x, y, z), k 2 and h 2 are scalars, j i 2 is a three-vector. We have also defined At this point, we can focus only on the four functions P n , n = {0, 1, 2, 3}, which will give us the four second-order coefficients, λ 0 ≡ ητ Π , λ 1 , λ 2 and λ 3 , respectively. Since the boundary theory is defined in flat space-time, this procedure will not allow us to find the coefficient κ. Furthermore, we know that in Landau frame there are no other transport coefficients coming from either the scalar or the vector sector. Still, we need to use the constraint equation r 2 f 0 (r)E rr + N GB E vr = 0 and the dynamical equation E rr = 0 to eliminate h 2 , k 2 and their derivatives from the dynamical equations for P n . The remaining differential equations for P n can be solved perturbatively to an arbitrarily high order in λ GB . Here we outline what we believe is the most efficient way to extract information from the functions P n necessary to recover the four transport coefficients λ 0,1,2,3 . First, the functions P n are expanded in series near the boundary as Then the metric (3.26) with P n expanded as in Eq. (3.31) is substituted into the full secondorder metric, and the energy-momentum tensor (3.3) is computed. The main observation is that in the limit r → ∞, finite contributions to T µν only depend on the coefficients of P n proportional to r −4 , i.e. T µν depends on p 3 and p 4 . In order to find the four coefficients, we use the fact that all four differential equations for P n (r) can be written in the form of Eq. (3.22), i.e. as ∂ r [Q(r)∂ r P n (r)] − R n (r) = 0, (3.32) where Q and R expanded to the desired order in λ GB , and the function Q is the same in all four cases. The differential equations can be formally solved, as in Eq. (2.37), by writing Fortunately, in Eq. (3.33) it is sufficient to take the inner integral over r whose integrand depends on F 0 (r) expanded to the desired order of λ GB . Integration constants C n are fixed by requiring regularity at the horizon. The coefficients D n may remain undetermined since we only need the specific terms in the r → ∞ expansion. Thus, using the expansion (3.31) in the differential equations (3.33) we find all p (4) n . For example, from the equation obeyed by P 0 we obtain p (1) 0 which allows us to find p (4) 0 to the desired order in λ GB by setting to zero the coefficient in front of r −5 .

Transport coefficients
Once the coefficients p (4) n are known, the full second-order metric can be used to determine the expansion of the energy-momentum tensor (3.3) near r → ∞ and read off the transport coefficients ητ Π , λ 1 , λ 2 and λ 3 from the coefficients of tensors (3.27) -(3.30). The results are in exact agreement with the corresponding terms of the λ GB -expansions of the four non-perturbative second-order transport coefficients given by Eqs. (1.20), (1.22), (1.23) and (1.24), as well as with those computed in ref. [49] to linear order 29 .
The conclusion of this section is that fluid-gravity duality applied to Gauss-Bonnet holographic fluid allows to determine all the transport coefficients of second-order hydrodynamics, except κ, but only the shear viscosity η is determined non-perturbatively in λ GB : the coefficients τ Π and λ 1,2,3 are found only as series in λ GB , due to technical problems related to evaluating integrals of Appell function. Finally, we note that within the fluidgravity approach one is able to check the Haack-Yarom relation order by order in λ GB and find that it is violated at quadratic order as shown in Eq. (1.26).

Gauss-Bonnet transport from three-point functions
The full non-perturbative expressions for the Gauss-Bonnet transport coefficients can be found by computing the three-point functions 30 of the energy-momentum tensor in the hydrodynamic approximation and using the Kubo-type formulae derived in refs. [29,126,127]. The retarded three-point functions are defined following the recipes of the Schwinger-Keldysh closed time-path formalism [128,129]. Part of the material in this section has some overlap with refs. [30,86] and is included here for convenience and continuity.

An overview of the method
In the Schwinger-Keldysh formalism, given a Lagrangian L [φ, h], where φ collectively denotes matter fields and h is a metric perturbation around a fixed background g, the degrees of freedom are doubled: φ → φ ± , g → g ± , h → h ± , where the index ± labels the fields defined either on a "+"-time contour running from t 0 towards the final time t f > t 0 or the "−"-time contour, where the time runs from t f backwards to t 0 . When the theory is considered at finite temperature T = 1/β, the two real-time contours can be joined by a third, imaginary time, contour running between t f and t f − iβ. Fields defined on this imaginary time contour will be denoted by ϕ. The generating functional of the energy-momentum tensor correlation functions is given by For all fields, it will be convenient to use Keldysh basis φ R = 1 2 (φ + + φ − ) and φ A = φ + −φ − . Upon computing the variation, classical expectation values obey φ + = φ − . Thus, all fields with an index A will vanish and one can define T ab ≡ T ab R : .
(4.2) 29 In matching those expressions, one should recall that the horizon scale r+ in the fluid/gravity calculation is promoted to a field b(r), with b0 fixed by Eq. (3.9). 30 In holography, the first equilibrium real-time three-point and four-point functions in strongly coupled N = 4 SYM at finite temperature were computed in refs. [123][124][125].
The expectation value of T R at x = 0 can be expanded as where G RAA... denote the fully retarded Green's functions [130] obtained by 31 where the ellipses indicate further insertions of ∂h R in the expression with the derivatives as well as the T ab A insertions into the n-point function on the right-hand side of Eq. (4.4). We follow refs. [29,126] and use Kubo formulae for pressure and transport coefficients of a conformal fluid derived by exciting fluctuations of the relevant metric components. Choosing the spatial momentum along the z direction, one turns on h xy , h xz and h yz perturbations to obtain (4.7) By turning on h xy , h tx and h ty components, we find ∂ 2 ∂p z ∂q z G xy,tx,ty RAA (p, q), (4.8) and, finally, by considering h xy , h ty and h xz perturbations, we obtain ∂ 2 ∂p 0 ∂q z G xy,ty,xz RAA (p, q). (4.10) A consistency check on our calculations is provided by the following two Kubo formulae which both give the expression for pressure: G xy,tx,ty RAA (p, q). Note that our definitions of transport coefficients are the same as in ref. [2] (see Appendix B for a digest of notations and conventions used in the literature). 31 See e.g. ref. [126].

The three-point functions in the hydrodynamic limit
The three-point functions are computed by solving the Einstein-Gauss-Bonnet equations of motion (2.1) to second order in relevant perturbations, where the book-keeping parameter is used to indicate the order of perturbation. The Dirichlet condition h (2) µν = 0 is imposed at the boundary [29]. Once the bulk solutions are found, one should take the triple variation of the on-shell action with respect to the boundary values h µν (r → ∞) to find the correlators. A simplifying feature of this procedure is that since equations of motion are solved to order 2 , only the boundary term contributes to the three-point function, and hence no bulk-to-bulk propagators appear in the calculation.
To compute the three-point functions used in the Kubo formulae above, we need to turn on the following sets of metric perturbations: 2) h xy = h xy (r)e i(p z +q z )z , h tx = h tx (r)e ip z z , h ty = h ty (r)e iq z z , (4.14) 3) Here, we outline the steps leading to obtaining the three-point function G xy,xz,yz RAA . First, we find the bulk solutions for h (1) xy , h (1) xz and h (1) yz imposing the standard incoming wave boundary condition at the horizon and the condition h µν at the boundary. As in section 2.1, it will be convenient to work with the radial variable v defined by Eq. (2.21).
Since the metric fluctuations in the set (4.13) are independent of the spatial momentum, all three of them obey the same 32 differential equation (2.20), and thus h (1) xy , h (1) xz and h (1) yz will have the same functional dependence on v. Moreover, we can use the solution to the equation already obtained in section 2.1, with q set to zero and the relevant frequencies, p 0 + q 0 , p 0 and q 0 , inserted instead of w, respectively. Thus, for h (1) xy we find the expression and similar formulas for h (1) xz and h (1) yz . We can deal with the remaining integral in Eq. (4.16) in the same way as in section 2.1, by integrating order-by-order in the near-boundary expansion.
Next, we need to look for the second-order solution h (2) xy , which includes the first-order metric back-reaction. The differential equation again has the form of Eq. (2.35) and can be solved using the same methods. The relevant part of the solution takes the following form: with a complicated and unilluminating expression for h(v) not shown here explicitly.
With the second-order solution in hand, we substitute the resulting formula for g µν + r 2 h The other three-point functions, G xy,tx,ty RAA and G xy,ty,xz RAA , are computed using the same procedure, with the differential equations always taking the form of (2.20). The only difference is that we cannot impose the in-falling boundary conditions on perturbations h tx or h ty in Eq. (4.14), and similarly on h ty in Eq. (4.15), because they only fluctuate in the z-direction and not time. Regularity then demands setting h tx = h ty = 0 at the horizon. Consequently, h xy in Eq. (4.14) also needs to vanish at the horizon. 33

Second-order transport coefficients and the zero-viscosity limit
Having computed in the hydrodynamic approximation the three-point functions G xy,xz,yz RAA , G xy,tx,ty RAA and G xy,ty,xz RAA , we can use the Kubo formulae to compute pressure (4.11), shear viscosity (4.5) and all second-order transport coefficients (4.6) -(4.10). The result for pressure coincides with the one in Eq. (2.42), and the shear viscosity is confirmed to be 33 The full expressions for the other two three-point functions are very cumbersome and will not be written here explicitly. For an example of a technically simpler but conceptually identical calculation in N = 4 SYM theory, see refs. [29,30].
The derivatives of the coefficients with respect to λ GB are shown in Fig. 3. The coefficients λ 3 and κ are monotonically decreasing functions of λ GB as can be seen from whereas the coefficients ητ Π , λ 1 and λ 2 are not.

Charge diffusion from higher-derivative Einstein-Maxwell-Gauss-Bonnet action
Can first-order transport coefficients other than shear viscosity be tuned to zero with a suitable choice of higher derivative bulk terms, and can this be done simultaneously with tuning to zero the viscosity? In this section, we compute non-perturbative corrections to the well known result for the U (1) charge diffusion constant at infinite coupling [100] in a hypothetical boundary theory dual to Einstein-Maxwell-Gauss-Bonnet gravity with the charge neutral black brane background (2.2).

The four-derivative action
We are interested in the four-derivative Einstein-Maxwell-Gauss-Bonnet action whose equations of motion involve at most second derivatives. Such theories were previously considered in refs. [131,132], and in the context of an effective target-space heterotic string theory action in [133]. 34 The higher-derivative Maxwell terms may appear as a result of compactification, e.g. of a higher-dimensional Gauss-Bonnet action. Here we construct the necessary action directly. We begin by considering the Einstein-Maxwell-Gauss-Bonnet theory with the most general four-derivative Maxwell field Lagrangian, where Λ = −6/L 2 , the Gauss-Bonnet Lagrangian L GB is given by Eq. (1.16), and The coupled equations of motion for g µν and A µ following from the action (5.1) are written in Appendix E. To make third-and fourth-order derivatives of the fields vanish in the equations of motion (E.1), we must impose the following constraints on the coefficients α n α 4 = α 6 , 8α 4 + α 5 − 4α 6 = 0, where we have defined the dimensionless couplings β 1 ≡ α 4 /L 2 , β 2 ≡ α 7 /L 2 , β 3 ≡ α 11 /L 2 and β 4 ≡ α 8 /L 2 . To simplify the Lagrangian further, we notice that the term proportional to β 4 can be rewritten as hence the entire expression vanishes due to the cyclic property of the Riemann tensor. Thus the Lagrangian L A leading to second-order equations of motion is given by Therefore, there are altogether four parameters, λ GB , β 1 , β 2 and β 3 , entering the secondorder equations of motion of the theory. One may wonder if a black hole (brane) solution with non-perturbative values of these parameters exists. The black brane metric (2.2) is automatically a solution of the theory when A µ = 0. It is also possible to find perturbative corrections in β 1 , β 2 and β 3 to the five-dimensional AdS-Reissner-Nordström metric. However, we were not able to find a generalization of the solution (2.2) with non-trivial A µ and fully non-perturbative non-vanishing β 1 , β 2 and β 3 . 35

The U (1) charge diffusion constant
To compute the charge diffusion constant in a hypothetical neutral liquid dual to the bulk action constructed in the previous section we follow the procedure outlined in [102]. We begin by perturbing the trivial A µ = 0 background vector field as A µ → A µ + a µ and writing the electromagnetic field strength corresponding to the linearized perturbation as F = da. Given the trivial background A µ = 0, the metric fluctuations decouple from a µ and can be set to zero. 36 In the equations of motion, the terms proportional to α 7 and α 11 (i.e. β 2 and β 3 ) only contribute to quadratic or higher orders in the expansion in . Hence, they will not contribute to the charge diffusion constant. The linearized equations obeyed by a µ read Vector field fluctuations can be decomposed into transverse and longitudinal modes, with charge diffusion coming from the low-energy hydrodynamic excitations in the longitudinal sector. By choosing the spatial momentum along the z-direction, the relevant gaugeinvariant variable in the longitudinal sector is We use the variable u = r 2 + /r 2 , with the boundary at u = 0 and horizon at u = 1. Then we impose the incoming wave boundary condition required for the calculation of retarded correlators [99] by writing 11) 35 An asymptotically AdS black hole solution to the theory considered in this section with β1 = 0 was found in an integral form and studied in [131]. Unfortunately, for β = 0 the equations are significantly more complicated. In particular, in the relevant metric ansatz, ds 2 = −e 2λ dt 2 + e 2ν dt 2 + . . ., the relation λ = −ν is no longer true. 36 Charge diffusion in a three dimensional boundary theory, including the β1 term, was computed in a neutral Einstein-Hilbert black brane background in [135].
where the function Z 4 (u) regular at the horizon can be found perturbatively in µ 1, with q and w scaling as w → µ 2 w and q → µq. We find it useful to introduce a new variable w , so that u = w 2 − γ 2 GB / 1 − γ 2 GB . The boundary is now at w = γ GB and horizon at w = 1. At order O(µ 0 ), the function Z 4 can be written as Z 4 = C 1 + C 2 z(w ), where z(w ) is a solution of the equation We solve for z(w ) and impose the boundary conditions z(γ GB ) = 1 and z(1) = 0. The constant C 2 can then be expressed as a function of C 1 , w, q and other parameters of the theory by substituting z(w ) into the original differential equation, expanding to order O(µ 2 ) and imposing regularity at the horizon. The hydrodynamic quasinormal mode can be found by solving the equation Z 4 (w, q) = 0 at the boundary for w. The dispersion relation has the form where D is the charge diffusion constant of the dual theory. For the Gauss-Bonnet coupling in the interval λ GB ∈ [0, 1/4] (1 ≥ γ GB ≥ 0) we find 37 where β ≡ 1 + 48β 1 and γ GB ≡ √ 1 − 4λ GB . We can now consider various limits. For the two-derivative Maxwell field in Gauss-Bonnet background, i.e. for β 1 = 0 (β = 1), we find the expression 37 It is also possible to write an explicit formula for D in the interval λGB < 0 (γGB > 1) but here we are mostly interested in the dissipatinless limit λGB → 1/4.

Conclusions
Together with refs. [12,30,86], the present paper is an attempt at a comprehensive investigation of the second-order transport properties, energy-momentum tensor correlation functions and quasinormal spectrum in the Gauss-Bonnet holographic fluid in D = 3 + 1 dimensions non-perturbatively in Gauss-Bonnet coupling λ GB . The existence of a strongly coupled CFT dual to classical non-Einsteinian gravity such as Gauss-Bonnet gravity at finite λ GB would be an interesting alternative to the standard scenario of gauge-gravity duality. However, the work of Camanho et al. [78] appears to cast a serious doubt on such a possibility, reducing the status of the curvature-squared terms to that of a perturbative correction. At the same time, we have not found any obvious pathology in hydrodynamic properties of the hypothetical dual field theory at finite λ GB . The curvature-squared terms are interesting even as corrections to the Einstein's gravity description of a dual field theory, the second-order nature of the Gauss-Bonnet equations of motion making it easier to search for the new features such as the extra poles of the correlators not seen at λ GB = 0. The analysis of gravitational quasinormal spectrum in ref. [12] and in the present paper shows that the analytic structure of dual thermal correlators is qualitatively different depending on the sign of λ GB (understood as inverse coupling), with the λ GB < 0 case showing "normal" (e.g. qualitatively similar to N = 4 SYM at finite 't Hooft coupling and infinite N c and having a potential to connect to the kinetic regime) features, and the λ GB > 0 case demonstrating various anomalies (whose precise meaning remains to be understood, possibly invoking various monotonicity arguments). 38 On the other hand, constraints on λ GB may come from different considerations such as the recent argument for λ GB > 0 in ref. [83] based on unitarity. Fortunately, corrections coming from R 2 and R 4 higher derivative terms seem to be very similar [12] in uncovering a qualitative picture of transport and other properties at large but finite coupling. (ERC Grant agreement 307955). A.O.S. would like to thank the Institute for Nuclear Theory at the University of Washington (Seattle, USA) and the Mainz Institute for Theoretical Physics (Mainz, FRG) for their kind hospitality and partial support during the completion of this work.
A Second-order transport coefficients of N = 4 SYM at weak and strong coupling For the finite-temperature N = 4 SU (N c ) supersymmetric Yang-Mills (SYM) theory in d = 3 + 1 dimensions in the limit of infinite N c and infinite 't Hooft coupling λ = g 2 Y M N c , first-and second-order transport coefficients were computed, correspondingly, in [137] and [2,116] using methods of gauge-gravity and fluid-gravity dualities: Coupling constant corrections to the coefficients (A.1), (A.2) can be determined from the higher-derivative terms in the low-energy effective action of type IIB string theory where γ = α 3 ζ(3)/8, L is the AdS curvature scale, and the ratio α /L 2 is related to the value of the 't Hooft coupling λ in N = 4 SYM via α /L 2 = λ −1/2 . The effective fivedimensional gravitational constant is connected to N c by κ 5 = 2π/N c . The term W is given in terms of the Weyl tensor C µνρσ by Corrections to all first and second-order transport coefficients are known [23][24][25][26][27][28][29][30]: Leading order results for the third order coefficients θ 1 and θ 2 entering the hydrodynamic dispersion relations are known as well [20]: Additional explicit results for the linear combinations of N = 4 SYM third order coefficients can be found in ref. [20]. Other coupling constant corrections to the results at infinitely strong t'Hooft coupling in finite temperature N = 4 SYM include corrections to the entropy [31,32], photon emission rate [33], and poles of the retarded correlator of the energymomentum tensor [34].
In N = 4 SYM at weak coupling, the shear viscosity has been computed in ref. [138]. The second-order transport coefficients in various theories at weak coupling (QCD with either 0 or 3 flavours, QED, λφ 4 ) were determined by York and Moore [139]. In conformal kinetic theory (at weak coupling) one finds 2ητ Π + λ 2 = 0 [2,139,140]. Curiously, in the theories considered in [139] the Haack-Yarom relation (1.7) at weak coupling can be expressed as 13) where and P are energy density and pressure, correspondingly, and C 1 and C 2 are theoryspecific constants (e.g. C 1 ≈ 6.10517, C 2 ≈ 6.13264 for λφ 4 theory). It appears that at weak coupling one finds H = 0 also in other (nearly conformal) examples [141]. It would be interesting to compute H(λ) directly in N = 4 SYM at weak coupling. Another interesting finding of ref. [139] is that at weak coupling the coefficients κ and λ 3 vanish to order ∝ T 2 /λ 4 (but may be non-zero at ∝ T 2 /λ 2 ). We note that λ 3 = 0 in the limit λ → ∞ but has a non-trivial coupling dependence as can be seen from (A.10).
The energy-momentum tensor of a neutral conformal relativistic fluid considered in the Landau frame is written as T ab = εu a u b + P ∆ ab + Π ab , (B.1) where ∆ ab ≡ g ab + u a u b , pressure P and energy density ε are related by the conformal fluid equation of state in four dimensions, P = ε/3, and where D ≡ u a ∇ a . We use the following definitions (in our case, d = 4) where by construction the resulting tensors are transverse, u a A ab = 0, traceless, g ab A ab = 0, and symmetric. The tensor σ ab is a symmetric, transverse and traceless tensor involving first derivatives of the velocity field The vorticity Ω µν is defined as an anti-symmetric, transverse and traceless one-derivative tensor The Haack-Yarom relation in our notations reads whereas the conformal kinetic theory result [139] is Notations and conventions used in refs. [116], [87] We label the objects used in refs. [87,116] with the letter "R", e.g.
We note that in the paper by Shaverin and Yarom [49], the notations for σ µν SY , vorticity ω µν SY and their coupling λ SY 2 σ µ α,SY ω αν SY are the same as in ref. [47]. The relations between our transport coefficients (i.e. the ones in [2], [139], [30]) and the ones used in [49] are C Boundary conditions at the horizon in the hydrodynamic regime In this Appendix, we clarify the procedure of imposing the incoming wave boundary condition at the horizon on a (gauge-invariant) fluctuation Z given by a perturbative series in the hydrodynamic regime (w 1 and q 1). Consider such a solution Z 1 near the horizon u = 1: where q is ignored for simplicity. Here, the function F (regular at u = 1 by Fröbenius construction) is found perturbatively as a series in w 1, F (u, w) = F 0 (u) + wF 1 (u) + w 2 F 2 (u) + · · · , (C.2) where F i (u) satisfy the equation of motion obeyed by Z to a given order in w, with F i (1) = S i for i ≥ 0, and S i are constants independent of u and w. Now consider another solution, Z 2 , near u = 1, where C(w) is a function of w only, and G is found perturbatively by solving the differential equation obeyed by Z by a series in w 1, where the coefficients are given by P 2 = rf λ GB f − r , (D.10) where f (r) is given by Eq. This coordinate is more convenient for taking the limit of zero temperature.

E Equations of motion of Einstein-Maxwell-Gauss-Bonnet gravity
The equations of motion of Einstein-Maxwell-Gauss-Bonnet gravity following from the action (5.1) form a system of two coupled PDEs: R µν − 1 2 g µν R + g µν Λ = T GB µν + 2κ 2 5 T A µν , (E.1) Here, the gravitational energy-momentum tensor term is given by and the Maxwell field contribution has the form