Hydrodynamic dispersion relations at finite coupling

By using holographic methods, the radii of convergence of the hydrodynamic shear and sound dispersion relations were previously computed in the ${\cal N} = 4$ supersymmetric Yang-Mills theory at infinite 't Hooft coupling and infinite number of colours. Here, we extend this analysis to the domain of large but finite 't Hooft coupling. To leading order in the perturbative expansion, we find that the radii grow with increasing inverse coupling, contrary to naive expectations. However, when the equations of motion are solved using a qualitative non-perturbative resummation, the dependence on the coupling becomes piecewise continuous and the initial growth is followed by a decrease. The piecewise nature of the dependence is related to the dynamics of branch point singularities of the energy-momentum tensor finite-temperature two-point functions in the complex plane of spatial momentum squared. We repeat the study using the Einstein-Gauss-Bonnet gravity as a model where the equations can be solved fully non-perturbatively, and find the expected decrease of the radii of convergence with the effective inverse coupling which is also piecewise continuous. Finally, we provide arguments in favour of the non-perturbative approach and show that the presence of non-perturbative modes in the quasinormal spectrum can be indirectly inferred from the analysis of perturbative critical points.


Introduction
In the hydrodynamic regime, quantum field theory is expected to contain collective excitations such as sound waves [1,2]. These hydrodynamic modes are characterised in momentum space by their gapless dispersion relations ω = ω(q), where ω is the frequency of the mode and q is its wave-vector. In the simplest case of a relativistic neutral isotropic fluid, two hydrodynamic modes known as shear and sound modes have dispersion relations ω shear (q 2 ) = −iDq 2 + · · · , (1.1) These modes arise as linearised fluctuations of an equilibrium state and describe transverse momentum (shear) and longitudinal energy-momentum (sound) transfer. The coefficients of the series such as the speed of sound v s , the transverse momentum diffusion constant D = η/sT and the sound attenuation constant Γ = (ζ + 4η/3)/sT , where η and ζ are, respectively, shear and bulk viscosities and s is the equilibrium entropy density at temperature T , are determined by the underlying microscopic quantum field theory [2]. In the following, it will be convenient to use the frequency w = ω/2πT and the spatial momentum q = q/2πT normalised by the Matsubara frequency.
Recently, in the context of exploring the domain of applicability of hydrodynamics, the radii of convergence of the series (1.1), (1.2) have been investigated in some strongly interacting quantum field theories by using their dual gravitational descriptions in refs. [3][4][5][6][7][8] and by using field theory methods in refs. [9,10]. In particular, by promoting q 2 to a complex variable and analysing critical points of the associated spectral curves, in refs. [4,5], it was found that for the N = 4 supersymmetric SU (N c ) Yang-Mills theory (SYM) in the limit of infinite number of colours N c → ∞ and infinite 't Hooft coupling λ = g 2 Y M N c → ∞, the radii of convergence R = |q 2 | of the hydrodynamic series w = w(q 2 ) in the complex q 2 -plane are given by The physical reason behind the breakdown of the convergence of hydrodynamic series is the presence of the gapped non-hydrodynamic degrees of freedom whose spectra "cross levels" with the hydrodynamic degrees of freedom at some (generically complex) value of q 2 . Our main goal in this paper is to find the 't Hooft coupling constant corrections to the infinite coupling results (1.3), (1.4), similar to the coupling constant corrections to the entropy [11,12], shear viscosity [13,14] and other transport coefficients (see ref. [15] and references therein) computed for the N = 4 SYM theory earlier. Our methods are discussed in detail in section 2, and in appendices A and B.
Naively, one may expect that the radius of convergence R(λ) decreases with the coupling decreasing from its infinite value. Indeed, schematically [16], while at infinite coupling the characteristic spectral distance ν (∞) (set by the location of quasinormal modes in the dual gravity theory [17][18][19]) is coupling-independent, its counterpart ν (0) at small coupling (set by the eigenvalues of a suitable linearised collision operator [20,21], [22]) is coupling-dependent and parametrically small, hence, (1.5) where eq. (1.5) is clearly consistent with the results (1.3), (1.4). However, these expectations are shattered by a concrete calculation. Using perturbative methods only, in section 3 we find instead that at large coupling the radius of convergence increases with the coupling decreasing from its infinite value, namely,  (1.4). This result is unexpected. Admittedly, the large numerical coefficients in eqs. (1.7) and (1.8) may reflect the necessity of a non-perturbative "resummation" along the lines of ref. [23]. Indeed, applying such a resummation (discussed in detail below), we find that R shear (λ) and R sound (λ) become decreasing functions of the decreasing λ after the initial growth that is well-approximated by eqs. (1.7) and (1.8) (see fig. 1). In consequence, for the N = 4 SYM theory, our analysis implies that R(λ) should be a non-monotonic and piecewise continuous function. The dependence of R shear and R sound on γ ∼ λ −3/2 shown in fig. 1 constitutes the main result of this paper. The non-perturbative part of the curve R shear (γ) (the red segment of the curve in the right panel of fig. 1) coincides with the boundary of validity of hydrodynamics previously discussed in ref. [16]. We note that the piecewise character of this dependence is similar to the one recently observed for infinitely strongly coupled theories with finite chemical potential [7,24] and the Sachdev-Ye-Kitaev chain at finite coupling [9].
To understand the non-perturbative aspects of the analysis better, in section 4, we compute the radii of convergence of hydrodynamic series using the Einstein-Gauss-Bonnet gravity in five dimensions as a theoretical laboratory. There, the second-order bulk equations of motion can be solved fully non-perturbatively in the Gauss-Bonnet coupling [16,25,26], and thus the outcome of relevant perturbative resummations can be compared with exact results. We find (see figs. 13 and 16) that the radius of convergence decreases (and the dependence is piecewise continuous) with what can be phenomenologically identified as the direction of decreasing CFT coupling [16,[26][27][28], which is qualitatively similar to the results obtained for the N = 4 SYM theory in section 3.
The question of the hydrodynamic series convergence at finite coupling was recently addressed in ref. [10] for experimentally realisable fluids, and in ref. [9] for the Sachdev-Ye-Kitaev chain. The calculations of ref. [10] are based on estimating the size of the k-gap [29]  and show an increasing R with increasing Coulomb coupling strength, while ref. [9] finds a non-monotonic dependence, with R growing towards weak coupling. Radii of convergence of hydrodynamic series in relativistic kinetic theory (in the relaxation time approximation) were recently studied in ref. [30].
This paper is structured as follows. In section 2, we briefly review the method of critical points of spectral curves introduced in refs. [4,5] to compute the radii of convergence, as well as the non-perturbative resummation approach for theories with higher-derivative equations of motion. In section 3, we use the dual higher-derivative gravity to compute the radii of convergence for the shear and sound modes in the N = 4 SYM theory at large but finite 't Hooft coupling. This analysis is done perturbatively and non-perturbatively by using the "resummed" version of the first-order theory. In section 4, we perform similar calculations in Einstein-Gauss-Bonnet gravity to check the validity of our approach. We also demonstrate level-crossings at higher momenta. Then, in section 5, we discuss the validity of non-perturbative "resummations" used in holography. We first consider a toy algebraic example and then study in detail the shear channel of the Einstein-Gauss-Bonnet theory. These examples are used to draw plausible conclusions about the N = 4 SYM theory and the emergence of purely relaxing gapped modes in that theory. We conclude with a discussion of open problems in section 6. Appendix A is a short review of the methods involving critical points and quasinormal level-crossing. Appendix B introduces a useful method for determining the Puiseux exponent at a critical point by analysing the coefficients of hydrodynamic series of a dispersion relation w = w(q 2 ). Finally, appendices C and D contain the coefficients of the differential equations used in the paper.

Critical points, quasinormal level-crossings and the "resummation"
In this section, we briefly review the methods used to obtain the main results of the paper. These methods were formulated in refs. [4,5], where the interested reader can find more details and examples.

Critical points, Puiseux series and quasinormal level-crossing
In momentum space, the hydrodynamic dispersion relations arise from the hydrodynamic spectral curve P H (q 2 , w) = 0 given by the zeros of the determinant of the matrix of linearised fluctuations around an equilibrium state [4,5]. Using the symmetry of the system and applying the Newton polygon method, one can write generic expressions for the dispersion relations describing transverse momentum and longitudinal energy-momentum fluctuations in terms of the converging Puiseux series centred at the origin: The modes (2.1) and (2.2) are gapless. The coefficients c n and a n of the series are real functions proportional to the transport coefficients. In particular, c 1 = 2πT D, a 1 = ±v s , a 2 = −ΓπT . In the underlying quantum field theory, the full spectral curve P (q 2 , w), which reduces to P H (q 2 , w) in the hydrodynamic limit, is proportional to the denominator of the twopoint retarded correlation function of the corresponding conserved current (here, the energymomentum tensor). The spectral curve equation P (q 2 , w) = 0 contains the full spectrum of modes w = w(q 2 ), gapless and gapped. Identifying gapless modes with (2.1) and (2.2), one reads off the transport coefficients in terms of the quantum field theory parameters. Each Puiseux series is an expansion around a critical point (q 2 c , w c ) of order p which is a solution of the following set of equations: The order p determines the number of branches of the curve at the critical point. The analytic properties of the branches can be found by using e.g. the Newton polygon method, as explained in ref. [5] and references therein. Accordingly, a critical point may constitute a branch point singularity (or worse) or be a regular point depending on the coefficients of the original complex curve. For example, a point with p = 1 is always a regular point, as guaranteeed by the implicit function theorem, in which case a Puiseux series is the ordinary Taylor series. In terms of eq. (2.3), the shear mode is a Puiseux series in q 2 of order p = 1 around the origin (q 2 c , w c ) = (0, 0) (i.e., a Taylor series around a regular point), whereas for the sound mode, the origin is a branch point singularity generating a Puiseux series of order p = 2 [4,5].
The radii of convergence of the series (2.1) and (2.2) are set by the locations of the closest to the origin singularities of the functions w shear (q 2 ) and w sound (q 2 ), respectively, in the complex q 2 -plane. Critical points of the spectral curve P (q 2 , w) = 0 having branch points at q 2 = q 2 c are the common source of such singularities. At the critical points with p > 1, the equation P (q 2 c , w c ) = 0 has multiple roots, and hence the hydrodynamic mode "collides" with one or more gapped modes in the complex w-plane. If the corresponding q 2 = q 2 c is a branch point, this is the level-crossing phenomenon 1 , albeit happening here at complex values of frequency and momentum squared. We illustrate this with simple examples in Appendix A. In Appendix B, we also introduce a method based on the Darboux theorem which allows one to compute the Puiseux exponent at a critical point closest to (but different from) the origin by analysing the coefficients of a series centred at the origin.
Computing a spectral curve in quantum field theory, even perturbatively, is a difficult problem. However, for strongly interacting theories with gravity dual descriptions, this task is rather straightforward. Indeed, the recipe for computing the two-point retarded correlators from dual gravity [31] implies that the spectral curve P (q 2 , w) = 0 is determined by the boundary value Z(u = 0, q 2 , w) of the solution Z(u, q 2 , w) to the bulk equations of motion for the fluctuations coupled to the relevant conserved current: P (q 2 , w) = Z(u = 0, q 2 , w) = 0 [4,5]. Analytic properties of the spectral curve such as the location of branch point singularities are thus inherited from the properties of the bulk ODEs. In practice, the ODEs are sufficiently complicated and have to be solved numerically. Having such a solution, one first solves eqs. (2.3) to find the critical points in the complex q 2 -plane, and then determines the degree of singularity at the critical points by considering the quasinormal mode behaviour in the complex w-plane under the monodromy q 2 = |q 2 |e iφ , where φ ∈ [0, 2π] (see Appendix A). The closest to the origin (in the complex q 2 -plane) critical point exhibiting a branch point singularity sets the radius of convergence of the series (2.1), (2.2). In the N = 4 SYM theory at infinite N c and infinite 't Hooft coupling, the critical points closest to the origin in the shear and sound channels are located at Shear : Sound : leading to the radii of convergence (1.3) and (1.4). In general, a multitude of critical points is expected to exist in the complex q 2 -plane, representing level-crossings among two or more branches of the spectrum. Moreover, at finite N c or at weak coupling described by kinetic theory, one may also expect other types of singularities to appear [30].
Here, we continue working in the limit N c → ∞, and extend the approach of refs. [4,5] to bulk gravity theories with higher derivative terms, i.e., to the domain of large but finite 't Hooft coupling.

Non-perturbative "resummation"
Inverse 't Hooft coupling corrections in the N = 4 SYM theory arise from higher-derivative terms in the dual type IIB string theory low energy effective action (see e.g. refs. [11][12][13][14][15]). In a holographic calculation of a quasinormal spectrum, the bulk equations of motion typically produce a differential equation for the background fluctuation Z = Z(u, w, q 2 ) of the form [13][14][15] where u is the radial coordinate in the bulk with u = 0 the location of the boundary, γ is a small parameter proportional to the inverse coupling (e.g. γ ∼ λ −3/2 in the N = 4 SYM with the 't Hooft coupling λ), and the right hand side comes from the leading higher-derivative correction to Einsten-Hilbert action (e.g. from the R 4 term in type IIB supergravity). To avoid issues such as Ostrogradsky instability (see e.g. ref. [26] and references therein), the higher-derivative terms in eq. (2.6) are usually treated as perturbations of the second-order ODE, and its left-hand-side is used to eliminate all derivatives higher than the first one from the right-hand-side, ignoring contributions of order γ 2 and higher. The resulting equation, is a homogeneous linear second-order ODE whose coefficientsĀ andB now depend on γ. As discussed above, the spectral curve is then determined by the boundary value of the solution Z(u, w, q 2 , γ) to that ODE, i.e. P (q 2 , w) ≡ Z(u → 0, w, q 2 , γ) = 0.
In the standard approach, one looks for a perturbative solution to eq. (2.7) in the form Z = Z 0 + γZ 1 . Similarly, the perturbative ansatz for the spectrum is w = w 0 + γw 1 , where w 0 is the quasinormal frequency at γ = 0. Alternatively, one can solve eq. (2.7) without assuming a perturbative ansatz. Such a non-perturbative solution, if it can be expanded in series in powers of γ ≪ 1, will not be fully quantitatively correct beyond linear order in γ, since both in the original equation (2.6) and in the steps leading to eq. (2.7) terms of order γ 2 and higher were ignored. Quantitatively, the solution only captures the non-perturbative effects in γ related to eq. (2.7) and in this sense only partially "resums" the contributions nonlinear in γ in the approximation to the full solution. However, such a solution may provide a more faithful qualitative approximation to the exact solution at finite γ. Moreover, if the exact solution is non-perturbative in γ ≪ 1, any perturbative ansatz would necessarily miss it completely, whereas the non-perturbative approach is capable of describing the situation qualitatively correctly. The choice of a correct ansatz is the crucial step in singular perturbation theory [32]. We discuss these issues in more detail in section 5 and illustrate them with simple examples. In the context of holography, partial "resummations" have been used in refs. [23], [16,33] and criticised in ref. [34]. A crucial feature of such a "resummation" in holography, first pointed out in [16], is that the quasinormal spectrum now contains new, non-pertubative gapped modes which seem to play an important role in describing physics at finite coupling qualitatively correctly [16,33,35,36]. We shall see in section 3 that the situation with the radii of convergence is similar: the non-perturbative "resummation" reverses the tendency seen in eqs. (1.7), (1.8), making the radii to decrease (after an initial rise) with the coupling decreasing. In section 4, we compare this behaviour with that in the Einstein-Gauss-Bonnet theory, where both perturbative and non-perturbative results are available, using it as a theoretical laboratory to test our methods, and find a qualitative agreement with the N = 4 SYM case. Curiously, in section 5 we find that it is in fact possible in some cases to infer the existence of non-perturbative critical points by using perturbative data. While we are able to explicitly demonstrate this in the Einstein-Gauss-Bonnet theory, for the N = 4 SYM theory, this may serve as an indicative argument that the same behaviour is plausible.

