Complexified quasinormal modes and the pole-skipping in a holographic system at finite chemical potential

We develop a method to study coupled dynamics of gauge-invariant variables, constructed out of metric and gauge field fluctuations on the background of a AdS$_5$ Reissner-Nordstr\"om black brane. Using this method, we compute the numerical spectrum of quasinormal modes associated with fluctuations of spin 0, 1 and 2, non-perturbatively in $\mu/T$. We also analytically compute the spectrum of hydrodynamic excitations in the small chemical potential limit. Then, by studying the spectral curve at complex momenta in every spin channel, we numerically find points at which hydrodynamic and non-hydrodynamic poles collide. We discuss the relation between such collision points and the convergence radius of the hydrodynamic derivative expansion. Specifically in the spin 0 channel, we find that within the range $1.1\lesssim \mu/T\lesssim 2$, the radius of convergence of the hydrodynamic sound mode is set by the absolute value of the complex momentum corresponding to the point at which the sound pole collides with the hydrodynamic diffusion pole. It shows that in holographic systems at finite chemical potential, the convergence of the hydrodynamic derivative expansion in the mentioned range is fully controlled by hydrodynamic information. As the last result, we explicitly show that the relevant information about quantum chaos in our system can be extracted from the pole-skipping points of energy density response function. We find a threshold value for $\mu/T$, lower than which the pole-skipping points can be computed perturbatively in a derivative expansion.


Introduction
Low energy dynamics of many body systems near thermal equilibrium can be effectively described in terms of IR variables. These variables themselves are of two types. First are those corresponding to fast thermalizing excitations. These are in correspondence with non-conserved quantities and swallowed by thermal vacuum, very soon after being excited. Second are those corresponding to conserved quantities whose excitations cannot be locally thermalized. They are relaxed just through the transport in long wavelength compared to local thermalization time scale.
In fact the first-type degrees of freedom are in local equilibrium defined by the secondtype ones. Hydrodynamics is an effective theory which describes the dynamics of latter, namely slow modes, in a many-body system. Hydrodynamic equations are then conservation equations of stress tensor T µν and charge current J µ in the system. The main idea of hydrodynamic is that, in local equilibrium, each of these quantities are given in a perturbative expansion of local thermodynamic variables and their derivatives. Such perturbative expressions are called hydrodynamic constitutive relations. Using symmetry arguments, the general form of constitutive relations can be fixed up to a set of unknown coefficients, which are referred to as transport coefficients. While hydrodynamics is a universal regime of thermal systems, the associated transport coefficients depend on the underlying microscopic theory by which the system is described.
Considering hydrodynamics as a classical field theory, one can find simply the response functions of the conserved quantities (see refs. [1,2]); the corresponding poles are then the so-called hydrodynamic modes. Hydrodynamic modes are the longest lived modes around thermal equilibrium in the system, so ω(q → 0) = 0. In contrast, the excitations corresponding to the first-type degrees of freedom, mentioned in the first paragraph, are short-lived in the sense that ω(q → 0) = 0. They are called non-hydrodynamic modes.
Hydrodynamic modes can be in general found from the linearized hydrodynamic equations perturbatively, order by order in a derivative expansion (see ref. [1] for a comprehensive review). However, without knowing the microscopic of system, the result will be restricted to general perturbative expressions containing unknown transport coefficients. Moreover, even by having the microscopic theory, computing the transport coefficients in general is limited to lower orders in derivative expansion and also to perturbative regime of microscopic field theory [3].
When the microscopic theory has a holographic dual [4], however, extracting information will be relatively straightforward not only in the hydrodynamic limit, but also beyond that. In fact holography allows to know the microscopic of certain strongly coupled field theories via studying their gravitational dual in a one-higher-dimensional AdS space time. In last two decades by analyzing perturbations around specific AdS black brane solutions, significant information about the strongly coupled N = 4 SYM plasma have been found: transport coefficients of first order hydrodynamics specifically η/s [5], location of sound poles at the same order [6], the leading 't Hooft coupling correction to η/s [7], the spectrum of quasinormal modes [8], second order transport coefficients [9,10], infinite 't Hooft coupling correction to second order transport coefficients [11], anomalous transport coefficients [12,13], third order transport coefficients for linear [14] and full non-linear hydro [15] and many other related issues.
At first sight higher order terms in hydrodynamic constitutive relations seem just to improve the accuracy of derivative expansion. However, as it was found in refs. [17], large order behavior of hydrodynamic expansion might have some information about non-hydrodynamic modes as well. By numerically computing energy density for a holographic boost-invariant flow, up to terms with 240 derivatives, authors of ref. [17] found the radius of convergence of hydrodynamic derivative expansion to be zero. Let us emphasize that this result is associated with an expansion over the inverse proper time in position space. Then by using Borel resummation techniques, they found the frequency corresponding with the leading singularity in Borel-transformed hydrodynamic stress tensor to be exactly the same as frequency of lowest-lying non-hydrodynamic mode, found in refs. [18,19] 1 .
In a recent paper [22], information about non-hydrodynamic modes has been extracted from large order behavior of derivative expansion via a new approach. Instead of working with full non-linear hydrodynamics to study a highly symmetric flow in position space [17], the choice of ref. [22] is to work with linear hydrodynamics in the complex q plane. The focus of latter paper is to study the dispersion relation of hydrodynamic shear mode at large orders in derivative expansion, in a holographic model in 2 + 1 dimensions with finite chemical potential. Considering the multi-sheet structure of exact ω shear (q), the closest non-analytic points on the longest-lived sheet are found to be as q = ±iq * [22], where q = |q|. The author of [22] then argues that q * sets a finite radius of convergence for the hydrodynamic expansion of the shear mode 2 . As explicitly mentioned in ref. [22], the physical origin of q * is the collision between a hydrodynamic and non-hydrodynamic mode on the imaginary q axis.
Much more recently, it was shown that the same result could be found from analytic structure of spectral curves in classical hydrodynamics 3 . In ref. [24], the above-mentioned collision points were found from the associated spectral curves in a holographic neutral fluid in 3+1 dimensions, for both sound and shear hydrodynamic modes. In the mentioned reference, such collision points are called level-crossing points. In fact the theory of [24] is the general theory for studying spectral curves and level-crossing.
The value of q * found in ref. [22] is inversely proportional to µ. This simply implies an infinite radius of convergence for the shear mode in the limit of vanishing µ. But it obviously contradicts with result of ref. [24]. However since the former was obtained at a specific fixed value of µ/T , one may expect it to change when the complex momentum behavior of other gapped modes is taken into account with T /µ varying [25].
To investigate how actually q * may vary with µ/T and also to study the convergence of other modes, in this paper we follow the issue in a holographic model at finite chemical potential in 3 + 1 dimension. In the gravity side we consider a AdS 5 Reissner-Nordström black brane. In the bulk of AdS 5 , such black brane is identified with the temperature T as well as a parameter Q, both related to N = 4 SYM boundary theory, where Q is a monotonically increasing function of µ/T .
The main difficulty in our problem comes from the fact that in our model parturbations of gravity and gauge field in the bulk are coupled. One may think of finding master fields and then deriving decoupled equations governing their dynamics [26] 4 . But we choose to work with coupled equations! The reason is related to the numerical method that we want to adopt. We construct the generalized version of Frobenius expansion used in ref. [8] to find the quasinoraml modes associated with coupled differential equations. The advantage of working with coupled equations then is that it lets each of our results at Q = 0 be comparable with some counterpart result within refs. [8] and [24]. Let us note that our numerical method, turns out to be reliable in the range 0 ≤ Q < 0.88, or equivalently 0 ≤ µ T 4. However, in our plots we will demonstrate the results within the range 0 ≤ Q < 0.8.
We first construct gauge-invariant variables out of bulk field perturbations, in the Spin 0, 1 and 2 representations of SO(2) group corresponding to isotropy of transverse plane perpendicular to q. This is actually the subject of § 2. In § 3 we develop a new method to find the quasinormal modes of coupled differential equations in the bulk. In § 4 we will derive coupled equations governing dynamics of gauge-invariant variables in each of Spin channels, on the AdS 5 RN black brane background. Our equations can be regarded as non-trivial generalizations of decoupled equations associated with an AdS 5 Schwarzschild black brane accompanied by a probe gauge field [8].
In each channel of Spin, we find the spectrum of quasinormal modes by considering variations of both q = q/(2πT ) and Q. Spin 0 and Spin 1 spectra, each turns out to be a superposition of two types of poles. In Spin 0 channel, poles are in correspondence with fluctuations of scalar components of stress tensor and charge current. The hydrodynamic modes associated with former fluctuations are sound modes while that of related to the latter fluctuations is the diffusion mode. In Spin 1 channel, poles correspond to fluctuations of vector components of stress tensor and charge current. In this case, only one hydrodynamic mode does exist which iz related to stress tensor fluctuations. That is actually the shear mode. Finally, from the single dynamical equation of Spin 2 channel, we find poles corresponding to fluctuations of tensor components of stress tensor. As expected, there is no any hydrodynamic mode in this channel.
The spectrum of quasinormal modes at finite chemical potential have been already computed in AdS 4 RN model [32][33][34]. But in holographic models in 3+1 dimensions, well-known results are just limited to Spin 1 and Spin 2 channels [35]. To best of our knowledge our study is the first computation of quasinormal modes in Spin 0 channel of AdS 5 RN model. In addition, by analytically solving the coupled dynamical equations in the hydrodynamic limit, we find the spectrum of hydrodynamic modes as well. Again, this is the first analytic computation of hydrodynamic modes from AdS 5 RN black brane. However it should be noted that we will be able to find analytic solutions in the bulk at small µ/T . We then use fluid/gravity [12,13] to confirm our results by explicit computation of hydrodynamic modes at finite µ/T .
In the second part of the paper, in § 5, we study quasinormal modes at complex momenta. We find radii of convergence of derivative expansion corresponding to dispersion relations of all three hydrodynamic modes, separately, within the range 0 ≤ Q ≤ 0.8. For the shear mode, level crossing between hydrodynamic and non-hydrodynamic poles in Spin 1 channel sets the radius of convergence of associated dispersion relation. It turns out that the radius of convergence monotonically increases when Q increases from 0 to 0.8.
Our surprising results are mostly related to Spin 0 channel. Firstly we motivate that for each of sound and diffusion modes, one has to find the critical points of spectral curve lying on their own branch of Puiseux series 5 . Doing so, for the diffusion case we show that within the whole range of Q, the level crossing happens between diffusion pole and gapped poles associated with charge current spectrum. But sound modes are found to collide with various types of poles, depending on the value of µ/T . Their first collisions are actually with either gapped poles of stress tesnor spectrum or gapped poles of charge current spectrum or even with diffusion pole! The latter is a novel aspect of level crossing phenomenon, specific to holographic systems at finite chemical potential. It can be regraded as a counterexample for the statement that finite radius of convergence of hydrodynamic derivative expansion is determined by the interplay between hydrodynamic and non-hydrodynamic modes. In other words, we numerically find that in a specific range of µ/T , the convergence radius of derivative expansion is fully determined by hydrodynamic information.
In the last part of the paper, in § 6, we study one another aspect of quasinormal modes. Following recent studies on the relation between hydrodynamics and quantum chaos in maximally chaotic systems [36][37][38], we will show that the pole-skipping point of energy density response function in our system precisely coincides with the chaos point of the system. The latter is found without using any assumption on the value of µ/T . This result indeed provides a new support for the hydrodynamic origin of quantum chaos [37]. Then we focus on Spin 0 channel and numerically find the dispersion relation of sound and diffusion poles at purely imaginary momenta. For a typical value of µ/T , we show that the above-mentioned pole-skipping points lie actually on the sound curve.
In the very last part, we discuss the possibility of finding the chaos point by using the derivative expansion. Exploiting the Spin 0 channel results of § 5, we will find that the chaos point does not always lie with in the domain of validity of hydrodynamics. It turns out that for values of Q smaller than a critical value Q c , derivative expansion is sufficient to extract the chaos information from the energy density response function. However as Q > Q c , chaos points lie outside of domain of convergence of hydrodynamic derivative expansion. It simply means that for this range of Q, the pole-skipping point of energy density response function has to be found non-perturbativley.
Finally in § 7, we briefly review our results and discuss possible follow-up directions.

