The complex life of hydrodynamic modes

We study analytic properties of the dispersion relations in classical hydrody- namics by treating them as Puiseux series in complex momentum. The radii of convergence of the series are determined by the critical points of the associated complex spectral curves. For theories that admit a dual gravitational description through holography, the critical points correspond to level-crossings in the quasinormal spectrum of the dual black hole. We illustrate these methods in N = 4 supersymmetric Yang-Mills theory in 3+1 dimensions, in a holographic model with broken translation symmetry in 2+1 dimensions, and in con- formal field theory in 1+1 dimensions. We comment on the pole-skipping phenomenon in thermal correlation functions, and show that it is not specific to energy density correlations.


Introduction
Hydrodynamics is an effective theory of fluids valid at sufficiently long times and sufficiently large distances. Classical hydrodynamics is formulated by combining conservation laws for energy, momentum, and other charges (such as the particle number in non-relativistic systems) with the constitutive relations expressing the conserved fluxes in terms of the JHEP11(2019)097 hydrodynamic variables: local temperature, fluid velocity, and charge density [1]. The constitutive relations are written down on the basis of symmetry, using the derivative expansion. The constitutive relations which contain no derivatives of the hydrodynamic variables are said to describe "perfect fluids", corresponding to "zeroth-order" hydrodynamics. The constitutive relations which contain up to one derivative of the hydrodynamic variables are commonly said to describe "viscous fluids", corresponding to "first-order" or "Navier-Stokes" hydrodynamics. One can proceed by adding terms with more derivatives of the hydrodynamic variables to the constitutive relations, hoping to improve the hydrodynamic description of the actual physical fluids. The constitutive relations which contain up to n derivatives of the hydrodynamic variables then give rise to n th -order hydrodynamics. In this paper, we will use simple facts from the theory of complex curves in order to study some aspects of convergence of the above derivative expansion. Our focus will be on relativistic fluids, and we shall illustrate general results with the examples of strongly interacting quantum field theories possessing a dual holographic description. A neutral homogeneous and isotropic relativistic fluid whose energy-momentum tensor is conserved supports collective excitations in the form of shear and sound hydrodynamic modes [2]. The collective modes arise from the analysis of linearised fluctuations of energy and momentum densities around the equilibrium state at temperature T . The mode's frequency ω is related to the wave vector q by the dispersion relation ω = ω(q). In hydrodynamics, the dispersion relations are given by Shear mode: ω shear (q) = −iDq 2 + · · · , (1.1) Sound mode: where v s is the speed of sound, and D and Γ s can be expressed through relevant transport coefficients. In d s spatial dimensions, we have represent the contributions from first-order hydrodynamics, while the ellipses denote terms with higher powers of q, which can be matched to transport coefficients in second-and higher-order hydrodynamics [3][4][5]. To all orders in the hydrodynamic derivative expansion, the dispersion relations (1.1), (1.2) are represented by infinite series in q. Are these hydrodynamic series convergent? If so, what are their radii of convergence and what determines them? If the series are only asymptotic, can they be resummed? The shear mode (1.1) appears to be an expansion in powers of q 2 , whereas the sound JHEP11(2019)097 mode (1.2) seems to contain odd and even powers of |q| = q 2 in its real and imaginary part, respectively. Can we prove that no other power of q 2 or non-analytic terms appears in the hydrodynamic expansions at any order? 1 These questions were considered in our recent short paper [8], and we shall provide more details here.
In general, proving convergence of a perturbative series and finding the corresponding radius of convergence may constitute rather difficult (and rather different) problems. For example, the convergence of the so called 1/Z series representing the lowest energy eigenvalue of the two-electron atom with the nucleus of charge Z was rigorously proven by Kato in 1951 [9] but reliably computing the actual value of the radius of convergence from the series coefficients and their Padé extensions remained a controversial problem for many decades [10,11]. Yet another example is the series solution to Kepler's equation whose "mysterious" convergence properties were discussed by the giants such as Laplace and Cauchy and were finally understood as arising from the critical points of the associated curve in the complex eccentricity plane (see appendix C).
The problems involving re-summing all-order hydrodynamic expansions and finding the radius of convergence of hydrodynamic series have been previously addressed in the literature. All-order linearised hydrodynamics was investigated by Bu and Lublinsky [12,13] using fluid-gravity correspondence. Re-summations of the hydrodynamic derivative expansion have been also considered in the framework of kinetic theory [14] and in cosmological models [15]. By far the largest body of work on the subject has been done in the setting of the boost-invariant flow, where the gradient expansion in the position space typically produces asymptotic rather than convergent series, and the Borel-Padé and "resurgence" methods [16] can be used to re-sum the series and extract information about gapped modes in the spectrum from the hydrodynamic series [17][18][19][20][21][22][23][24][25][26]. In ref. [27], Withers studied the convergence properties and analytic continuation of a dispersion relation for the sheardiffusion mode in a holographic model involving a dual Reissner-Nordstrom-AdS 4 black brane. There, a finite radius of convergence resulted from a branch point singularity at a certain value of the purely imaginary momentum. From the point of view of the quasinormal spectrum, this point corresponds to the collision of the modes or level-crossing, very similar to the discussion in ref. [8] and in the present paper.
The prediction of classical hydrodynamics is that the above frequencies ω(q) appear as poles of the retarded two-point functions 2 of the energy-momentum tensor in thermal equilibrium, as q → 0 [2,28]. Assuming that the prediction of classical hydrodynamics is correct and that the actual response functions computed from quantum field theory (viewed as functions of ω) indeed have isolated poles at ω = ω(q) with ω(q → 0) → 0, we can define the function ω(q) as the location of the relevant pole. In what follows, we shall discuss several models where the poles of the full response functions can be readily analysed, both analytically and numerically, for generic values of q, either real or complex, 1 In this paper, we consider classical hydrodynamics only, ignoring the effects of statistical fluctuations such as those discussed for example in refs. [6,7]. Such fluctuation effects are suppressed in the holographic models we study below. 2 We shall often call the retarded two-point functions "response functions", as they form the basis of the linear response theory.

JHEP11(2019)097
small or large. This will allow us to study the analytic properties of the derivative expansion in hydrodynamics by studying the dispersion relations ω(q) obtained from the poles of the relevant retarded functions in thermal equilibrium. In general, the hydrodynamic dispersion relations arise as solutions to P (q 2 , ω) = 0 , (1.5) where P is proportional to the inverse of the corresponding retarded two-point function. 3 The hydrodynamic dispersion relations ω(q) are solutions to (1.5) obeying ω(q → 0) → 0, where q 2 and ω are treated as complex variables. We shall refer to P (q 2 , ω) as the hydrodynamic spectral curve. In order to obtain P (q 2 , ω) from classical hydrodynamics, we choose a set of hydrodynamic variables ϕ a (such as the fluid velocity and temperature), and linearise the hydrodynamic equations about the equilibrium state, ϕ a = ϕ (eq) a +δϕ a . In the absence of external sources, the hydrodynamic equations are homogeneous and, upon Fourier transforming, can be written as a set of linear algebraic equations, K ab δϕ b = 0. The hydrodynamic spectral curve is then simply P = det K.
In n th -order (classical) hydrodynamics, P (q 2 , ω) is a polynomial of a finite order, and eq. (1.5) defines a complex algebraic curve. The theorems of analysis such as the implicit function theorem and the theorem of Puiseux then determine the structure and properties of ω(q). In particular, these theorems guarantee that for any finite order of the derivative expansion, the dispersion relations ω(q) are given by series converging in some vicinity of the origin (q 2 , ω) = (0, 0), and the same is true as n → ∞, provided P (q 2 , ω) is an analytic function at (0, 0). We discuss the spectral curve and the associated dispersion relations of the hydrodynamic modes in section 2 of the paper.
In addition to the spectral curve, we shall also study the retarded functions of the energy-momentum tensor. As a simple example, the prediction of 1 st -order hydrodynamics for the retarded function of the transverse momentum density is [2] where D = η/( + p), as before. When viewed as a function of ω, there is a simple pole given by the shear mode dispersion relation (1.1). When viewed as a function of two variables q 2 and ω, one can immediately see that the point (q 2 * , ω * ) = (0, 0) is a singular point of G R (q 2 , ω), such that the value of the response function at (q 2 * , ω * ) is undefined. This is commonly expressed by saying that the limits ω → 0 and q → 0 do not commute. Such indeterminacy-type singularities in physical response functions can also exist at finite non-zero (q 2 * , ω * ). This was recently explored for the sound mode (retarded function of the energy density) in refs. [29][30][31], where the phenomenon of the indeterminacy-type singularities at non-zero (q 2 * , ω * ) was called "pole-skipping". Put differently, pole-skipping happens when P (q 2 * , ω * ) = 0, and the residue of the corresponding pole of G R (q 2 , ω) vanishes at (q 2 * , ω * ). In the example of eq. (1.6), the "skipping" of the shear pole happens at the origin. 3 Here, the proportionality is assumed to be trivial in the sense that the equations G −1 R = 0 and (1.5) are equivalent.