Convergence of hydrodynamic series in the N = 4 SYM theory
We begin by studying the coupling dependence of the radii of convergence of the hydrodynamic shear and sound modes of the N = 4 SU (N c ) SYM theory in the N c → ∞ limit and at large but finite 't Hooft coupling λ. Our analysis uses its gravitational dual, namely, the type IIB supergravity with higher-derivative terms in the action. For the N = 4 SYM theory, 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, with α ′ set by the length of the fundamental string, and the term W 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 [37]. Moreover, as discussed in [38], for the purposes of computing the corrected quasinormal spectrum, one can use the Kaluza-Klein reduced five-dimensional action where W is now given by eq. (3.2) in 5d. The effective five-dimensional gravitational constant is related to the rank of the gauge group SU (N c ) by the expression κ 5 = 2π/N c . The parameter γ is related to the value of the 't Hooft coupling constant λ in the N = 4 SYM theory via γ = λ −3/2 ζ(3)L 6 /8. This parameter is dimensionless in units of L. Higher derivative terms in the equations of motion are treated as perturbations in γ. In the following, we shall use λ and interchangeably.
The black brane solution to the equations of motion following from the action (3.3), which is dual to an equilibrium thermal state of the CFT at temperature T , is given by [11,12] where f (u) = 1 − u 2 . The radial coordinate is denoted by u, with the boundary located at u = 0 and the horizon at u = 1. To leading order in γ, the functions A(u) and B(u) were found to be The correction to the N = 4 SYM entropy density is then given by [11] where s 0 = 2π 2 N 2 c T 3 /3 is the Stefan-Boltzmann entropy density of the ideal gas of particles in the N = 4 SYM theory (i.e., in the theory at λ = 0). The metric (3.5) was also used to compute 't Hooft coupling constant corrections to the ratio of shear viscosity η to entropy density [13,14], 8) and to all of the second-order transport coefficients of the N = 4 SYM plasma. 2 Computing transport coefficients and, in general, correlation functions of the energy-momentum tensor in the N = 4 SYM plasma requires considering small fluctuations of the metric g µν = g µν is the background (3.5). Due to translational invariance and spatial isotropy of the background, we can Fourier transform the fluctuations and choose the direction of spatial momentum along z, so that Following the recipes of ref. [19] and choosing the radial gauge h uν = 0, one can write down the linearised equations of motion for the three gauge-invariant linear combinations Z i , i = 1, 2, 3, of the modes h µν (u, ω, q) in the scalar, shear and sound channels, respectively [16,19,39]. The linearised equations of motion obtained following the procedure outlined in section 2 are given in the three channels by [16] . The coefficients are given explicitly in Appendix C. As discussed in section 2, Figure 2: Quasinormal spectrum in the shear channel of the N = 4 SYM for γ = 1 · 10 −5 and q 2 = 1 (left panel). The hydrodynamic shear mode at w ≈ −0.60064i is shown by the red circle. The new feature, not seen in a perturbative calculation, is the appearance of an extra mode on the imaginary axis (shown by the red square), ascending from complex infinity with γ increasing [16]. With real q 2 increasing, the hydrodynamic mode moves down the imaginary axis, while the new non-perturbative mode moves up. They collide at q 2 * ≈ 2.72 and move off the axis for q 2 > q 2 * (right panel).
As an example, the shear channel quasinormal spectrum for γ = 1 · 10 −5 is shown in fig. 2. Its novel feature, discussed in detail in ref. [16], is the existence of the non-perturbative (in γ) gapped modes on the imaginary axis. The highest (closest to the real axis) of those modes is shown in fig. 2 by the red square: with real q 2 increasing, this mode moves up the axis and at q 2 = q 2 * it collides with the hydrodynamic shear mode (shown in fig. 2 by the red circle). For q 2 > q 2 * , the two modes move off the imaginary axis, effectively destroying the diffusive pole of the correlator. In ref. [16], this was interpreted as the end of the hydrodynamic regime at sufficiently large spatial momentum (small wavelength), where the microscopic effects prevail over the collective ones. The dependence q 2 * = q 2 * (γ), shown in fig. 3, suggests that the domain of applicability of the hydrodynamic description is smaller at larger γ (i.e., at smaller 't Hooft coupling), but it seems to extend to arbitrarily large momentum in the limit of infinite coupling (at γ → 0) [16].
To see how this qualitative picture is amended at very large but finite 't Hooft coupling, we now consider the radius of convergence of the hydrodynamic shear and sound dispersion series in this theory, which requires us to solve eq. (2.3) and look for critical points with p = 2. With P (q 2 , w) given by P (q 2 , w) = Z(w, q 2 ) ≡ Z(u = 0, w, q 2 ) for any channel (we omit the index "i" labeling the channel), we are therefore looking for solutions w = w c and q 2 = q 2 c to the system with ∂ 2 w Z(w, q 2 ) ̸ = 0. In the non-perturbative approach (in γ), critical points follow directly   Figure 4: The closest to the origin branch points (see Table 1) of the N = 4 SYM shear mode dispersion relation w = w(q 2 ) in the complex q 2 -plane at γ = 1 · 10 −5 . The critical points seen in perturbation theory and their branch cuts are shown by blue colour, the two non-perturbative critical points on the real axis are shown in red. The radius of convergence at γ = 1 · 10 −5 , R shear ≈ 2.27, is determined by the pair of branch points at q 2 ≈ 1.93 ± 1.19i. We note that it is the sole inclusion of the non-perturbative red critical point on the positive real q 2 -axis that leads to the non-perturbative result of fig. 3. from eq. (3.11), where Z(w, q 2 ) can be found by constructing a Frobenius series solution Z(u, w, q 2 ) to eq. (3.10) around the horizon in the standard way [17][18][19], and setting u = 0. For γ = 0, this is the same procedure as the one used in refs. [4,5].
To find the critical points perturbatively, we expand equations (3.11) as and also expand Together, these expansions yield a system of equations at O(γ 0 ): (3.14) and a system at O(γ): Eqs. (3.14) and (3.15) are sufficient to find w c,0 , w c,1 , q 2 c,0 and q 2 c,1 . More explicitly, eq. (3.15) allows us to express As before, the function Z 0 (w, q 2 ) can be obtained as the boundary value of the Frobenius solution to eq. (3.10) with γ = 0 and Z 1 (w, q 2 ) as the boundary value of the Frobenius solution to the corresponding inhomogeneous equation.

Shear channel
At infinite 't Hooft coupling, the shear mode dispersion relation w = w(q 2 ) has numerous branch point singularities [4,5]. At finite coupling, we expect those singularities, now parametrised by γ ∝ λ −3/2 , to move in the complex q 2 -plane with γ varying. As discussed in section 2.2, one can compute relevant corrections by using either the "conservative" perturbative or the non-perturbative approach. # w c (non-perturbative) Table 2: The six closest to the origin (in the complex w-plane) critical points for γ = 1 · 10 −5 .

Perturbative calculation
For a perturbative calculation of the coupling constant correction to the radius of convergence in the shear channel of the N = 4 SYM theory, we use eq. (3.10) with i = 2. To first order in γ ∝ λ −3/2 , from eqs. (3.16) we find the following first set of critical points closest to the origin in the complex q 2 -plane: The value of q 2 c in eq. (3.17) gives the convergence radius |q 2 c | quoted in eq. (1.7). The radius increases with γ increasing (i.e. with the coupling λ decreasing from its infinite value). The numerical coefficients multiplying the parameter γ in eqs. The next closest to the origin critical point (i.e. the critical point with larger value of |q 2 | than (3.17)) is located on the negative real axis of q 2 :  Notice again the large numerical coefficients multiplying the perturbative parameter γ in eqs. (3.21)-(3.24). For illustration, several perturbative closest to the origin critical points in the complex q 2 -plane for γ = 1 · 10 −5 are shown in fig. 4 in blue colour.

Non-perturbative calculation
For a non-perturbative calculation, we solve eqs. (3.10) and (3.11) numerically without assuming γ to be small. We observe three qualitatively different scenarios of quasinormal modes' behaviour, and illustrate them by showing the modes at γ = 1 · 10 −5 , γ = 2 · 10 −5 and γ = 3 · 10 −5 , respectively (see figs.  Tables 1 and 2, where a comparison between perturbative and non-perturbative results is also made. For γ = 1 · 10 −5 , the location of the critical points is in a reasonably good agreement with the perturbative results (3.17), (3.18) and (3.19), (3.20), except when the collision of modes involves the mode on the imaginary axis which has no perturbative analogue. b) At γ = 2 · 10 −5 , the top modes in the spectrum are shown in fig. 7. Here, the first level-crossing (i.e., the level-crossing with the minimal value of |q 2 |) occurs between the two top gapped modes and the non-perturbative mode on the imaginary axis, as shown in the top left panel of fig. 7. However, the shear mode is not affected by this crossing: its first non-analyticity still arises as a result of the collision with the top two gapped modes as shown in the top right panel of fig. 7. This collision sets the radius of convergence of the shear mode at R shear (γ) = |q 2 c | ≈ 2.31. In the complex q 2 -plane, the position of the branch point singularities is thus qualitatively the same as at γ = 1 · 10 −5 (see fig. 6, left panel), but the radius of convergence R shear (γ) increases with γ increasing (i.e. with the coupling λ decreasing). c) Finally, at γ = 3 · 10 −5 , the top modes in the spectrum are shown in fig. 8. In this case, the situation is qualitatively different. Now the first level-crossing occurs at a real value of q 2 c ≈ 2.12157 and at w c ≈ −2.096i, which is a result of the collision between the shear mode and the non-perturbative mode on the imaginary axis (top left panel in fig. 8). This is the collision of the type shown in fig. 2 and interpreted in ref. [16], where it was discovered, as the end point q 2 * (γ) of the hydrodynamic regime (in the sense that for real q 2 > q 2 * the hydrodynamic purely imaginary shear mode does not exist). The radius of convergence in the q 2 -plane is set by the corresponding value of |q 2 * | = |q 2 c | ≈ 2.12157 (see fig. 6, right panel, where this situation is shown schematically). Thus, at γ = 3 · 10 −5 , the radius of convergence is determined by the non-perturbative mode. In this regime, the convergence radius decreases with γ increasing.
The three examples considered above illustrate the general situation fully. For infinitesi- mally small γ, the radius of convergence of the shear mode dispersion relation is an increasing function of γ. In the complex q 2 -plane, the obstacle to convergence is the pair of critical points, as shown in fig. 6 (left panel). These points move away from the origin with γ increasing. At γ = γ * ≈ 2.05 · 10 −5 , the picture changes qualitatively, as the transition between the regimes a) and c) occurs. Now the new critical point arising from the level-crossing with the non-perturbative mode is located closer to the origin in the complex q 2 -plane than the previous pair of critical points ( fig. 6, right panel). This new critical point is located on the positive real axis of q 2 and corresponds to the end point of hydrodynamic regime as discussed in ref. [16]. This point moves closer to the origin with γ increasing. The dependence of the radius of convergence on γ is thus a piecewise continuous function 3 , shown in fig. 1 (left panel), which is the main result of this section. The results of the perturbative and non-perturbative calculations for the closest to the origin critical points (the ones seen in perturbation theory) are in good agreement, as can be observed from Table 1.
We emphasise that the non-perturbative effects discussed in this section are at best qualitative, since terms of order γ 2 and higher will inevitably modify them. We believe these effects are qualitatively correct as they fit well with various physical expectations [16,23,33,36]. In particular, the existence of the non-perturbative critical point makes the radius of convergence decrease with the 't Hooft coupling decreasing from its infinite value. Admittedly, an alternative conservative point of view would simply limit any considerations by the range of γ up to γ ≲ 2 · 10 −5 , beyond which the perturbative treatment becomes unreliable.

