From strong to weak coupling in holographic models of thermalization

We investigate the analytic structure of thermal energy-momentum tensor correlators at large but finite coupling in quantum field theories with gravity duals. We compute corrections to the quasinormal spectra of black branes due to the presence of higher derivative $R^2$ and $R^4$ terms in the action, focusing on the dual to $\mathcal{N}=4$ SYM theory and Gauss-Bonnet gravity. We observe the appearance of new poles in the complex frequency plane at finite coupling. The new poles interfere with hydrodynamic poles of the correlators leading to the breakdown of hydrodynamic description at a coupling-dependent critical value of the wave-vector. The dependence of the critical wave vector on the coupling implies that the range of validity of the hydrodynamic description increases monotonically with the coupling. The behavior of the quasinormal spectrum at large but finite coupling may be contrasted with the known properties of the hierarchy of relaxation times determined by the spectrum of a linearized kinetic operator at weak coupling. We find that the ratio of a transport coefficient such as viscosity to the relaxation time determined by the fundamental non-hydrodynamic quasinormal frequency changes rapidly in the vicinity of infinite coupling but flattens out for weaker coupling, suggesting an extrapolation from strong coupling to the kinetic theory result. We note that the behavior of the quasinormal spectrum is qualitatively different depending on whether the ratio of shear viscosity to entropy density is greater or less than the universal, infinite coupling value of $\hbar/4\pi k_B$. In the former case, the density of poles increases, indicating a formation of branch cuts in the weak coupling limit, and the spectral function shows the appearance of narrow peaks. We also discuss the relation of the viscosity-entropy ratio to conjectured bounds on relaxation time in quantum systems.


Introduction
Nuclear matter produced in heavy ion collisions at RHIC and LHC appears to be well described by relativistic fluid dynamics at the time shortly after the collision, i.e. for t > τ H , where the "hydrodynamization" time τ H is of the order of 1 − 2 fm/c [1][2][3][4][5][6]. The hydrodynamic description fits the available experimental data well provided the shear viscosity -entropy density ratio of the resulting nuclear fluid is low, η/s ∼ /4πk B . An interesting and not fully understood question is how the matter reaches the hydrodynamic stage of its evolution so quickly and which physical mechanisms are responsible for such a rapid thermalization at intermediate values of QCD coupling. The regime of intermediate coupling can in principle be approached from either the weak or the strong coupling side and accordingly, issues related to thermalization have been studied in kinetic theory at weak coupling and in gauge-string duality (holography) at strong coupling. While the kinetic theory approach and the holographic methods are very different, it is clear that in one and the same theory (e.g. in N = 4 supersymmetric SU (N c ) Yang-Mills (SYM) theory at infinite N c ) one should expect an interpolation between strong and weak coupling results for observables describing thermalization, similar to the coupling constant dependence of the shear viscosity -entropy density ratio [7,8] or pressure [9,10]. The goal of this paper is to investigate such a dependence for a number of models where corrections to known holographic results at infinitely strong coupling can be computed by using higher derivative terms in the dual gravity action.
Among relevant observables, we focus on the hierarchy of times characterizing the approach to thermal equilibrium. In simple models of kinetic theory, the appropriate time scales emerge as eigenvalues of the linearized collision operator, with the largest eigenvalue, τ R , essentially (within a specified approximation scheme) setting the time scale for transport phenomena [11][12][13][14] (see Section 2 for details). In particular, for the shear viscosity in the non-relativistic kinetic theory one typically obtains [15] η = τ R n k B T , (1.1) where n is the particle density. The relativistic analogue of Eq. (1.1) is where s is the volume entropy density 1 . In kinetic theory, the relaxation time τ R is simply proportional to the (equilibrium) mean free time for corresponding particles or quasiparticles and thus the internal time scale associated with the kinetic operator acquires a transparent physical meaning. In the regime of validity of Eq. (1.2), the dependence of η/s on e.g. the coupling is the same as the dependence on the coupling of τ R T and thus we expect the ratio η/sτ R T to be (approximately) constant in that regime. Another interesting feature of kinetic theory models is the breakdown of the hydrodynamic description for sufficiently large values of the wave vector q > q c and the appearance of the strongly damped Knudsen modes [16]. We shall see that these phenomena have their counterparts in the regime of strong coupling despite the fact that kinetic theory is not applicable in that regime.
It is believed that the quark-gluon plasma created in heavy ion collisions at energies available at RHIC or LHC is a strongly interacting system, for which a direct or effective (via a suitable quasiparticle picture) application of kinetic theory is difficult to justify. Instead, insights into the time-dependent processes at strong coupling are obtained by studying qualitatively similar strongly coupled theories having a dual holographic description in terms of higher-dimensional semiclassical gravity. Holography [17][18][19][20][21] provides a convenient framework for studying non-equilibrium phenomena in strongly interacting systems. The dynamics and evolution of non-equilibrium states in a strongly interacting quantum many-body system is mapped (in the appropriate limit) into the dynamics and evolution of gravitational and other fields of a dual theory. Holography should in principle be capable of encoding all types of non-equilibrium behavior. In particular, evolution of the system towards thermal equilibrium is expected to be described by the dynamics of gravitational collapse. Numerical and analytical studies of processes involving strong gravitational fields including black holes and neutron stars mergers resulting in black hole formation and particles falling into black holes show a characteristic scenario in which a primary signal (strongly dependent on the initial conditions) is followed by the quasinormal ringdown (dependent on the final state parameters only) and then a late-time tail (see e.g. [22], [23]). A holographic description of fully non-equilibrium quantum field theory states via dual gravity has been developed over the last several years and the results suggest that the quasinormal spectrum (i.e. the eigenvalues of the linearized Einstein's equations of the dual black brane background) and in particular the fundamental (the least damped non-hydrodynamic) quasinormal frequency play a significant role in the description of relaxation phenomena. Recent studies (including sophisticated numerical general relativity approaches) of equilibration processes in the dual gravity models [24][25][26][27][28][28][29][30][31][32][33] reveal that the hydrodynamic stage of evolution is reached by a strongly coupled system long before the pressure gradients become small and that the relevant time scales are essentially determined by the lowest quasinormal frequency, even for non-conformal backgrounds [31,[34][35][36][37][38]. The characteristic time scale here is set by the inverse Hawking temperature of the dual equilibrium black hole.
A seemingly natural question to ask is whether the relation between transport phenomena and the relaxation time(s) familiar from kinetic theory exists also at strong coupling and if yes, how it changes as a function of coupling. Is there a limiting value of the wave vector beyond which hydrodynamic description breaks down at large but finite coupling? Extrapolating kinetic theory results to the regime of intermediate coupling was the subject of recent investigation by Romatschke [39]. In holography, these questions can be studied by computing coupling constant corrections to the full quasinormal spectra using the appropriate higher derivative terms in dual gravity. Recently, such corrections have been studied in Refs. [40], [41].
In this paper, we compute the quasinormal spectra of metric perturbations of the gravitational background with R 4 higher derivative term (dual to N = 4 SYM at finite temperature and large but finite 't Hooft coupling), and for the background with R 2 terms including Gauss-Bonnet gravity in d = 5 dimensions. Normally, higher derivative terms are treated as infinitesimally small corrections to the second order equations of motion of Einstein gravity, otherwise one is doomed to encounter the Ostrogradsky instability and related problems. Accordingly, extrapolating results from infinitesimal to finite values of the corresponding parameters requires caution. Gauss-Bonnet and more generally Lovelock gravity are good laboratories since their equations of motion are of second order and thus can handle finite values of the parameters multiplying higher derivative terms. However, such theories appear to suffer from internal inconsistencies for any finite value of the parameters [42] (for an apparently dissenting view, see [43]). The passage between Scylla and Charybdis of those two difficulties may be hard to find, if it exists at all. We find some solace in the fact that our results show a qualitatively similar picture regardless of the exact form of higher derivative terms used.
The paper is organized as follows: Our main results are summarized in Section 2, where we also review some facts about the relaxation times in quantum critical, kinetic and gravitational systems, adding a number of new observations along the way. In Section 3, we compute the (inverse) 't Hooft coupling corrections to the quasinormal spectrum of gravitational fluctuations in AdS-Schwarzschild black brane background modified by the higher derivative terms and discuss the relaxation time behavior, the density of poles and the inflow of extra poles from infinity. In Sections 4 and 5, correspondingly, a similar procedure is applied to Gauss-Bonnet gravity and to the background with generic curvature squared terms. We briefly discuss the results in the concluding Section 6. Some technical issues and comments about our numerical procedures appear in the Appendices.

Relaxation times at weak and strong coupling
In this Section, we briefly review the appearance of the hierarchy of relaxation times in kinetic theory, holography and some models of condensed matter physics, emphasizing their similarities and adding some new observations. In this context, at the end of the Section, we list the main results of the present paper.
In kinetic theory, transport coefficients and relaxation time(s) are intimately related. To be clear, by the relaxation time we mean the characteristic time interval during which a local thermal equilibrium (e.g. a local Maxwell-Boltzmann equilibrium) is formed everywhere in the system. We are not interested in the momentum-dependent equilibration time-scales of the densities of conserved charges (these densities always relax hydrodynamically) which are, strictly speaking, infinite in the limit of vanishing spatial momentum. Consider, for illustration, non-relativistic Boltzmann equation obeyed by the one-particle distribution function F (t, r, p) where U (r) is the external potential and C[F ] is the Boltzmann collision operator containing details of the interactions. For small deviations from the local thermal equilibrium described by the distribution function F 0 (r, p), the kinetic equation can be linearized by the ansatz where ϕ 1. The ansatz (2.2) leads to the evolution equation where L 0 is a linear integral operator resulting from linearization of C[F ]. Formal solution to Eq. (2.3) with the initial condition ϕ(0, r, p) = ϕ 0 (r, p) can be written in the form [44] ϕ(t, r, p) = e tL ϕ 0 (r, p) = 1 2πi where R s = (sI − L) −1 is the resolvent whose analytical structure in the complex s-plane determines the relaxation properties. In some simple cases, such as e.g. the relaxation of a low-density gas of light particles in a gas of heavy particles, the resolvent can be constructed explicitly and the time dependence fully analyzed [45]. Generically, however, the time evolution is not known explicitly. For spatially homogeneous equilibrium distributions and perturbations, a simple ansatz ϕ(t, p) = e −νt h(p) reduces the linearized kinetic equation to the eigenvalue problem for the linear collision operator: The eigenvalues of L 0 determine the spectrum of (inverse) relaxation times in the system. One can then write a general solution of the linearized kinetic equation in the form where the coefficients C n are determined by the initial conditions and the sum should be replaced by an integral if the spectrum turns out to be continuous. The hierarchy {ν n } in Eq. (2.6) is clearly reminiscent of the hierarchy of imaginary times of the quasinormal modes in the dual gravity treatment of near-equilibrium processes at strong coupling. The spectrum of the operator L 0 for (classical) particles interacting via the potential U (r) = α/r n has been investigated by Wang Chang and Uhlenbeck [11] and by Grad [13]. The spectrum consists of a five-fold degenerate null eigenvalue, corresponding to conserved quantities and the rest of the spectrum which can be discrete (for n = 4) or continuous [11,13,14], with or without a gap, depending on n (see Fig. 1). The time dependence is obviously sensitive to the type of the spectrum: discrete spectrum leads to a clear exponential relaxation, whereas continuous spectrum implies a more complicated pattern including a pure power-law fall-off in the gapless case. Assuming the spectrum is discrete and denoting τ R = 1/ν min , in the relaxation time approximation, when the sum in (2.6) is dominated by a single term with ν n = ν min , we find Generalization to weakly inhomogeneous systems gives [14,46] ∂F ∂t The equation (2.8) has been remarkably successful in describing transport phenomena in systems with a kinetic regime 2 [12], [44]. In particular, assuming τ R = const, for the shear Figure 1: The spectrum of a linear collision operator: a) discrete spectrum, b) continuous spectrum with a gap, realized for the interaction potential U = α/r n , n > 4, c) gapless continuous spectrum, realized for the interaction potential U = α/r n , n < 4, d) Hod spectrum (see text): 0 ≤ ν min ≤ ν c . In all cases, ν = 0 is a degenerate eigenvalue corresponding to hydrodynamic modes (at zero spatial momentum).
viscosity one obtains the result (1.1). Estimates of τ R based on Ritz variational method relate the relaxation time to the mean free time: where σ is the interaction cross-section. The account above may look too schematic but a more detailed treatment is available in the standard kinetic theory [44] (including relativistic and quantum cases [14], [48]), in the mathematical theory of Boltzmann equation [49] and in thermal gauge theory [50,51].
Do the relations between transport coefficients and relaxation time(s) similar or identical to the ones in Eqs. (1.1) and (1.2) hold beyond the regime of applicability of kinetic theory and in the absence of quasiparticles? One may appeal to dimensional analysis and the uncertainty principle [52] or "general wisdoms" [21] when arguing for an affirmative answer 3 but in all cases the concept of weakly interacting quasiparticles seems to be lurking behind such reasoning. At the same time, the concepts of relaxation time and transport are meaningful irrespective of whether or not the kinetic theory arguments are applicable. In particular, in condensed matter physics, considerable attention has been drawn to the studies of quantum critical regions [53], where the characteristic time scales of strongly interacting theories at finite temperature are of the order of τ ∼ /k B T (see Fig. 2). Moreover, estimates of thermal equilibration time τ R in relevant models suggest that [53] where C is a constant of order one, with the inequality saturated in the quantum critical region. In some models, the constant C can be computed analytically. For the quantum 3 Indeed, the characteristic time scale in the kinetic regime is the mean free path τ ∼ t mf p and in the regime of strong coupling it is the inverse temperature of a dual black hole, τ ∼ /kBT . Assuming η/s ∼ τ kBT , we have η/s ∼ √ mT /nσ in the first case and η/s ∼ /kB in the second.