JHEP11(2019)097
A study of pole-skipping at non-zero (q 2 * , ω * ) will be another focus of our paper. 4 In fact, the original motivation for our work was to see whether the pole-skipping singularities in response functions and the non-zero radius of convergence of the hydrodynamic derivative expansion might be related to each other.
To illustrate our approach in a simple exactly solvable example, in section 3 we discuss the holographic bottom-up model studied, in particular, by Davison and Goutéraux in ref. [33]. The advantage of the model is that the effects of translation symmetry breaking on hydrodynamics can be studied in a controlled manner, and that the hydrodynamic and non-hydrodynamic degrees of freedom can be easily separated. The explicit breaking of the translation symmetry means that momentum is no longer conserved, and the sound mode is absent from the spectrum as q → 0. Thus, the only modes with ω(q → 0) = 0 are the diffusive modes of the energy and charge densities. For a certain special value of the translation symmetry breaking parameter in the model, the response functions of the energy and charge densities can be found exactly for all (not just small) momenta [33]. One then finds the following dispersion relation for the diffusive modes: It is clear that the corresponding hydrodynamic series converges in the circle |q| < |q c | = 1/2 due to the branch point singularities at q c = ±1/2. We shall study the exact and approximate spectral curves and obstruction to convergence in this model in section 3.
Our main example, considered in section 4, is the N = 4 supersymmetric SU(N c ) Yang-Mills theory at infinite N c and infinite 't Hooft coupling, which we will abbreviate as "strongly coupled N = 4 SYM theory". In this theory, real-time equilibrium correlation functions can be calculated by the methods of holographic duality [34,35] (see for example refs. [36][37][38][39] for an introduction to the holographic duality and applications). The dispersion relations of the shear and sound modes in the strongly coupled N = 4 SYM theory obtained by holographic methods are shown in figure 1. Using the units = c = 1, we plot the dispersion relations in terms of dimensionless variables w ≡ ω/2πT and q ≡ |q|/2πT . The function w(q) for the shear mode is purely imaginary for real q, while w(q) for the sound mode has both real and imaginary parts for real q. The functions w(q) in figure 1 appear to be smooth and generally unremarkable functions of real positive q. Their behaviour at q 1 has a clear hydrodynamic interpretation [40,41] fully compatible with eqs. (1.1), (1.2), and their asymptotics as q → ∞ were studied in refs. [42,43]. Thus, at least in the case of the N = 4 SYM theory, if the series (1.1), (1.2) have finite radii of convergence, the obstacle to convergence must lie either at negative q = q 2 , or more generally, at complex momenta. By studying complex q, one indeed finds that the functions w(q) in the N = 4 SYM theory have finite non-zero radii of convergence: |q sound  Figure 1. Dispersion relations of the hydrodynamic modes in the strongly coupled N = 4 SYM theory, obtained using the dual holographic description. The dispersion relations are plotted in terms of dimensionless w ≡ ω/2πT and q ≡ |q|/2πT , with complex w as functions of real positive q. The left panel shows w shear (q) for the shear mode, the right panel shows w + sound (q) for one of the two sound modes. In the left panel, the actual −Im w shear (q) for the shear mode is shown by the solid red curve, and the analytic hydrodynamic approximation to O(q 8 ) (computed in section 4.1) is shown by the dashed blue curve. The blue dot indicates the pole-skipping point at q * = 3/2, w * = −i, discussed in section 4.7. The right panel shows Re w + sound (q) (solid red curve) and −Im w + sound (q) (dashed red curve) for the "+" sound mode. The straight dotted line indicates the light cone Re w = q. and demonstrate that the level-crossing phenomenon is also observed in the scalar channel of the correlators.
One of our main results concerns the response functions of the energy-momentum tensor in the strongly coupled N = 4 SYM theory. It was shown in refs. [29,31] that there is a pole-skipping singularity in the retarded two-point function of the energy density at (q 2 * , w * ) = (−3/2, i), at which point the sound pole is "skipped". The sound dispersion curves pass through the point (q 2 * , w * ), as illustrated in figure 2. We observe that in close analogy with the sound mode, the shear mode dispersion relation (1.1) passes through the point q 2 * = 3/2, w * = −i (see figure 1). It turns out that this is not an accident: we will show that the pole-skipping phenomenon in the strongly coupled N = 4 SYM theory is exhibited not only by the response functions which give rise to the sound mode (energy density correlations), but also by the response functions which give rise to the shear mode (transverse momentum density correlations). Moreover, we find that the pole-skipping at non-zero complex momentum also occurs in response functions of those components of the energy-momentum tensor that are not related to hydrodynamic modes. For example, for q along the z direction, the fluctuations of T xy are gapped, and the response function of T xy has no hydrodynamic singularities. Nevertheless, the gapped singularities of the retarded function of T xy pass through (q 2 * , w * ) = (−3/2, −i), at which point one of the gapped poles is "skipped". This is illustrated in figure 3. In other words, all retarded functions of T µν in strongly coupled N = 4 SYM theory exhibit pole-skipping at |q * | = 3/2 and |w * | = 1. The retarded functions of the energy-momentum tensor and the convergence of the derivative expansion in the N = 4 SYM theory are discussed in section 4.
A natural question to ask is whether pole-skipping happens within the domain of validity of the hydrodynamic approximation, as far as the convergence of the hydrodynamic JHEP11(2019)097  Figure 2. The analytically continued sound mode frequencies in the strongly coupled N = 4 SYM theory, obtained using the dual holographic description. The dimensionless frequencies w ± sound of the two sound modes are plotted for purely imaginary dimensionless spatial momentum q, with the "+" branch in red and the "−" branch in blue. The frequencies w ± sound are purely imaginary at imaginary q. At small momenta, the curves are linear with slopes ±v s , with v s = 1/ √ 3. The curves pass through pole-skipping points (q * , w * ) = (±i 3/2, i) indicated by the blue dots. derivative series is concerned. In other words, if the hydrodynamic dispersion relation w i (q) has a finite radius of convergence |q i c | and pole-skipping in the corresponding response function happens at |q i * |, how does |q i c | compare with |q i * |? In the strongly coupled N = 4 SYM theory we have |q sound c | = √ 2, |q shear c | ≈ 1.49 [8], and |q * | = 3/2. Thus |q * | < |q c |, and therefore pole-skipping in correlation functions takes place within the convergence domain of the hydrodynamic derivative expansion. On the other hand, in the model of ref. [33], we have |q c | = 1/2, |q * | = √ 2, hence |q * | > |q c |, so pole-skipping occurs outside the convergence domain of the hydrodynamic derivative expansion. This indicates that the "skipping" of hydrodynamic poles in retarded functions of energy and momentum JHEP11(2019)097 densities is not directly related to the convergence radius of the derivative expansion in hydrodynamics.
More generally, pole-skipping singularities in real-time response functions at non-zero momentum do not have to have any relation to hydrodynamics at all. As an example, we consider response functions of spin-zero operators in 1+1 dimensional conformal field theory (CFT). For a primary operator, the Euclidean correlation function on R 2 is fixed by conformal symmetry. Performing a conformal transformation to the cylinder R × S 1 gives Euclidean thermal correlations functions [44], which can be Fourier transformed and analytically continued to produce exact real-time retarded functions G R (w, q) [34]. These functions have no hydrodynamic poles, yet we will see that there is an infinite number of pole-skippings at non-zero values of (q * , w * ). We discuss this in detail in section 5. Our conclusions and discussion of the issues raised in the paper appear in section 6.

Hydrodynamic spectral curves
We start with a brief review of how the hydrodynamic dispersion relations are derived. Consider hydrodynamics of a neutral homogeneous and isotropic relativistic fluid in flat space in d s spatial dimensions. We are interested in linearised fluctuations in a homogeneous and isotropic equilibrium state, T µν = T µν eq. + δT µν , where T µν denotes the expectation value of the symmetric energy-momentum tensor operator, and the equilibrium state is characterised by T 00 eq. = , T ij eq. = pδ ij , T 0i eq. = 0, where and p are the equilibrium energy density and pressure. The equations of hydrodynamics follow from the conservation of the energy-momentum tensor, ∂ µ T µν = 0. Translation invariance of the equilibrium state implies that we can Fourier transform the fluctuations and take all variables to be proportional to exp (−iωt + iq·x). Furthermore, rotation invariance allows us to choose the direction of the z axis along q. We then have the following system of conservation equations for the linearised fluctuations: −ω δT 0a + q z δT za = 0 , (2.1a) −ω δT 00 + q z δT z0 = 0 , where we use the index a and subsequent Latin indices to denote any of the d s − 1 spatial directions orthogonal to z. The above conservation equations need to be supplemented by the constitutive relations which express δT µν in terms of the hydrodynamic degrees of freedom. For linearised hydrodynamics, a convenient choice of the degrees of freedom is the energy density δT 00 and momentum density δT 0i . This choice implies that we only need the constitutive relations for the spatial stress, δT ij = δT ij (δT 00 , δT 0k ). The constitutive relations will contain derivatives of δT 00 and δT 0k , as is needed for example to describe the viscosity of fluids. We will organise the constitutive relations according to the number of derivatives of the JHEP11(2019)097 hydrodynamic variables. The hydrodynamics of k-th order is determined by the constitutive relations in which δT ij contains up to k derivatives of δT 00 and δT 0i . It is then straightforward to write down the linearised constitutive relations at any order, by noting that under the spatial SO(d s ), the stress fluctuation δT ij is a rank-two tensor, momentum density δT 0i is a vector, and the energy density δT 00 is a scalar. For example, in the first-order hydrodynamics of ref. [1], we have where η is the shear viscosity, ζ is the bulk viscosity, and the ellipses denote terms with more than one derivative of δT 00 , δT 0i . Combining the constitutive relations (2.2) with the conservation equations (2.1) gives a system of linear equations for the fluctuations δT 00 and δT 0i . The equations have non-trivial solutions provided the corresponding determinant vanishes: where v 2 s = ∂p/∂ is the speed of sound squared, and D, Γ s are defined by eqs. (1.3), (1.4). In fact, rotation invariance implies that the most general linearised constitutive relations in momentum space take the following form: δT nm = − iA q n δT 0m + q m δT 0n + δT 00 (Bq n q m + Cδ nm ) + iq l δT 0l (Dq n q m + Eδ nm ) , (2.4) where A, B, C, D, E are scalar functions of ω and q 2 . Substituting the constitutive relations (2.4) into the conservation equations (2.1), we find a system of d s +1 linear equations for d s +1 hydrodynamic variables. This system has non-trivial solutions provided the determinant of the corresponding matrix vanishes. The vanishing of the determinant is equivalent to the vanishing of Here the coefficients are γ η ≡ A, γ s ≡ 2A − E − Dq 2 , H = Bq 2 + C. Thus, the shear and the sound modes decouple as a consequence of rotation invariance. 5