Sound channel
For the sound channel, the analysis follows the strategy used in the previous section and in refs. [4,5] very closely. The relevant equation of motion is eq. (3.10) with i = 3.

Perturbative calculation
Solving the equations perturbatively to linear order γ, we find the following correction to the location of the closest to the origin critical point  Figure 8: Quasinormal spectrum in the shear channel at γ = 3 · 10 −5 , computed non-perturbatively in γ. The radius of convergence is set by the level-crossing between the shear mode and the nonperturbative mode on the imaginary axis (left panel). This critical point coincides with the endpoint of the hydrodynamic regime discussed in ref. [16]. An example of critical points with a larger |q 2 | is shown in the right panel.
Eq. (3.25) gives the radius of convergence R sound = |q 2 c | stated in eq. (1.8). The next critical point is given by As in the case of the shear mode, the coefficients of the perturbative expansion are quite large: the correction becomes comparable to the γ = 0 result already for γ ∼ 10 −4 − 10 −5 .

Non-perturbative calculation
Non-perturbative treatment implies solving the equation of motion (3.10) without assuming γ to be small. As in the shear mode case, there are essentially two qualitatively different scenarios of the distribution of critical points. We illustrate them by showing trajectories of quasinormal modes in the complex w-plane at complex q 2 = |q 2 |e iφ , φ ∈ [0, 2π], for γ = 2·10 −5 and γ = 4.5 · 10 −5 in figs. 9 and 10, respectively. In fig. 9, left panel, the critical point limiting the radius of convergence of the sound mode's dispersion relation arises from the level-crossing between that mode and the top gapped modes (this situation is qualitatively the same as at γ = 0 [4,5]). The pair of critical points closest to the origin in the complex q 2 -plane is well approximated by the perturbative expression (3.27) for γ ≲ 3 · 10 −5 ( fig. 1, right panel). The radius of convergence, R sound = |q 2 c |, increases with γ (see the blue curve in fig. 1, right panel).
For γ > γ * ≈ 3.225 · 10 −5 , the situation changes qualitatively, as illustrated in fig. 10. Now the sound mode first collides with the non-perturbative mode (in fig. 10, i.e. at γ = 4.5 · 10 −5 , this happens at q 2 c = 0.9568 ± 1.7083i, w c ≈ ∓0.16513 − 1.8681i, corresponding to R sound = |q 2 c | ≈ 1.958). In this regime, the radius of convergence, with γ (see the red curve in fig. 1, right panel). This dependence is very similar, although not identical, to the one observed in fig. 1 (left panel) for the shear channel. As before, in a conservative approach we would limit ourselves to the blue part of the curve in the right panel of fig. 1, which ends in the region where perturbation theory becomes unreliable. We believe, however, that the red part of the curve, although not quantitatively precise, reflects the dependence of the radius of convergence on coupling qualitatively correctly. The piecewise smooth dependence on the coupling, shown in fig. 1 (right panel), is likely to persist even with γ 2 and higher terms in the action taken into account quantitatively correctly, since it corresponds to the discrete change of "status" of the closest to the origin critical point, even though the critical points move continuously in the complex plane with varying coupling.