Gauge invariant variables
The back ground solution on which we would like to find the quasinormal modes can be written in general as whereq = q 0 /r p h is the electric charge density on the horizon. r h is the radial location of the horizon: f (r h ) = 0. We take the fluctuations of metric and gauge field to be h µν (r)e −iωt+iqz and a µ (r)e −iωt+iqz , respectively, where z = x p . In the following we mostly focus on p = 3 case. We then find the specific combinations of perturbations which are invariant under the simultaneous general diffeomorphism, denoted by ξ µ , and gauge transformations φ in the bulk. Under such transformations, It is convenient to arrange perturbations in the gauge invariant representations of SO(2) group corresponding to the rotational symmetry of the two dimensional plane perpendicular to q. The corresponding transformations have been given in the Appendix A. In the following, we specify the general form of probable gauge invariant quantities in different channels of SO(2).

Spin 0 channel
In this channel six perturbations G i ∈ {a t , a z , h tt , h tz , h zz , h}, where h = 1 2 (h xx + h yy ), are coupled to each other. Any first order gauge invariant variable, say G inv , is naturally a linear combination of G i 's: On the other hand, there are in general four diff and gauge parameter functions in this channel: ζ i ∈ {φ, ξ t , ξ z , ξ r }. The task is to find a i 's such that G inv to be independent of ζ i 's. Under the diff and gauge transformations (2.2), G i 's are transformed as with β i being a linear combination of α i 's. Gauge invariance demands β i 's must vanish. So we find four equations among six parameters α i . As a result, two of the parameters α i remain free; it simply shows that in the sound channel we deal with two gauge invariant variables. Let us now find them in details. From (A.1), the above mentioned two gauge invariant variables can be found either by • finding an appropriate combination of its first four lines; or by • finding an appropriate combination of its last three lines.
In the first case, we take h zz + 1 h tt + 2 h tz + 3 h and then demand its diff-gauge transformed be independent of ξ t , ξ z and ξ r . We find In the second case, we take a z + e 1 a t + e 2 h combination. Again, demanding its diff-gauge transformed be independent of ξ r and ξ t fixes the coefficients e 1 and e 2 : In summary, the gauge-invariant variables in the Spin 0 channel can be written as where A t = a t /c h , A z = a z /c h , H tt = h tt /af , H tz = h tz /a, H ij = h ij /a (i, j = t) and H = h/(p − 1)a. We have also defined c h = c(r h ). In the above equations, 0 denotes the Spin 0.

