Many-body chaos and energy dynamics in holography

Recent developments have indicated that in addition to out-of-time ordered correlation functions (OTOCs), quantum chaos also has a sharp manifestation in the thermal energy density two-point functions, at least for maximally chaotic systems. The manifestation, referred to as pole-skipping, concerns the analytic behaviour of energy density two-point functions around a special point $\omega = i \lambda$, $k = i \lambda/v_B$ in the complex frequency and momentum plane. Here $\lambda$ and $v_B$ are the Lyapunov exponent and butterfly velocity characterising quantum chaos. In this paper we provide an argument that the phenomenon of pole-skipping is universal for general finite temperature systems dual to Einstein gravity coupled to matter. In doing so we uncover a surprising universal feature of the linearised Einstein equations around a static black hole geometry. We also study analytically a holographic axion model where all of the features of our general argument as well as the pole-skipping phenomenon can be verified in detail.


I. INTRODUCTION
Over the last few years there has been exciting progress in characterizing chaotic behavior in quantum many-body systems using out-of-time ordered correlation functions (OTOCs) [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15][16][17][18][19] For instance, for a large class of systems, it has been observed that 1 (1.1) where V and W are generic few-body operators, β 0 = 1/T is the inverse temperature, and is a small parameter inversely proportional to the number of degrees of freedom. The exponential growth of (1.1) is reminiscent of the diverging trajectories of two initially infinitesimally separated particles in classical chaotic systems. Thus λ is often referred to as the quantum Lyapunov exponent, and v B , which describes the speed at which the growth propagates in space, as the butterfly velocity. For later purposes, let us note that the equation (1.1) has the form of a plane wave with purely imaginary values of both frequency ω and momentum k: (1. 3) The quantum nature of (1.1) is highlighted by the fact that λ has an upper bound λ ≤ λ max = 2πk B β 0 [5] (henceforth = k B = 1), which is saturated by a variety of systems including holographic theories, two-dimensional conformal field theories (CFTs) in the large central charge limit, and strongly coupled SYK models. Below we will refer to systems which saturate the bound as maximally chaotic systems.
In particular, two recent developments showed that in addition to (1.1), quantum chaos also has a sharp manifestation in the thermal energy density two-point functions [26,27]. The manifestation, referred to as pole-skipping in [27], was first observed numerically in a holographic system in [26], and then derived as a general prediction of an effective field theory (EFT) for chaos [27]. More explicitly, it was proposed in [27] that the chaotic behavior (1.1) can be captured by the propagation of an effective chaos mode, which, at least for maximally chaotic systems, coincides precisely with the hydrodynamic mode for energy conservation. The phenomenon of pole-skipping is a direct consequence of the same mode playing the dual role of capturing both energy conservation and chaos. The EFT description of chaos also provides a simple explanation for connections between (1.1) and energy diffusion observed in previous literature [13,14,[22][23][24][25][26][27]; they correspond to the behaviour of a single effective mode, albeit at different scales: the OTOC (1.1) at the scale given by (1. 3), and transport at scales of ω, k 1/β 0 . They are related by an O(1) extrapolation.
Pole-skipping can be checked to hold for the SYK chains studied in [13]. And it has also been verified in two-dimensional CFTs in the large central charge limit [32].
The purpose of this paper is to provide further support for the generality of the manifestation of quantum chaos in energy density correlation functions. We show that the pole-skipping phenomenon is universal in general holographic systems described by the Einstein gravity coupled to matter fields, thereby generalising the analysis of [26] for the Schwarzschild black hole in AdS 5 . We will provide a general analytic derivation of this phenomenon and elucidate its gravitational origin. We will also make a connection with the gravitational shock wave analysis of (1.1), thus directly establishing that the behaviour of an OTOC and pole-skipping have the same gravitational origin. In doing so we will uncover a surprising universal feature of the linearised Einstein's equations coupled to matter around a static geometry. A hint of this universal feature came from the analysis of [26] which initiated the study of linearised gravitational perturbations at the special point (1.3). We emphasise that for a general gravity system there is no near-horizon SL(2, R) symmetry (as there is in AdS 2 or the BTZ black hole). This makes it clear that pole-skipping is a phenomenon that occurs more generally than just in highly symmetric cases.
Before summarising the main contents of the paper, let us first review the phenomenon of pole-skipping in the energy density two-point functions. The simplest example which exhibits this phenomenon is the SYK chain of [13] for which the energy density retarded two-point function has the form (1. 4) where C is some constant and D E is the energy diffusion constant. zero, so at the special point the would-be pole skips. It has been argued in [27] that this is a generic phenomenon for maximally chaotic systems. More generally, writing (1. 5) it follows from energy conservation that G R should exhibit hydrodynamic poles at small ω, k, i.e. a(ω, k) should have a line of zeros at (1. 6) where in the second equality we have performed a small k expansion with the upper (lower) line for a system with (without) momentum conservation. Now the statement of pole-skipping is: 1. When analytically continued to imaginary ω and k, (1.6) passes precisely through the special point (1.3); 2. b(ω, k) has a line of zeros which also passes through the special point (1.3).
The same phenomenon is also present in G R T 00 T 0i , G R T 0i T 0i as they are related by Ward identities. For the rest of paper we will simply focus on G R T 00 T 00 . An immediate consequence of having zeros of a(ω, k) and b(ω, k) passing through the same point (1.3) is that there, G R T 00 T 00 = 0 0 , and thus the energy density two-point does not have a well defined value at that specific ω and k! More explicitly, let us consider ω = iλ + δω, k = ik 0 + δk . (1.7) Then, as δω, δk → 0, we find G R T 00 T 00 (ω, k) = ∂ ω b(iλ, ik 0 ) δω δk + ∂ k b(iλ, ik 0 ) ∂ ω a(iλ, ik 0 ) δω δk + ∂ k a(iλ, ik 0 ) .
(1. 8) In other words, the function becomes infinitely multiple-valued at (1. 3), depending on the slope δω/δk at which the point is approached. Thus if the gravity analysis is to reproduce the pole-skipping phenomenon, something special must be happening to the linearised Einstein's equations at (1.3). We will see this indeed to be the case.
We now summarise the main results of the paper on holographic systems. To illustrate the general discussion of Section II, in Section III we proceed to study this phenomenon in detail in a specific holographic model which describes a class of three-dimensional strongly coupled field theories with broken translational symmetries [33].

II. NEAR-HORIZON EINSTEIN'S EQUATIONS AND POLE-SKIPPING
In this section, we will demonstrate a remarkably universal property of the linearised Einstein's equations coupled to general matter content. We will then use this property to argue that in general, the retarded energy density two-point correlator of holographic theories with Einstein gravity exhibits pole-skipping at the location of Eq. (1.3). This property also enables us to make a connection with the gravitational shock wave analysis [2][3][4] of (1.1), thus establishing that the behaviour of the OTOC and pole-skipping have the same gravitational origin.
The discussion of this section is general, thus somewhat abstract. In the next section we will then examine a family of examples very explicitly.

A. Setup
Let us start by considering Einstein's equations that arise from a bulk action of the form where Λ = −d(d+1)/2L 2 , L M is the matter Lagrangian and L is the AdS radius that henceforth we set to 1. We will allow for a rather general matter Lagrangian L M , but will assume that the theory admits an equilibrium configuration corresponding to a homogeneous and isotropic black brane. We write the background metric as where f (r) is the emblackening factor which vanishes at the horizon f (r 0 ) = 0. The Hawking temperature is given by T = r 2 0 f (r 0 )/(4π).
Note that for a metric of this form the special frequency and momentum in (1.3) are extracted by constructing the Dray-'t Hooft shock wave describing a gravitational perturbation where β 0 = 1/T is the inverse temperature. Solving this equation gives an exponential profile for c(x) that leads to the form (1.1) with the Lyapunov exponent and butterfly velocity given Here we wish to study the retarded energy density correlation function G R T 00 T 00 (ω, k) near (2.4). This correlation function can be extracted from solving the linearised gravitational perturbation equations around (2.2) subject to ingoing boundary conditions at the horizon. It is therefore convenient to introduce ingoing Eddington-Finkelstein (EF) coordinates (v, r, x i ) with , (2.5) in terms of which the background metric (2.2) is To calculate the energy density correlation function we then need to study the perturbation equations of the metric δg vv (r, v, x) = δg vv (r)e −iωv+ikx , together with all the other metric perturbations and matter fields that couple to this mode. The fields that couple in this channel are δg vv , δg rr , δg vx , δg vr , δg x i x i , δg rx and δϕ, where δϕ schematically represents any matter fields that couple to these perturbations. 2

B. Near-horizon expansion
The retarded Green's function of energy density G R T 00 T 00 is governed by solutions of the gravitational and matter perturbation equations that are regular at the horizon in ingoing EF coordinates. It is therefore useful to expand these perturbations near the horizon as where T µν (r 0 ) is the bulk stress-energy tensor of the background matter fields supporting the black brane (2.2) and δT µν (r 0 ) describes the matter perturbations. A priori, (2.8) therefore depends on the matter content of the theory, i.e. the precise form of L M in (2.1). However, for a large class of black brane solutions that we examined, including AdS gravity coupled to scalar and gauge fields, we find that the identity holds automatically for any value of ω and k as a consequence of regularity at the horizon.
In particular, in Appendix A we show that the identity (2.9) holds for the commonly studied Einstein-Maxwell-Dilaton-Axion gravity theories. As such, in all these theories the vv component of Einstein's equations at the horizon takes a remarkably universal form that depends only on the metric perturbations x i x i of the near-horizon solution to each other. However, when ω = 2πT i then all other fields decouple at the horizon from δg    x i x i . In other words, precisely at (1.3) there is one fewer equation to solve at the horizon than at a generic value of ω and k.
As a consequence we can deduce that at the location (1.3) there exists an extra linearly independent ingoing solution to Einstein's equations. For the specific cases of AdS 4 gravity or the axion model that we study in detail in Section III we can check this by explicitly constructing the solutions (2.7) by solving the equations of motion order-by-order in the expansion (2.7). For instance, in AdS 4 gravity then, after fixing radial gauge, we find that the most general ingoing solution to Einstein's equations is specified at generic ω and k by 4 parameters in this nearhorizon expansion. These can be chosen to be the parameters δg  With a slight risk of oversimplifying the problem, but in aid of conceptual clarity, let us use a single scalar field as an illustration of the implications this extra solution. Recall that a scalar field φ has the following expansion near the boundary r → ∞: where the A-term (B-term) is the non-normalisable (normalisable) term and corresponds to the external source (expectation value). The dual retarded Green's function is then given by G R = B/A with φ obeying ingoing boundary conditions at the horizon [36]. For generic ω and k, at the horizon, φ has an ingoing and an outgoing mode. When we choose the ingoing mode, the ratio of B over A becomes completely fixed, and so does the retarded Green's function.
Suppose now that at the special point (1.3), we have an extra ingoing mode at the horizon, which means that both independent solutions of φ are now allowed at the horizon. Then, there is no constraints on A or B, and the retarded Green's function becomes infinitely multiplevalued, depending on one free parameter which determines the linear combination of two nearhorizon solutions. In particular, we can always choose a combination so that there is only the normalisable term near the boundary (i.e. A = 0), which then leads to a pole in G R . Similarly, there exists another combination such that the normalisable term is absent (i.e. B = 0), which then corresponds to a zero of G R .
The present situation with metric perturbations is a bit more complicated, but the essence is the same: the extra freedom in δg vv at the horizon should generically lead to one free parameter in G R T 00 T 00 at the special point. We will shortly see that when we move slightly away from the special point, the free parameter can be taken to be the slope δω/δk approaching the point, thereby precisely mirroring the situation of equation (1.8). In fact, after fixing all the bulk gauge freedom and solving the constraints, the sector of metric perturbations associated with δg vv always reduces to a single scalar degree of freedom [37]. On the boundary theory side this can be understood as a result of Ward identities for G R T 00 T 00 , G R T 00 T 0x , G R T 0x T 0x , whereby all of these two-point functions are proportional to a single scalar function. We will see this explicitly in the example studied in the next section.
To close this subsection, we note that the fact that at this special point the metric component vv can be tuned independently of the other fields at the horizon resonates with the analysis of [26] which found a special null ansatz solution to Einstein's equations at (1.3). The solution constructed in [26] is indeed a special case of our general ingoing solutions corresponding to the case with δg (0) Note that even though we have only moved slightly away from the special point, Einstein's equation (2.13) now imposes a non-trivial constraint on the near-horizon expansion that must be satisfied by the ingoing mode (in addition to those coming from the smooth limits of X = 0).
As such for a fixed δω/δk there is no longer an extra solution. However, the singular nature of (1.3) is reflected in the fact that the constraint (2.13) depends on the slope δω/δk with which we move away from (1.3). This means that near (1.3) then, even after fixing the usual independent parameters in the near-horizon expansion, we have a family of different ingoing modes parameterised by δω/δk. The slope δω/δk now acts as an extra parameter in the nearhorizon solution and can be tuned to match the ingoing solution onto different asymptotic solutions at the boundary.
In particular, by choosing δω/δk appropriately, we can ensure that we have an ingoing mode that matches continuously onto the normalisable solution at the boundary. This can be achieved by choosing δω/δk to satisfy where δg (0) µν corresponds to the near-horizon expansion of the normalisable solution. We therefore deduce that if we move away from the special point (1.3) along the slope (2.14), then we will see a line of poles in the energy density correlation function that passes through (1.3) (as well as in the correlation functions of the operators dual to the fields that couple to δg vv ). Furthermore, given knowledge of the normalisable solution of the perturbation equations then (2.14) gives a prediction for the slope of this line of poles in terms of the metric components of this mode. Note that we could also choose a different slope so that the ingoing mode matches onto a solution with different asymptotics at the boundary. In particular one can choose the slope to instead match to a solution with the expectation value T 00 (ω, k) β 0 = 0, from which we also deduce the existence of a line of zeroes in the dual Green's functions passing through (1.3), with a different slope that is again related to the metric components of this solution through the same equation (2.14).
Whilst the above discussion gives a simple gravitational argument for the phenomenon of pole-skipping for a broad class of black brane solutions, it is somewhat abstract. In particular, although the existence of an extra parameter in the ingoing solution should generically result in an extra parameter in the asymptotic behaviour of δg vv , the precise dependence of the asymptotic δg vv on the extra parameter (i.e. the slope δω/δk) depends on the details of the specific gravitational theory. To show how the argument works more directly, in the next section, we examine in detail a holographic model in which we can analytically study the energy density two-point Green's function. We will see explicitly all the features mentioned above, and confirm the prediction (2.14). Note that the above discussion applies in the vicinity of Eq. (1.3), thus only showing that a line of poles passes through (1.3). It is not clear that the line of poles actually goes over to those that correspond to hydrodynamic excitations at small ω and k. This will be checked by using explicit examples in the next section.

III. ENERGY DENSITY GREEN'S FUNCTION NEAR SPECIAL POINT
In this section, we study the connection between chaos and the energy density correlator in a specific holographic model [33,34], which has a free parameter m controlling the strength of momentum dissipation. We choose to study this example because it provides a simple and soluble model, which nevertheless exhibits a wealth of transport behaviour as one varies m, thus allowing us to demonstrate explicitly the insensitivity of the pole-skipping phenomenon to long distance transport properties of a system. Remarkably, in this model, we will be able to analytically compute the retarded energy density correlator near the special point (1.3) in Fourier space. This will allow us to demonstrate explicitly how the general considerations of Section II are borne out in a specific example. In particular, we will be able to prove that both the line of poles and the line of zeroes of G R T 00 T 00 (ω, k) always pass through (1.3) in these models. Moreover, we will be able to derive an analytic expression for the slope of the line of poles (of the dispersion relation ω(k)) as it passes through the point (1.3), which we find to be in perfect agreement with a numerical calculation of the same quantity. In Appendix B, we show that the same analytic discussion can also be applied to the Schwarzschild black hole of AdS 5 , dual to a four-dimensional CFT at a finite temperature and studied initially in [26].
Here, we focus on the 2 + 1 dimensional boundary theories dual to 3 + 1 dimensional gravity coupled to massless scalar ('axion') fields with action In ingoing EF coordinates, this action has the black brane solution [33] When m = 0, (3.2) is simply the Schwarzschild-AdS 4 solution and our results in this limit are identical to those that would be found in the absence of matter. Increasing m causes the scalar fields to backreact on the metric, and in the extreme limit m/T → ∞ the metric near the horizon becomes AdS 2 × R 2 .
The parameter m sets the strength of the source of a scalar operator φ (0) i = mx i that explicitly breaks the translational symmetry of the dual field theory. This symmetry breaking radically alters how energy is transported over long distances, resulting in a much richer variety of energy dynamics than in the translationally invariant case [34].
As we have explained, our goal is to study the energy density correlator of these theories near and depends explicitly on the dimensionless parameter m/T .
We emphasise that all of our results hold even in the translationally invariant m = 0 case.
We allow for m = 0 simply to illustrate the generality of the pole-skipping phenomenon, and its insensitivity to the long distance transport properties of the theory. Analogous results to those we will present can also be found in higher dimensions, and in Appendix B we outline how to obtain these for the Schwarzschild-AdS 5 solution. This solution was studied numerically in [26], and our analysis provides an explanation of the results found there.

A. Master field perturbation equations
In order to calculate the energy density correlation function we need to study the linearised gravitational perturbation equations about the spacetime (3.2). In a general theory, studying these linearised perturbation equations is rather complicated, since at non-zero frequency ω and momentum k there are many coupled fields. For instance, if one aligns the momentum in the x direction then for this axion model, the correlator G R T 00 T 00 (ω, k) can be extracted from solving the equations of motion for the 'longitudinal' perturbations {δg vv , δg xx , δg yy , δg rr , δg vx , δg xr , δg vr , δϕ 1 } after Fourier transforming δg µν = δg µν (r)e −iωv+ikx .
Fortunately, in the simple model (3.1) the dynamical equations for the perturbations can be decoupled by working with suitable gauge-invariant 'master field' variables [34]. In particular, the energy density correlation function is simply related to the dynamics of the following variable which obeys the equation of motion d dr where primes denote derivatives with respect to r and we have defined The retarded two-point Green's function of the energy density T 00 is determined by solving (3.5) subject to ingoing boundary conditions at the black hole horizon. From the near-boundary expansion of the solution ψ = ψ (0) + ψ (1) r −1 + . . . one can then read off the retarded Green's function using the usual holographic dictionary as up to contact terms.
As motivated earlier, we are particularly interested in demonstrating that there is always a line of poles ω(k) in (3.7) that passes through (1.3). Since we are working in terms of the master field (3.4), which is related to metric components in a complicated way, the usual notion of normalisable and non-normalisable modes is a bit subtle in terms of ψ. However, we can see that (3.7) will have a pole when there is an ingoing solution with ψ (0) (ω, k) = 0 and ψ (1) (ω, k) + iωψ (0) (ω, k) = 0. It can be checked that a solution with these properties corresponds to a normalisable metric perturbation, and we will refer to it as a normalisable mode ψ n .
Likewise, there will be a line of zeroes us consider the near-horizon behaviour of (3.5) at generic k. Then taking the near-horizon limit of (3.5) one finds the following power law solutions to (3.5): At a generic value of ω, k there is therefore a single regular solution in ingoing coordinates, is the corresponding outgoing solution.
However, we can see from (3.5) that at the special value k 2 = −k 2 0 the near-horizon behaviour of (3.5) changes since (k 2 + r 3 f ) now vanishes at the horizon. Taking the near-horizon limit of (1.3) one now instead finds power law solutions of the form The indices are shifted by one at this special value of k. For a generic ω it is clear that there is still one regular solution in ingoing coordinates corresponding to the (r − r 0 ) power law.
However, one can now see that something special occurs when we also set iω = −2πT in (3.9).
In this case we have power law solutions As such, precisely at the location (1. 3) it appears that both near-horizon solutions for ψ are regular in ingoing coordinates.
The existence of two regular solutions for ψ at this point can be seen more directly from the fact that in this simple example we can exactly solve (3.5) at the special point (1.3). At this location Ω(ω, k) = 0 and hence (3.5) simplifies considerably to d dr This is a first-order equation for ψ and so it is straightforward to find the general solution, Note that the integral in (3.12) is regular as r → r 0 . We therefore manifestly have regular solutions near the horizon for any values of c 1 , c 2 , and these solutions can be expanded as We note that the considerable simplification Ω(ω, k) = 0 that occurs at the location (1.3) does not generalise to higher dimensions. It is this simplification that allowed us to write down the exact expression (3.12) for the general solution for ψ. However, we emphasize that a simplification like this is not necessary to realise the key property that there are two independent solutions for ψ that are regular at the horizon. In Appendix B, we discuss the case of Schwarzschild-AdS 5 and illustrate this property by constructing regular series solutions for ψ near the horizon.

C. Solutions near special point
In order for the choice of ingoing boundary conditions on ψ to be non-trivial, it is necessary to move slightly away from the location (1.3) in Fourier space. To do this we will consider perturbing a small distance from (1.3) by taking where we have introduced q as a convenient parameterisation of the direction δω/δk in which we move away from (1.3) Note that near the horizon, the function appearing in the equation (3.5) for ψ takes the form Our goal now is to construct the form of this ingoing solution as we take → 0 and approach the point (1.3). We can do this by dividing the radial direction into two regions and performing a matching calculation (see Figure 1). In particular for (r − r 0 ) /b we are free to safely 3 Note that an exception to our analysis is provided by the special case m 2 = 2r 2 0 . In this axion model one finds ignore in (3.5). In this 'UV regime' the solutions are then given by solving (3.5) at ω = iλ and k = ik 0 and simply take the form (3.12).
However, close to the horizon (r − r 0 ) ∼ /b the approximation of ignoring breaks down.
We therefore need to solve separately for the solution in this 'IR regime'. To do this we can consider the scaling and then construct the ingoing solution perturbatively in ψ r (y) = ψ 0 (y) + ψ 1 (y) + . . . .
The regular solution to (3.19) simply corresponds to a constant which we can set to unity c 1 = 1 to normalise our solution. At this order in our expansion the regular solution therefore only has one of the power laws in (3.10). To match to the UV we also need to determine the coefficient of the y 1 term. This requires us to go to next order in our expansion (3.18) and solve for ψ 1 (y). At O( ) we find that the equation of motion takes the form where the forcing term f 1 (y, m 2 , q, r 0 ) is determined by expanding the equation of motion and then inserting the zeroth-order solution ψ 0 = 1. We can then again solve this equation subject to regularity at the horizon to find ψ 1 (y). Expanding this solution at large y we find that the leading behaviour of ψ 1 is where b + (q) is explicitly given by the formula .

(3.23)
After transforming back to the radial coordinate r we then find that the ingoing solution takes the form ψ r (r) = 1 + b + (q)b(r − r 0 ) + . . . , (3.24) in the IR region, where the . . . indicate terms that vanish as → 0. 4 The key point is that as we take → 0 the ingoing mode in Explicitly performing this matching we find that the ingoing solution (3.24) matches to the UV .
whereÑ can be written in terms of the integral of a known function and is defined through (C4) and (C6).
We have focused on describing the matching to ψ n since this establishes the existence of a pole passing through (1.3). However, by moving away from the point (1.3) along a different direction q we could find an ingoing solution that matches to any UV solution. In particular we could also pick a (different) q = q z to use (3.27) to match to the coefficients a 1,nn and a 2,nn of a solution which has no normalisable component. As such there will always also be a line of zeros of the Green's function passing through (1.3) at a slope q = q z . For a general q the UV solution is the linear combination of these modes given by (3.25). In Appendix C we extract the full Green's function from (3.25) and verify it indeed has the form  [26,27] and is referred to as pole-skipping. To compare these discussions we can simply insert the near-horizon expansion where a 1 and a 2 are related to the coefficients of the metric and δϕ 1 in the near-horizon expansion (2.7). The coefficient a 1 has a simple form and can immediately be read off from the near-horizon expansion of the metric as Determining the coefficient a 2 is more complicated as after inserting our expansion (2.7) into  . From this we investigate the behaviour of this same pole at small ω, k and explore the connection with hydrodynamic modes seen in [26,27].
In particular, let us recall that in the translationally invariant (pure gravity) case m = 0, the hydrodynamic limit of the energy density Green's function (3.7) is dominated by gapless sound modes. That is, at small ω and k it exhibits hydrodynamic sound poles with dispersion relations [38] ω(k) = ±v s k − i k 2 8πT + . . . , (3.36) where v s = 1/ √ 2 is the speed of sound and . . . denote higher-order corrections in k that in principle can be computed using higher-order hydrodynamics (see e.g. [39]). As for the analogous AdS 5 theory studied numerically in [26], we observe that the full dispersion relation ω(k) of the mode which approaches ω(k) = −iD E k 2 + . . . , (3.37) where D E is the energy diffusion constant. Numerically tracking the dispersion relation of this pole, we find that for any value of m/T that ω(k) continues to pass through the location (1.3).
We emphasise that this is very non-trivial. There are a family of different dispersion relations ω(k), parameterised by m/T , whose qualitative features change dramatically as m/T is varied (see Figure 2). However, in all cases we find that regardless of the details of the full dispersion relation it always satisfies the constraint ω(ik 0 ) = iλ as can be seen in the left hand panel of and exactly 2 for m/T → ∞. In the right hand panel of Figure 3 we plot this formula (red line) and also the numerical values of the slope (black dots) extracted from the dispersion relation (e.g. those shown in Figure 2). The numerical results agree perfectly with the analytic expression (3.28), validating the details of the matching argument described in Section III.
As can be seen in Figure 2, as m/T is increased, the leading-order hydrodynamic approximations of the dispersion relations ω(k) become a better and better approximation to the exact dispersion relation near the special point (1.3). Using the fact that D E (m/T → ∞) → v 2 B /λ [23], our result that δω δk /v B → 2 in this limit indicates that the hydrodynamic approximation to the dispersion relation (3.37) is exact in the vicinity of the point (1.3) when m/T → ∞. This is consistent with the observation that hydrodynamic collective modes exist even over time scales much shorter than T −1 in holographic theories with an AdS 2 factor in the near-horizon metric [40].  of m 2 = 2r 2 0 the above discussion in terms of the gauge invariant mode ψ is somewhat subtle, as a result of the fact the parameter b controlling our matching calculation vanishes for this choice of m. We perform a detailed analysis of this special case in Appendix D. For the purposes of our discussion in the main text, we simply note that at this value of m/r 0 , the theory described by (3.1) dramatically simplifies and has an enhanced SL(2, R) × SL(2, R) symmetry [34]. This enhanced symmetry allows one to obtain an analytic expression for the Green's function G R T 00 T 00 for any ω, k. Here we will use this expression to demonstrate the pole-skipping phenomenon very explicitly. In particular, for this choice of m 2 = 2r 2 0 the retarded Green's function takes the form 5 (see Ref. [34]) From which we see that the Green's function has an infinite family of poles ω(k) satisfying where n, p are non-negative integers. Expanding the pole with n = 0 gives the hydrodynamic mode ω = −iD E k 2 + . . . , (3.40) with an energy diffusion constant D E = r −1 0 . As in Section III E we can consider tracking this pole as we increase imaginary k. The full dispersion relation can be readily read from (3.39):  Whilst for the purposes of the main text we have simply focused on using (3.38) to illustrate this phenomenon, the enhanced symmetry of the theory at m 2 = 2r 2 0 means we can perform a very detailed analysis of the origin of pole-skipping in this example. As we explain in Appendix D, there are subtleties involved in discussing pole-skipping in terms of the gauge invariant field ψ at this point. Fortunately however this special case is sufficiently simple that we are able to also explicitly discuss this phenomenon in terms of the family of ingoing metric solutions discussed in Section II. In particular it is possible for this choice of m to find exact analytic expressions for the metric solutions at (1.3) everywhere in the bulk, and hence study their UV asymptotics. As we discuss in Appendix D this allows us to confirm the existence of an extra parameter in the ingoing metric solution, and explicitly see that this extra parameter gives rise to the phenomenon of pole-skipping according to the discussion in Section II.

IV. DISCUSSION
In this paper we have shown that in general holographic models dual to Einstein gravity coupled to matter, remarkable signatures of many-body chaos exist in the energy density twopoint correlation functions. A key element for the discussion is the observation that one of the Einstein's equations becomes trivial at the horizon at the special point (1.3) determined by chaos parameters, which leads to a general argument for the phenomenon of pole-skipping [26,27]. We then illustrate how the general argument works in a specific holographic model in great detail.
As emphasised in [27], the phenomenon of pole-skipping can be considered as a "smokinggun" for the fact quantum many-body chaos is tied with energy conservation. The results of this paper give further strong support for this surprising, but likely profound connection (at least for maximally chaotic systems). In particular, from the perspective of the hydrodynamic dispersion relation (1.6), the statement that ω h (k) must pass through (1.3) is a highly nontrivial one. In the Einstein-Axion model we explicitly saw that this happened regardless of the detailed behavior of the full dispersion relation, which could vary dramatically as we changed the ratio m/T . The special point (1.3) lies outside the range of the small ω, k expansion and thus the full dispersion relation ω h (k) is needed to extrapolate to the point. This means that in the small k expansion of the second equality of (1.6) an infinite number of terms (which arise from higher-order hydrodynamics) must conspire for ω h (k) to pass through (1.3). In the proposal of [27], this conspiracy is ensured by an emergent shift symmetry in the hydrodynamic EFT for chaos. Thus, the results of this paper can also be considered as support for this "hidden" shift symmetry in hydrodynamics.
We end with a discussion of various questions for future research.

Immediate generalizations
The general argument we presented in Section II for pole-skipping implies that this phenomenon also occurs for charged black holes and/or for those with additional scalar fields, and it would be interesting to test this explicitly. Furthermore, whilst our analytic arguments show that there is a pole that always passes through (1.3), it was only through the help of numerics that were able to see this pole to always be connected to the dispersion relation of the hydrodynamic mode. It would be interesting to understand if one can argue that this should always be the case from gravity, perhaps by tracking the holographic diffusion poles that were recently constructed using horizon data in [44].

Shift symmetry in hydrodynamics from gravity
As commented above, the confirmation of pole-skipping in general holographic theories (which all have maximal chaos) can be considered as support for the shift symmetry in an all-order quantum hydrodynamics proposed in [27]. It would be extremely interesting to In Section II we saw that the reason that pole-skipping occurred at (1.3) was a consequence of the fact that at ω = iλ the Einstein's equation (2.10) reduced to a decoupled equation vv , that took the same form as the equation determining the shock wave profile (1.3). However, as noted in [26], there is not an immediate connection between such a perturbation and the delta function shock-wave in Kruskal-Szekeres coordinates. We hope to explore this question further in the future.

Higher-derivative gravity theories and stringy corrections
So far, all examples of pole-skipping are for maximally chaotic systems. It is clearly of crucial importance to understand what happens to non-maximally chaotic systems.
For holographic systems, stringy corrections decrease the Lyapunov exponent from being maximal [4], so it is very interesting to see how stringy corrections affect the pole-skipping.
Perhaps a more immediate question is to see whether the phenomenon persists or gets modified when higher-derivative corrections are included. For example, Gauss-Bonnet gravity (see [42,45,46]) may be a good testing ground for this purpose.
More generally, it is also important to study this phenomenon in weakly coupled systems.
A recent discussion [28] of chaos using kinetic theory gives some encouraging indications on the connection between chaos and energy dynamics even at weak coupling.

Theories with weak energy dissipation
Finally, it would be interesting to explore what happens if a system no longer has exact energy conservation. In the holographic context, weak energy dissipation can be introduced by working with massive gravity with an appropriate choice of graviton mass [47].
This, or axion models with time-dependent fields, could provide a nice laboratory to study how the Lyapunov exponent and pole-skipping are affected by energy dissipation.
"Entanglement in Quantum Systems" and the Simons Foundation for partial support during the completion of this work.
Appendix A: Stress-energy tensor in Einstein-Axion-Dilaton Theories Here we wish to demonstrate that the condition (2.9) that we needed to claim we could ignore the matter terms in (2.8) holds in a large class of matter theories. We will consider a general Einstein-Axion-Dilaton matter theory with an action of the form with a matter Lagrangian which has an equilibrium black hole solution (2.6) sourced by the fields with other components of the Maxwell field A µ vanishing.
We wish to demonstrate that for this general class of matter fields (2.9) holds identically for any ω, k. That is For the Lagrangian (A2) the stress-energy tensor is from which we can read off the relevant component T vr of the background stress-energy tensor as We also need to work out δT vv , and hence need to vary the stress-energy tensor (A5) once more with respect to the fields g µν , A µ , φ, ϕ. At leading order this gives where ψ i = {A µ , φ, ϕ i } denotes the matter fields.
Evaluating the vv-component of this for the black hole solution with the matter fields in (A3) gives where we have used that g rv = 1, g rr = −g vv . Now since g vv (r 0 ) = −r 2 0 f (r 0 ) = 0 then assuming the quantity in square brackets is regular at the horizon we immediately see by comparing to (A6) that indeed holds for this class of black holes. In all these theories the vv component of the Einstein equations therefore reduces at the horizon to the universal form (2.10) that depends only on the near-horizon expansion of the metric perturbations.
Appendix B: Generalisation to AdS 5 In the main text we studied the AdS 4 axion model. However, as we have noted, the phenomenon of pole-skipping was first noticed in a gravitational setting in [26], which studied Einstein gravity in AdS 5 . Here we will show that the matching argument presented in Section III can easily be generalised to explain the pole-skipping observed in [26]. The action we consider is and the AdS 5 -Schwarzschild solution of interest to us can be written in ingoing EF coordinates. Note that for this theory the special location (1.3) in Fourier space corresponds to The generalisation of the 'master field' (3.4) to Schwarzschild-AdS 5 is which obeys the equation of motion d dr where We again find that when k 2 = −k 2 0 the two solutions near the horizon are of the form (3.9) and thus at the special location (1.3) the general solution for ψ that is regular at the horizon is of the form (3.10).
Infinitesimally away from the location (1.3), the solution continues to take the form (3.10) away from the horizon with the ratio a 2 /a 1 fixed by the direction in Fourier space in which one moves. To quantify this, we again parameterise the perturbation away from (1.3) by (3.14) and note that near the horizon We first solve in the IR regime r − r 0 ∼ /b of the spacetime by scaling the radial coordinate as in (3.17) and then solving perturbatively for ψ, as in Section III C. After demanding regularity at the horizon, the solution is where . . . indicate terms that vanish as → 0, and b + (q) = q − 5 24r 2 0 . (B9) In the UV region r − r 0 /b, the solution has the form (3.10). Matching (3.10) to (B8) where the solutions overlap tells us that imposing ingoing boundary conditions enforces the relation is equivalent to the condition where δg µν are again the coefficients in the near-horizon expansion of the metric components. Converting q back to δω/δk then gives perfect agreement with the equation (2.14) that relates these coefficients to the direction q in which one moves away from the location (1.3).

Appendix C: Explicit matching to normalisable mode
In Section III C we determined the form of the ingoing solution for ψ (3.24) near (1.3). Here we will extract the form of the Green's function (3.7) and a prediction for the slope of the line of poles passing through (1.3) by explicitly matching to the UV solution (3.12). In particular we saw in (3.25) that after matching the ingoing solution to the UV we are left with a solution of the form from which the Green's function near (1.3) can be extracted by expanding (C1) near the UV boundary as ψ UV = ψ (0) + ψ (1) /r + . . . and with k 2 0 given by (3.3) and the . . . in (C2) refer to terms that vanish as → 0. 6 From expanding (C1) one can explicitly read off where have defined the integral and F (r 0 ) was defined in (3.13).
From (C3) this implies we need to chose q = q p by imposing withÑ which is equivalent to our prediction (3.27). Using the explicit forms of b, T, b + (q) we can use (C5) to solve for q p and deduce there is a line of poles passing through (1.3) with a slope (C7) 6 We will again assume we are not at the special value m 2 = 2r 2 0 so that k 2 0 − m 2 = 0.
Likewise there will be a line of zeroes when (C1) corresponds to the non-normalisable solution with ψ (0) = 0. From (C3) this corresponds to choosing q = q z such that Indeed from (C1) it is straightforward to read off the Green's function itself as Here we wish to discuss in detail the special case of m 2 = 2r 2 0 , in order to both highlight various subtleties with the discussion in terms of the gauge invariant mode ψ and also to use this case to elaborate on the general argument for pole-skipping we provided in Section II.
Pole-skipping in terms of ψ In particular, whilst our arguments in Sections III B and III C are generically valid, the description of pole skipping in terms of the variable ψ is somewhat different for the choice of m 2 = 2r 2 0 . To see why m 2 = 2r 2 0 is special it is useful to note that the function k 2 + r 3 f (r) = k 2 + m 2 that appears in the equation (3.5) for ψ becomes a constant for this value of m. At the special point (1.3) k 2 = −k 2 0 = −m 2 it therefore vanishes identically, and we must be more careful in describing the behaviour of (3.5) at (1.3). For the choice m 2 = 2r 2 0 the equation of motion (3.5) dramatically simplifies to from which the retarded energy density Green's function can be extracted by solving (D1) subject to ingoing boundary conditions and using G R T 00 T 00 (ω, k) = k 2 (k 2 + 2r 2 0 ) Note that upon setting ω = iλ and k = ik 0 in (D1) this no longer reduces to the equation (3.11) which holds for any other m 2 = 2r 2 0 . Solving (D1) at (1.3) one now finds that the two linearly independent solutions for ψ are from which we see that there is only one solution ψ 1 that is regular at the horizon.
Therefore whilst for any other m 2 = 2r 2 0 there are two linearly independent regular solutions for ψ at (1.3), for the specific choice m 2 = 2r 2 0 there is instead a unique regular solution. Note that this observation is perfectly consistent with smoothly taking the limit m 2 → 2r 2 0 in the solution (3.25) we derived by matching, even though the matching procedure is formally invalid at this point as b → 0. In that limit we find bb + (q) = −1/2r 0 and hence the ingoing solution in Note that although there is a unique solution for ψ 1 , the Green's function (D2) still cannot be defined at this point. To see this note that the ingoing solution ψ 1 is normalisable in the UV, that is the denominator in (D2) vanishes for ψ 1 . However the numerator in (D2) vanishes at k 2 = −k 2 0 because of the explicit factor of k 2 + 2r 2 0 in the relation between G R T 00 T 00 (ω, k) and the asymptotics of ψ.
In contrast to our explanation in Section II, at this value of m it appears that the poleskipping phenomenon is unrelated to the existence of an extra ingoing solution. As we will illustrate shortly, this is not the case. The aforementioned subtleties associated with ψ at the point m 2 = 2r 2 0 , ω = iλ, k = ik 0 are in fact indicators that ψ is not a suitable variable for capturing the general gauge-invariant solution to the Einstein equations at this point. A more careful analysis reveals that at this point ψ in fact obeys a first order equation of motion (which enforces c 2 = 0), as does the other decoupled gauge-invariant variable (see (E2) below). To illustrate the origin of pole-skipping at this value of m, we will therefore examine the solutions for the metric perturbations directly. 7 Furthermore, this also is consistent our discussion in Section II that there should be an extra free parameter for the metric solutions at (1.3) corresponding to δg

Discussion in terms of metric solutions
The enhanced symmetry at m 2 = 2r 2 0 allows the linearised gravitational equations of motion to be solved analytically. The following is an exact solution to these equations at m = 2r 2 0 , ω = iλ, k = ik 0 in ingoing EF coordinates 60r 3 60δg (0) yy r 5 + 60δg (0) yy r 4 r 0 + 8c 5r 2 + 5rr 0 + 6r 2 0 + 3r 5 0 22δg We are working here in a different gauge than in Section II. This solution is regular at the horizon r = r 0 and depends on six constants: the five sources of the dual energy-momentum tensor and scalar operator δg Calculating the expectation value of the dual energy-momentum tensor of this solution using [33,48] T µν = lim r→∞ The full dispersion relations ω(k) of the poles of G R T 00 T 00 can be found by solving Einstein's equations numerically. In this appendix we briefly describe how we did this to produce the numerical results shown in Figures 2 and 3 of section III E.