Convergence of hydrodynamic series in the Einstein-Gauss-Bonnet theory
We now consider the radii of convergence of gapless quasinormal modes in five-dimensional Einstein-Gauss-Bonnet gravity. Although this theory may not have a healthy QFT dual [41], it is nevertheless a very useful theoretical laboratory for relevant bulk calculations in higherderivative gravity: by design, its equations of motion are second-order in derivatives and can thus be solved fully non-perturbatively in terms of the higher-derivative coupling. The Einstein-Gauss-Bonnet action in 5d is 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 a dimensionless parameter. The black brane metric solution of Einstein-Gauss-Bonnet equations of motion can be found analytically and is given by 4 where and the constant N GB can be chosen to normalise the speed of light at the boundary to c = 1: The position of the horizon is at r = r 0 . The Hawking temperature corresponding to the solution (4.2) 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 onto γ GB ∈ [0, 1), with λ GB = 0 corresponding to γ GB = 1.
To compute the quasinormal mode spectrum, we again consider linearised metric perturbations. In analogy with our discussion in section 3, we arrive at three gauge-invariant equations of motion for the scalar, shear and sound channels (i = 1, 2, 3, respectively): where the coefficients A (i) and B (i) are given in Appendix D. All relevant details regarding the theory and the derivation of these equations can be found in ref. [26].