Spin 1 channel
From (A.2), it is obvious that E α = iωA α (α = x, y) are the two gauge invariant objects. In addition, two another gauge-invariant variables can be constructed by taking h zα + κh tα and searching for appropriate κ . Demanding this combination be independent of diffgauge parameter functions, one obtains κ = q/ω. Therefore the gauge invariant objects in this channel are (due to the SO(2) symmetry, we just work with two of those four) 3 Quasinormal modes from coupled differential equations As we will see later, except for the Spin 2 channel, the dynamics in each of the other two channels is governed by a specific set of coupled differential equations. In this section, before finding the equations themselves, we argue how to find, in general, the associated quasinoaml modes from such coupled equations 6 . Following our notation in the previous section, let us take the two gauge invariant variables in these channels as Z and E. In one specific channel, say Spin 0, upon imposing the ingoing boundary condition at the horizon, the solution to these fields can be formally written as The pair of coefficients C Z and C E are fixed by the value of Z and E at the horizon, so y Z (1) = g E (1) = 0. Moreover, the functions y Z (u) and g E (u) vanish at Q = 0 as well because we know at Q = 0 the two equations decouple.
In order to find associated quasinormal modes, we need the near boundary expansion of these fields. Considering the conformal dimension of Z and E, one writes with constants c 1 and c 2 being unimportant for our next discussions. Using these solutions, the bilinear boundary bulk action takes the following form: with k = (ω, k). By making an appropriate unitary transformation U , the middle matrix in above is simply diagnolized: The transformed variables are then given by (3.5) The holographic AdS/CFT duality [4] implies that variables Z and E couple to specific operators O 1 and O 2 in the boundary theory. From the near boundary expansion of Z and E, we can then find the Green's function of these operators. Substituting (3.2) into (3.5), we find: with the coefficients given by Note that the forms of B Z and B E are unimportant for our following discussion. Now, the Green's functions are simply given by Then the quasinormaol modes of perturbations Z and E are defined by the following equations [40][41][42] A Z = 0, Substituting (3.7) into (3.9) and considering the fact that det U = 0, we arrive at the following set of equations These equations leads to nontrivial solutions if This is our first result in this paper; the equation from which, we can find the quasinormal modes of coupled perturbations Z and E. Let us summarize. In order to find the quasinoraml modes of coupled perturbations Z and E in either spin 0 or spin 1 channels, we first consider the corresponding coupled differential equations. We should find the solutions to them which are ingoing at the horizon. The result can be formally written in the form of (3.1), up to two arbitrary normalization constants C Z and C E . Then equation (3.11) determines the spectrum of quasinormal modes in the associated channel.
4 Quasinormal spectrum and hydrodynamic modes in N = 4 SYM theory at finite chemical potential Our system of interest is described by dynamics of metric and a U (1) gauge field in the bulk of AdS. The corresponding action is given by where S bdy is the boundary counter term. The equations of motion are given by: We work in the unite where L = 1. The solution in the Poincare patch is written as it follows r 6 and m and q being two constants. It is convenient to re-scale the quantities with the radius of the outer horizon, R, namely the largest positive root of f (r) = 0. We may write In the rescaled coordinates, the outer horizon locates at ρ = 1 while the boundary is identified with ρ → ∞. For further requirements, we need to work in a system with finite domain for the radial coordinate; for this purpose, we make the change ρ = 1/ √ u. The function f then takes the following form (4.5) Using this together with (4.4), the Hawking temperature T and the chemical potential µ of the boundary theory are found to be From the expressions (4.6), one can derive Q as a function of µ/T : This relation simply shows that there is a one-to-one map between µ/T and Q. Thus in the following, we take T and Q as the two independent thermodynamic variables. Then the bulk metric and gauge field can be finally written in (t, x, y, z, u) coordinates as the following: (4.8)

Spin 0 Channel
We turn on the set of perturbations {h tt , h tz , h zz , h, A t , A z } in the radial gauge. On the RN background solution (4.8), the gauge invariant variables (2.7) then take the following form Note that on the background solution (4.8),q = Q. To find the dynamical equations governing dynamics of the above two variables, we have to replace in spin 0 equations of (4.2) (see (2.7) and explanations given below that). Doing so, we firstly find dynamical equations of the latter "six" perturbations. The difficult task is to eliminate all these perturbations in favor of the two gauge invariant variables Z 0 and E z . After long computations, which are not shown here, we arrive at the two following coupled differential equations where we have put a i coefficients in front of Z 0 and Z 0 and have done the same for b i coefficients with E z and E z . Considering w = ω/2πT and q = k/2πT , the coefficients are found to be as At this point it should be noted that when Q = 0, the above equations reduce exactly to the pair of decoupled equations (4.5b) and (4.35) in [8]. Finding the quasinormal modes as well as hydrodynamic modes from these equations is the subject of following subsections.

Quasinormal modes
The analytic solution to equations (4.10) is unknown; so to find the associated spectrum of quasinormal modes, one may proceed with numerically solving them. To this end, we combine the method developed in § 3 together with the Frobenius expansions of Z 0 and E z in the bulk (see Appendix B for more details). In what follows, the corresponding numerical results will be given. In Fig.1, the typical arrangement of poles has been demonstrated for two cases in this channel. In the left panel, we have just shown the quasinormal modes associated with q = 1 at Q = 0.5. As can be seen, we split them into two sets, denoted by dots and stars. The idea for such spitting comes from our previous knowledge about the arrangement of poles in sound channel as well as in the R−current channel on the AdS-Schwarzschild background [8]. While in the latter case R−current is decoupled from the spin 0 fluctuations of stress tensor, on the AdS RN background they both together constitute the spin 0 channel. Therefore by comparing the result given in this Figure with the associated spectra of AdS-Schwarzschild background, we have been able to distinguish between stars and dots. In the right panel of Fig.1, we have compared the left panel spectrum with the one corresponding to the same Q but a smaller q. One clearly notices that as q decreases, all complex poles Figure 1: Left panel: Stars represent poles of the charge current Green's function. Dots correspond with poles of energy density Green's function. Right panel: As q decreases, all poles stay at a finite distance from the real axis, except for the one marked with a large star and the two marked with large dots. Large dots therefore manifest the existence of two sound modes and the large dot corresponds with existence of a diffusive U (1) charge mode in the boundary N = 4 SYM theory. Figure 2: The spectrum of eight lowest complex quasinormal modes and the lowest purely imaginary one in spin 0 channel, at q = 1. Every colorful path-like set of points, starting in purple and ending in red, shows the change of one specific mode when Q discretely increases from 0 to 0.8. We have considered 16 regular steps by increments of ∆Q = 0.05. move away from the horizontal axis, but large dots and large star move towards the origin. This simply shows that the spectral curve of spin 0 fluctuations includes three branches of Puiseux series passing through the origin of complex plane. In the next subsection, we explicitly derive the equation of these branches in the vicinity of origin. These three branches correspond to three hydrodynamic modes: two sound modes together with one diffusion mode. It should be also noted that the poles other than these three are all gaped modes. So far we have just talked about the quasinormal modes at a fixed finite value of Q. In Fig.2, we have depicted part of the spectrum of quasinormal modes at the fixed momentum q = 1 for several values of Q within the range 0 ≤ Q ≤ 0.8. One observes that by approaching towards extremal limit 7 , even at a finite fixed value of momentum, e.g. at q = 1, the non-hydrodynamic (gaped) modes move away from the origin. As mentioned in the Introduction, our numerical method does not cover the whole range 0 ≤ Q ≤ √ 2.
It would be interesting to try other numerical methods to investigate the extrapolation of colorful paths depicted in Fig.2.