Figure 2:
Phase diagram of the d = 1 + 1 quantum Ising model [53]. The relaxation time in the quantum critical region is determined by the lowest quasinormal frequency of the BTZ black hole.
Ising model in d = 1+1 dimension serving as one of the main examples illustrating quantum critical behavior in [53], the relaxation time of the order parameterσ z having the anomalous dimension ∆ = 1/8 in the quantum critical region is determined by the correlation function of a 1 + 1-dimensional CFT at finite temperature. The (equilibrium) retarded two-point correlation function of an operator of (non-integer) conformal dimension ∆ in momentum space is given by [54] where C ∆ is the normalization constant and we put T L = T R = T . The correlator has a sequence of poles at where n = 0, 1, 2, .... Note that these are precisely the quasinormal frequencies of the dual BTZ black hole [55], [54]. At zero spatial momentum, the lowest quasinormal frequency determines the relaxation time  [53]. Curiously, inserting ∆ = 2 (the scaling dimension of the energy-momentum tensor) into Eq. (2.12) and using Eq. (1.2), one formally 5 finds η/s = 1/4π. In holography, the importance of the quasinormal spectrum as the fundamental characteristic feature of near-equilibrium phenomena in a dual field theory has been recognized early on [56], [57], [58] and later it was observed [55] and shown [54], [59] that the quasinormal frequencies correspond to poles of the dual retarded correlators. A typical distribution of poles in the complex frequency ω plane at fixed spatial momentum q of an equilibrium retarded correlator computed via holography in the supergravity approximation (e.g. at infinite 't Hooft coupling and infinite N c in N = 4 SYM) is shown in Fig. 3 (right panel), where the spectrum of a scalar fluctuation is shown [60]. For correlators of conserved  quantities such as the energy-momentum tensor, the spectrum, in addition to an infinite tower of gapped strongly damped modes ω n = ω n (q), contains also a sector of gapless hydrodynamic modes ω = ω(q) with the property ω(q) → 0 for q → 0 [59,60,62]. Asymp-totics of these spectra were computed in Refs. [63,64] (for large n) and in Ref. [65] (for large q). Curiously, at weak coupling the correlators at finite spatial momentum q seem to have branch cuts stretching from −q to q rather than poles [61] (see left panel of Fig. 3). At zero spatial momentum, the branch cuts reduce to a sequence of poles on the imaginary axis [61]. These issues are further discussed in [39] and in the present paper.
Finite coupling corrections to the quasinormal spectrum can be computed by using higher derivative terms in the appropriate supergravity action. Such corrections for gravitational backgrounds involving R 4 higher derivative term were recently computed in Refs. [40,41]. In this paper, we consider R 4 and R 2 terms, including Gauss-Bonnet gravity. We find a number of novel features in addition to those reported in Refs. [40,41]. Our observations can be summarized as follows (see Sections 3, 4, 5 for full details): • The positions of all poles change with the coupling. In the shear channel in particular, two qualitatively different trends are seen depending on whether η/s > /4πk B or η/s < /4πk B (see Fig. 4). In the first case (realized, for example, in N = 4 SYM), the symmetric branches of non-hydrodynamic poles lift up towards the real axis 6 and the diffusion pole moves deeper down the imaginary axis. In the second case (corresponding to known examples of the dual gravity actions with curvature squared corrections, in particular, Gauss-Bonnet gravity with positive coupling), the branches move up only very slightly and the diffusion pole comes closer to the origin.
• For η/s > /4πk B , the density of poles in the symmetric branches increases monotonically with the coupling changing from strong to weak values as shown schematically in Fig. 4. Qualitatively, this seems to be compatible with the poles merging and eventually forming branch cuts (−∞, q] and [q, ∞), where q is the spatial momentum, in the complex frequency plane at vanishing coupling. For η/s < /4πk B , however, the density of poles decreases and they seem to disappear from the finite complex plane completely in the limit of vanishing viscosity.
• In the holographic models we considered, the function η/s τ R T is a slowly varying function of the coupling, with an appreciable change in the vicinity of infinite coupling only, suggesting that approximations of the type η/s ∼ const τ R T are not unreasonable in the strongly coupled regime even though they cannot possibly follow from kinetic theory arguments.
• In view of the relation between η/s and relaxation time, a bound on quasinormal frequencies of black branes similar to the one proposed by Hod for black holes [66] may imply a bound on η/s. This is further discussed in Section 6.
• As η/s increases well beyond /4πk B and the poles approach the real axis, we expect them to be visible as clear quasiparticle-like excitations (i.e. well-defined, high in amplitude and very narrow peaks) in the appropriate spectral function of the dual field theory, well known from weakly coupled theories. This is indeed the case (see Section 4.2.6 for a calculation of the shear channel spectral function in the Gauss-Bonnet theory where these feature can be seen explicitly).
• An inflow of new poles from complex infinity is observed at finite coupling. The new poles ascend from the negative infinity towards the origin along the imaginary axis as the coupling changes. The behavior of these new poles as a function of coupling also depends on whether η/s > /4πk B or η/s < /4πk B . In the Gauss-Bonnet model with η/s < /4πk B (i.e. with positive values of the Gauss-Bonnet coupling), the poles reach the asymptotic values known analytically [67], without interfering with the hydrodynamic pole. However, in models with η/s > /4πk B (N = 4 SYM or Gauss-Bonnet holographic liquid with negative coupling), a qualitatively different picture is observed. In this case, in the shear channel, the least damped new pole reaches the hydrodynamic pole at a certain value of the coupling (for each fixed q), the two poles merge and then move off the imaginary axis. Furthermore, as the coupling constant varies at fixed q, the poles previously describing the hydrodynamic excitations (diffusion and sound) become the leading (i.e. having the smallest Im|ω|) poles of the two symmetric branches. We interpret these phenomena as the breakdown of the hydrodynamic gradient expansion at some value of the coupling (for each q). Phrased differently, at each value of the coupling λ, there exists a critical value of the wave vector q c (λ) such that for q > q c (λ) the hydrodynamic description becomes inadequate. In the holographic models we considered, the function q c (λ) is a monotonically increasing function of the coupling suggesting that the range of validity of the hydrodynamic description is larger at strong coupling. Details are reported in Sections 3 and 4. This is reminiscent of the weak coupling kinetic theory behavior mentioned earlier [16] and also the one described in [39], although our interpretation is somewhat different from the one in Ref. [39].
The reported observations (admittedly, made only for a few holographic models and suffering from various limitations mentioned above) seem to suggest the following picture: First, the relations such as (1.2) may still hold in the regime of the coupling where the kinetic theory approach used to derive them can no longer be justified. This may explain why using the kinetic theory formally outside its regime of applicability can still give results compatible with experimental data. Second, it seems that for a fixed value of the coupling, there exist critical length-and time-scales beyond which the hydrodynamic approximation fails. The dependence of these critical scales on coupling extracted from the holographic models suggests that hydrodynamics has a wider range of applicability at strong coupling in comparison to weaker coupling. This appears to be compatible with the widely reported "unreasonable effectiveness of hydrodynamics" in models of strongly coupled plasma.

Coupling constant corrections to equilibrium energy-momentum tensor correlators in strongly interacting N = 4 SYM theory
For N = 4 supersymmetric SU (N c ) Yang-Mills (SYM) theory in d = 3 + 1 (flat) dimensions, corrections in inverse powers of the 't Hooft coupling λ = g 2 Y M N c at infinite N c to thermodynamics [9,68] and transport [8,[69][70][71][72][73][74][75] have been computed using the higher derivative R 4 term [76,77] in the effective low-energy type IIB string theory action 7 . In particular, for the shear viscosity to entropy density ratio the coupling constant correction to the universal infinite coupling result is positive [8,69]: The result (3.1) can be found, in particular, by computing the λ −3/2 correction to the hydrodynamic (gapless) quasinormal frequency in the shear channel of gravitational perturbations of the appropriate background. Coupling constant corrections to the full quasinormal spectrum of gravitational perturbations of the AdS-Schwarzschild black brane background, dual to finite-temperature N = 4 SYM were previously computed by Stricker [40] (see also [41]). In this Section, we reproduce those results and find some new features focusing on the relaxation time and the behavior of the old and new poles.

Equations of motion
The source of finite 't Hooft coupling corrections is the ten-dimensional low-energy effective action of type IIB string theory where γ = α 3 ζ(3)/8 and the term W is proportional to the contractions of the four copies of the Weyl tensor Considering corrections to the AdS-Schwarzschild black brane background and its fluctuations, potential α corrections to supergravity fields other than the metric and the five-form field have been argued to be irrelevant [82]. Moreover, as discussed in [83], for the purposes of computing the corrected quasinormal spectrum one can use the Kaluza-Klein reduced five-dimensional action The full set of α 3 terms in the ten-dimensional effective action is currently unknown. Corrections involving the self-dual Ramond-Ramond five-form were considered in Refs. [78][79][80][81]. Following the arguments in [82], in this paper we assume that the (unknown) corrections to fields whose background values vanish to leading order in α 3 for a given supergravity solution will not modify the quasinormal spectrum at order α 3 and thus can be neglected. We thank A. Buchel and K. Skenderis for discussing these issues with us.
where W is given by Eq. (3.3) in 5d. The parameter γ is related to the value of the 't Hooft coupling λ in N = 4 SYM via γ = λ −3/2 ζ(3)L 6 /8 (we set L = 1 in the rest of this Section). Higher derivative terms in the equations of motion are treated as perturbations and thus any reliable results are restricted to small values of the parameter γ. The effective five-dimensional gravitational constant is connected to the rank of the gauge group by the expression κ 5 = 2π/N c .
To leading order in γ, the black brane solution to the equations of motion following from (3.4) is given by [9,68] where f (u) = 1 − u 2 , r 0 is the parameter of non-extremality of the black brane geometry and the functions Z t and Z u are given by The γ-corrected Hawking temperature corresponding to the solution (3.5) is T = r 0 (1 + 15γ)/π. For the isotropic N = 4 SYM medium, we now consider fluctuations of the metric of the form g µν = g µν is the background (3.5). We Fourier transform the fluctuations with respect to t and z to introduce h µν (u, ω, q), choose the radial gauge with h uν = 0 and follow the recipes in [59] to write down the equations of motion for the three gauge-invariant modes Z i = Z (0) i + γZ (1) i , i = 1, 2, 3, in the scalar, shear and sound channels, respectively. Explicitly, the three modes and the corresponding equations of motion are given by the following expressions 8 :

Scalar channel
Shear channel Sound channel The functions G 1 , G 2 and G 3 appearing on the right hand side of the equations can be found in Appendix A. Here and in the rest of the paper we use the dimensionless variables The equations of motion are solved numerically and the quasinormal spectrum is extracted using the standard recipes [8,59,60,62,69,70]. Our numerical approach is described in Appendix C.

The spectrum of the metric fluctuations
Given the smooth dependence of the equations of motion on the parameter γ, we may expect the eigenvalues to shift somewhat in the complex frequency plane with respect to their γ = 0 positions. This is indeed the case, as noted previously in Refs. [40,41] and the details of this shift are interesting. In addition to this, we observe an inflow of new poles from complex infinity along the imaginary axis. These poles are non-perturbative in γ (the relevant quasinormal frequencies scale as 1/γ) but under certain conditions they are visible in the finite complex frequency plane and can even be approximated by analytic expressions. The new poles appear in all three channels of perturbations. In the shear and sound channels, they interfere with the hydrodynamic poles and effectively destroy them at still sufficiently small, q-dependent values of γ. A qualitatively similar phenomenon is observed in Gauss-Bonnet gravity where the equations of motion are second order and fully non-perturbative (see Section 4).

Scalar channel
The scalar equation of motion (3.8) is solved numerically with the incoming wave boundary condition at u = 1 and Dirichlet condition at u = 0 for fixed small values of γ > 0. A typical distribution of the quasinormal frequencies (poles of the scalar components of the energy-momentum retarded two-point function of N = 4 SYM) in the complex frequency plane is shown in Fig. 5.
The two symmetric branches of the modes move up towards the real axis relative to their γ = 0 position. Here and in all subsequent calculations, we do not observe the bending of the quasinormal modes with large real and imaginary parts towards the real axis reported earlier in Ref. [40]. Rather, our findings agree with the results of Ref. [41], where the two branches lift up without bending. The two branches become more and more horizontal with the 't Hooft coupling decreasing and move closer to the real axis.
At the same time, the distance between the poles in the branches decreases: in a sense, there is an inflow of new poles from complex infinity along the branches. This last effect is too small to be noticeable e.g. in Fig. 5 because in N = 4 SYM we are restricted to the γ 1 regime. Extrapolating to larger values of γ (smaller values of 't Hooft coupling) would not be legitimate with the R 4 corrections treated perturbatively but it is conceivable that in the limit of vanishing 't Hooft coupling the poles in the two branches merge forming two symmetric branch cuts (−∞, −q] and [q, ∞). We shall see more evidence for this behavior in Gauss-Bonnet gravity, where the equations of motion are second-order and the coupling dependence is fully non-perturbative (see Section 4). The closeness of the two branches of poles to the real axis at intermediate and small values of the 't Hooft coupling raises the question of the behavior of the spectral function and the appearance of quasiparticles. Again, this is investigated in detail in Gauss-Bonnet gravity in Section 4, where we are not constrained by the smallness of the perturbation theory parameter.
We also observe a novel phenomenon: a sequence of new poles ascends along the imaginary axis towards the origin as γ increases from zero to small finite values. The first of these poles reaches the vicinity of the origin at γ ∼ 0.01. One can find a crude analytic approximation for this top pole by solving the equation in the regime |w| 1 (for simplicity, we also take |q| 1). We assume the scaling w → w and q → q for 1, so that to first order in , the function are found perturbatively to first order in γ. To find the quasinormal frequency, we solve the polynomial equation Z 1 (u = 0, w, q) = 0, looking for a solution of the form w(q). To leading order in q, we find a gapped pole on the imaginary axis with the dispersion relation . (3.14) As shown in Fig. 6, the analytic approximation (3.14) works better for larger values of γ.
For γ → 0, the pole recedes deep into the complex plane along the negative imaginary axis (the approximate formula (3.14) is compatible with this observation but breaks down when |w| becomes large).

Shear channel
In the shear channel, the distribution of poles at finite coupling is similar to the one in the scalar channel. The exception is the gapless hydrodynamic pole on the imaginary axis responsible for the momentum diffusion. The poles are shown in Fig. 7 for several values of γ. General properties of non-hydrodynamic poles described in detail for the scalar channel are observed here as well. The new feature is the interaction between the diffusion pole and the first of the new poles rising up from complex infinity along the imaginary axis with increasing γ. The dispersion relation for the diffusion pole is given by the formula [84][85][86] where in the absence of the chemical potential ε + P = sT . In N = 4 SYM theory, one has [8,69,70,72,72,73,[85][86][87][88] The coupling constant correction to the coefficient θ 1 of the third-order hydrodynamics is currently unknown. However, for q 1, the q 2 term in Eq. (3.15) dominates and the pole moves down the imaginary axis with γ increasing, in agreement with our numerical findings.
For certain values of γ 1, the leading new pole ascending the imaginary axis approaches the hydrodynamic pole. The two poles collide on the imaginary axis at some critical value of γ at fixed q (or equivalently, at some q = q c (γ) at fixed γ) and then for larger γ they symmetrically move off the imaginary axis, both having acquired non-zero real parts (see Fig. 8). At this point, the hydrodynamic pole (3.15) ceases to exists and for q > q c the hydrodynamic description appears to be invalid. We interpret this as the breakdown of hydrodynamics at sufficiently large, coupling-dependent value of the wave-vector. The function q c (γ) is shown in Fig. 9. It is monotonically decreasing with γ suggesting that hydrodynamics has a wider applicability range at larger 't Hooft coupling as far as the spatial momentum dependence is concerned.
The phenomenon just described can be approximated analytically in the region of small w and q (although this approximation is not very precise quantitatively, it captures the behavior of the poles correctly). Indeed, solving the equation (3.10) perturbatively in w 1 and q 1 (still with γ 1) and imposing the Dirichlet condition Z 2 (u = 0, w, q) = 0, we find a quadratic equation This equation has two roots parametrized by γ and q, At fixed q and sufficiently small γ, the roots are purely imaginary, moving closer to each other with increasing γ. Finally, the two roots merge and then acquire non-zero real parts for larger γ. The physical meaning of the solutions (3.20) and (3.21) becomes transparent from their small q expansions: where w g is given by Eq. (3.14). Here, the mode (3.22) is the standard hydrodynamic momentum diffusion pole predicted by Eq. (3.15), whereas the mode (3.23) approximates the new gapped pole moving up the imaginary axis. Note that the gap w g in the mode (3.23) is the same as in the scalar channel. Using Eqs. (3.22) and (3.23), we can find an approximate analytic expression for the function q c (γ) plotted in Fig. 9: As is evident from Fig. 9, the analytic approximation becomes more precise with larger γ.

Sound channel
The quasinormal spectrum in the sound channel is found by solving Eq. (3.12) and imposing the Dirichlet condition Z 3 (u = 0, w, q) = 0. The distribution of poles in the complex frequency plane at various values of the coupling is shown in Fig. 10. The movement of the poles with varying coupling is qualitatively similar to the one observed in the scalar and shear channels. The two gapless sound poles symmetric with respect to the imaginary axis have the dispersion relation predicted by hydrodynamics [84][85][86] where c s = 1/ √ 3 for conformal fluids in d = 3 + 1 dimensions, Γ = 2η/3(ε + P ) and ε + P = sT at zero chemical potential. For N = 4 SYM theory, the coefficients η/s, τ Π and θ 1 are given in Eq. (3.18) and The full γ-dependence of the quartic term in Eq. (3.25) is currently unknown. Figure 9: Critical value of the spatial momentum q c , limiting the hydrodynamic regime, as a function of higher derivative coupling γ in the shear channel of N = 4 SYM. Hydrodynamics has a wider range of applicability in q at smaller γ (larger 't Hooft coupling).
With γ increasing, the leading new gapless pole rising along the imaginary axis approaches the region of the sound poles (see Fig. 11). For w 1 and q 1, the equation (3.12) can be solved perturbatively and from the Dirichlet condition Z 3 (u = 0, w, q) = 0 one finds a quintic equation (3.27) Expanding further in γ 1 and q 1, we obtain the following analytic expressions for the three closely located modes of interest: Here, the Eq.

Coupling constant dependence of the shear viscosity -relaxation time ratio
The dependence of real and imaginary parts of the smallest in magnitude quasinormal frequencies in the symmetric branches on γ (at fixed q) in the scalar, shear and sound In all three channels, a relatively strong dependence of the spectrum on γ in the vicinity of γ = 0 changes to a nearly flat behavior at larger values of γ. As discussed in the Introduction, these data can be used to test whether the relations between transport coefficients and the relaxation time typical for a kinetic regime of the theory may still hold at strong coupling. In kinetic theory, the hierarchy of relaxation times arises as the non-hydrodynamic part of the spectrum of a linearized Boltzmann operator (see Section 2). At strong coupling, it seems natural to associate this hierarchy with the (inverse) imaginary parts of the quasinormal spectrum frequencies. In particular, the relaxation time τ R can be defined as where ω F is the fundamental (lowest in magnitude) quasinormal frequency. The prediction of kinetic theory is that Eq. (1.2) holds at least at weak coupling, i.e. that the ratio η/s τ R T is approximately constant. In Fig. 15, we plot the ratios η/s τ k T , k = 1, 2, 3, 4, as functions of γ using the data for τ k = 1/Im w k of the leading four non-hydrodynamic quasinormal frequencies (including the fundamental one) in the shear channel at q = 0. Curiously, although rapid decrease of all four functions is seen in the vicinity of γ = 0, the dependence changes to a nearly flat one very quickly, already at γ ≈ 2 × 10 −3 (corresponding to the 't Hooft coupling λ ∼ 18), which is well within the regime of small γ. Note that the 't Hooft coupling correction to η/s for γ ≈ 2 × 10 −3 is approximately 25%. Thus, the naive use of kinetic theory expressions such as Eq. (1.2) may not be so disastrous at moderate or even strong coupling. We shall see in the next Sections that the features discussed here for the specific gravity dual with the higher derivative term of the type R 4 are also observed for gravity backgrounds with R 2 terms, in particular Gauss-Bonnet gravity.    Figure 15: The ratios η/sτ k T , k = {1, 2, 3, 4}, as functions of γ in N = 4 SYM.

Relaxation time and poles of energy-momentum tensor correlators in a theory dual to Gauss-Bonnet gravity
The action of Einstein-Gauss-Bonnet gravity in five space-time dimensions is given by where the scale l 2 GB of the higher derivative term can be chosen to be set by a cosmological constant, l 2 GB = λ GB L 2 , where λ GB is the dimensionless parameter. The coefficients of the curvature-squared terms ensure that the equations of motion following from the action (4.1) are second order in derivatives. Thus, in the absence of Ostrogradsky instability and other difficulties usually induced by the dynamics with higher derivatives, Gauss-Bonnet and more generally Lovelock theories, are popular theoretical laboratories for studying non-perturbative effects of higher-derivative couplings. For example, the shear viscosityentropy ratio in a (hypothetical) conformal fluid dual to five-dimensional Gauss-Bonnet gravity turns out to be [89] η and this result is obtained without the assumption |λ GB | 1, i.e. non-perturbatively in the coupling. However, as pointed out and investigated in detail in Refs. [89][90][91][92][93][94], for λ GB outside of a certain interval, the dual theory suffers from pathologies associated with superluminal propagation of high momentum modes. More recently, Camanho et al. [42] argued that Gauss-Bonnet theory suffers from causality problems in the bulk that can only be cured by adding higher spin fields. This would effectively imply that Gauss-Bonnet and, most likely, general Lovelock theories 9 should loose their privileged non-perturbative status and be treated as any other theory with higher derivative terms, i.e. the coupling λ GB in, for example, Eq. (4.2) should be seen as an infinitesimally small parameter (see, however, Refs. [43,98]). We note those difficulties but will not constrain λ GB beyond its natural (here, limited by the existence of the black brane solution) domain λ GB ∈ (−∞, 1/4] in the following.
Our goal is to compute the quasinormal spectrum of gravitational fluctuations of the Gauss-Bonnet black brane metric 10 dual to a thermal state of a boundary CFT. Here 9 See Refs. [95][96][97] for relevant work in Lovelock theories. 10 Exact solutions and thermodynamics of black branes and black holes in Gauss-Bonnet gravity were considered in [99] (see also [100][101][102][103]). and the constant N GB can be chosen to normalize the speed of light at the boundary to c = 1: The Hawking temperature corresponding to the solution (4.3) is given by where we introduced the notation γ GB ≡ √ 1 − 4λ GB . We shall use λ GB and γ GB interchangeably in the following. The range λ GB < 0 corresponds to γ GB ∈ [1, ∞) and the interval λ GB ∈ [0, 1/4] maps into γ GB ∈ [0, 1], with λ GB = 0 corresponding to γ GB = 1.

Equations of motion
Fluctuations h µν (r, t, z) of the Gauss-Bonnet black brane metric (4.3) can be decomposed into the scalar, shear and sound channels in the standard way [54,59]. The corresponding gauge-invariant combinations Z 1 , Z 2 , Z 3 of the metric fluctuations h µν (r, ω, q) (Fourier transformed in the variables along the brane directions) in the three channels are given by

7)
Shear: 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 B. To find the quasinormal spectrum in the three channels, we impose the "incoming wave" boundary conditions at the horizon at u = 1 [54], where the functions Z i are regular at u = 1. The quasinormal spectra w = w(q) are then solutions to the equations Z i (u = 0, w, q) = 0. They can be found numerically. In addition, in the regime w 1 and q 1, some frequencies are determined analytically. In all three channels, it will be convenient to use a new variable 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) and the results for that point, which are identical to those of N = 4 SYM theory at infinite 't Hooft coupling, have to be obtained independently.

The spectrum of the metric fluctuations
The quasinormal spectra in Einstein-Gauss-Bonnet theory obtained non-perturbatively in λ GB show the properties qualitatively similar to the ones discussed in Section 3.2 for the AdS-Schwarzschild background corrected by the R 4 term. In this Section we show the numerical results and analytic approximations for the spectra in the three channels, including the details of the breakdown of the hydrodynamic regime. There are some novelties in the Gauss-Bonnet case. First, not being restricted by the perturbative nature of the higherderivative coupling, we are able to explore the coupling dependence to a fuller extent than in N = 4 SYM. In particular, we are able to say more about the spectral function and the density of poles in the complex plane than we could in N = 4 SYM owing to the restriction γ 1. Second, Gauss-Bonnet gravity (and gravity with generic R 2 terms) provides an example of a holographic model, where the shear viscosity -entropy density ratio can be greater or less than 1/4π, depending on the sign of λ GB . We find qualitatively different patterns in the behavior of relaxation time and other quantities in those two regimes.

Scalar channel
The spectrum of gravitational perturbations in the scalar channel is shown in Fig. 16. Two different regimes are observed depending on the value of η/s.
For η/s > 1/4π (corresponding to λ GB < 0), the behavior of the poles is qualitatively the same as in N = 4 SYM: the two symmetric branches of gapped poles lift up towards the real axis monotonically with |λ GB | increasing, the distance between the poles decreases suggesting a formation of branch cuts (−∞, −q] and [q, ∞) in the limit |λ GB | → ∞. Observing the motion of individual poles in the symmetric branches, one can say that there is an inflow of new poles from complex infinity along the branches with |λ GB | increasing. The dependence of real and imaginary parts of the top three poles in the symmetric branches on λ GB at q = 0.5 is shown in Fig. 17. Within the limits of numerical accuracy, this dependence is monotonic. One may notice that the functions become flat for large negative λ GB .
When the poles in the two branches are sufficiently close to the real axis, we expect the spectral function to show the distinct quasiparticle peaks. We shall discuss this in detail for the shear channel, see subsection 4.2.6. There is a new pole rising up the imaginary axis from complex infinity 11 . The position of the new pole in the regime w 1, q 1, γ GB 1 can be determined analytically by solving the equation for Z 1 perturbatively and imposing the condition Z 1 (u = 0, w, q) = 0: The mode remains purely on the negative imaginary axis and approaches the origin as . This result is confirmed numerically. The residue vanishes in the limit γ GB → ∞, and so the pole disappears in that limit [67]. For η/s < 1/4π (corresponding to λ GB > 0), the poles in the two branches become more sparse relative to their λ GB = 0 distribution (see Figs. 16 and 18). In sharp contrast with the η/s > 1/4π case, here the branches lift up very slightly, almost infinitesimally, relative to their λ GB = 0 positions. As shown in Figs. 16 and 18, an outflow of poles along the branches to complex infinity is observed and it is conceivable that the poles of the two branches are eventually completely pushed out of the finite complex plane. At the same time, there are still new poles rising up the imaginary axis. In the limit λ GB → 1/4 (γ GB → 0) they are seen numerically to approach the positions (known exactly [67]) where n 1 and n 2 are non-negative integers. The limit of vanishing shear viscosity λ GB → 1/4 is difficult to explore numerically. However, the observed behavior is consistent with analytic results available at λ GB = 1/4. Indeed, exactly at λ GB = 1/4 the equations of motion can be solved in terms of hypergeometric functions and the quasinormal spectrum is determined exactly [67]. The only quasinormal frequencies at λ GB = 1/4 are the ones given by Eq. (4.14). This is consistent with the picture we observe numerically for 0 < λ GB < 1/4.

Shear channel
The distribution of the poles in the shear channel is shown in Fig. 20 and the coupling dependence of the real and imaginary parts of the top three poles in the symmetric branches can be seen in Fig. 21. The behavior of the poles in the symmetric branches is qualitatively similar to the one observed in the scalar channel. In the limit λ GB → 1/4, the new poles Figure 21: Real (left panel) and imaginary (right panel) parts of the top three quasinormal frequencies in the symmetric branches in the shear channel of Gauss-Bonnet at q = 0.5.
moving up the imaginary axis approach the q-independent positions known analytically [67], where n 1 and n 2 are non-negative integers. Figure  currently known analytically to quartic order in q [86] and is given by where the transport coefficients were defined in Section 3. For Gauss-Bonnet gravity, solving the equation for the shear mode analytically, perturbatively in w 1, q 1 and nonperturbatively in γ GB and imposing the Dirichlet condition Z 2 (0) = 0, we find (4.17) The full set of non-perturbative first-and second-order hydrodynamic transport coefficients in Gauss-Bonnet theory was computed in [104]. The coefficients relevant for the dispersion relation (4.16) are given by Thus, the value of the third-order coefficient θ 1 in the Gauss-Bonnet theory can now be read off Eq. (4.17): (4.20) In the limit of λ GB → 0 (γ GB → 1), this reproduces the corresponding result for N = 4 SYM theory found in Ref. [86], (4.21) The behavior of the momentum diffusion pole depends on whether η/s is greater or less than 1/4π. For η/s < 1/4π (0 < λ GB < 1/4), the pole moves up the imaginary axis relative to its λ GB = 0 position and approaches the origin. It completely disappears from the spectrum at λ GB = 1/4 [67]. For η/s > 1/4π (−∞ < λ GB < 0), its behavior is qualitatively similar to the one observed in N = 4 SYM: it moves down the imaginary axis and collides with the top new pole moving up the axis from complex infinity at which point the hydrodynamic description seemingly fails. Then the two poles move off the imaginary axis into the complex plane. For sufficiently large values of |λ GB |, this phenomenon happens in the range of small w, q and thus can be approximated analytically (e.g. for λ GB ∼ −3, the merger of the poles occurs at |w| ∼ 0.1, q ∼ 0.1). Solving the shear mode equation of motion perturbatively in w 1, q 1, we find a pair of quasinormal frequencies (4.22) (4.23) whose motion in the complex plane approximates the numerical observations quite well (see Fig. 22). Expanding the above expressions for w 1 and w 2 to second order in q, we find the standard hydrodynamic pole of Eq. (4.17) where the gap w GB g is identical to the one in Eq. (4.13). The behavior of the poles is qualitatively the same as in the N = 4 SYM theory. The diffusion pole moves down the imaginary axis while the new gapped pole moves up as λ GB decreases from 0 towards negative values. Then the two poles collide at some q-dependent value of λ c GB = λ c GB (q) and move off the axis. An analytical approximation for this dependence (or more conveniently, for q c = q c (γ GB )) can be found from the condition w 1 (q c ) = w 2 (q c ). We interpret this condition as the condition indicating inadequacy of hydrodynamic description for q > q c (λ GB ). Equating the expressions (4.24) and (4.24), we obtain Note that if we used instead the un-expanded Eqs. (4.22), (4.23), we would find q The discrepancy is due to the additional q corrections not captured by Eqs. (4.24), (4.24). The dependence q c = q c (λ GB ) obtained numerically as well as the analytic approximation (4.26) are shown in Fig. 23.

Sound channel
The poles in the sound channel are shown in Fig. 24 and the behavior of the real and imaginary parts of the three leading non-hydrodynamic poles in the symmetric branches is demonstrated in Fig. 25. We observe the same features of the coupling dependence of the spectrum as in the other channels. The two symmetric branches lift up from their λ GB = 0 positions, moving swiftly towards the real axis and becoming more dense in the case of η/s > 1/4π (λ GB < 0) and moving only slightly, becoming more sparse and apparently disappearing from the finite complex plane for η/s < 1/4π (0 < λ GB < 1/4). There are new ○ gapped poles rising up the imaginary axis regardless of the sign of λ GB . For η/s < 1/4π (0 < λ GB < 1/4), they reach the asymptotic values where n 1 and n 2 are non-negative integers, in the zero viscosity limit λ GB → 1/4. Here the modes (4.27) are the exact quasinormal frequencies at λ GB = 1/4 [67]. In the regime η/s > 1/4π (λ GB < 0), the top new gapped pole moving up the imaginary axis with |λ GB | increasing gradually approaches the level of the two symmetric sound mode poles and becomes aligned with them (see Fig. 26). For larger values of |λ GB |, all three poles move closer to the real axis, with the sound poles now becoming parts of the symmetric branches.
When the three poles are close to the origin, one can try to build an analytic approximation Figure 25: Real (left panel) and imaginary (right panel) parts of the top three quasinormal frequencies in the symmetric branches in the sound channel of Gauss-Bonnet at q = 0.5. Figure  by solving the equation for Z 3 perturbatively in w 1, q 1. The Dirichlet condition (Z 3 (0) = 0) then gives the equation (4.28) The three roots, w 1,2,3 , can be found analytically, but the expressions are too cumbersome to present here. Their expansions in q to quadratic order are given by where the gap w GB g is the same as in the scalar and shear channels (Eqs. (4.13) and (4.25), respectively). The poles (4.29) correspond to the sound wave modes.
Defining the critical momentum q = q c (γ GB ) as the one at which the hydrodynamic expansion no longer serves as an adequate description of the low-energy limit of the theory, we may choose the equation Im[w 1 (q c )] = Im[w 2 (q c )] = Im[w 3 (q c )] to represent such a condition. Solving this for q c (γ GB ), we find exactly the same function (4.26) as in the shear channel. Note, however, that the agreement between our numerical results and the analytic approximation is less satisfactory than in the shear channel (see Fig. 26), apparently due to a stronger q dependence in the sound channel.

The density of poles and the appearance of branch cuts
In Section 3 we observed that for N = 4 SYM theory correlators, the density of poles in the two symmetric branches increases with 't Hooft coupling decreasing. In Gauss-Bonnet theory, the same phenomenon can be investigated in more detail since we are not constrained by infinitesimally small values of the higher derivative coupling. In all channels, the density of non-hydrodynamic poles in the two branches monotonically increases for η/s > /4πk B and decreases for η/s < /4πk B . Although the trend is apparent already from Figs. 16, 20, 24, here we show the density of poles as a function of the coupling constant in Fig. 27. The density is determined by selecting a region of the complex w plane, counting the number of poles in the symmetric branches inside that region and computing the resulting number density. The dependence in Fig. 27 is monotonic within the bounds of our numerical accuracy. The situation for λ GB > 0 (η/s < /4πk B ) is clear: as λ GB → 1/4 (η/s → 0), the poles in the symmetric branches become less and less dense and in the limit they disappear from the finite complex plane altogether, as confirmed by analytic calculation at λ GB = 1/4. For λ GB < 0 (η/s > /4πk B ), the poles become more and more dense, the symmetric branches lift up toward the real axis with |λ GB | increasing, and one may conjecture that in the limit λ GB → −∞ they merge to form branch cuts in the complex plane of frequency along (−∞, −q] and [q, ∞). Numerically, we observe that Re[w] of the leading quasinormal mode in the (right) branch of poles monotonically approaches the line w = q for large |λ GB | (see Fig. 28) which supports the conjecture that ±q are the branch points of the correlator in the limit λ GB → −∞. Note that all other poles (the ones not belonging to the symmetric branches at finite λ GB ) in all channels either join the branches (in the sound and shear channels) or disappear due to vanishing residues (scalar and sound channels) in the limit λ GB → −∞. Thus, in that limit, the analytic structure of the correlator is represented by the branch cuts (−∞, −q] ∪ [q, ∞) (see Fig. 29). This resembles the zero temperature limit of the thermal correlator in a CFT dual to Einstein gravity with no higher derivative corrections. Such correlators are known analytically only for BTZ background. For example, the ∆ = 2 thermal correlator has the form [54] where ψ(z) is the logarithmic derivative of the Gamma-function with poles at z = −n, n = 0, 1, 2, .... In the zero-temperature limit w 1, q 1, the poles merge forming two branch cuts running from the branch points ω = ±q to infinity parallel to the imaginary axis. For large z, Binet's formula implies ψ(z) ∼ ln z and thus in the limit of zero temperature the correlator (4.31) becomes G R ∼ k 2 ln k 2 , where k µ = (−ω, q). Similarly, the zerotemperature limit of the energy-momentum correlator in a 4d CFT dual to Einstein gravity is G R ∼ (−ω 2 + q 2 ) 2 ln (−ω 2 + q 2 ). This function has branch points at ω = ±q, ∞ joined by the branch cuts (−∞, −q] ∪ [q, ∞).

Figure 29:
The conjectured analytic structure of thermal correlators in holographic Gauss-Bonnet theory in the limit λ GB → −∞.

Coupling constant dependence of the shear viscosity -relaxation time ratio in Gauss-Bonnet theory
The coupling constant dependence of the ratio η/s τ R T in Gauss-Bonnet theory shows the same qualitative features as in N = 4 SYM discussed in Section 3.2.4. In Fig. 30, we plot the ratios η/(s τ k T ), where τ k , k = 1, 2, are defined as τ k = 1/|Im ω k | for the two smallest in magnitude non-hydrodynamical quasinormal frequencies ω k at q = 0. We identify τ R with τ 1 , ω 1 being the fundamental frequency. The functions are monotonic, changing rapidly in the vicinity of λ GB = 0 and flattening out in the region |λ GB | ≈ 3 − 6. As in N = 4 SYM theory, the kinetic theory result (1.2) seems to hold at intermediate coupling. Figure 30: The ratios η/sT τ k , for k = {1, 2}, as functions of λ GB in the shear channel of the Gauss-Bonnet theory. Here, τ k are defined as τ k = 1/|Im ω k | for the two smallestin-magnitude non-hydrodynamical quasinormal frequencies ω k at q = 0. The error bars correspond to resolution errors in the w-plane.

Shear channel spectral function and quasiparticles at "weak coupling"
Since the non-hydrodynamic poles in the symmetric branches approach the real axis with |λ GB | increasing (i.e. at weaker coupling), one may expect the corresponding spectral function to develop a structure resembling quasiparticle peaks. We check this by computing the spectral function in the shear channel. Choosing the spatial momentum along the z axis, the shear channel retarded energy-momentum tensor correlator G R xz,xz (w, q, λ GB ) in Gauss-Bonnet theory can be computed as follows [67]: where the function C is given by and Z 2 (u) is the solution of the shear channel equation of motion obeying the incoming wave boundary condition at the horizon and normalized to one at the same u = ε → 0. The solution Z 2 (u) can be written as Z 2 (u) = A 2 Z I 2 (u) + B 2 Z II 2 (u), where Z I 2 (u) and Z II 2 (u) are the two local Frobenius expansions at the boundary (see e.g. [105]). In terms of A 2 and B 2 , the retarded Green's function (4.32) is given by 12 The spectral function is then computed as In Fig. 31, we plot the dimensionless spectral function where κ 2 5 is the Newton's constant from the Gauss-Bonnet action (4.1) and T the Hawking temperature (4.6). As |λ GB | increases and the symmetric branches of poles approach the real w axis, the appearance of quasiparticle-like peaks in the spectral function is clearly seen. As a result of the quasinormal modes now having |Im[w]| |Re[w]| at large |λ GB |, the peaks become sharp and very narrow. Since the density of poles increases with |λ GB |, the density of peaks increases as well. In the limit |λ GB | → ∞, they presumably form a continuum.

Generic curvature squared corrections to quasinormal spectra of metric perturbations
In this Section, we comment on the quasinormal spectrum in the theory with general curvature squared terms in the action, multiplying ln u, it is numerically more convenient to compute B2 by subtracting off the logarithmic term, as was done in [105]. We where the cosmological constant is Λ = −6/L 2 . For the special choice of the parameters α 1 , α 2 , α 3 given by the action (5.1) coincides with the Gauss-Bonnet action (4.1). Generically, however, the action (5.1) leads to the equations of motion involving derivatives up to the fourth order. In this case the higher derivative terms in (5.1) are treated perturbatively and the parameters α i are assumed to be infinitesimally small. We can find the corresponding quasinormal spectra by using a field redefinition and the known results for Gauss-Bonnet and N = 4 SYM theories. One may notice [89] that the action (5.1) with α 3 = 0 is equivalent via a field redefinition and an additional rescaling to the Einstein-Hilbert action with the same cosmological constant and modified Newton's constant (which does not enter the vacuum equations of motion). 13 Consider now a gauge-invariant (with respect to infinitesimal metric perturbations) mode Z (h µν ) that is linear in metric perturbations, δg µν = h µν . To linear order, the Ricci and Einstein tensors are invariant under diffeomorphisms, hence so is g µν R. It therefore follows that g µν andḡ µν transform identically under the diffeomorphisms and so Hence, when α 3 = 0, the quasinormal modes of Z h µν are also those of Z (h µν ), which means that the quasinormal mode spectrum of the AdS-Schwarzschild black brane (dual to thermal N = 4 SYM theory at infinite 't Hooft coupling) is exactly the spectrum of the theory defined by (5.1) with α 3 = 0.
To include the α 3 contributions, we can use the fact that the perturbative (in α i ) quasinormal spectrum generically has the form where ω * 0 are the quasinormal frequencies of the AdS-Schwarzschild black brane. Moreover, the above discussion shows thatω * 1 = 0,ω * 2 = 0. Keeping in mind the identification (5.2) and considering the linearized quasinormal spectrum in Gauss-Bonnet theory, we conclude that λ GBω * 3 /2 = λ GBω * GB . Hence, the quasinormal spectrum of a background defined by the action (5.1) has the form where ω * 0 is the corresponding frequency in the spectrum of AdS-Schwarzschild black brane with no higher derivative corrections included andω * GB is the coefficient of the term linear in λ GB in the corresponding spectrum of the Gauss-Bonnet theory. Thus, the coupling dependence of the relaxation time and other properties of the spectrum described in previous sections are qualitatively the same as the ones in the Gauss-Bonnet theory (and N = 4 SYM theory with large but finite 't Hooft coupling). In particular, one observes a qualitative difference between the regimes with η/s > /4πk B and η/s < /4πk B similar to the one described in the previous Section.

Discussion
In this paper, we have studied the influence of higher derivative R 2 and R 4 terms on the quasinormal spectra of gravitational perturbations of black branes. In a dual QFT, this corresponds to changing the 't Hooft coupling or its analogue from infinite to large but finite value. Understanding, even qualitatively, how the physical quantities responsible for thermalization change from strong to weak coupling would be important both from a conceptual and a phenomenological point of view. We were looking for robust, model-independent qualitative features the higher derivative terms may bring about. Vulnerabilities of this approach are quite obvious. While N = 4 SU (N c ) SYM is a well defined unitary theory, higher derivative corrections in its dual gravity description are only partially known even to leading order in γ ∼ λ −3/2 at infinite N c and those terms must be treated perturbatively. Moreover, as emphasized recently in [41], different physical quantities may have very different sensitivity to coupling corrections, and the smallness of the perturbative parameter γ may not necessarily be a good indicator of the size of corrections. In contrast, the second order equations of motion of Gauss-Bonnet gravity can be treated fully non-perturbatively. However, the (hypothetical) dual field theory suffers from causality violation and even the bulk theory may need higher spin fields to mend the problems (the latter would imply that higher derivative corrections can only be treated perturbatively, i.e. the theory loses its special status with respect to Ostrogradsky instability). Unphazed by these uncertainties, we proceed to investigate coupling corrections in both theories and are encouraged to observe qualitatively similar results in both cases. Our findings are summarized at the end of Section 1.
One curious feature we find is the behavior of quasinormal spectrum leading to a breakdown of the hydrodynamic description at a coupling-dependent critical value q c of the spatial momentum. In both N = 4 SYM and Gauss-Bonnet theories, the dependence on coupling implies that hydrodynamics has a wider applicability range at strong coupling. It may be interesting to investigate the convergence properties of the hydrodynamic derivative expansion at finite coupling, possibly along the lines of Refs. [106,107].
Another qualitatively similar feature for both theories is the coupling dependence of the ratio of the shear viscosity to the product of relaxation time, entropy density and temperature. This quantity is (approximately) constant in kinetic theory at weak coupling. From the dual gravity with higher derivative corrections we find that this ratio changes rapidly in the vicinity of infinite coupling and then shows a very weak (essentially flat) dependence on coupling when the coupling is further decreased to large but finite values. Similar behavior is expected for other transport coefficients. Admittedly, corrections from the unknown higher derivative terms may influence the dependence at intermediate coupling. Yet, if correct, the observed tendency may help to explain certain phenomenological success and "unreasonable effectiveness" of kinetic theory methods far beyond their justified domain.
We also found that the behavior of coupling corrections to quasinormal spectrum and related quantities depends strongly on whether η/s > /4πk B or η/s < /4πk B . In the regime of η/s < /4πk B , the symmetric branches of quasinormal modes exhibit monotonically increasing |Im ω|. Since this could lead to the relaxation time of the system τ R decreasing below any possible lower bound (see Eq. (2.9)), it is conceivable that this regime is pathological. Earlier work was focused on looking for possible pathologies (e.g. causality violation) in the ultraviolet sector of the theories having the regime η/s < /4πk B , and constraining higher derivative couplings accordingly. However, inconsistencies in this regime may exist in the infrared sector as well. As the qualitative behavior of the spectra critically depends on the sign of the correction to η/s = /4πk B , we note that the relation between η/s and the relaxation time τ R raises the possibility that the bound on η/s speculated upon 14 in Ref. [52] is related to a bound on relaxation time. In the holographic models considered in this paper, both η/s and τ R T are monotonic functions of the coupling. For η/s decreasing below /4πk B , the relaxation time also decreases below its value at infinite coupling. Is there a minimal relaxation time possibly correlated with the viscosity bound? Are there any universal constraints on the constant C in Eq. (2.9)? Curiously, in 2006 Hod [66] suggested a universal bound on relaxation time in any system: For the black hole quasinormal spectrum, the inequality (6.1) means that there exists at least one quasinormal frequency whose imaginary part lies in the strip 0 > Im ω ≥ −πk B T / in the complex frequency plane or, in terms of w = ω/2πk B T , in the strip In the language of the kinetic theory linear collision operator spectrum, the bound implies see Fig. 1d. Apparently, the inequality (6.2) holds for black holes (for black holes, the bound suggests that a black hole has (at least) one channel of slowly decaying perturbation modes which respect (6.2)). At first glance, however, the relaxation time bound is void of meaning since one expects the hydrodynamic modes to be always present in any system in the thermodynamic limit and they may relax arbitrarily slowly for sufficiently long wavelengths (in other words, gapless quasinormal frequencies corresponding to hydrodynamic modes are 14 A number of strongly interacting many-body systems -from quark-gluon plasma and cold atoms [3,108,109] to dusty plasmas [110,111] and rare gases and molecules in the vicinity of the critical point [112] have η/s /4πkB. always present in the strip (6.2) for sufficiently small spatial momentum q). Moreover, even if we regard the bound (6.1) as the bound obeyed by the (non-hydrodynamic) relaxation time τ R (and correspondingly, the inequality (6.2) as the one for the fundamental nonhydrodynamic quasinormal frequency), it appears to be violated in all black brane channels (see e.g. tables of quasinormal frequencies in Refs. [59,60,62]). Nevertheless, we believe the question of whether holography or black hole physics implies an inequality of the type (6.3) or (6.1) is an interesting one in view of its apparent validity for black holes and its possible connection to viscosity bound.

C Numerical methods used
In this Appendix, we briefly review the numerical methods used to compute quasinormal modes and spectral functions. The relevant equations of motion are (3.8), (3.10), (3.12) for the N = 4 and (4.10) for the Gauss-Bonnet case. In all approaches, the first step involves computing the index of the corresponding regular singular point at the horizon, i.e. factoring out the singular part of the solution: where ν ± = ± i 2 w in all channels for N = 4 SYM and Gauss-Bonnet theories, andZ i (u) is the Frobenius solution regular at the horizon. A similar factorization is done at the boundary.
Leaver's method In this method, originally introduced in [113] and then used in [114][115][116][117], the regular function is represented by a series whose convergence radius reaches both the horizon and the boundary,Z = Nmax n=0 a n (w, q) (z − z 0 ) n . (C.1) In principle, z 0 could be chosen arbitrarily, however, the optimal value is z 0 = 1/2. By inserting Eq. (C.1) into the equations of motion one gets a set of N max + 1 equations for N max + 1 unknowns. Written in a linear algebra language, For a given value of q, the quasinormal modes correspond to frequencies w for which the system (C.2) admits a non-trivial solution, i.e. the quasinormal spectrum is determined by the equation det M nm (w, q) = 0. (C. 3) In this work, N max has been chosen to be 150. This method determines the quasinormal modes, i.e. the poles of the corresponding Green's functions but not the residues. The method is efficient in determining highly damped modes.
Integration When one is interested in the region of the complex frequency plane close to the origin, e.g. in hydrodynamic poles, or when the residue of the poles is needed, one could directly numerically integrate equations of motion, e.g. via the Runge-Kutta method [118,119]. Having already imposed the ingoing boundary condition at the horizon, there is still one boundary condition left. The numerical integration starts at the horizon (more precisely, at fixed small distance away from the horizon). In order to determine the starting values for the function and its derivative, one needs to iteratively solve the equations of motion by expanding the solution in series around the horizon. The freedom coming from the remaining boundary condition is encoded in a constant undetermined by this process which can be set to one without loss of generality. The quasinormal modes w are then determined by solving numerically the equation where Z * is a numerical solution and is a boundary cut-off. This method was used in this work to determine the spectral functions in Section 4.2.6.
Spectral method In Ref. [41], another numerical approach was used to compute quasinormal modes. It is based on the spectral method of solving ODEs. In this method, the regular functionZ is expanded in Chebyshev polynomials where c is the vector consisting of the N + 1 coefficients {c i } andM , an (N + 1) × (N + 1) matrix, consists of the evaluated equation of motion on the i-th grid point. A non-trivial solution exists only when detM = 0. This determinant is a polynomial in w and its roots correspond to the quasinormal modes, easily computed numerically. One of the advantages of this method is the rapid convergence of the modes with N , requiring a relatively small matrixM .
In this paper, we have used all three described methods, where convenient, finding consistent results (within numerical resolution).