Shear channel
The shear channel spectrum is determined by eq. (4.6) with i = 2. As in the case of the N = 4 theory, we compare calculations done perturbatively and non-perturbatively in the Gauss-Bonnet coupling λ GB .

Perturbative calculation
Solving eqs.  The radius of convergence for the shear hydrodynamic mode is then given by R shear = |q 2 c |: The coefficient in front of λ GB is a positive number, which is similar to the N = 4 case. However, the "physical" regime of the hypothetical field theory dual to Gauss-Bonnet gravity corresponds to λ GB < 0 [16,[26][27][28]47]. In this case, the radius of convergence decreases with |λ GB | increasing, the trend opposite to the one found for the N = 4 SYM theory in section 3.1.
We also compute perturbative λ GB corrections to the next three higher critical points, finding   fig. 11 and fig. 12. Their characteristic feature, explored in refs. [16,26], is the presence of non-perturbative modes located (for real q 2 ) on the imaginary axis in the complex w-plane. At sufficiently small values of |λ GB |, these modes lead to new critical points only for large values of |q 2 |: the closest to the origin critical points setting the radius of convergence of the shear mode are not affected by them (see fig. 11, where the spectrum is shown for λ GB = −0.01 and various values of complex q 2 ). At larger values of |λ GB |, however, the situation changes qualitatively. Now the closest to the origin critical point and thus the radius of convergence are set by the non-perturbative mode (this is illustrated by fig. 12, where the spectrum is shown for λ GB = −0.03). The transition between the two regimes occurs at λ GB ≈ −0.0198. The dependence of the radius of convergence of the shear mode's dispersion relation on Gauss-Bonnet coupling is shown in fig. 13.

Sound channel
Finally, we repeat the same analysis for the sound channel of the Einstein-Gauss-Bonnet theory using eqs.

Perturbative calculation
To linear order in the perturbative expansion in λ GB , we find the closest to the origin pair of critical points at  Hence, the radius of convergence in the sound channel is given (perturbatively) by For "physical" values of λ GB (λ GB < 0), this is a decreasing function of |λ GB | (which is different from the N = 4 SYM theory).

Non-perturbative calculation
Computing critical points in the sound channel of the Einstein-Gauss-Bonnet theory nonperturbatively in λ GB , we find that for sufficiently small |λ GB |, the situation remains qualitatively the same as in the λ GB = 0 case: the radius of convergence is determined by the level-crossing between sound modes and the top gapped modes in the complex w-plane, as shown in fig. 14 (left panel). The non-perturbative mode also leads to critical points, but this happens at larger value of |q 2 | (see the right panel of fig. 14). At larger value of the coupling |λ GB |, the picture changes qualitatively and the radius of convergence of the sound dispersion relation is now determined by the level-crossing with the non-perturbative quasinormal mode, as shown in the left panel of fig. 15. The transition between the two regimes happens at λ GB ≈ −0.0338. The dependence of R sound (including both perturbative and non-perturbative results) is shown in fig. 16.