Hydrodynamic limit
Although the complete analytic solution to equations (4.10) is unknown, in the hydrodynamics limit, w(q → 0) = 0, we are able to find analytical behavior of the solution. At first sight, equations (4.10) may seem impossible to become decoupled. But as we will show, at small charge limit, namely when Q 1, they decouple in the hydrodynamic limit. What we are going to do is to solve (4.10) in the hydro limit and up to second order in the perturbative expansion over Q.
Demanding the solutions be ingoing at the horizon, simply fixes the near horizon behavior of Z 0 and E z : In order to enter the hydrodynamic expansion, we apply the rescaling w → w and q → q to the dynamical equations as well as to (4.13). A quick look at the the coefficients (4.11) and (4.12) reveals that Such scaling simply means that in each of the equations (4.10), the function Z 0 must be one order higher in expansion, than the function E z . Thus the appropriate ansatz for the functions Z 0 (u) and E z (u) in the double expansion over and Q is given by (4.16) Substituting the above expressions into equations (4.10) and expanding over and Q, we obtain a set of second order ordinary differential equations for the functions Z n,m 0 and E n,m z . Starting from the lowest order in and Q, one firstly finds Z 0,0 0 (u) from the first line of (4.10) up to an unknown coefficient. Regularity at u = 1 together with fixing its value at the same point, namely Z 0,0 0 (1) = C 1 , picks out a unique regular solution to Z 0,0 0 (u). Using this solution, then the next function that can be found is E 0,0 z , from the second equation of (4.10). Again, regularity at u = 1 and demanding E 0,0 z (1) = C 2 fix the solution. In the Appendix D, we have listed the corresponding perturbative solutions at higher orders, according to the ordering we have found them through the perturbation theory. It should be noted that at every step of perturbation, one of the two boundary conditions is regularity at u = 1 and the second one is the the normalization of the solution at the same point: Quasinormal modes are by definition [40], the modes obtained upon applying the Dirichlet boundary condition to (4.16). One thus writes According to our earlier discussions these two equations are coupled. By explicit computations (Appendix D) we find 8 where C 1 and C 2 are the normalization coefficients defined in (4.17). The m ij coefficients, up to first order in derivatives, are found to be where y = w/q. Obviously, in the limit Q = 0, the two non-diagonal coefficients m 12 and m 21 vanish and one is left with two decoupled equations in (4.19); C 1 m 11 = 0 and C 2 m 22 = 0. Requiring C 1 , C 2 = 0, these equations then give the well-known dispersion Solving the recent equation to second order in q and Q, we find the dispersion of spin 0 hydrodynamic excitations in a holographic charged fluid as it follows The expression in the first line of (4.22) is the dispersion relation of sound modes in the charged fluid. The second line is showing the hydrodynamic diffusion of the U (1) charge. Finally the expression given in the third line does obviously correspond to a gaped mode which is beyond the regime of hydrodynamics. To best of our knowledge, this is the first computation of hydrodynamic modes in a AdS RN background. 9

Spin 1 Channel
We turn on the set of perturbations {h tx , h zx , h, A x } in the radial gauge. As mentioned in § 2, the gauge invariant quantities associated with this channel are identified with Z 1 and E x (2.8). Just like in the spin 0 channel, the difficult part of the computation here is to find dynamical equations governing dynamics of these quantities. Combining all spin 1 components of (4.2) and after performing long computations, which are not shown here, we have eliminated {h tx , h zx , h, A x } in favor of Z 1 and E x . Eventually we have arrived at the following two coupled dynamical equations (4.23) 9 In [28], using the master field method, only the first term of sound mode, namely w = 1/ √ 3q has been found. Authors of [28] have also computed the first order hydrodynamic transport coefficients via using the relevant Kubo formulas. Figure  As it must be, one can readily check that when Q = 0, the above equations reduce exactly to the pair of decoupled equations (4.5a) and (4.26) in [8]. In the following two subsections, we proceed with numerically and analytically solving the above equations, respectively.

Quasinormal modes
The analytic solution to coupled equations (4.23) is unknown. However, just like what was done in the spin 0 channel, we can use the Frobenius expansions of Z 1 and E x to find the corresponding quasinormal modes via the method developed in § 3. Firstly in order to get familiar with the typical arrangement of quasinormal modes in the spin 1 channel, in Fig.3 and in the left panel we have shown the spectrum at q = 1 for Q = 0.5. We have split the poles into two sets, denoted by stars and dots. The idea of such splitting originates from the corresponding spectra of quasinormal modes associated with decoupled variables Z 1 and E x on an AdS-Schwarzschild background. One naturally concludes that dots in the above figures correspond to the poles of spin 1 fluctuations of stress tensor, T xz . The stars then identify the poles of spin 1 fluctuations of boundary charge current, J x .
In the right panel of Fig.3 we have compared the left panel spectrum with the spectrum associated for the same Q but at a smaller q. As it can be seen, when q decreases, all complex modes move away from the real axis. At the same time, the large dot mode becomes closer and closer to the origin. It simply shows that there is only one branch of Puiseux series associated with the spectral function of Spin1 channel that passes through the origin. In other words, there is only one gapless mode in this channel which is actually the hydrodynamic shear mode. In next subsection we explicitly derive the dispersion relation of this mode. Figure 4: The eight lowest complex quasinormal modes and the lowest purely imaginary one in spin 1 channel, at q = 1. Every colorful path-like set of points, starting in purple and ending in red, shows the change of one specific mode when Q discretely increases from 0 to 0.8. We have considered 16 regular steps by increments of ∆Q = 0.05.
Our discussion on quasinormal modes in this channel has been so far limited to the case with a fixed value of Q. In Fig.4, we have demonstrated part of the spectrum of quasinormal modes at the fixed momentum q = 1 for several values of Q within the range 0 ≤ Q ≤ 0.8. One observes that by approaching towards extremality, the gaped modes move away from the real axis, while the gapless mode becomes close to the origin. It would be interesting to study the spectrum in the extrapolation region 0.9 < Q ≤ √ 2 to see what the fate of this mode in the extremal limit will be [43].

Hydrodynamic limit
In this subsection we are going to find an analytic solution to the dynamical equations (4.23) in the hydrodynamic limit w(q → 0) = 0. As mentioned in [8], it turns out that the appropriate rescaling in this channel is w → 2 w and q → q with 1. But even in this limit, equations (4.23) are coupled; then as was the case in the spin 0 channel, we proceed with perturbatively expanding equations over Q as well. This expansion for Q 1 together with the derivative expansion over make it possible to find an analytic solution for Z 1 and E x .
At this point we have to investigate the behavior of coefficients in (4.23) under the above-mentioned hydro rescaling. Doing so we find that the appropriate ansatz for the functions Z 1 and E x is given by: To leading order in and second order in Q, we find (4.27) It is obvious that at Q = 0 limit, equations (4.26) decouple and one finds w = −iq 2 /2, the well-known shear mode in a holographic neutral fluid [8]. When Q = 0, equations are coupled and have non-trivial solutions if and only if det(s ij ) = 0; this gives which shows how the finite density of U (1) charge modifies the dispersion of the shear mode in the system.