JHEP11(2019)097
If the constitutive relations (2.4) are given by a local derivative expansion, then the functions γ η (q 2 , ω), γ s (q 2 , ω) and H(q 2 , ω) are power series in ω and q 2 with finite values at ω = 0, q 2 = 0: with v s , D and Γ s as above. Truncating the derivative expansion at order k then gives a sequence of algebraic equations defined by finite polynomials in q 2 and ω, For general complex ω and q 2 , eqs. (2.6), (2.7), or (2.9), (2.10) define complex algebraic curves 6,7 which we will call hydrodynamic spectral curves. The complete dispersion relations of the i-th mode, ω i = ω i (q 2 ), can be obtained by solving eqs. (2.6), (2.7) for ω, while the corresponding approximate expressions arising in k-th order hydrodynamics are found from eqs. (2.9), (2.10). The hydrodynamic dispersion relations are the solutions which satisfy ω i (q 2 → 0) = 0. They correspond to infinite relaxation times for infinite-wavelength perturbations of conserved densities, i.e. to the conservation of energy and momentum. Note that the polynomials F sound (q 2 , ω) are not defined uniquely because of the freedom to organise the derivative expansion in hydrodynamics, such as the choice of "frame" and the use of on-shell conditions [2,46]. As an example, the conservation equations (2.1) imply that the factors of ω in the constitutive relations (2.4) can be eliminated at each order in the derivative expansion. Thus one can organise the derivative expansion in such a way that the functions γ η (q 2 , ω), γ s (q 2 , ω) and H(q 2 , ω) are all ω-independent at each given order in the expansion. Then eqs. (2.9), (2.10) give simple explicit expressions for ω shear (q 2 ) and ω sound (q 2 ) in terms of three scalar functions γ η (q 2 ), γ s (q 2 ) and H(q 2 ), whose small-q limits are given by eq. (2.8). In this way of implementing the derivative expansion, the hydrodynamic dispersion relations are the only solutions to (2.9), (2.10). Of course, other choices of organising the derivative expansion are possible where the ω-dependence in γ η (q 2 , ω), γ s (q 2 , ω) and H(q 2 , ω) is retained, and non-hydrodynamic (gapped) modes appear in addition to the hydrodynamic modes.
The above discussion was in the context of classical hydrodynamics. A similar factorisation of the shear and sound modes happens in the full response functions of the energy-momentum tensor, without any hydrodynamic assumptions [41]. For example, for the wave vector along z, the shear mode is described by fluctuations of δT 0a , where the direction a is orthogonal to z. The condition that the inverse of the equilibrium response function of the T 0a operator vanishes can be written as P shear (q 2 , ω) = 0. In general, P shear (q 2 , ω) is a complicated function which describes both hydrodynamic (long-distance, long-time) and non-hydrodynamic (short-distance, short-time) physics. For small (appropriately defined) q 2 and ω, the exact P shear (q 2 , ω) will reduce to the above F shear (q 2 , ω),

JHEP11(2019)097
provided hydrodynamics is a valid effective description of the system. The same applies to the sound mode: the condition that the inverse of the equilibrium response function of the T 00 operator vanishes can be written as P sound (q 2 , ω) = 0. For small (relative to the appropriately defined scale) q 2 and ω, the exact P sound (q 2 , ω) will reduce to the above F sound (q 2 , ω), provided hydrodynamics is a valid effective description of the system. In this paper, we will only be studying physical systems in which the near-equilibrium physics governed by the conserved densities can be described by classical hydrodynamics. In other words, we will assume that the functions γ η (q 2 , ω), γ s (q 2 , ω), H(q 2 , ω) are defined by the exact response functions of the T µν operator.

Small-momentum expansions
For a given spectral curve, the small-q 2 expansion of the hydrodynamic dispersion relation ω i (q 2 ) can be found using the theorem of Puiseux (see appendix A and refs. [45,47,48]).
Starting with the shear mode, let us assume that γ η (q 2 , ω) is analytic at (q 2 , ω)=(0, 0) and thus eq. (2.6) defines an analytic curve for complex (q 2 , ω) ∈ C 2 . This is of course not guaranteed a priori and should be established by independent methods, e.g. by finding the exact response function. 8 The analyticity of γ η (q 2 , ω) implies that F shear (q 2 , ω) is analytic at the origin as well, and, since ∂F shear /∂ω = 1 = 0 at (0, 0), the origin is a regular point of the analytic curve (2.6). Then the analytic implicit function theorem (see appendix A) guarantees that for sufficiently small q 2 and ω, there exists a unique solution of eq. (2.6) of the form convergent in a neighbourhood of q 2 = 0. The radius of convergence is determined by the location of the nearest to the origin critical point of the curve (2.6).
Continuing with the sound mode, let us again assume that γ s (q 2 , ω) and H(q 2 , ω) are analytic functions at (0, 0). The function F sound (q 2 , ω) is then analytic at the origin as well. Now we have ∂F sound /∂ω = 0 at (0, 0), and thus the origin is a critical point of the spectral curve. On the other hand, ∂ 2 F sound /∂ω 2 = 2 = 0 at (0, 0), thus we expect the sound dispersion relation to have p = 2 branches. The Puiseux series expansions are then given by eqs. (A.6), (A.7), in other words ω (j) sound (q 2 ) can be represented by series in non-negative powers of (q 2 ) 1/m j , where m j are positive integers, and j = 1, 2 labels the two branches corresponding to the two sound modes. Following the general analysis of algebraic curves, the integers m j may be found using the Newton's polygon method (see refs. [45,47,48] for details). For analytic γ s and H, we have the expansions H nm ω n q 2m . The coefficients in front of various powers of ω in the expression (2.7) are then given by H 0k q 2k+2 , (2.14) The vertices of the Newton polygon are thus given by (0, 1+k 0 ), (1, 1+k 1 ), (2, 0), (3, 1+k 3 ), (4, 1 + k 4 ), . . ., where k 0 , k 1 , . . . are the smallest indices such that H 0k 0 = 0 as well as H 0k 1 = 0 or/and γ s 0k 1 = 0, etc. The Newton polygon for the sound mode is shown in figure 4, where it is assumed that H 00 = 0 and that either H n,0 = 0 or γ s n−1,0 = 0 (or both) are non-zero for n = 3, 4, . . .. The exponents of x ≡ q 2 in the Puiseux series are given by the negative slopes of the polygon's lines, i.e. by 1/2 for H 00 = 0. Thus m j = 2, and the lowest order term in the two branches is (2.20) From the Newton polygon, it is clear that H 00 = 0 is the necessary and sufficient condition for the fractional powers of q 2 to appear in the dispersion relation. In the hydrodynamic derivative expansion, H 00 = ∂p/∂ = v 2 s is the speed of sound squared. One may have H 00 = 0 at the point of a phase transition (see e.g. ref. [50]) in which case the dispersion relation contains only positive powers of q 2 .
Generically, for H 00 = v 2 s = 0, the sound mode dispersion relation will be given by the two branches of Puiseux series in (q 2 ) 1/2 converging in some neighbourhood of the point q 2 = 0, where a n ∈ R and a 1 = c s is the speed of sound. In particular, for q 2 = e ±iπ |q 2 |, the functions ω ± sound (q 2 ) are purely imaginary as anticipated in ref. [29].