Non-perturbative quasinormal modes and singular perturbation theory
In our analysis of the radii of convergence in the N = 4 SYM theory, as well as in previous works on higher-derivative holography [23], [16,33,35,36], the non-perturbative "resummation" and the appearance of the non-perturbative quasinormal modes played a rather prominent role. In particular, the existence of these non-perturbative features seems to be fully consistent with the requirement of a physically reasonable interpolation between strongly coupled (holographic) regime and weakly coupled (e.g. kinetic) regime in the same theory. This interpolation, even for simplest theories such as CFTs considered at finite temperature, is not fully understood [16,33,36,[48][49][50][51][52]. Admittedly, relying -even only qualitatively -on the non-perturbative treatment in models arising as truncations of a perturbative expansion may seem to be unwarranted [34]. However, before dismissing such an approach as ineffable nonsense, one may wish to consider examples where it is known to be successful. In this section, we first outline the problem as we see it, and then discuss in detail a simple example of an algebraic equation containing a small parameter, where similar issues arise.
Consider the full set of solutions to an equation (algebraic or differential), which we 2 (ϵ), . . .} and so on. A natural question to ask is in what sense the solutions X (0) , X (1) , X (2) , . . . approximate the true solution X. Note that even the number of roots in each "truncated" set X (n) depends on n. For example, the equation has a single solution X (0) = {x 2 (ϵ)}, one of which is perturbative and another one is non-perturbative in ϵ: The solution x (1) 1 (ϵ) can be constructed order by order in ϵ via standard perturbation theory (i.e., assuming a perturbative ansatz for the solution), whereas the x (1) 2 (ϵ) is "invisible" in the standard perturbative approach yet it can be built consistently using singular perturbation theory [32]. Now imagine that eq. (5.4) is a truncation of the equation At any finite N , among the N + 1 roots of the equation (5.7), one is perturbative in ϵ (it is a regular perturbation of the solution x = 1 of the equation (5.3)) and the remaining N roots are non-perturbative: they disappear to infinity in the limit ϵ → 0. These extra N roots are located (approximately) along the circle |x| = 1/ϵ in the complex x-plane. Finally, we note that eq. (5.7) can be regarded as a truncation at order ϵ N of an exact function Note that eq. (5.8) has only one solution, whose small ϵ expansion coincides (for |ϵ| < 1) with the perturbative solution of the equation (5.7) at the appropriate order. Thus, in this example, truncating the small ϵ expansion of eq. (5.8) at order N produces one correct and N spurious roots located approximately at the boundary of analyticity of the function L (∞) [x, ϵ], i.e. at |x| = 1/ϵ. This simple example seems to reinforce the "conservative" approach to the quasinormal spectrum in higher-derivative gravity suggesting that only perturbative solutions can be trusted. Non-perturbative solutions exist (and can be constructed via singular perturbation theory) at each given order of the expansion in a small parameter, but these solutions appear to be artefacts of the expansion and disappear when the full function is considered. However, such a verdict might be too quick, as the example in the next subsection shows.

An algebraic equation example
Consider the algebraic equation where ϵ is a parameter. At ϵ = 0, eq. (5.10) has a single root, x = x 0 = −i. For ϵ > 0, however, there are infinitely many solutions, parametrised by ϵ (see fig. 17, left panel, where the roots of eq. (5.10) closest to the origin in the complex x-plane are shown for ϵ = 0.5). Note that all these solutions but one are non-perturbative in ϵ, since they must disappear from the set of solutions in the limit ϵ → 0 leaving the single root x = x 0 at ϵ = 0.
The solutions to eq. (5.10) can be constructed as series in ϵ ≪ 1. One such solution is the finite ϵ correction to the solution x 0 = −i of the equation at ϵ = 0. Using the standard perturbation theory, we find The non-perturbative roots can be found analytically by using the methods of singular perturbation theory [32]. Introducing a new variable x =x/ϵ and taking the limit of ϵ → 0 in eq. (5.10) while keepingx fixed, we find the equation The infinite set of solutions to this equation,x n = iπ(1 + 4n)/2, where n ∈ Z, then gives all the non-perturbative roots of the original eq. (5.10) as The series in eqs. (5.11) and (5.13) converge 5 for ϵ < |ϵ c | ≈ 0.2625. We now replace L[x, ϵ] by a finite polynomial (a truncated Taylor expansion of eq. (5.10)), and ask whether the 2N + 2 roots of eq. (5.14) approximate the exact solutions of eq. (5.10) as we increase N . Our previous discussion implies that at least one such approximation, a regular perturbative correction to the zeroth-order root x 0 = −i, should exist. Indeed, solving the corresponding equations L (2) [x, ϵ] = 0, L (4) [x, ϵ] = 0, etc., perturbatively in ϵ, we reproduce the series (5.11) term by term. All others roots are non-perturbative in ϵ, disappearing to complex infinity in the limit ϵ → 0.
To mimic our approach to the quasinormal spectra in higher-derivative gravity, we now solve eq. (5.14) numerically at each order of N , without assuming ϵ to be small. A generic result is illustrated in the right panel of fig. 17, where ϵ = 0.5, and both the exact solutions to eq. (5.10) and the roots of the polynomial (5.14) with N = 50 are shown. It is clear that in a (large) region of the complex plane containing the origin, the solutions of eq. (5.10) (both perturbative and non-perturbative) are well approximated by some of the roots of the polynomial (5.14). There are also "spurious roots", forming the top and the bottom red "arcs" in fig. 17, right panel: they do not approximate any solution. The exact solutions of eq. (5.10) located outside of the domain bounded by the red arcs are not approximated by any of the roots of the polynomial (5.14) and remain "invisible". The domain bounded by the red arcs of spurious roots increases with N increasing. Our main conclusion is that it is possible in principle to approximate at least some of the exact non-perturbative solutions by the non-perturbative roots of the perturbative truncation of the exact equation. Of course, this example is not a justification of our "resummation" of quasinormal spectra but we hope it shows that such an approach has the right to exist. We now ask a different question. As noted in footnote 5, the radius of convergence of the The solution of eqs. (5.15) with the smallest |ϵ| is x ≈ −2.1176i, ϵ ≈ −0.2625. In the complex x-plane, the critical point corresponds to the level-crossing between the perturbative mode x 0 and the nearest non-perturbative solution of eq. (5.10), as shown in fig. 18. Now suppose that we only have access to the successive polynomial approximations (5.14) to the full equation (5.10). Can we detect the existence of the correct critical point and hence the non-perturbative solution by using only these polynomial approximations and the regular perturbation theory? The answer is affirmative: solving perturbatively eqs. (5.15) with L[x, ϵ] replaced by the polynomial approximation (5.14) of order 2N + 2, we find that the (smallest in magnitude) critical value ϵ c (N ) converges to the "exact" result ϵ c ≈ −0.2625 with N increasing (see fig. 19, left panel). Since the critical point corresponds to the level-crossing (see fig. 19, right panel), i.e., a collision between perturbative and non-perturbative roots in the complex x-plane (see fig. 18), the above observation suggests that the perturbative calculation "knows" about this collision, even though the non-perturbative root responsible for it is inaccessible in perturbation theory and remains "invisible". This raises a question of whether the non-perturbative quasinormal modes arising in higher-derivative gravity can be indirectly detected by a perturbative analysis of the relevant critical points. We investigate this question in the next section.