Spin 2 Channel
As discussed in § 2, by turning on metric perturbations {h xy , h xx , h yy } in the radial gauge, two spin 2 gauge invariant variables Z 2 and W 2 are excited (see (2.9)). We find that on the RN background solution (4.8), these quantities commonly obey the following equation (4.29) With no need to the follow § 3, we take the following Frobenius expansion (4.30) Substituting it into (4.29), we find coefficients c n all in terms of c 0 . Then by applying the Dirichlet boundary condition we arrive at the spectral curve of fluctuations in this channel: By keeping sufficient number of terms in the sum, one can numerically find the spectrum of quasinormal modes. In Fig.5, we have demonstrated the spectra associated with Q = 0.5 for two values of q. For both cases, we have shown five lowest-lying quasinormal modes. When q decreases, all complex modes move away from the real axis. At the same time, the purely imaginary mode becomes closer to the origin, however, it never reaches the origin. This observations is enough to conclude that there is no any branch of poles passing through the origin when q → 0. In other words, spin 2 channel of fluctuations are always gapped.
Our discussion on quasinormal modes in this channel has been so far limited to the case with a fixed value of Q. In Fig.6, we have demonstrated part of the spectrum of quasinormal modes at the fixed momentum q = 1 for various values of Q within the range 0.45 ≤ Q ≤ 0.8. One observes that by approaching towards extremality, the gapped modes move away from the real axis, while the purely imaginary mode becomes close to the origin. It would be interesting to study the spectrum extrapolating region 0.9 < Q ≤ √ 2 to see what the fate of this mode in the extremal limit will be. It is clear that both non-hydrodynamic and possible hydrodynamic modes are encoded in this equation. Hydrodynamic modes, w(q 2 → 0) = 0, are expected just to live in the vicinity of origin (0, 0). The small-q expansion of w(q 2 ) then can be found by using the theorem of Puiseux. The Puiseux analysis implies that the domain of convergence of Puiseux series centered 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 [24]. The critical points of the spectral curve, themselves, can be found by solving the following set of equations: In [24], it has been shown that the critical points obtained from these equations are exactly the level-crossing points of the quasinormal modes in the associated channel. The focus of [24] is to study the quasinormal modes of a holographic neutral fluid. However, when fluid carries a conserved charge as well, the arrangement of complexified quasinormal modes may change significantly. Specifically in the spin 0 channel, in addition to two sound branches of the Puiseux series passing through the origin, the diffusion branch passes through the same point too. It may possibly give rise to emergence of new critical points, due to crossing between sound and diffusion branches. As we explicitly show in next subsection, such critical points will appear in the spin 0 channel when the parameter Q exceeds a specific threshold. 10 It should be noted that we did not explicitly write down the spectral curve equations anywhere in this paper. However it would be useful to note that in the hydrodynamic limit, the spectral curve of the spin 0 channel is obtained by substituting (4.20) into (4.21). In the spin 0 channel, it is simply given by the determinant of (4.27).

Spin 0 channel
In § 2 we formally argued how to find the quasinormal modes of coupled variables in the bulk. For the present spin channel, the spectral equation, namely F 0 (q 2 , w; Q) = 0, is nothing other than (3.11), obtaining from (4.10). We already solved this equation numerically in previous section and the corresponding spectrum of quasinormal modes typically was shown in Fig.1. From the latter one notices that F 0 (q 2 , w; Q) = 0 has three eigen frequencies in the vicinity of (0, 0). Assuming the analyticity of F 0 at (0, 0), then the implicit function theorem gives these three branches by Puiseux series as it follows Following [24], what we refer to as the hydrodynamic derivative expansion is the above series in the momentum space. Let us recall that the two sound branches are located symmetrically with respect to imaginary axis in the complex w plane. As a result, in half of the w plane, say for instance within Re w > 0, there are exactly two branches passing through the origin; w + sound together with w diffusion . By naively applying the statement of [24] to the present case, one may conclude that the distance between (0, 0) and the critical point of F 0 , the nearest to origin, identifies the radius of convergence of the derivative expansion in this channel. But an immediate follow-up question is: to which series given in (5.3) such radius corresponds?
Needless to say that w + sound and w diffusion may have different radii of convergence. In fact the radius of convergence of w + sound (w diffusion ) is identified with the distance between origin and the nearest critical point located on the sound (diffusion) branch of Puiseux series. Thus in addition to find the nearest critical points, via solving (5.2), one has to be careful about positioning of them on the considered branch of Puiseux series.
We have numerically found the nearest critical points to the origin, on both sound and diffusion branches of Puiseux series, for several values of Q, within the range 0 ≤ Q ≤ 0.8. The result has been shown in Fig.7. Typical behavior of diffusion mode, shown by the red curve, seems to be qualitatively the same for the whole range of Q. The situation for sound mode, the blue curve, however, is more complicated. It should be studied in three different intervals; (i) 0 ≤ Q ≤ 0.386, (ii) 0.386 ≤ Q ≤ 0.610 and (iii) 0.610 ≤ Q ≤ 0.8.
To explore more on the relation between radius of convergence of the derivative expansion and Q we now start to study the complex life of quasinormal modes. We assume q 2 to be a complex number q 2 = |q 2 |e iθ . Then θ = 0 simply corresponds to poles with real q 2 , already shown in Fig.1. As before, we show them by (large) dots and (large) stars in upcoming figures. Then we let the phase θ change from 0 to 2π. For a given |q| 2 , such change of θ corresponds to moving dots and stars along some trajectories in the complex w plane. The interaction of poles via crossing of their trajectories is the main issue that we will discuss in details for each of the three intervals mentioned above, separately.
(i) 0 ≤ Q ≤ 0.386: In this interval we choose to show the results associated with Q = 0.3. See Fig.8. It is clear that by increasing |q| 2 , trajectories of poles become more complicated. Let us just firstly consider the highest large dots. These pols actually lie on the sound branch of the Puiseux series near the origin. When |q c | 2 ≈ 2.32, their trajectory collides with that of (dot) gapped mode poles, at two points (marked by black dots in the bottom row plots). One then concludes that at Q = 0.3 the radius of convergence of the derivative expansion for the w ± sound is |q sound c | ≈ (2.32) 1/2 ≈ 1.52.
Since the above-mentioned collision happens between two (dot) poles which both belong to the stress tensor spectrum, we call the crossing of the associated trajectories the T T −crossing. As a result, within the range 0 ≤ Q ≤ 0.386, this is T T −crossing which determines the radius of convergence of the derivative expansion for the sound branch. As pointed out in Fig.7, the AdS-Schwarzschild case with Q = 0, studied in [24], falls into the same range.
(ii) 0.386 ≤ Q ≤ 0.610: In this interval we choose to show the results associated with Q = 0.5. See Fig.9. By increasing |q| 2 , the poles tend to collide. Interestingly, it turns out that the first collision of the sound poles (the highest dots) is with the other  Fig.7. After the collision, for instance at |q 2 | = 2.34, the orbits of sound pole and the two nearest dot gapped poles are no longer closed: four of them exchange their positions as the phase θ increases from 0 to 2π. This is the manifestation of the T T −crossing.
hydro pole in this channel, namely the diffusion pole denoted by a star on imaginary axis. The collision point is identified with critical value of momentum |q c | 2 ≈ 2.04. One then concludes that at Q = 0.5 the radius of convergence of the derivative expansion for the sound branch of Puiseux series is |q sound c | ≈ (2.04) 1/2 ≈ 1.43. Let us emphasize that while one of the two colliding poles comes from the stress tensor spectrum, the other one belongs to the spectrum of charge current. Therefore it is reasonable to call such crossing Figure 9: Poles of the retarded two-point function in the spin 0 channel at Q = 0.5, in the complex w−plane, at various values of the complexified momentum q 2 = |q 2 |e iθ . Large dots and large stars correspond to the poles with purely real momentum (i.e. at θ = 0). As θ increases from 0 to 2π, each pole moves counter-clockwise following the trajectory whose color changes continuously from blue to red.  Fig.7. After the collision, for instance at |q 2 | = 2.06, the orbits of sound poles and the diffusion pole are no longer closed. They exchange their positions cyclically as the phase increases from 0 to 2π. This is the manifestation of the T J −crossing.
of trajectories the T J −crossing.
(iii) 0.610 ≤ Q ≤ 0.8: In this interval we choose to show the results associated with Q = 0.7. See Fig.10. At this value of Q, specifically, we illustrate crossings associated with critical points of both sound and diffusion branches of Puiseux series. It turns out Figure 10: that by increasing |q 2 |, the first collision of diffusion pole would be with the nearest star gapped poles (see the black dots in top middle and top right panels). This occurs at |q 2 c | ≈ 1.50. Since both colliding poles are star poles, we call such crossing of trajectories the J J −crossing. The first collision of the sound poles, however, will occur at a larger value of momentum. As can be seen in the bottom panels of Fig.10, at |q 2 c | ≈ 1.95 sound poles collide with the nearest purely imaginary gapped pole. Since a dot pole collides with a star pole, such crossing of trajectories is a T J−crossing.
Before ending this subsection let us give two comments regarding our results. Firstly, as mentioned implicitly earlier and also can be seen by the behavior of red curve in Fig.7, the critical point of diffusion branch, not only at Q = 0.7 but also in the whole range of 0 ≤ Q ≤ 0.8 is of JJ−crossing type. Secondly, one may wonder if T J−crossings associated with intervals (ii) and (iii) are the same. In fact within (ii) the sound pole collides with the other hydro pole, namely the diffusion pole; while within (iii) it collides with a gapped pole belonging to the spectrum of charge current. The T J−crossing within the interval (ii) is then quite novel in the sense that it implies convergence of the derivative expansion for 0.386 ≤ Q ≤ 0.610 is controlled totally by the hydrodynamic information. This is a new aspect of level-crossing specific to systems at finite chemical potential. In previous studies at vanishing µ, this level-crossing was found as the result of interplay between hydro and non-hydrodynamic poles.