Convergence of the hydrodynamic series
The shear and sound dispersion relation series (2.11) and (2.21) have non-zero radii of convergence as long as the corresponding spectral curves (2.6) and (2.7) are given by the functions of q 2 and ω analytic at the origin (0, 0). One way to find the radius of convergence of the series is to analyse the behaviour of the coefficients a n , c n at large n. This behaviour will of course depend on the microscopic details of the particular physical system, and may be difficult to study in practice (see, however, ref. [27]). Instead, here we use the spectral curves to determine the radii of convergence. The Puiseux analysis implies that the domain of convergence of Puiseux series centred e.g. at the origin is the circle whose radius is set by the distance from the origin to the nearest critical point of the associated spectral curve. Critical points of the spectral curve F (q 2 , ω) = 0, where q 2 , ω ∈ C, are determined by the conditions There are p > 1 branches of the solution ω = ω(q 2 ) in the vicinity of the critical point, provided that For example, the origin (0, 0) is the critical point (with p = 2) of the sound hydrodynamic spectral curve (2.7), as discussed in section 2.2. If the spectral curves happen to be known exactly or approximately (as will be the case in the holographic models we study below), eqs. (2.22) provide an efficient method to find the radii of convergence, without performing the large-n analysis of the expansion coefficients. When the function F (q 2 , ω) is a polynomial, the condition (2.22) means that the equation F (q 2 , ω) = 0, where F (q 2 , ω) is considered as a polynomial in ω with q 2 -dependent coefficients, has multiple roots at ω = ω c . This is equivalent to the condition that the discriminant of the polynomial F (q 2 , ω) vanishes. As an example, consider the first-order hydrodynamics of [1], where the spectral curves following from eq. (2.3) are 9 As mentioned in section 2.1, the expressions for the truncated spectral curves F (k) at each order k of the hydrodynamic derivative expansion are not unique due to the freedom allowed by the "frame" choice. Correspondingly, the critical points determined by the approximate spectral curves F (k) depend on the "frame" choice. This dependence becomes less and less pronounced with k increasing and disappears in the limit k → ∞. Thus, although the critical points of the exact spectral curve are "frame"-independent, the rate of convergence of the approximate location of the critical points to the exact values can be affected by the "frame" choice.

JHEP11(2019)097
These equations are simple enough to be solved explicitly: eq. (2.24) is solved by ω = −iDq 2 , whereas the solutions of eq. (2.25) are where in the expansion we only kept terms quadratic in q ≡ q 2 , since we expect the coefficients in front of the higher powers of q to be modified by higher-derivative terms in the hydrodynamic expansion. The series in q on the right-hand side of eq. (2.26) has the radius of convergence which coincides with (2.27). In what follows, we will be studying models where the convergence radii of the small-q expansions can be determined from the exact response functions of the theory, without resorting to the derivative expansion of hydrodynamics.

A holographic model with translation symmetry breaking
To illustrate the methods discussed in section 2.3, we consider the holographic model with translation symmetry breaking [52], studied, in particular, in ref. [33]. The model is a bottom-up gravity construction in 4d describing a hypothetical dual 2+1-dimensional QFT with broken translational invariance. In the context of the present paper, the significance of the construction discussed in ref. [33] lies in the fact that it provides exact analytic formulae for the current and energy-momentum correlator two-point functions at a special, self-dual symmetry point (see section 4 of ref. [33] for details). In particular, among the poles of the correlation function, one finds a gapless excitation whose dispersion relation is known analytically, possibly a unique example in holography.
The bulk action of the model is given by [33,52] (see also ref. [31]) where Λ = −3/L 2 (we set L = 1 in the following). The background solution of interest involves the AdS-Schwarzschild black brane with translational invariance broken by the linear dilaton fields, and a vanishing Maxwell field.

JHEP11(2019)097
First, we recast the setup of ref. [33] in the form more convenient for our purposes, having in mind the variables used in ref. [41]. The metric at the special symmetry point has the form (see section 4 of ref. [33]) where f = 1 − u, and we have used the coordinate u = r 2 0 /r 2 . The horizon is located at u = 1 and the boundary at u = 0. The Hawking temperature of the background is T = r 0 /2π. The dilaton fields are given by φ 1 = √ 2r 0 x and φ 2 = √ 2r 0 y. They will not play any role in the following.
Considering the Maxwell field fluctuations in the background (3.2), coupled to the current operator on the boundary, one finds the equations of motion , the equation of motion reads The equation (3.6) has 4 regular singular points (located at u = 0, 1, 1 − w 2 /q 2 , ∞) and thus is of the Heun type. The indices of this equation are, respectively, Note that the local solution at u = 0 does not contain logarithms. The exact solution to eq. (3.6) can be written as where G ≡ a t is the solution of the hypergeometric equation 11

The exact spectral curve
For the boundary value of the electric field, we find where The condition E x (u = 0) = 0 determines the quasinormal modes and thus the poles of the corresponding current-current correlators [34,41]. One can also write down the explicit analytic expressions for the correlators themselves [33], but this will not be necessary: the expression F (q 2 , w) = 0, where F (q 2 , w) is given by eq. (3.14), is the exact spectral curve containing all information about the poles of the two-point function. There are two sequences of quasinormal frequencies The solutions E ± x,n ± (u) themselves (the quasinormal modes) have the form where P ± n ± (u, q 2 ) are polynomials of degree n ± in u with q 2 −dependent coefficients. In particular, Note that apart from the prefactors determined by the indices, the solutions (3.17) are polynomials, with P ± 0 = 1. We are especially interested in the gapless mode The power series in the second line of eq. (3.19) converges in the circle |q 2 | < 1/4, due to the branch point singularity of the function at q 2 = 1/4 evident from eq. (3.19). The same conclusion can be obtained by analysing critical points of the spectral curve (3.14).

JHEP11(2019)097
Indeed, the critical points are determined by the equations (2.22) whose solutions are It is also clear that ∂ 2 w F (q 2 c , w c ) = 0 and thus there are two branches of the spectral curve emerging at each critical point. Put simply, the critical points occur when the two quasinormal frequencies in eq. (3.16) collide. This happens for real q in case of n + = n − , and for purely imaginary q for n + = n − . The critical point closest to the origin q 2 = 0 is at q 2 = 1/4 (it corresponds to the collision of the modes w + 0 and w − 0 ). This value determines the radius of convergence of the series in eq. (3.19).
Exact spectral curves are rare: in addition to (3.14), we are aware of only one example (involving the exact R-current two-point correlators in N = 4 SYM theory in appropriate limit [53]) for a QFT in dimension higher than 1 + 1, and even in that case it is only known exactly for q 2 = 0: (3.21) One can use the Weierstrass decomposition where γ is the Euler-Mascheroni constant, to write eqs. (3.14), (3.21) as infinite products: this ilustrates explicitly why the critical points determined by the condition (2.22) correspond to multiple roots of (infinite order) polynomials.