A signature of non-perturbative modes from the perturbative analysis of the Einstein-Gauss-Bonnet critical points
The quasinormal spectrum in the shear channel of Einstein-Gauss-Bonnet gravity contains modes located on the imaginary axis which are non-perturbative in the Gauss-Bonnet coupling λ GB (see fig. 11 and refs. [16,26]). As discussed in section 4, for some range of parameters, the top (closest to the origin) non-perturbative mode collides with the hydrodynamic shear mode at real q 2 c and purely imaginary w c , as illustrated in fig. 12. This collision is followed by another collision between the emergent pair of propagating modes with the pair of gapped modes from the standard "Christmas tree" sequence. For a fixed value of λ GB , the second collision again occurs at real q 2 c . We illustrate this in fig. 20, where the shear channel Einstein-Gauss-Bonnet spectrum is shown at λ GB ≈ −0.04956 for several values of real q 2 : as q 2 is increased from q 2 = 0.81 to q 2 = 1.17, the hydrodynamic mode and the non-perturbative mode on the imaginary axis approach each other, colliding at q 2 c ≈ 1.18544 and forming a pair of propagating modes, which move off the axis and at q 2 c ≈ 1.87785 collide again, now with the pair of the "Christmas tree" gapped modes. The two (sets of) critical points resulting from these collisions are given by We emphasise that the existence of the critical points (5.17) and (5.18) is the direct consequence of the existence of the non-perturbative mode on the imaginary axis, absent at λ GB = 0. Indeed, for λ GB = 0 (which corresponds to the N = 4 SYM theory at infinite 't Hooft coupling), the hydrodynamic shear mode travels unobstructedly along the imaginary w axis towards negative infinity as the real q 2 is increased: no critical point exists in the theory for purely imaginary w c at purely real and positive q 2 c . In contrast, at finite λ GB , the shear mode collides with the non-perturbative mode on the imaginary axis at real and positive q 2 c given by eq. (5.17), and then the resulting two propagating modes give rise to the pair of critical points at another real value of q 2 c (eq. (5.18)). The pair of the propagating modes leading to the pair of critical points in eq. (5.18) does not exist in the perturbative spectrum: it is created by the collision of the shear mode and the non-perturbative mode on the imaginary axis.
Can the existence of non-perturbative critical points such as the ones in eq. (5.18) be inferred from the perturbative data (i.e., from eqs. (2.3), solved perturbatively in λ GB ), in analogy with what has been observed in section 5.1? The goal is to find a perturbative approximation to the pair of critical points (5.18) with real q 2 .
Perturbatively, to leading order in λ GB , the closest to the origin critical point is given by eq. (4.7), reproduced here for convenience: Eq. (5.19) implies that the critical value of q 2 is purely real for λ GB ≈ −0.0519780, and is given by We conjecture that the critical point (5.22) at λ GB ≈ −0.0519780 is the perturbative approximation to the "exact" critical point (5.18) at λ GB ≈ −0.04956. To see that this is indeed the case, we extend the perturbative expansion in eqs. (5.19), (5.20) to order O(λ N GB ). At each N , we then numerically compute the value of λ GB (denoted by λ GB (N )) from the algebraic condition that sets Im q 2 c (λ GB ) = 0 (among multiple roots λ GB (N ), we choose the solution which is real and closest to the non-perturbative value λ GB ≈ −0.04956). Then, we use thus obtained λ GB (N ) to evaluate w c (N ) and q 2 c (N ) from the perturbative series in λ GB . The numerical sequences w c (N ), q 2 c (N ) do not converge and have to be Padé-resummed. The Padé-resummed values then converge to the corresponding values for the non-perturbative critical point (5.18) with N increasing (see fig. 21).
Hence, a perturbative analysis of eqs. (2.3) (admittedly, aided by a Padé resummation) seems to be capable of reproducing at least one of the non-perturbative critical points of the full quasinormal spectrum. Since such a point can only occur due to the presence of a nonperturbative mode in the spectrum, we conclude that the calculation indirectly confirms the existence of the non-perturbative mode itself.
Another example of an inherently non-perturbative critical point that can be reproduced this way is the critical point located at λ GB ≈ −0.02327 : The appropriate perturbative critical point is given by eq. (4.14). Requiring Im q 2 c (λ GB ) = 0, to first order in λ GB we have λ GB ≈ −0.024486 : Continuing the expansion (5.24) to higher orders of λ GB as described above, we find that the sequence of order-N critical points λ GB (N ), w c (N ), q 2 c (N ) converges to the non-perturbative result (5.23) with N increasing.

A signature of non-perturbative modes from the perturbative analysis of critical points in the N = 4 SYM theory
We now turn to the perturbative analysis of critical points in the N = 4 SYM theory. Can we detect the existence of non-perturbative modes in the quasinormal spectrum by using the approach of the previous section? The recipe is to look for a pair of critical points with real value of q 2 . Considering the lowest perturbative critical point (3.17) and finding γ from the condition Im q 2 c (γ) = 0, we find to leading order in γ: γ ≈ −0.0006290 : Such a pair of critical points at a single real q 2 c does not exist in an infinitely strongly coupled theory (at γ = 0). At finite γ, it can arise as a result of the collision on the imaginary axis between the shear mode and the non-perturbative mode, followed by the second collision between the pair of the resulting propagating modes and the two gapped "Christmas tree" modes, similarly to what happens in the Einstein-Gauss-Bonnet case as described in section 5.2. We conjecture that eqs. (5.25) are approximations to the non-perturbative critical point appearing in the right panel of fig. 8. Of course, unlike in the Einstein-Gauss-Bonnet theory, here we have no access to γ 2 and higher terms in the action and thus cannot verify this conjecture. Nevertheless, we can interpret the existence of the point (5.25) as an indirect evidence for the existence of the non-perturbative modes in the full quasinormal spectrum of the gravitational background dual to the N = 4 SYM theory.

Discussion
In this paper, we have determined the dependence on the coupling of the radii of convergence of the shear and sound hydrodynamic dispersion relations in the strongly coupled N = 4 SYM theory. This dependence is shown in fig. 1. Limiting ourselves to the results obtained using the standard perturbation theory only, the coupling constant dependence of the radii in the shear and sound channels is given by eqs. (1.7), (1.8) (the blue lines in fig. 1), respectively. These perturbative results suggest that the radii of convergence increase with the 't Hooft coupling decreasing from its infinite value.
However, we argued that the presence of non-perturbative modes in the quasinormal spectrum of the dual black brane background at large but finite coupling will modify the perturbative result by making the radii's dependence on the coupling piecewise continuous, as shown in fig. 1. This type of dependence is familiar from the finite density examples of holographic theories at infinite coupling [7,24] and from the Sachdev-Ye-Kitaev (SYK) chain at finite coupling [9]. The origin of this similarity is not entirely clear to us. 6 In the SYK case, it is the coupling dependence that gives rise to the piecewise behaviour, in parallel to what we observe in this paper. The strong-weak dependence in that case is, however, reversed. In the charged case, the chemical potential normalised by temperature plays the role of the coupling in the present paper being responsible both for the extra modes on the imaginary axis and the piecewise continuous dependence of the radius of convergence. One plausible guess is the emergence of approximate symmetries in all these cases in line with what was discussed in ref. [35].
Curiously, in the regime where the non-perturbative mode becomes relevant, the radius of convergence of the shear mode dispersion relation coincides with the endpoint q 2 c (λ) of the hydrodynamic regime introduced earlier in ref. [16]. We then repeated our analysis for the case of the Einstein-Gauss-Bonnet theory using it as a theoretical laboratory to test our methods. Again, due to the presence of non-perturbative quasinormal modes in the spectrum, we found the piecewise dependence of the radii of convergence on the coupling (see figs. 13 and fig. 16 for the shear and sound channel results, respectively). The role of the non-perturbative modes in a quasinormal spectrum thus appears to be significant. Their properties are similar to the properties of non-perturbative roots of algebraic equations with a small parameter, where singular perturbation theory is often useful. We have discussed in detail a simple example of an algebraic equation whose perturbative and non-perturbative roots can be consistently found using singular perturbation theory. Although at present we cannot offer a generalisation of such a discussion to differential operators, we believe a qualitatively similar picture should hold there, too, and therefore the non-perturbative quasinormal modes can represent the qualitatively correct feature of the full theory. Finally, we show that the presence of the nonperturbative modes in the spectrum can be indirectly inferred by the analysis of critical points using the perturbative data only. From the work in refs. [4,5] it is clear that thermal two-point functions of the energymomentum tensor in the N = 4 SYM theory at infinite 't Hooft coupling contain multiple branch point singularities in the complex plane of q 2 . At large but finite coupling, the location of these singularities acquires a dependence on the coupling. Our analysis in the present paper suggests that the property of being a singularity closest to the origin (and thus setting the radius of convergence) may change discretely as a function of the coupling, leading to the piecewise nature of the dependence of the radii of convergence on coupling, and that such a change is induced by the presence of non-perturbative quasinormal modes in the spectrum of a dual gravitational theory. It would be interesting to investigate whether a similar phenomenon is observed at small but finite coupling. Of special interest are the studies of the hydrodynamic series convergence in strongly coupled theories with non-zero chemical potential [6,7], [53], where the coupling dependence has not yet been explored.
is contained in the quasinormal spectrum of the bulk equations of motion. In terms of a gauge-invariant variable Z, a typical bulk equation of motion is a second-order ODE where u is the bulk radial coordinate with u = 0 corresponding to the boundary of the dual gravity background. The dependence of the coefficients of the equation on q 2 reflects the rotation invariance of the theory, whereas the dependence on w is a consequence of the choice of the boundary condition at the horizon. The quasinormal spectrum w i = w i (q 2 ) is determined by the equation which also defines the spectral curve. The critical point condition (2.3) means that p ≥ 2 branches w i = w i (q 2 ) collide at (w c , q 2 c ) ∈ C 2 . The branches are locally represented by Puiseux series whose analytic structure is determined by the coefficients of the spectral curve via, e.g., the Newton polygon method. The situations when w i ∼ (q 2 − q 2 c ) −ν , where ν is a fractional power, e.g. ν = −1/2, are common. The closest to the origin (in the complex q 2 -plane) critical point with such a branch point singularity limits the convergence of the series w i = w i (q 2 ) centered at q 2 = 0 and thus sets its radius of convergence, R = |q 2 c |. This is the phenomenon of quasinormal level-crossing, analogous to the quantum-mechanical level-crossing [4,5]. When ν is a negative integer or zero, the branches are locally analytic, and we have "level-touching" rather than "level-crossing" [5].
In practice, the bulk ODEs are sufficiently complicated and have to be solved numerically. With such a solution in hand, one first solves eqs. (2.3) to find the critical points in the complex q 2 -plane, and then determines the degree of the singularity at the critical points by considering the quasinormal mode behaviour in the complex w-plane under the monodromy q 2 = |q 2 |e iφ , where φ ∈ [0, 2π]. We illustrate the difference between "level-crossing" and "level-touching" by the following simple example.
Consider the complex curves where x 2 ∈ C, y ∈ C, and the coefficients a, b, c are some fixed complex numbers. Applying the criterium (2.3) to P 1 and P 2 , we find that both curves have the p = 2 type crtitical point at (x 2 c , y c ) = (b/c, a). The two branches of the curve P 1 are given by whereas the two branches of P 2 are y (1) Figure 22: Level-touching: the two branches y (1,2) 1 (x 2 ) of the spectral curve P 1 at complex x 2 = |x 2 |e iφ , with φ varying from 0 to 2π, at fixed |x 2 1 | < |x 2 c | (left panel) and |x 2 2 | > |x 2 c | (right panel). 2 (x 2 ) of the spectral curve P 2 at complex x 2 = |x 2 |e iφ , with φ varying from 0 to 2π, at fixed |x 2 1 | < |x 2 c | (left panel) and |x 2 2 | > |x 2 c | (right panel).
For the curve P 1 , the two branches at the critical point are analytic functions of x 2 . For the curve P 2 , the critical point is a branch point singularity. Numerically, one can distinguish between the regular and the singular behaviour by finding local solutions y i = y i (x 2 ) to eqs. (A.3), (A.4) at complex x 2 = |x 2 |e iφ , with fixed |x 2 | and φ ∈ [0, 2π]. The "trajectories" traced by the branches as the phase φ varies from 0 to 2π are shown in fig. 22 for the curve P 1 and in fig. 23 for the curve P 2 at fixed |x 2 1 | < |x 2 c | < |x 2 2 |. The regular branches of P 1 approach and touch each other at x 2 = x 2 c , with individual trajectories remaining closed under monodromy at |x 2 | > |x 2 c |, as shown in fig. 22. This is "level-touching", observed e.g. in the case of the BTZ quasinormal spectrum [5]. In case of a branch point singularity at the critical point, the two branches merge into a single trajectory for |x 2 | > |x 2 c |, as they are mapped into each other under the monodromy ( fig. 23). This is the level-crossing phenomenon [4,5].
B How to reconstruct the Puiseux exponent from the power series Given the hydrodynamic expansion (e.g. (2.1) or (2.2)) of the dispersion relation w = w(q 2 ) around q 2 = 0, one can determine, at least in principle, the exponent of the Puiseux expansion w ∼ (q 2 − q 2 c ) −ν at the nearest to the origin critical point q 2 c , given by eq. (2.3). This can be done by using the Darboux theorem (see e.g. [54], Theorem 11.10b) which states that if a function p(t) has a singularity at t = t 0 of the form where r(t) is an analytic function and ν is not a negative integer or zero, then the coefficients a n of the Taylor expansion of p(t) at the origin, p(t) = ∞ n=0 a n t n , have the following asymptotic form as n → ∞ a n ∼ Γ(n + ν) Γ(ν)t n 0 n! r(t 0 ) − (ν − 1)t 0 r ′ (t 0 ) (n + ν − 1) + (ν − 1)(ν − 2)t 2 0 r ′′ (t 0 ) 2!(n + ν − 1)(n + ν − 2) + · · · . (B.2) Keeping only the leading term in eq. (B.2), we find ν = lim n→∞ t 0 (n + 1) a n+1 a n − n .

(B.3)
Then, using the subdominant terms in eq. (B.2), one can in principle reconstruct the function r(t) by recovering its derivatives. In practice, a finite (and often relatively small) number of the coefficients a n (e.g. computed numerically) is sufficient to determine ν with a good precision. Thus, if t 0 is the critical point closest to the origin, one can determine (or, strictly speaking, conjecture) whether this point is a singularity of p(t) (i.e., whether ν is a fractional or positive integer number) and therefore sets the radius of convergence and reconstruct the Puiseux expansion around t 0 by analysing the coefficients a n . A complication arises when there are two or more critical points located at the same distance from the origin. This is the case for the hydrodynamic dispersion relations in the N = 4 SYM theory at infinite 't Hooft coupling, where for the shear and sound modes we have a pair of complex conjugate critical points in the complex q 2 -plane [4]. Such cases were studied for example in ref. [55]. Instead, we find it more convenient to reduce the case with two critical points to the previous one with the help of a conformal transformation. Indeed, let t = t 1 and t = t 2 , where |t 1 | = |t 2 |, be the critical points located at the same distance from the origin t = 0. By performing a Möbius transformation and requiring that under the map, 0 → 0, t 1 → 1 and t 2 → z 2 , where |z 2 | > 1, we reduce the situation to the one of the Darboux theorem: e.g., the exponent ν 1 of the branch point t = t 1 is determined by applying the Darboux procedure to the point z = 1. The inverse transformation is t = dz − b a − cz . (B.5) We must also require that the singularity at z ≡ z s = a/c stays outside of the unit circle in the complex z-plane, so that any analytic part of p(t) remains analytic after the transformation.
Applying this procedure to the hydrodynamic series of the N = 4 SYM theory using the data of ref. [4], we find ν = −1/2 to a good precision, confirming that the dispersion relations have branch point singularities of the square root type at the closest to the origin critical points. This is fully consistent with the characteristic quasinormal mode behaviour described in Appendix A.