Spin 1 channel
The typical spectrum of quasinormal modes in this channel was already shown in Fig.3. The spectrum includes poles associated with transverse component of stress tensor (represented by dots) as well as those associated with transverse component of charge current (represented by stars). We argued that there would exist only one branch of Puiseux series passing through (0, 0), namely the shear mode.
In this subsection we investigate how the radius of convergence of the derivative expansion for shear mode changes with Q. To this end, we numerically find the spin 1 spectrum of quasinormal modes at complex momenta. The result has been given in Fig.12. The smoothness of plot suggests that for the whole range of 0 ≤ Q ≤ 0.8, only one specific type of level-crossing corresponds to the critical point of w shear (q 2 ).
To be more precise, in Fig.12, we have shown the situation of complexified poles at several values of |q 2 | associated with Q = 0.3. As it is seen, by increasing |q 2 | firstly the highest star gapped poles collide. It does actually occur at some critical value of 0.25 < |q 2 c | < 0.75. By further increasing (the top right panel and bottom left panel), trajectories of highest dot gapped poles, which are still closed, come close to that of the diffusion pole, namely the single dot pole on imaginary axis. The first critical point of spectral curve on the shear branch of Puiseux series, or equivalently the first collision of dot poles, turns out to be at |q 2 c | ∼ 2.70. We have demonstrated the situation of poles just before and just after the collision in the bottom middle and bottom right panels,  Figure 11: Radius of convergence of the derivative expansion versus Q in the spin 1 channel. As Q increases, the domain of convergence of w shear monotonically increases. The intersection point with the vertical axis, related to the neutral fluid case, was found in [24]. Figure 12: Poles of the retarded two-point function in the spin 0 channel at Q = 0.3, in the complex w−plane, at various values of the complexified momentum q 2 = |q 2 |e iθ . Large dots and large stars correspond to the poles with purely real momentum (i.e. at θ = 0). As θ increases from 0 to 2π, each pole moves counter-clockwise following the trajectory whose color changes continuously from blue to red. At respectively. We then conclude that at Q = 0.3, the radius of convergence of the derivative expansion for the dispersion relation of shear mode is |q c | ∼ (2.70) 1/2 ∼ 1.64.
We have checked that the spectrum of poles and specifically the collision points for all other values of Q, in the range 0 ≤ Q ≤ 0.8, are qualitatively the same as what we saw in the case of Q = 0.3. Thus in contrast to spin 0 channel, in the current channel the level-crossing corresponding to radius of convergence of the derivative expansion is only of one type. Since colliding poles on the shear branch both belong to spectrum of (vector components of) stress tensor fluctuations, we call such type of crossing point as the T T −crossing.
Let us briefly compare our result about the shear mode with that of associated with AdS 4 RN model [22]. As mentioned in the Introduction, the latter reference finds the radius of convergence of the shear mode to be inversely proportional with µ, resulting in an infinite radius of convergence at µ = 0. This is, however, in contradiction with explicit computations of [24], as well as with our finding in Fig. 12. Our result shows that at least for the range of µ/T studied in this paper, radius of convergence of the shear mode increases with increase of µ/T . Additionally as has been pointed out to in the figure, at µ = 0 we find a finite radius of convergence in full agreement with [24].
Before ending this section we briefly review the results in the following table.
Spin Type of mode range of Q type of level-crossing type of colliding poles  Table 1: We have classified all different types of level-crossing associated with critical points of hydrodynamic spectral curves, in terms of spin of mode, value of Q and type of colliding poles. For example, the first row of the table is saying that sound mode is an eigen frequency for the spectral curve of spin 0 fluctuations. When 0 ≤ Q ≤ 0.386, derivative expansion associated with dispersion relation of sound converges; the radius of convergence is identified with collision of a hydro mode, actually the sound pole, belonging to spectrum of stress tensor T , and a non-hydro mode, which is actually gapped, belonging to the same spectrum, namely T . It should be noted that by use of (4.7), one can rewrite the third column data in terms of µ/T as well., Let us also note that we did not study spin 2 channel in this section. The reason is that its corresponding spectral curve does not have any branch passing through the origin (see Fig.5). The latter is another way of saying that there is no any spin 2 hydrodynamic mode.