The hydrodynamic approximation to the spectral curve
The expansion of the expression wF (q 2 , w) for small w, q (assuming the scaling w → w, and so on. Using this expansion to solve the equation wF (q 2 , w) = 0 for w perturbatively in q 2 , one reproduces the series in eq. (3.19). The equation F k (q 2 , w) = 0 defines the hydrodynamic spectral curve of order k as discussed in section 2. At each order, the critical points are determined by the condition (2.22). In figure 5, we plot the corresponding value |q 2 c | for each of the spectral curves F k (q 2 , w) = 0 for k = 2, 3, . . . , 13. The resulting points converge rapidly to the exact value |q 2 c | = 1/4.

Pole-skipping in the full response functions
Finally, we comment on the relationship between the critical point defining the hydrodynamic series radius of convergence and the pole-skipping point in the holographic model with broken translation symmetry, which occurs at |q 2 * | = 2 > |q 2 c | [31]. Because gapless excitations in the current and the energy density correlators have the same dispersion relation, we can directly use w(q) from eq. (3.19) to discuss both the charge and the energy sectors. At the pole-skipping point, the hydrodynamic series stated in the second line of eq. (3.19) diverges but can be resummed using the Borel transform 12 where I n (x) is the modified Bessel function. Since this is a series in q 2 and not q, the corresponding Borel integral has the form where

JHEP11(2019)097
each with their respective region of q for which the integral converges. Together, we find that the Borel integral representation (3.28) of the series converges for q ∈ C in the region defined by the function C(q): This is a significant improvement over the convergence region of the hydrodynamic power series, |q 2 | < |q 2 c | = 1/4. In particular, the Borel series is well suited for studying purely imaginary q. By writing q = i , with ∈ R, we see that Hence, the Borel representation of the series converges for all values of imaginary q, including w(q) at the chaos point w * (q * = √ 2i) = i. Finally, for q such that C(q) < 1/4, it is easy to check that the sum of four terms in (3.28) indeed reproduces the full dispersion relation (3.19).

Response functions in strongly coupled N = SYM theory
In this section, we use holographic duality to find the spectral curves, determine the radii of convergence of hydrodynamic series and analyse the pole-skipping phenomenon in the three channels of the response function of the energy-momentum tensor in the N = 4 SU(N c ) SYM theory at infinite 't Hooft coupling and infinite N c . The details of the duality are well known, and the necessary ingredients can be found e.g. in refs. [35,41,55]. In short, holography reduces the study of the response functions to the analysis of the fluctuations of the dual gravitational background involving a black hole with translationally invariant horizon -the AdS-Schwarzschild black brane.
The equations of motion describing fluctuations of the gravitational background dual to finite-temperature N = 4 SYM theory are of the Heun type [56], and the exact analytic solution for the spectral curve similar to eq. (3.14) is not available. The equations can be solved perturbatively and analytically in w 1, q 1, as was done in refs. [35,40], thereby giving a hydrodynamic approximation to the spectral curve, or numerically, for arbitrary w and q, along the lines of ref. [41]. In this section, we consider and compare these two approaches.

Shear mode: hydrodynamic approximation to the spectral curve
We start with the analysis of the shear channel of the energy-momentum tensor response function [41]. In the hydrodynamic approximation, the spectral curve can be found analytically from the boundary value of the solution to the ODE obeyed by one of the components of the shear perturbation in dual gravity (see section 6.2 of ref. [35] for details):

JHEP11(2019)097
where G is regular at u = 1. Rescaling w → λ 2 w and q 2 → λ 2 q 2 , assuming λ 1, and looking for a perturbative solution of the form we find the following equation for coupled components G n , G n−1 , G n−2 : where G n = 0 for n < 0, and we can set G 0 (1) = 1, G i (1) = 0, i ≥ 1, without loss of generality. 13 The explicit formulae for G 0 (u), G 1 (u) and G 2 (u) obeying the boundary conditions at u = 1 are written in appendix B. The solution of the homogeneous equation is given by g = C 1 g 1 (u) + C 2 g 2 (u), where Note that the Wronskian is W (g 1 , g 2 ) = 1/f . Then one can write the following expression for G n , n ≥ 2: where Boundary conditions at u = 1 (regularity of G n (u) at u = 1 and G n (1) = 0 for n > 0) require C 1 = 0, C 2 = 0. Thus, we have the equation determining G n from G n−1 and G n−2 , 13 See appendix C of ref. [57].

JHEP11(2019)097
where F n is given by eq. (4.8). This can be written as Note that F n (t) ∼ 1/(1 − t) at t → 1, and so all integrals converge. In particular, The explicit expressions for G i (0), with i = 1, 2, 3, 4, are written in appendix B. The results for G 0 (0), G 1 (0) and G 2 (0) coincide with those in ref. [3]. The results for G 3 (0) and G 4 (0) are new. The condition G(0) = 0 at this order in the hydrodynamic expansion defines the algebraic curve Here, ζ(z) is the Riemann zeta function. Solving eq. (4.13) perturbatively in q 2 , one finds the dispersion relation for the shear mode (4.14) The first two terms in (4.14) agree with the ones obtained in refs. [35] and [3] correspondingly, we obtain the algebraic curves F k (q 2 , w) = 0 for k = 1, 2, 3, 4. Using the condition (2.22), we find that there are no solutions at k = 1, whereas for k = 2, 3, 4, we have In figure 6, the k = 2, 3, 4 approximations are shown in the complex plane of q 2 c together with the "exact" value q 2 c ≈ 1.8906469 ± 1.1711505i obtained via the quasinormal levelcrossing method (see section 4.2). In contrast with the holographic model with broken translation symmetry (see figure 5), the convergence is slow (admittedly, we only have 3 points in the present case). However, we learn an important lesson, namely, that the critical point can be located at a generic complex value of q 2 . In the next section, we use a more efficient method of the quasinormal modes level-crossing to determine the critical points.

Shear mode: full spectral curve
To compute the spectral curve numerically for all values of q 2 and w, it is more convenient to use the ODE obeyed by the gauge-invariant perturbations. The gauge-invariant shear mode gravitational perturbations of the AdS-Schwarzschild black brane are described by the function Z 1 (u), where u is the radial coordinate ranging from u = 0 (asymptotic boundary) to u = 1 (event horizon) [41]. The function Z 1 (u) obeys the equation

JHEP11(2019)097
where f (u) = 1 − u 2 , and the solution must obey the incoming wave boundary condition at the horizon, Z 1 (u) ∼ (1 − u) −iw/2 as u → 1 [34,41]. The full shear spectral curve is then given by The spectral curve (4.19) can be determined numerically by e.g. using N terms of the Frobenius series solution at u = 1 [41]. The stability of the numerical procedure is ensured by checking that adding several more terms to the series does not appreciably change the result. Applying (numerically) the criterion (2.22) to (4.19), we find the critical point closest to the origin, In figures 7 and 8, we show the corresponding quasinormal spectrum (solutions w = w(q 2 ) of eq. (4.19)) in the complex plane of frequency w, parametrised by q 2 . The difference with previous works is that now, we treat q 2 as a generic complex variable, q 2 = |q 2 |e iθ , and vary its phase, θ ∈ [0, 2π], at specific fixed values of the magnitude |q 2 |. From eq. (4.18), it is clear that if Z 1 (u; q, w) is a solution satisfying the incoming wave boundary condition at the horizon, then Z 1 (u; q, −w) is also a solution. The spectrum in figures 7 and 8 therefore appears to be symmetric with respect to the imaginary axis. For real q 2 , the spectrum coincides with the one originally found in ref. [58].
At small |q 2 |, the poles follow closed orbits as the phase θ varies from 0 to 2π (figure 7, top panels). With the parameter |q 2 | increasing, the poles start feeling the presence of each other, and their orbits become more complicated. Finally, at |q 2 c | ≈ 2.224, the hydrodynamic shear pole collides with one of the two closest non-hydrodynamic poles (see figure 7, bottom panels). For |q 2 | > |q 2 c |, those three poles no longer have closed orbits: as the phase θ varies from 0 to 2π, they interchange their positions cyclically (figure 7, bottom right panel). For even larger |q 2 |, other gapped poles become involved in this collective motion of quasinormal modes in a similar manner (figure 8).

Sound mode: hydrodynamic approximation to the spectral curve
The gauge-invariant perturbation corresponding to the sound mode obeys the equation [41] The full sound spectral curve is constructed from the solution Z 2 (u; q 2 , w) obeying the incoming wave boundary conditions at the horizon and is given by   Figure 9. Quasinormal spectrum level-crossing: the real (blue curves) and the imaginary (red curves) parts of the hydrodynamic shear mode and the closest gapped quasinormal mode dispersion relations plotted as functions of |q 2 | at the fixed phase θ ≈ 0.338858π of the complex momentum q 2 = |q 2 |e iθ . At |q 2 | = |q 2 c | ≈ 2.224, the level-crossing occurs.

JHEP11(2019)097
Similarly to the shear mode case, one can first find a hydrodynamic approximation to eq. (4.25) analytically. For w 1 and q 1, eq. (4.24) can be solved perturbatively. Writing to enforce the boundary condition at the horizon, rescaling w → λw and q → λq, and assuming λ 1, we find the following recurrence relation for the functions z n (u): with z −1 = z −2 = 0. To fourth order in λ, we find the hydrodynamic algebraic curve to be (see also eq. (4.18) in ref. [3]) The form of F (q 2 , w) for general w and q is at present not known to O(λ 5 ). However, assuming that w can be expanded in a series in powers of q, in ref. [5], F (q 2 ) was computed to order O(λ 5 ) = O(q 5 ), which was sufficient to find the coefficient in the sound mode dispersion relation multiplying q 4 (i.e. to third order in the gradient expansion): The first two terms in the dispersion relation (4.29) coincide with those obtained in ref. [40]. The third term was found in ref. [3] (see also ref. [69]). The fourth term was computed in ref. [5]. It appeared earlier in ref. [13] (with an incorrect coefficient in front of ln 2) and (correctly) in ref. [70] in the context of the fluid-gravity correspondence.
Here, we use the form of the algebraic curve (4.28) and apply the condition (2.22) to show that in the sound channel, the small w and q expansion of the spectral curve F (q 2 , w) qualitatively correctly accounts for two sets of critical points:

Sound mode: full spectral curve
As discussed in section 4.3, the origin (4.30) is a critical point of the sound mode dispersion relation w = w(q 2 ), as predicted by hydrodynamics (see section 2.2). Proceeding as in section 4.2, we find the first set of critical points nearest to the origin at within the limits of our numerical accuracy. Curiously, although eq. (4.24) looks rather complicated, one can check that with w and q 2 given by (4.32), a simple analytic solution satisfying the correct boundary conditions at u = 1 and u = 0 is available. Explicitly, the JHEP11(2019)097 Figure 11. Quasinormal spectrum (poles of the retarded energy-momentum tensor two-point function in the N = 4 SYM theory) in the sound channel plotted in the complex w-plane for different values of the complex momentum q 2 = |q 2 |e iθ . Large dots in all four plots correspond to the location of the poles for purely real momentum, q 2 (i.e. at θ = 0) [41]. What is shown is the level-crossing phenomenon for the critical points in eq. (4.35).
two solutions corresponding to the pair of critical points in (4.32) are where C ± 2 are arbitrary constants. Although we were not able to show analytically that ∂ w Z(u = 0) = 0 at (4.32) as well, we have verified this numerically to high precision. The existence of the critical point (4.32) implies that the convergence radius of the sound mode dispersion relation is given by |q c sound | = √ 2 ≈ 1.41421. The next set of critical points is located at The behaviour of poles in the complex frequency plane is shown in figures 10 and 11, and the quasinormal level-crossing phenomenon is presented in figure 12.

Scalar mode: full spectral curve
In the scalar channel, the gauge-invariant metric perturbation Z 3 (u) obeys the equation [41] The full spectral curve is constructed from the solution Z 3 (u; q 2 , w) obeying the incoming wave boundary conditions at the horizon and is given by The quasinormal spectrum in the scalar channel has no hydrodynamic modes, but it exhibits the phenomenon of level-crossing for the gapped modes, as shown in figure 13. The  Figure 12. Quasinormal spectrum level-crossing in the sound channel: the real (blue curves) and the imaginary (red curves) parts of the hydrodynamic sound mode and the closest gapped quasinormal mode dispersion relations plotted as functions of |q 2 | at the fixed phase θ = π/2 of the complex momentum q 2 = |q 2 |e iθ . At |q 2 | = |q 2 c | = 2, the level-crossing occurs. Figure 13. Quasinormal spectrum (poles of the retarded energy-momentum tensor two-point function in the N = 4 SYM theory) in the scalar channel plotted in the complex w-plane for different values of the complex momentum q 2 = |q 2 |e iθ . Large dots in all plots correspond to the location of the poles for purely real momentum, q 2 (i.e. at θ = 0) [41]. There are no hydrodynamic modes in this channel, but the gapped modes exhibit the level-crossing phenomena at complex values of momenta given by eq. (4.38).
first two sets of critical points nearest to the origin are given by These are branch points in the scalar (gapped) dispersion relation w = w(q 2 ).

Analytic structure of the hydrodynamic dispersion relations
From the discussion above, it is clear that the shear mode dispersion relation w shear (q 2 ) is an analytic function of complex q 2 in the circle |q 2 | < |q 2 c | ≈ 2.224. Since the appropriate second derivative of the spectral curve at the critical point is non-zero (i.e. p = 2 in the Puiseux language of section 2.3; see eq. (2.23)) which corresponds to a collision of two quasinormal modes, the critical point is the branch point singularity of the square root type, with the Puiseux series in powers of ± q 2 − q 2 c providing the extension beyond the radius of convergence. We show the critical points, the radius of convergence and the appropriate branch cuts in the complex plane of q 2 in figure 14 (left panel).
For the sound mode dispersion relation, considered as a function of q 2 ∈ C, the origin is a branch point, and the corresponding Puiseux series is given by eq. (2.21). It will be more convenient to consider the sound dispersion relation w sound (q) as a function of complexified magnitude q ∈ C of the wave-vector q. The critical points, the radius of convergence and the appropriate branch cuts in the complex plane of q are shown in figure 14 (right panel).

Pole-skipping in the full response functions
As already discussed in the Introduction, analytically continued hydrodynamic modes appear to be connected to the parameters of an OTOC related to the microscopic manybody quantum chaos. The apparent connection is provided by the phenomenon of poleskipping [29][30][31][32], whereby a pole and a zero of a two-point correlation function collide for some w, q ∈ C. In the sound channel of the energy-momentum tensor retarded two-point function, the pole-skipping has been studied in the context of holography [29,31,32], effective field theory [30] and two-dimensional conformal field theory in the limit of large central JHEP11(2019)097 charge [71]. Here, we extend the discussion in refs. [29,31,32] to show that pole-skipping also occurs in other channels.
Consider a retarded two-point function G R (w, q) of the energy-momentum tensor components at finite temperature. Schematically, and modulo tensor structure, the correlator can be written as where a(w, q) necessarily contains a gapless hydrodynamic mode w = w(q) (either shear or sound) in the appropriate channels [41]. More generally, let Z d = {w = w(q) : a(w(q), q) = 0} and Z n = {w = w(q) : b(w(q), q) = 0}, where we assume for simplicity that all zeros are simple. Then, pole-skipping occurs at generically complex (q * , w * ) ∈ P = Z n ∩ Z d .
For theories with available gravity dual descriptions, the set P can be determined directly either by computing Z n and Z d (the set Z d is nothing but the quasinormal spectrum) or from the dual gravity equations of motion (see below). In the case of energy-momentum tensor correlators of the N = 4 SYM theory in the limit of infinite N c and infinite 't Hooft coupling λ, pole-skipping in the three channels occurs at points (q * , w * ) given by Sound channel : Shear channel : Scalar channel : We observe that |q * | = 3/2, |w * | = 1 in all three channels. The connection to the Lyapunov exponent λ L and the butterfly velocity v B is given by the formulae Sound channel : Shear channel : Scalar channel : where * = 3/2, O * = 1, and v B = O * / * . It is clear from figures 2 and 1 that the sound and the shear dispersion relations pass through their respective "chaos" points (4.40) or (4.41). In the scalar channel, which has no hydrodynamic modes, pole-skipping is exhibited by (one of the pair of) the lowest-lying gapped modes in the spectrum. This can be seen from figure 3. Thus, in the N = 4 SYM theory at infinite 't Hooft coupling, the values of λ L and v B defined by pole-skipping in the complexified dispersion relations of the lowest-lying modes (either hydrodynamic or gapped) coincide with those obtained from the appropriate limit of the OTOC: 14 (4.46) 14 To subleading order in the inverse 't Hooft coupling expansion, the butterfly velocity is vB = 2/3 1 + 23ζ(3) 16 λ −3/2 while the relevant (long-distance) Lyapunov exponent remains uncorrected [32].

JHEP11(2019)097
In fact, irrespectively of the channel in question, we can define the (maximal holographic) Lyapunov exponent and the butterfly velocity through the pole-skipping location exhibited by the mode closest to the origin in the complex plane of w: Pole-skipping points (q * , w * ) can be found directly from the dual gravity equations of motion [29,31]. To show this explicitly for the N = 4 SYM theory, we follow the arguments of ref. [31] and examine the horizon behaviour of Einstein's equations in the infalling Eddington-Finkelstein (EF) coordinates (v, r, . (4.49) We perturb the 5d AdS-Schwarzschild metric ds 2 = g µν dx µ dx ν = −r 2 f (r)dv 2 + 2dvdr + r 2 d x 2 to first order, g µν → g µν + δg µν (r)e −iωv+ikz , and expand the (regular) metric fluctuations around the horizon r = r 0 as Similarly to what was observed in ref. [31] for the sound channel, we find that in any channel, there exists a linear combination of the components E µν , which vanishes identically at the horizon r → r 0 at the pole-skipping values of the parameters (q * , w * ). Explicitly, Sound channel : lim Shear channel : lim In other words, pole-skipping occurs at values of the parameters (q * , w * ) for which the rank of the matrix E µν decreases at the horizon. We note also that for the N = 4 SYM theory, the chaos point |q 2 * | = 3/2 lies within the radius of convergence of the hydrodynamic series in both the shear (|q 2 c | ≈ 2.224) and the sound (|q 2 c | = 2) channels.

Pole-skipping and level-crossing in 2d thermal CFT correlators
In a 2d CFT, the (equilibrium) retarded finite-temperature two-point correlation function of an operator of non-integer scaling dimension ∆ and spin zero in momentum space is JHEP11(2019)097 given by the expression 15,16 [34] cosh (πq)−cos (π∆) cosh (πw)+i sin (π∆) sinh (πw) , where C ∆ is the normalisation constant, and we put T L = T R = T . Note also that here, in 1 + 1 dimensions, the symbol q denotes q/2πT , rather than |q|/2πT . Similar formulae can be written for integer ∆ [34]. The correlator (5.1) has a sequence of poles at where n = 0, 1, 2, . . .. These are precisely the quasinormal frequencies of the dual BTZ black hole [34,72]. In this section, we shall examine these correlators for their pole-skipping and level-crossing properties. 17,18

Pole-skipping in the full response functions
The zeros of the correlator (5.1) come from the zeros of the expression in the square brackets, cosh (πq) − cos (π∆) cosh (πw) + i sin (π∆) sinh (πw) and are given by where n 1 , n 2 = 1, 2, 3, . . .. Note that the zeros of eq. (5.3) with n 1 , n 2 = 0, −1, −2, . . . are not among the zeros of the correlator since they are identically (i.e. for arbitrary w, q) cancelled by the poles of the first two Gamma-functions in eq. (5.1). 15 In ref. [34], the expression for G R (ω, q) was derived holographically from dual gravity. For integer ∆, it was further checked in ref. [34] that thus obtained formula coincides (up to normalisation) with the retarded correlator obtained from 2d CFT. 16 The expression for G R (ω, q) in the form (4.16) of ref. [34] assumes w, q ∈ R. To be valid for generic w, q ∈ C, it has to be rewritten in the form (5.1). We would like to thank D. Vegh for pointing this out. 17 Similar issues have been recently independently studied in ref. [73]. The results of ref. [73] agree with ours whenever they overlap. 18 Here, we only consider correlation functions of 2d CFT operators with scaling dimension ∆ and spin s = 0. The energy-momentum, having ∆ = 2 and s = 2, is not of this type. Its finite-temperature twopoint function (see e.g. refs. [71,74]) in momentum space has a pole corresponding to a mode propagating on the light cone. The corresponding dispersion relation line passes through the "chaos" point of that correaltor [71], just as it does in 4d.

(5.11)
Therefore, the pole-skipping points for ∆ = 2 are given by q * = ±i(n + 1) , (5.12) where n = 0, 1, 2, . . .. In particular, for the pair of poles that lies closest to the origin in the complex w plane (ones with n = 0 among those in eq. (5.10)), the branch w + 0 passes through the (lowest-lying) pole-skipping point q * = i, w * = −i, whereas the branch w − 0 does not. The branch w − 0 passes through the pole-skipping point q * = −i, w * = −i, whereas the branch w + 0 does not (see figure 15). Finally, we note that for ∆ = 1, the correlator is directly proportional to the sum of two ψ−functions [34], and there is no pole-skipping. 19 Introducing Q = n * and N = n+n * , the pole-skipping condition can be written as q * = ±i(∆+N −2Q), w * = −iN , where N = 1, 2, . . . and Q = 1, . . . , N . This coincides with the results in ref. [73]. We would like to thank R. Davison for pointing out an error in eqs. (5.6), (5.7) in the first version of this paper. Figure 15. Pole-skipping and level-crossing points in a 2d CFT correlator of operators with conformal dimension ∆ = 2, for the smallest in magnitude poles w ± 0 and w ± 1 . The red stars at Im q = ±1 correspond to the pole-skipping points and the red dots label the critical points of levelcrossings of w ± 0 . The blue stars at Im q = ±2 correspond to the pole-skipping points and the blue dots label the critical points of level-crossings of w ± 1 .

BTZ spectrum level-crossing
Since the quasinormal spectrum w ± n (q) is known explicitly (see eq. (5.2)), the level-crossing points can be found directly. Such level-crossing points were used above in theories with gapless excitations to determine the radius of convergence of their hydrodynamic series. Considering complex q = |q|e iθ , we have Re w ± n = ±|q| cos θ ≡ X , (5.15) and thus the orbits followed by the poles in the complex w plane when the phase θ changes from 0 to 2π are circles X 2 + (Y + 2n + ∆) 2 = |q| 2 , n = 0, 1, 2, . . . . (5.17) The poles move counter-clockwise and collide on the imaginary axis of w (moreover, at integer values of |w| if ∆ is an integer). More precisely, two poles collide when w − n = w + m , m = n, i.e. when (the case n = 0 or m = 0 is treated separately below) The motion of poles in the complex frequency plane and their level-crossings are illustrated for ∆ = 1/8 and |q| = 1 in figure 16. For larger |q|, the trajectories intersect but the poles miss each other, so there is no phenomenon of one trajectory crossing into and continuing as the other. In a sense, here, we have "level-touching" rather than "level-crossing". The nearest critical points are thus In particular, for ∆ = 2, we have 20

JHEP11(2019)097
We observe that for non-integer ∆, the values of q corresponding to the pole-skipping and the level-crossing do not coincide. The same is (trivially) true for ∆ = 1 as well, where there is no pole-skipping in the correlator at all. For ∆ = 2, however, a curious picture emerges. Consider again the "sound" mode w ± 0 (5.14). The imaginary parts of the two branches obey The pole-skipping points found in section 5.1 are q * = i, w * = −i and q * = −i, w * = −i, and the critical points are q c = ±in, w c = −i(n + 2), n = 1, 2, . . .. For ∆ = 2, and the mode with n = 1,

Discussion
In this paper, we introduced spectral curves as a useful tool for investigating analytic properties of gapless collective excitations in classical hydrodynamics. 21 We showed that the dispersion relations of hydrodynamic modes, such as shear and sound modes, are generically given by Puiseux series expansions in rational powers of the spatial momentum squared. These series are guaranteed to converge in some vicinity of the origin (the point with zero frequency ω = 0 and zero spatial momentum q 2 = 0 in the (ω, q 2 ) ∈ C 2 space), so long as the analyticity of the spectral curve at the origin can be proven (e.g. by holographic or other means). Thus, given the analyticity of the spectral curve, the asymptotic nature of the series for hydrodynamic modes in momentum space can be automatically ruled out. The radius of convergence of the series is given by the distance from the origin to the critical point of the spectral curve nearest to the origin. After developing the general theory, we then used holography as a theoretical laboratory where all these of features can be seen and analysed explicitly. Before focusing on the main example of the strongly coupled N = 4 supersymmetric Yang-Mills theory at finite temperature, to illustrate our method, we first considered the holographic model with broken translation symmetry, where an exact spectral curve is available. We have shown that the critical points of the spectral curves can be found by studying quasinormal spectra at complex values of spatial momentum: the critical points correspond to the collisions of quasinormal frequencies (poles of dual correlators) in the complex frequency plane at critical (generically, complex) values of spatial momentum. These values also set the radii of convergence for the dispersion relation. We call this phenomenon the quasinormal mode level-crossing, in analogy with the wellknown phenomenon of level-crossing for eigenvalues of Hermitian operators. Applying these methods to the strongly coupled N = 4 supersymmetric Yang-Mills theory, we found that the gradient expansions for the hydrodynamic shear and sound modes have finite radii of convergence given by q c sound = √ 2 ω 0 for the sound mode and 21 See footnote 1.

JHEP11(2019)097
by q c shear ≈ 1.49 ω 0 for the shear mode, where ω 0 = 2πT is the fundamental Matsubara frequency. Thus, in both channels, the hydrodynamic series converge up to the order of |q|/T ∼ O(10), which is a vast improvement over the naive expectation that |q|/T 1 provides a natural expansion parameter for hydrodynamic dispersion relations. The obstruction to convergence in the example of the N = 4 SYM theory comes from the collision of poles of the two-point correlation function of the energy-momentum tensor at complex q 2 .
As mentioned in the Introduction, the problem of extending the hydrodynamic modes in the complex momentum plane beyond the branch point singularity was recently investigated by Withers [27] in the context of a holographic model in 2 + 1 dimensions with finite chemical potential. The shear-diffusion mode series could be Padé-resummed and extended beyond the branch point singularity, which was in that case located at an imaginary value of momentum. The main focus of ref. [27] was on the possibility of recovering the full spectrum from the hydrodynamic derivative expansion, similar to recovering the non-hydrodynamic modes from asymptotic series via Borel resummation and resurgence [17,18,25]. The quasinormal spectrum in the holographic models with finite temperature T and non-vanishing chemical potential µ such as the one considered in ref. [27] is rather complicated and changes qualitatively as the parameter T /µ is varied (in particular, it involves pole collisions even at real values of the momentum) [59,66,[75][76][77][78][79]. In ref. [27], the shear-diffusion mode was found to have a radius of convergence inversely proportional to the chemical potential. Naively, this would imply infinite radius of convergence in the limit of vanishing µ, in apparent contradiction with our results. However, the result of ref. [27] was obtained at a specific fixed value of T /µ, and we expect it to change when the complex momentum behaviour of other gapped modes in the model is taken into account with T /µ increasing. This will require further study. It would be also interesting to extend the results of the present work to the sound channel (not considered in ref. [27]) as well as to other holographic models with finite chemical potential such as the STU model [80], and other models [81][82][83], including those in the large D limit [84].
Pole collisions in the correlation functions appear in holographic models in different contexts [33,57,[59][60][61][62][63][64][65][66][67][68]. No less interesting are collisions among poles and zeros of the correlators known as pole-skipping [29][30][31][32]. What we have shown here is that this phenomenon, known to exist in the sound channel of strongly coupled N = 4 SYM theory [29], exists also in the shear and scalar channels of the energy-momentum correlators. The conjectured connection to the OTOC thus allows one to determine the parameters quantifying microscopic many-body chaos (scrambling time and butterfly velocity) by considering the complexified behaviour of the lowest-lying modes (those with the smallest |ω| in the spectrum) in any channel, be it a channel with or without hydrodynamic modes. In general, the critical points and the pole-skipping points are different. We have analysed the 2d CFT finite-temperature correlators and the spectra of the dual BTZ black hole to demonstrate this explicitly. What this implies for the relation between chaos and hydrodynamics is that the "chaos" (or pole-skipping) points can lie within or outside of the radius of convergence of the hydrodynamic series. In particular, while this is not the case in the holographic model with broken translation symmetry considered in section 3, we found that in the JHEP11(2019)097 N = 4 SYM theory, pole-skipping points for both of the two hydrodynamic modes lie within the radius of convergence of the corresponding dispersion relations. This provides an explanation for the observation of the fast convergence of the hydrodynamic series to the exact chaos point in ref. [29].
Can finiteness of the convergence radius of the hydrodynamic modes dispersion relations expansion be taken as a criterion for validity of hydrodynamics? By analogy, one may think of a free particle whose dispersion relation ω = p 2 + m 2 − m = p 2 /2m + . . . has branch points located at p = ±im, and for which the failure of the convergence of the gradient expansion corresponds to the breakdown of the non-relativistic approximation. We hope our results may be of interest for studies of higher-order hydrodynamics necessary for improving the precision of hydrodynamic predictions and also for justifying the construction of the effective field theory of hydrodynamics formulated as a gradient expansion [85][86][87][88][89][90][91][92][93][94][95]. As already mentioned in ref. [8] in the context of the discussion of the "unreasonable effectiveness" of hydrodynamics as an effective theory, many previous studies have reported the divergence of the derivative expansion in relativistic hydrodynamics [15,17,18,18,25]. Possibly, the asymptotic nature of the expansion appearing in those publications should be viewed as a reflection of the singular nature of the state about which this expansion is performed, rather than a generic property of the hydrodynamic gradient expansion itself. On the other hand, even for a free particle, the momentum space and position space pictures look different in this respect: the small-momentum expansion of the corresponding dispersion relation has a finite radius of convergence, whereas e.g. for the position space propagator, the large-time expansion is only asymptotic. 22 This issue needs further investigation. The role of the non-hydrodynamic degrees of freedom is the common feature of the mentioned works and the present paper.
Of special interest is the dependence of the radii of convergence on coupling. In ref. [8], using eq. (2.27) as a crude approximation and the results of refs. [96,97], we argued that in the N = 4 SYM theory, the radius of convergence is smaller at weaker coupling. This, of course, requires the actual study of the spectrum at finite coupling. More generally, in the context of the problem of interpolating between weak and strong coupling regimes of the same theory at finite temperature [60,98], one may note 23 that the problem of convergence of hydrodynamic series has been raised and partially investigated in the 1960s in kinetic theory [99]. This approach, together with recent studies of relevant issues at weak coupling [14,[100][101][102], may deserve more attention in the context of the problem of the validity of the hydrodynamic description at finite coupling.

JHEP11(2019)097
for stimulating and illuminating discussions, and for suggestions on the relevant algebraic geometry literature. P.K. would like to thank the Rudolph Peierls Centre for Theoretical Physics at the University of Oxford for hospitality during the initial stage of this work, and the organisers of the KITP programme "Chaos and Order", where part of this work was completed. P.K.'s work was supported in part by NSERC of Canada. A.O.S. would like to thank Moscow State University and especially A.V. Borisov, as well as the Kadanoff Center for Theoretical Physics at the University of Chicago for hospitality during the final stage of this work. He also thanks P. Glorioso, D.T. Son, M. Stephanov, P.B. Wiegmann and F. Essler, J. March-Russell, S. Parameswaran for discussions in Chicago and Oxford, respectively, and the participants of the seminars at the Institute for Nuclear Research, Steklov Mathematical Institute and Lebedev Physical Institute of the Russian Academy of Sciences for critical questions and useful suggestions. The work of P.T. is supported by an Ussher Fellowship from Trinity College Dublin. We would like to thank A. Buchel, M. Heller, A. Kurkela and J. Noronha for illuminating correspondence. We also would like to thank B. Withers for correcting our inadvertent partial misrepresentation of his results in the first version of our preprint [8].

A Analytic implicit function theorem and Puiseux series
Here, we collect the necessary information from complex analysis regarding the following problem. Given an implicit function f (x, y) = 0, where x, y ∈ C, we would like to find explicit solution(s) in the form y = y(x), at least locally in the vicinity of some point (x 0 , y 0 ), where y(x) may be represented by a finite or infinite series in x. We would like to determine, furthermore, what sets the radius of convergence of such series.
A simple example is provided by the function f (x, y) = x 2 + y 2 − 1 = 0. Since f (x, y) is a polynomial, it determines a complex algebraic curve. Singular points of f (x, y) are determined by the simultaneous solution of the equations f (x, y) = 0, f ,x (x, y) = 0, f ,y (x, y) = 0, where the comma subscript denotes the partial derivative with respect to the argument after the comma. Clearly, this particular curve has no singular points. It does, however, have the so-called "points of multiplicity 1" or "one-fold points", where f ,x (x, y) = 0 or f ,y (x, y) = 0 (but not both simultaneously). These are sometimes called critical points. We are interested in the local behaviour of y = y(x) near a critical point defined by the conditions f (x, y) = 0, f ,y (x, y) = 0. In our example, there are two such points: (x, y) = (±1, 0). The series representation y = y(x) in the vicinity of e.g. (x, y) = (1, 0) has two branches: This is an example of the Puiseux series, i.e. the power series with fractional exponents. These series converge in the circle with the centre at (x, y) = (1, 0) and radius R = 2 which is the distance from (1, 0) to the nearest critical point, (x, y) = (−1, 0). One may be interest in the behaviour y = y(x) in the vicinity of a regular point, where f ,y (x, y) = 0, for example, near (x, y) = (0, 1) in our example. Here, since f ,y (x, y) = 0, JHEP11(2019)097 the implicit function theorem guarantees that we can compute the derivatives y (x), y (x) and so on, and represent y(x) by the Taylor series in the vicinity of x = 0, This series is convergent in the circle of radius R = 1, determined by the distance from the point x = 0 to the nearest critical point(s) at x = ±1. In general, for an implicit function given by the equation f (x, y) = 0, where f (x, y) is either a finite polynomial in x and y, or an analytic function at a point (x, y) (i.e. a polynomial of an infinite degree), the behaviour at a regular point is governed by the analytic implicit function theorem [103], and the behaviour in the vicinity of a critical point is determined by the Puiseux theorem. In the former case, y = y(x) is represented by a Taylor series converging in some vicinity of a regular point. In the latter case, it is represented by a Puiseux series converging in some vicinity of a critical point. We now recall some facts from complex analysis [104] and explain the Puiseux construction [45,47].

Definition.
A function, f (x, y), from a neighbourhood of (x 0 , y 0 ) ∈ C 2 to C is called analytic at (x 0 , y 0 ) if near (x 0 , y 0 ) it is given by the uniformly convergent power series Theorem (Analytic implicit function). If f (x 0 , y 0 ) = 0 and f ,y (x 0 , y 0 ) = 0, there exist > 0 and δ > 0 so that D (x 0 ) × D δ (y 0 ) is in the neighbourhood where f is defined, and g is a map of D (x 0 ) into D δ (y 0 ) so that f (x, g(x)) = 0 and for each x ∈ D (x 0 ), g(x) is the unique solution of f (x, g(x)) = 0 with g(x) ∈ D δ (y 0 ). Moreover, g(x) is analytic in D (x 0 ) and g (x) = − ∂f ∂x (x, g(x)) ∂f ∂y (x, g(x)) . (A.5) Similarly, one can compute higher derivatives of g(x) and represent it by a Taylor series around x = x 0 convergent in D (x 0 ). Note that the statements of the theorem are local, e.g. the size of the domain D (x 0 ) is unspecified, it is only known that it exists for some > 0. In other words, we know that the radius of convergence of the series of g(x) around x = x 0 is non-zero but its value is left unspecified. In the example above, we saw that the value of the radius of convergence is determined by the distance from the centre of the expansion x 0 to the nearest critical point of f (x, y). Note also that the statements of the theorem depend crucially on f (x, y) being an analytic function at (x 0 , y 0 ) (in particular, a finite order polynomial in x and y).
Now we return to the original problem: for f (x, y) = 0, where x, y ∈ C, find explicit solution(s) in the form y = y(x), at least locally in the vicinity of some point (x 0 , y 0 ), where y(x) may be represented by series (possibly infinite) at x = x 0 . For simplicity, we set x 0 = 0. First, we check f (0, y). If this is a polynomial in y of degree n, then the equation f (0, y) = 0 has n roots y i , i = 1, . . . , n. There are two possibilities.

JHEP11(2019)097
Local behaviour at regular points: all the roots y i , i = 1, . . . , n of the equation f (0, y) = 0 are distinct. Then f ,y (0, y i ) = 0, and the analytic implicit function theorem guarantees the existence of a unique Taylor expansion y = y(x) at x = 0.
the radius of convergence of the series solving Kepler's equation is determined by the critical points in the complex eccentricity plane. Kepler's law of motion of a planet in an elliptical orbit with eccentricity e, 0 < e < 1, is usually expressed in a parametric form [106] r = a (1 − e cos ψ) ,

JHEP11(2019)097
where r is the magnitude of the radius-vector from the centre of the force to the planet, a is the major semi-axis of the ellipse, T is the period of revolution and ψ ∈ [0, 2π] is the parameter known in astronomy as the eccentric anomaly. Knowing ψ(t), one can find the position of the planet in the polar coordinates (r(t), ϕ(t)) as a function of time using eq. (C.1) and the equation A formal series solution of Kepler's equation was found by Lagrange [107] who apparently was not concerned with the series convergence (more details can be found in the book [108]) ψ(τ, e) = τ + ∞ n=1 a n (τ ) e n n! , (C.5) where a n = d n−1 (sin n τ ) dτ n−1 . (C.6) As pointed out by Laplace [109], the series (C.5) converges for all τ ∈ [0, 2π] as long as |e| ≤ e L ≈ 0.662743 . . .. For e > e L , the series diverges for some values of τ , in a rather peculiar manner (see figure 17, left panel). What determines the radius of convergence e c (τ ) of the series (C.5)? This problem was investigated by Cauchy, Puiseux and Serret in a series of papers in 1849-1859 [108]. In modern language, the answer is the following. Treat Kepler's equation (C.4) as a complex analytic curve in the space of ψ ∈ C, e ∈ C, with τ remaining a real parameter. The critical points of the curve K = 0 obey eq. (C.4) as well as the equation

JHEP11(2019)097
The critical points closest to the origin in the complex eccentricity plane are shown in figure 18. Their location is parametrised by τ . The radius of convergence e c (τ ) is given by the distance from the origin to the nearest singularity. This distance is a monotonic function of τ in the interval [0, π] (half a period), with the minimum at τ = π/2 given by e c ( π 2 ) = e L . Thus, for 0 < e < e L , the series (C.5) converges for all τ ∈ [0, 2π]. The dependence of the radius of convergence on τ is shown in figure 17 (right panel).
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.