Quasinormal modes and the chaos point
Earlier than establishing the relation of quasinormal modes at complex momenta with convergence of the derivative expansion, they were found to be related with quantum chaos. The first observation of such relation turns to ref. [36]. Before recalling the main idea of the latter reference, let us recall that the exponential decrease of the out-of-timeorder-correlator (OTOC) in a thermal large-N system is the manifestation of quantum chaos [44][45][46][47][48][49][50][51][52][53]. Here V and W are two generic few-body operators, λ is the quantum Lyapunov exponent and t * is scrambling time. v B is the butterfly velocity. Equation (6.1) can be also written in the form of a plane wave with purely imaginary values for both momentum and frequency The point (ω ch , q ch ) is called the chaos point 11 . In holography, (6.1) is mapped onto the the back reaction of an small amount of energy thrown towards the horizon of a two-sided black hole [45]. The resultant deformed geometry is described by a Dray-'t Hooft shock wave [54]. But as has been explained in [55], the same geometry can be found from linearized Einstein equations too. In addition, shock wave solution deals with dynamics of energy-momentum in the boundary theory. Then one may conclude that the same information about quantum chaos and scrambling can be extracted from studying the linear perturbations of metric around the horizon. It was actually shown to be the case firstly in [36] and then in [38]. The argument is that at point (ω * , q * ) = (iλ, iλ/v B ), one specific component of Einstein equations at the horizon becomes trivial. In the dual theory it corresponds to multi-valudeness of energy density response function at (ω * , q * ). This point is called the pole-skipping point [38]. More precisely, this point lies on the analytically continued dispersion relation of sound mode in the boundary theory [36]. The coincidence of (ω ch , q ch ) with (ω * , q * ) can be regraded as the existence of a direct link between hydrodynamics and quantum chaos, at least in holographic systems 12 .
In the following two subsections we investigate on the relation between hydrodynamics and quantum chaos in our holographic model.

Chaos point from shock wave computations
In this subsection we analytically compute the chaos point in a system dual to a AdS 5 RN black brane. We exploit the result of ref. [53]. In [53], based on the shock wave propagation picture in the bulk of AdS 5 [45], the butterfly velocities for an anisotropic Q-lattice has been computed. Considering the background as it was shown that butterfly velocities in longitudinal L and transverses directions T are given as the following In our case h T (r) = h L (r) = r 2 and so m = 6πT /r h . The isotropic butterfly velocity then reads This result is exact in the whole range of Q, including Q = 0 case corresponding to AdS 5 Schwarzschild black brane [45], as well as Q = √ 2 case corresponding to the extremal AdS 5 RN black brane. Considering the fact that all systems with gravity dual saturate the chaos bound [49], chaos point in our system is found to be given by In next subsections, we find pole-skipping points and discuss their relation with (6.7). 12 Pole-skipping has been also derived as a general prediction of effective field theory in maximally chaotic systems in ref. [37]. The pole-skipping has been explicitly shown to happen in 2-dim CFT at large central charge [56] and recently in higher dimensions [57]. See [58] for considering the stringy corrections.

Pole-skipping
As discussed earlier, pole-skipping points of energy density response function can be found from a near horizon analysis of perturbations associated with bulk fields. To this end, it is convenient to work with the ingoing Eddington-Finkelstein coordinates. In this system of coordinates, bulk solutions (4.3) take the following general form where v is the ingoing time coordinate. Because of the relation between energy dynamics and quantum chaos in maximally chaotic systems [36,38], it is natural to look for poleskipping points in spin 0 channel of fluctuations. But as was recently shown, the same information about chaos can be extracted from both spin 1 and spin 2 channels too at least in a holographic system with vanishing µ [25]. Thus for completeness, we study pole-skipping in spin 1 channel, as well.

Spin 0 Channel
By taking the momentum along the third boundary direction, spin 0 perturbations in the bulk are δg vv , δg vr , δg rr , δg vz , δg rz , δg xx , δg yy , δg zz , δA r , δA v and δA z . Following [38], we expand E vv component of Einstein equations around the horizon r = r h ; to linear order in perturbations, we arrive at Using metric functions given around (4.3), equation (6.11) at the horizon r h = R becomes It is clear that this equation trivially holds at (6.13) Therefore at exactly the above points the rank of matrix E (0) µν decreases by one. The latter is equivalent to say that line of poles of energy density response function in the boundary suddenly skips at (6.13). Thus the point (6.13) is nothing but the pole-skipping point of energy density response function 13 . We now confirm it numerically.
More concretely, according to [37], the above-mentioned pole-skipping points lie on the lines of hydrodynamic poles associated with energy density response function. Let us recall that energy dynamics is related to spin 0 channel. Hydrodynamic poles in this channel correspond to sound and diffusion modes. However, according to (5.3), these modes become purely imaginary at purely imaginary values of q. One then takes a fixed value of Q and numerically finds frequencies of these modes at several purely imaginary momenta other than the momentum of pole-skipping points. Then interpolation between the resultant points gives the dispersion relations we look for.
In Fig.13, we have shown the results related to Q = 0.4. In the left panel, we have shown that how pole-skipping points (6.14) lie on the sound branches. The same can be done for other values of Q, with numerical difficulties being more or less than the shown case. As a result, we understand that in our system each of the chaos points (6.7) does always coincide with a point at which pole of hydrodynamic sound mode is skipped. This result completes our discussion about the hydrodynamic origin of quantum chaos in a holographic system at finite chemical potential. In right panel of Fig.13, we have compared our results with those associated with AdS 5 Schwarzschild case with Q = 0 [25]. At this point one may ask: Does pole-skipping happen within the domain of validity of the hydrodynamic approximation? To answer this question, in Fig.14, we compare the absolute value of momentum at pole-skipping point, namely |q * |, with the radius of convergence of the hydrodynamic dispersion relation |q c |. It is seen that when Q exceeds the critical value Q c ≈ 0.668, |q * | lies outside of the domain of convergence of the derivative expansion. Is it physically reasonable?
Recalling |q * | = |q ch |, one notices that when Q > Q c the chaos point lies outside the regime of validity of hydrodynamic derivative expansion too. This is exactly one of reasons for which a quantum theory of hydrodynamics was constructed in ref. [37]. As explicitly mentioned by its authors, to capture the exponential behavior of (6.1), they constructed an effective hydrodynamic theory non-perturbativley in derivatives 14 . In this sense, our result shown in Fig.14 is reliable. Then the simple understanding of it is that when Q < Q c derivative expansion is sufficient to see pole-skipping while for the range Q > Q c a nonperturbative treatment is required to find it from energy density response function [37].
A very explicit example of pole-skipping outside the regime of convergence of hydrodynamics was found in the self-dual graviton-axion model [25]. While the derivative expansion associated with dispersion relation of gapless mode in this theory diverges at the pole-skipping point, it was shown that the Borel resummation of the series converges at the same point. Such Borel resummation is in fact one kind of non-perturbative methods which we already argued should be used in the range Q > Q c 15 .

Spin 1 Channel
Spin 1 perturbations in the bulk are δg vx , δg rx , δg rr , δg zx , δA r and δA x . Now the question is if there exists a combination of spin 1 components of dynamical equations at the horizon, namely E x = 0, which trivially holds at (6.13). In the vanishing µ limit, Maxwell's equations decouple and one finds E xz to be the expected combination [25]. This expression becomes trivial at (−w * , ∓iq * ).
14 At this point it is worth mentioning that by "hydrodynamic origin of quantum chaos" posed in [37], they actually mean a quantum hydrodynamic origin. However, what we have dealt with in this paper have been about classical hydrodynamics. 15 We thank Sašo Grozdanov for bringing this example to our attention.  Figure 14: Comparison between the convergence radius of the derivative expansion associated with hydrodynamic sound mode (denoted by blue curve, given first time in Fig.7) and the absolute value of momentum corresponding to chaos points (denoted by green curve). When Q < Q c , pole-skipping and consequently chaos points lie within the domain of convergence of classical hydrodynamics. Thus they can be found perturbatively in a derivative expansion from energy density response function. However, when Q > Q c , extracting any information about quantum chaos from energy density response function requires non-perturbative computations [37].
In our case with µ = 0, however, Maxwell's equations are coupled to Einstein equations. One finds that the simplest possible combination close to what we want is (6.15) Although at Q = 0 this gives the decoupled equation E (0) vx + i 2/3E (0) xz = 0, at a general value of Q the left hand side combination is still coupled to metric perturbations of the right hand side 16 . Thus at Q = 0 we do not expect to see any pole-skipping in this channe, at least in the upper half of the complex frequency plane 17 . We do not prove this statement explicitly but we just show that (−w * , ∓iq * ) do not correspond to pole-skipping of shear mode. In Fig.15 it has been demonstrated for Q = 0.5. The point does not lie on the dispersion relation of shear mode associated with Q = 0.5. Again we emphasize that it is not in contradiction with arguments of [37]. Even if the 16 One another decoupling case is Q = √ 2, i.e. the extremal case, in which E (0) xz = 0 decouples from the rest of equations. 17 We thank Sašo Grozdanov for pointing this out to us.  Figure 15: In a neutral fluid, Q = 0, hydrodynamic branch of poles in spin 1 channel has a pole-skipping point denoted by blue dot. This point turns out to be (−w * , ∓iq * ) and so includes the same information about the chaos point (w * , q * ). However, when Q = 0, it will no longer be the case. As an example, black dot denotes the point (−w * , ∓iq * ) evaluated at Q = 0.5; it obviously does not lie on the line of poles associated with Q = 0.5. point (6.16) was really a pole-skipping point in spin 1 channel, it would not have anything to do to with "hydrodynamic origin of quantum chaos". The reason is that the points (6.16) would lie in the lower half of complex frequency plane while those being responsible for exponential growth of OTOC have Im w > 0 and are located in the upper half plane. Before ending this section let us note that we intentionally do not study possible poleskipping points in spin 2 channel. Looking at eq.(4.53) in ref. [25], we understand that poles-skipping in this channel cannot happen in the upper half plane. As a result, existence or non-existence of such points is not related to quantum chaos at all and so will be unrelated to scope of our discussion.

Review, Conclusion and Outlook
In this paper we explored three aspects of quasinormal modes in a holographic system at finite chemical potential: 1 The dependence of the quasinormal mode spectra and the hydrodynamic excitations on the value of µ/T (or equivalently on Q), at different channels of spin. 2 The collisions of hydrodynamic and non-hydrodynamic poles at complex momenta, at different channels of spin, and specifically changes in pattern of collisions when µ/T varies. 3 The relational between pole-skipping of energy density response function and hydrodynamics in the system and also extraction of chaos information from energy dynamics by use of derivative expansion.
Let us emphasize that the necessary tool for studying all above aspects was initially provided in § 3. In fact besides all analytical and numerical methods and techniques used in the paper, it was equation (3.11) by use of which we were able to produce most of the results, i.e. those associated with spin 0 and 1 channels.
In part 1 , among all other results we think that finding analytic expressions for hydrodynamic excitations is of more importance. To best of our knowledge, hydrodynamic modes (4.22) and (4.28) were never found in the literature.
In part 2 , our main results are definitely related to spin 0 channel. We have shown that convergence radius of the derivative expansion for sound mode non-trivially depend on the value of µ/T . In vanishing µ limit, it has been shown that convergence of the derivative expansion both in momentum [22,24] and in position space [23] corresponds to the collision of hydrodynamic and non-hydrodynamic modes. However, for µ = 0 we have found that at a specific range of µ/T , this is the collision between two hydrodynamic modes, namely sound and diffusion modes, which determines the radius of convergence for the sound dispersion relation. Then it would be interesting to investigate how to extract information about the non-hydrodynamic sectors from the large order behavior of derivative expansion, in the mentioned range. We leave study on this issue to a future work.
In part 3 , we provided a new evidence for the hydrodynamic origin of quantum chaos. As we showed, for all values of µ/T , the chaos points found from shock wave computations precisely coincide with the pole-skipping points of energy density response function. However comparison between momentum of pole-skipping points with the radius of converegnece of sound mode leads to an interesting result. Let us recall that quantum chaos seems not to be captured by classical hydrodynamics. The reason is that hydrodynamics is applicable for time intervals ∆t much larger than 1/T , while in a maximally chaotic system the exponential behavior in (6.2) [49] can be captured when ∆t ∼ 1/T [37]. It shows why a to-all-order non-perturbative theory of hydrodynamics is needed to describe the chaos. On the other hand, we have found that in the range 0 ≤ µ/T 2.1, chaos information can be extracted perturbatively in the standard hydrodynamic derivative expansion. This results is simply the consequence of our another result saying that in the mentioned range, pole-skipping happens in the domain of convergence of the hydrodynamic derivative expansion.
It should be noted that the discussion of previous paragraph associated with part 3 , implicitly confirms the result shown in Fig.7 of part 2 . Let us assume that the blue curve in that figure was monotonically increasing by Q; analogous of what actually was found to be the case for diffusion and shear modes. As a result, in Fig.14, the green curve would totally be located below the blue one. Since increase in Q at a fixed µ is equivalent to a decrease in T , then one would conclude that even as T → 0, where the derivative expansion definitely breaks down, the chaos point could be still found perturbatively in a derivative expansion. This is absolutely wrong. Consequently, quantum nature of chaos forces the blue curve in Fig.7 not to be always upper than the green one.
Based on the above discussion, it becomes more important to study the quasinormal modes at higher values of Q than those studied in this paper [43], specifically in low temperature regime which is in correspondence with a near extremal black hole [71].
And finally in the tensor channel:

B Frobenius solution
In each channel, to find the quasinormal modes, one has to search for the solutions for the functions Z and E obeying two specific coupled second order differential equations. Following the explanations given in [8], the solutions obeying the ingoing boundary condition at the horizon can be represented as Substituting the above two solutions into the corresponding coupled differential equations, one finds the coefficients a n and b n of the series expansion, all in terms of a 0 and b 0 . Then by imposing the Dirichlet boundary condition at u = 0, one arrives at a system of coupled Algebraic equations which can be formally written as C 11 (w, q)a 0 + C 12 (w, q)b 0 =0 Quasinormal spectrum then can be determined by demanding the above equations have non-trivial solutions for (w, q): C 11 (w, q) C 22 (w, q) − C 12 (w, q), C 21 (w, q) = 0 (B.4) and solving Eq. (B.4) numerically; the latter is done by taking a sufficiently large but finite number of terms, n, in the sum.

C Comparison with explicit hydrodynamic computations
In this section we reproduce the hydrodynamic excitations found from the study of perturbations, given in (4.22) and (4.28), via explicit hydrodynamic computations on the boundary. The hydrodynamic regime of a holographic charged fluid has been studied in the context of fluid/gravity duality [12,13] and also in [16]. Following [13] we write the associated constitutive relations up to first order in derivative expansion, in 4-dimensions, as the following T µν = p(η µν + 4u µ u ν ) − 2ησ µν j µ = n u µ − D P ν µ ∂ ν n + 3(u λ ∂ λ u µ )n (C.1) with P µν = η µν + u µ u ν . The coefficients found in [12,13] can be written as