Hydrodynamic diffusion and its breakdown near AdS$_2$ quantum critical points

Hydrodynamics provides a universal description of interacting quantum field theories at sufficiently long times and wavelengths, but breaks down at scales dependent on microscopic details of the theory. In the vicinity of a quantum critical point, it is expected that some aspects of the dynamics are universal and dictated by properties of the critical point. We use gauge-gravity duality to investigate the breakdown of diffusive hydrodynamics in two low temperature states dual to black holes with AdS$_2$ horizons, which exhibit quantum critical dynamics with an emergent scaling symmetry in time. We find that the breakdown is characterized by a collision between the diffusive pole of the retarded Green's function with a pole associated to the AdS$_2$ region of the geometry, such that the local equilibration time is set by infra-red properties of the theory. The absolute values of the frequency and wavevector at the collision ($\omega_{eq}$ and $k_{eq}$) provide a natural characterization of all the low temperature diffusivities $D$ of the states via $D=\omega_{eq}/k_{eq}^2$ where $\omega_{eq}=2\pi\Delta T$ is set by the temperature $T$ and the scaling dimension $\Delta$ of an operator of the infra-red quantum critical theory. We confirm that these relations are also satisfied in an SYK chain model in the limit of strong interactions. Our work paves the way towards a deeper understanding of transport in quantum critical phases.


50
C. SYK chain model 54 Interacting quantum field theories are notoriously challenging, especially when there is no quasiparticle-based description of the state. To describe the late time, long wavelength dynamics of these states, one can instead rely on effective approaches such as hydrodynamics.
At long times and wavelengths, hydrodynamics provides an effective description of a system in terms of a few conserved quantities dictated by symmetries [13][14][15]. Their evolution is governed by local conservation equations for the densities and associated currents, along with constitutive relations expressing the currents in terms of the densities in a gradient expansion. The late time relaxation of the system back to equilibrium is governed by the hydrodynamic modes: poles of the retarded Green's functions of the densities with gapless dispersion relations [16].
While extremely powerful, hydrodynamics breaks down at sufficiently short scales set by the local equilibration time and length. At such scales, the dynamics of the system can no longer be truncated to just the evolution of the conserved densities. Additional degrees of freedom play a significant role, and appear as additional poles of the retarded Green's function with lifetimes comparable to those of the hydrodynamic modes. 1 In cases where the density response exhibits a parametrically slow mode arising due to a weakly broken symmetry, the breakdown of hydrodynamics manifests itself as a collision in the complex frequency plane at nonzero wavevector k eq between the hydrodynamic mode and the slow mode. In this case, it is often possible to augment the hydrodynamic description to incorporate this slow mode [6,[20][21][22][23][24][25][26][27][28]. But typically the modes relevant for the breakdown of hydrodynamics are not of this nature, and a more complete knowledge of a system's microscopic details is required to understand them.
In contrast to this, in the non-zero temperature quantum critical phases found near quantum phase transitions, transport properties are typically universal and are governed by scaling properties of the critical point [29]. We might then expect that the breakdown of hydrodynamics in a quantum critical phase displays a greater degree of universality.
However, it is often challenging to calculate concrete observables beyond thermodynamics due to the lack of analytically tractable models without any simplifying limits, such as a large number of degrees of freedom.
In this work, we exploit gauge-gravity duality to address hydrodynamic transport and its breakdown in quantum critical phases. Gauge-gravity duality maps the late time dynamics of certain large N c quantum field theories (where N c is the rank of the gauge group) to theories of gravity with a negative cosmological constant [15,30,31]. The relaxation of conserved densities back to equilibrium is captured exactly by the evolution of perturbations of asymptotically anti de Sitter (AdS) black holes, which can be studied to obtain a precise understanding of the breakdown of hydrodynamics and the modes responsible for it [32,33].
Even in the absence of a weakly broken symmetry, the breakdown of hydrodynamics can be characterized by an energy scale ω eq and wavenumber k eq , which are sensitive to the system's microscopic details. ω eq and k eq are defined as the absolute values of the complex frequency ω and complex wavenumber k at which the hydrodynamic pole of the retarded Green's function first collides with a non-hydrodynamic pole or branch point [34][35][36][37]. 2 The convergence properties of the real-space hydrodynamic gradient expansion in the linear regime are governed by k eq [39], 3 which also coincides with the radius of convergence of the small-k expansion of the hydrodynamic dispersion relation ω hydro (k). See [41][42][43][44] for recent applications of this.
In this work, we study the breakdown of hydrodynamics in certain low temperature (T ) states dual to black holes with nearly-extremal AdS 2 ×R 2 near-horizon metrics. These are examples of quantum critical phases with an emergent scaling symmetry in time [45][46][47][48], which has been dubbed 'semi-local quantum critical'. Formally, this scaling symmetry corresponds to an infinite Lifshitz scaling exponent, z = +∞: time scales but space does not. Such states are closely related to the Sachdev-Ye-Kitaev (SYK)-like models of electrons in strange metals, which are governed by the same type of infra-red fixed point in the 2 A different mechanism for the breakdown of hydrodynamics arises when interactions between hydrodynamic modes are included, in the guise of a non-analytic frequency dependence of the retarded Green's functions (see e.g. [14] for a review). This mechanism is expected to be suppressed in the large N limit [38]. 3 See [40] for a study of the convergence of nonlinear hydrodynamics in real time. 5 limit of large number of fermions and strong interactions [49][50][51][52][53][54][55][56][57][58][59][60]. 4 Specifically, we study the AdS 4 neutral, translation-breaking black brane of [61,62] and the AdS 4 -Reissner-Nordström (AdS 4 -RN) black brane, and the breakdown of the hydrodynamics governing the diffusive transport of energy, charge and momentum in their dual states. The states we are interested in do not include any slow modes in the sense described above. Instead, local equilibration is controlled by the intrinsic dynamics of the quantum critical degrees of freedom of AdS 2 ×R 2 .
We identify simple, general results for the local equilibration scales ω eq and k eq and confirm that these also apply to the SYK chain model studied in [44] in the limit of strong interactions.
Our first result is that the breakdown is caused by modes associated to the AdS 2 region of the geometry, and as a consequence ω eq is set by universal (i.e. infra-red) data via ω eq → 2π∆T as T → 0, (1) where ∆ is the infra-red scaling dimension of the least irrelevant operator that couples to the diffusion mode. This is in contrast to systems with a weakly broken symmetry, for which ω eq T ,but is in line with the expectation that the quantum critical dynamics is controlled by a 'Planckian' timescale τ eq ∼ 1/T [29,63]. More precisely, we find that at small k and T the Fourier space locations of the longest-lived non-hydrodynamic poles are inherited from infra-red Green's functions, and are approximately located at ω n ≡ −i(n + ∆)2πT for non-negative integers n. The breakdown is characterized by a collision, parametrically close to the imaginary ω axis, between the n = 0 mode (which has a weak k-dependence) and the hydrodynamic mode. This collision manifests itself as a branch-point singularity in the dispersion relation of the mode.
Secondly, we find that at low temperatures the corrections to the quadratic approximation −iDk 2 to the exact hydrodynamic dispersion relation are parametrically small such that the collision occurs when k is almost real and In other words, the scales k eq and ω eq governing the regime of validity of hydrodynamics are set simply by the diffusivity D and the scaling dimension ∆. In some of the examples we study (those involving diffusion of energy), the relevant diffusivity is controlled by an irrelevant deformation of the AdS 2 fixed point and in these cases the result (2) indicates that k eq is controlled by the same irrelevant deformation. A priori, the result (2) is quite surprising: it relates the radius of convergence to just the leading order term in the hydrodynamic expansion. This is a consequence of the AdS 2 fixed point.
By rearranging equation (2) we obtain an answer to the question raised in [7] of what the underlying velocity and time scales are that govern the diffusivity in non-quasiparticle systems. In all our examples they are set by the local equilibration scales where v eq ≡ ω eq /k eq and τ eq ≡ ω −1 eq are the velocity and timescale associated to local equilibration.
In the cases where diffusive hydrodynamics breaks down due to a parametrically slow mode protected by a weakly broken symmetry, D is typically set by τ eq and the speed of the propagating mode that dominates following the breakdown. 5 We emphasize that the breakdown of hydrodynamics is qualitatively different in the cases we study: there is not a single slow mode but a tower of AdS 2 modes with parametrically similar lifetimes set by ω n , and the breakdown does not produce a propagating mode with velocity v v eq .
More generally, the local equilibration time has been argued to set an upper bound on the diffusivity in [64,65]. All examples that we study are consistent with a bound of the form D v 2 eq τ eq for the range of parameters we have investigated. In the absence of a slow mode, it was proposed that low temperature diffusivities are set by the butterfly velocity v B and Lyapunov time τ L that characterize the onset of scrambling following thermalization of the system [66]. This was shown to robustly apply to the diffusivity of energy density D ε in holographic theories and SYK-like models [58,59,[67][68][69][70]. For the examples we study, τ −1 L = 2πT and which is furthermore true in general for states governed by an infra-red AdS 2 with the universal deformation [68].
As for our result (3), equation (4) can be viewed as a consequence of the excellent applicability of the quadratic approximation to the exact hydrodynamic dispersion relation up to the relevant scale. Specifically, pole-skipping analysis suggests that the energy diffusion , from which (4) follows assuming corrections to the quadratic, diffusive form −iD ε k 2 at k = iv −1 B τ −1 L are parametrically small as T → 0.
Our result (3) is more general than (4) in that it is true for all diffusivities in the examples we study, not just the diffusion of energy. Fundamentally this is because, by definition, all diffusive modes pass through the location set by (ω eq , k eq ), while only the energy diffusion mode satisfies the pole-skipping constraint above [74]. As a consequence we provide a new perspective on, and generalization of, the relations between equilibration, transport and scrambling and their applications in AdS 2 /SYK-like models of electrons in strange metals.
In the remainder of this work, we explain how we arrive at equations (1) and (2) before closing with comments on implications and the more general applicability of our results.

II. DIFFUSIVE HYDRODYNAMICS
The spectrum of hydrodynamic modes is dependent on the system under consideration, and by our definition each hydrodynamic mode has its own associated local equilibration scales ω eq and k eq . We will focus on hydrodynamic diffusion modes, which arise when a system has a current density j with constitutive relation where ρ is the corresponding conserved density. ρ will then obey the diffusion equation with diffusivity D ρ at leading order in the derivative expansion, and has a retarded Green's function [14,16] where χ ρρ ≡ lim ω→0 G ρρ (ω, k) is the static susceptibility of ρ, and ellipses denote terms with higher powers of ω and k. The dispersion relation of the hydrodynamic diffusion mode is 8 Hydrodynamic diffusion is a very general phenomenon. Even within the restricted class of systems that we study, the set of conserved densities that exhibit hydrodynamic diffusion varies. In this work, we will be interested in the diffusion of energy ε and of transverse momentum Π.

III. DIFFUSION IN A NEUTRAL HOLOGRAPHIC STATE
We begin with the AdS 4 neutral translation-breaking model [61,62] which is a classical solution of the action with spacetime metric supported by two scalar fields ϕ i = mx i (i = 1, 2) that break translational symmetry. The emblackening factor of the solution is where r 0 denotes the location of the horizon with associated temperature T . The linear perturbations can be written in terms of four decoupled variables and we will focus on the one exhibiting the single hydrodynamic mode of the system. See Appendix A for further details on this spacetime, and on the calculations leading to the results below.

A. Hydrodynamic mode
The hydrodynamic mode corresponds to diffusion of energy, with the small k dispersion relation (7) and diffusivity D ε → 3/2m −1 in the low T limit [20].
First we quantify corrections to this result that will enable us to understand the breakdown of hydrodynamics. In the low temperature limit T ∼ k 2 ∼ 1, the retarded Green's function of energy density G εε exhibits a pole located at where we have explicitly written all dependence. For suitably small k, this is an approximation to the dispersion relation of the hydrodynamic mode as we show in Figure 1. It becomes invalid near specific wavenumbers k 2 = k 2 n related to the breakdown of hydrodynamics, which will be addressed shortly.
It is important to note that (11) is different than the hydrodynamic expansion: corrections to the quadratic k 2 term are not being neglected as in the usual gradient expansion, but are parametrically small in this limit under consideration. One consequence of this is that if we define any wavenumber k * with k 2 * ∼ T at low T , and define ω * to be the location of the hydrodynamic pole at this wavenumber, then (11) implies D ε → iω * /k 2 * as T → 0. For example, choosing k * = iv −1 B τ −1 L results in the chaos relation (4) as described in the Introduction. 6 We will soon show that the breakdown of hydrodynamics at low T is characterized by a pole collision at k 2 eq ∼ T, ω eq ∼ T and thus the diffusivity can alternatively be expressed simply in terms of these scales by (2). But prior to exploring the pole collision that characterizes the breakdown of hydrodynamics, it is instructive to first understand the origin of the non-hydrodynamic mode responsible.

B. Infra-red modes
At low T the state is governed by an infra-red fixed point manifest in the emergence of a near-horizon AdS 2 ×R 2 metric with SL(2,R) symmetry (see Appendix A). Each linear perturbation of the spacetime can be characterized by ∆(k), a wavenumber-dependent scaling dimension of the corresponding operator with respect to this infra-red fixed point, and a corresponding infra-red Green's function [47,75] For the spacetime perturbation that exhibits a diffusive mode, ∆(k) = (1+ 9 + 8k 2 /m 2 )/2.
Although analytically reconstructing G εε from G IR is not easy, for our purposes it is enough to observe that G εε exhibits poles whose locations approach those of the poles of G IR as k, T → 0. Specifically, this means that in this limit G εε exhibits poles at ω →ω n = −i2πT (n + ∆(0)), n = 0, 1, 2, . . . ,  with ∆(0) = 2. The dispersion relation of the n = 0 pole is shown in Figure 1.
At low T , the infra-red modes (13) have a parametrically longer lifetime than the other non-hydrodynamic poles of G εε , and are responsible for the breakdown of hydrodynamics.
The wavenumbers at which our calculation of the low T dispersion relation (11) of the hydrodynamic mode is invalid are k 2 n = 8/3(2 + n)πmT + O(T 2 ), 7 for which the mode would be located at precisely ω(k n ) = ω n in the limit of low T . The natural interpretation would therefore be that the invalidity of the calculation at these values of k 2 can be traced to the nearby presence of the infra-red mode (assuming that the location of the infra-red mode has a weak k-dependence). A more refined calculation below confirms this, as well as the existence of a collision between these modes for complex k that signals the breakdown of hydrodynamics.
7 See (A43) for a precise expression for k n .

C. Breakdown of hydrodynamics
In order to extract the existence of the pole collision, a more refined perturbative computation of G εε at the points ω = ω n + δω, k 2 = k 2 n + δ(k 2 ) is required. This yields where we show only terms relevant for understanding the collision. The low T limit of each coefficient is while D n → D ε in the same limit. The comparable size of the δω and τ n (δω) 2 terms at frequencies δω ∼ T 2 indicates that for such frequencies G εε is dominated by two poles, whose dispersion relations are given by solving the quadratic equation (14) for δω(δk). In  Figure   3.
The collision closest to the origin of k-space (n = 0) signals the breakdown of hydrodynamics. The absolute value of k and ω at this collision are (as T → 0) from which our main results (1) and (2) follow.
The collision location asymptotically approaches real (imaginary) values of k (ω) as T → 0. More precisely, as T → 0 the phases of k and ω at the collision point are where k collision = k eq e iφ k and ω collision = ω eq e iφω .
Note that even after the pole collision formally indicating the breakdown of hydrodynamics, Figure 4 illustrates that the system continues to exhibit a diffusion-like mode described extremely well by the dispersion relation (11). In the limit of zero temperature, the tower of infra-red poles in AdS 2 ×R 2 coalesces in a branch cut along the imaginary axis [47,76].
As the branch cut passes through k = 0, we expect that a hydrodynamic-like series for the dispersion relation at T = 0 would contain non-analytic terms, as was found recently in [77] for a similar state.

IV. DIFFUSION IN A CHARGED HOLOGRAPHIC STATE
The AdS 4 -RN solution to Einstein-Maxwell gravity has a metric of the form (9) but supported by a radial electric field such that the emblackening factor is  The state exhibits two independent diffusive hydrodynamic modes, each associated to a different such variable. The first, corresponding to diffusion of energy and U (1) charge with diffusivity D ε , is analogous to the diffusive mode of the previous section. 8 The second corresponds to the transverse diffusion of momentum with diffusivity D Π . In the limit of low temperature, the diffusivities are 9 The variables exhibiting each of these modes have k → 0 infra-red scaling dimensions 8 In the low T limit both such modes have D ε = κ/c ρ with κ the open circuit thermal conductivity and c ρ the heat capacity [28,68]. 9 See Appendix B for the full, T -dependent expressions. [76, 78, 79] and numerical calculations confirm that at small k and T each corresponding retarded Green's function exhibits non-hydrodynamic poles at the locations (13). As before, a collision close to the imaginary ω axis between the longest-lived such pole and the hydrodynamic pole signifies the independent breakdown of hydrodynamics in each case.
At low T (and until the collision occurs) both hydrodynamic modes are described extremely well by the quadratic approximation to the hydrodynamic dispersion relation, while the locations of the longest-lived non-hydrodynamic poles depend only very weakly on k.
As a consequence, the equilibration scales in both cases are set by the simple formulae (1) and (2) as shown in Figure 5.
The phase of the collision wavenumber φ k in each case is shown in Figure 6, illustrating that the collision point asymptotically approaches real values of k as T → 0. Both cases here, and the result (17) in the previous example, are consistent with ∆ controlling the low T scaling of the phase via φ k ∼ T ∆−1/2 . As in the previous example, the system continues to exhibit diffusion-like modes even after the collision formally indicating the breakdown of hydrodynamics. The dispersion relations of these modes are extremely well approximated by the quadratic approximation to diffusive hydrodynamics, as shown in Figure 10 of Appendix B.

V. COMPARISON WITH SYK CHAIN
The SYK model is a (0+1)-dimensional theory of N interacting fermions that, in the limit of large N and strong interactions, is governed by the same effective action as a theory of gravity in a nearly-AdS 2 spacetime [49][50][51][52][53][54][55][56][57]. The SYK chain [58] is a higher-dimensional generalisation of this, which has served as a very useful toy model for studying diffusive energy transport in strange metal states of matter. As it exhibits the local quantum criticality characteristic of AdS 2 ×R 2 fixed points, it is natural to ask whether local equilibration in this explicit microscopic model is governed by our general results (1) and (2).
In [44], an SYK chain model with N Majorana fermions per site χ i,x and Hamiltonian was studied. The two terms represent q-body on-site and nearest neighbour interactions respectively, where the couplings J i 1 ...iq,x and J i 1 ...i q/2 j 1 ...j q/2 ,x are Gaussian random variables with zero mean. Remarkably, in the limit N q 2 1 an exact analytic expression for G εε was found for all values of the effective interaction strength 0 < v < 1 and the relative strength of on-site and inter-site interactions 0 < γ ≤ 1 [44]. The analytic expression is given in Appendix C, where the details of the model are also summarised.
Taking advantage of this result, a detailed study of the breakdown of diffusive hydrodynamics as a function of interaction strength v was performed in [44]. Here we will focus on the limit of strong interactions v → 1, which is equivalent to T → 0. In this limit, the longest-lived non-hydrodynamic modes are a series of infra-red modes located at precisely the frequencies ω n of equation (13) with ∆ = 2 as k → 0. Hydrodynamics breaks down due to a collision between the hydrodynamic mode and the longest-lived of these infra-red modes, and in Figure 7 we confirm that the local equilibration scales are given simply by analogously to (1)  interactions here happens for real k (i.e. φ k = 0). It would be interesting to determine whether finite 1/q corrections generate a small non-zero phase φ k . Consistently with the other examples we have presented, following the formal breakdown of hydrodynamics the spectrum still contains a mode whose dispersion relation is very well approximated by the quadratic approximation to diffusive hydrodynamics.

VI. OUTLOOK
There is good reason to expect that at least some of our results will generalise beyond the specific examples studied here to other states governed by AdS 2 infra-red fixed points.
While our key observation that the quadratic approximation to the hydrodynamic dispersion relation works parametrically well even for wavenumbers k 2 ∼ T seems unusual, it is nontrivially consistent with the result (4) that is indeed true for holographic AdS 2 ×R 2 fixed points with a universal deformation [68] as well as in related SYK chain models [58,59].
There are more general holographic and SYK-like systems governed by AdS 2 fixed points that exhibit additional diffusive modes beyond the two types we have studied. Of particular interest are non-translationally invariant systems with a U (1) symmetry, for which an Einstein relation relates the electrical resistivity to a diffusivity [7,59,68]. If our results (1) and (2) extend to such modes, they will therefore also provide a simple relation between the phenomenologically important electrical resistivity and the local equilibration scales of such strongly correlated systems.
Confirmation of the broader applicability of our result (3) for AdS 2 ×R 2 solutions would be an important step for quantifying diffusivities near general infra-red fixed points. One way to do this would be to identify a speed u and timescale τ such that in general D ∼ u 2 τ with the coefficient being T -independent. This is difficult even for the relatively simple with the identification u = v eq and τ = τ eq . 10,11 Provided that naive T -scaling holds for the equilibration scales in the other cases with finite Lifshitz exponent z (τ eq ∼ 1/T and v eq ∼ T 1−1/z ), which seems likely, this identification will then work for all fixed points.
For cases where the breakdown of hydrodynamics is due to a slow mode, it is the separation of scales between the decay rate of the slow mode Γ and that of typical nonhydrodynamic excitations T that allows one to augment the hydrodynamic description to incorporate the slow mode. Mathematically, the separation allows one to resum the hydrodynamic expansion into a square root form, valid at scales ω ∼ k ∼ Γ (see e.g. [20,25,27]).
This makes manifest the convergence properties of the hydrodynamic expansion, k eq ∼ Γ.
It would be very interesting if one could extract an analogous effective theory for the cases described here, taking advantage of the separation of scales between the decay rate of the infra-red modes (set by T ) and that of the other non-hydrodynamic excitations (set by the curvature of AdS 2 ). Such an effective theory would need to resum the effects of the entire tower of infra-red modes and would give a greater understanding of why the quadratic approximation to diffusive hydrodynamics is valid up to (and indeed beyond) the wavenumber It would also be interesting to study whether our results continue to hold for AdS 2 fixed points supported by a different hierarchy of scales (such as large angular momentum or magnetic field compared to temperature), and whether analogous results hold for other types of hydrodynamic modes near these fixed points.
It would be very useful to have a semi-holographic description of our results (along the lines of [82,83]), in which we couple a gapless hydrodynamic diffusion mode to the tower of infra-red modes associated with the AdS 2 region of the spacetime. Our results rely on the fact that the the two types of mode have very little effect on one another (see e.g. Figure   1) and a semi-holographic description may clarify exactly under what conditions this is the case. A related avenue would be to adapt the recently constructed holographic effective Schwinger-Keldysh action for diffusion to AdS 2 horizons [83][84][85][86].
Extending our work to Schrödinger z = 2 IR geometries would allow to make contact with ultracold atomic systems [87], which realize a strongly-interacting Fermi gas near unitarity [5].
In [88], the charge diffusivity of a cold atomic system coupled to an optical lattice was measured. In [89,90], the thermal diffusivity of high T c superconductors in the strange metallic regime was also reported. It would be interesting to investigate the deviations from diffusive hydrodynamics in these systems.
The examples we studied are all consistent with bounds on D of the type proposed in [64, 65] but with the velocity in the bound given by v eq (rather than the operator growth velocity, characteristic velocity of low energy excitations, or butterfly velocity). 12 Indeed, the arguments in [64, 65] assume that v eq is set by one of these velocities and thus imply the bound D v 2 eq τ eq . In this sense our results support the assumption of [64, 65] that the local equilibration is controlled by an underlying effective lightcone, even though the systems are non-relativistic. It would be very worthwhile to determine D, v eq and τ eq for other states, and over a wider parameter range, in order to establish the robustness of these observations and to determine whether v eq is set by a speed such as the butterfly velocity in general. [63] Jan Zaanen, "Superconductivity: Why the temperature is high," Nature 430, 512-513 (2004). In this appendix we present in more detail the calculations leading to the results presented in the main text for the neutral translation-breaking model of [61,62]. The action of this model is and it has the classical solution with metric and scalar field profiles ϕ i = mx i (i = 1, 2) that break translational symmetry. The horizon is located at r 0 and the asymptotically AdS boundary at r → ∞. The Hawking temperature of the solution is and further details of its thermodynamic properties can be found in [62].
At zero temperature, m 2 = 6r 2 0 and an AdS 2 ×R 2 metric emerges near the extremal horizon. To see this explicitly, one should change coordinates from (t, r) to (u, ζ) where and then take the → 0 limit of (A2) to obtain where the AdS 2 radius of curvature is L 2 = 1/3.
The hydrodynamics of this model were studied and explained in [20]. Over the longest distance and time scales, this system supports a single hydrodynamic mode corresponding to the diffusion of energy. The dispersion relation of this mode is Green's function is [20] ψ(r, ω, k) = where ω and k here are the frequency and wavenumber after a Fourier transform with respect to the t and x 1 coordinates in which the metric has the form (A2). This variable obeys the equation where primes denote derivatives with respect to r and For some calculations it will be convenient for us to convert to an ingoing Eddington-Finkelstein like coordinate system with time coordinate v = t + r * , where the tortoise coordinate r * is In this coordinate system the relevant gauge-invariant variable may be written as where ω here denotes the frequency after a Fourier transform with respect to the v coordinate. This variable obeys the equation of motion d dr where Ω(ω, k) = (2r 2 0 − m 2 )3ir 0 ω + k 2 (k 2 + m 2 ). (A13) Up to contact terms, the energy density Green's function can be expressed as [73] where in this expression ψ denotes the solution of (A12) that is regular at the horizon r = r 0 in ingoing coordinates. We will fix the overall normalisation of the solution as ψ(r 0 , ω, k) = 1, without loss of generality.

Near-horizon perturbation equations
We begin by analysing the relevant perturbation equation in the AdS 2 ×R 2 spacetime that emerges near the horizon at low temperatures. To do this, we follow [47] in introducing the location of the extremal horizon r e via 6r 2 e = m 2 and then defining This is similar to the Fourier space version of the limit (A4), but now keeping ω/T fixed as ω → 0. The temperature is proportional to ζ 0 . At leading order in the near-horizon, low temperature limit (i.e. → 0 limit), the equation (A8) becomes This is the equation of a scalar field in non-zero temperature AdS 2 with effective mass L 2 m 2 eff = 2 (1 + k 2 /m 2 ). As such, we can associate it with an operator of dimension ∆(k) = (1 + 9 + 8k 2 /m 2 )/2 at the infra-red fixed point of the field theory. The infra-red Green's functions G IR of such operators can be found explicitly by solving the equation (A16) and imposing the usual AdS/CFT rules at the AdS 2 boundary ζ → ∞ [47, 75] At any non-zero temperature the infra-red Green's function has a series of poles along the negative imaginary frequency axis at the locations To obtain G εε , one must extend the near-horizon solutions through the rest of the spacetime.
In practice this is very difficult to do, but schematically one expects a result of the form [47,75]

Low temperature dispersion relation
In this section we will derive the dispersion relation (11) for a pole of G εε in the limit of low T (with k 2 /T fixed). The key observation that will allow us to make progress is that the equation (A12) simplifies considerably when Ω = 0 i.e. when ω = ik 2 (k 2 + m 2 )/(3r 0 (2r 2 0 − m 2 )), allowing us to formally obtain G εε along this line in Fourier space. We will use this as a starting point for a perturbation theory, allowing us to determine G εε close to this line in Fourier space. From this we will find that in the low T limit (with k 2 /T fixed), the Green's function near this line typically exhibits a single pole with dispersion relation (11).
Exceptional cases arise near particular points on this line, for which a more sophisticated perturbation theory will be presented in the following section.
First we choose a point (ω * , k * ) in Fourier space that lies on the line Ω = 0, around which we will perturbatively evaluate G εε . Note that although we will often write ω * and k * independently in the equations that follow, they are not really independent but are constrained by Ω(ω * , k * ) = 0. It will be simplest to think of us choosing a wavenumber k * which then fixes ω * (k * ).
We now take frequencies and wavenumbers that are close to this point in Fourier space and look perturbatively for a solution to the equation (A12) of the form ψ(r, ω, k) = ψ(r, ω * , k * ) + δψ(r, ω * , k * ), where we will evaluate δψ(r, ω * , k * ) to first order in δω and δ(k 2 ). It will sometimes be convenient for us to repackage the two parameters δω and δ(k 2 ) as δω and δΩ where follows from (A13).

a. Leading order solution
The leading order equation (A23) can be integrated trivially for any k * , giving ψ(r, ω * , k * ) = C 0 + C 1 dr e 2iω * r * (k 2 where C 0 and C 1 are constants. Imposing the ingoing boundary conditions means that we require that ψ(r, ω * , k * ) has a Taylor series expansion near the horizon.
Near the horizon, the exponential term in the integrand e 2iω * r * ∼ (r − r 0 ) iω * 2πT (which for generic ω * corresponds to an outgoing mode). Thus for generic ω * , the ingoing solution is Substituting this into the expression (A14) we obtain the leading order Green's function This leading order result is exact for points lying exactly on the line Ω = 0 and the fact that it is finite mean that there are no poles lying exactly on this line when ω * = 0.
Although the ingoing solution for generic choices of k * is given by (A26), there are important exceptional cases. Assuming that k 2 * + r 3 0 f (r 0 ) = 0, the integrand of (A25) is of the form ∼ (r − r 0 ) iω * 2πT −1 near the horizon. Therefore for the exceptional choices of wavenumber the integral term in (A25) will also obey the appropriate ingoing boundary condition. For these cases, the leading order ingoing solution is not simply that in (A26), and we will address them separately in the next section. There are two additional exceptional cases that we will not discuss as they have been studied extensively elsewhere: , and the usual hydrodynamic expansion k 2 * = ω * = 0 [20].

(A30)
Each integral in this expression gives rise to an integration constant that we must fix to maintain the appropriate ingoing boundary conditions.
Consider first the inner integrand. Recalling that near the horizon e −2iω * r * ∼ (r −r 0 ) − iω * 2πT , the inner integral is only finite at the horizon for Im(ω * ) > −2πT . We will assume this condition for now, and will show in the next section that it can be relaxed without affecting the results. With this condition, we write the inner integral as where C 2 is a constant. Substituting this into (A30), we deduce that C 2 = 0 in order for δψ to obey ingoing boundary conditions at the horizon. Therefore the ingoing solution is given by δψ(r, ω * , k * ) = δΩ where we have also imposed the second boundary condition ψ(r 0 , ω, k) = 1.
We can now evaluate the energy density Green's function using equation (A14). From our point of view the important piece is the term (−r 2 ψ + iωψ) r→∞ in the denominator, which contains the information about poles of G εε . At ω = ω * + δω and k = k * + δk, this term has the form While the perturbative corrections in (A33) are formally small, the key point is that in an appropriate low temperature limit one of these corrections becomes anomalously large and thus is important. This allows us to deduce that at low T the Green's function near the Ω = 0 line in Fourier space is dominated by a single pole whose dispersion relation we will shortly compute.
The relevant low T limit is T → 0, keeping k 2 * /T fixed (and therefore ω * /T fixed), for which the integral on the first line of (A33) diverges like 1/T 2 . The divergence arises from the near-horizon region of the integral, and to isolate it we define the AdS 2 radial coordinate ζ as in equation (A4), rescale T → T , k 2 * → k 2 * , ω * → ω, and then expand the integrand in the → 0 limit to obtain Evaluating this integral up to a UV cutoff Λ of the near-horizon region (with Λ T ), we find the following divergent, cutoff-independent contribution to the integral on the first line of (A33) 35 The ellipsis denotes terms that are subleading in the → 0 limit. The double integral in equation (A33) also diverges like 1/ 2 in the same limit.
In the expression (A33) for the denominator of G εε , we observe the first and third terms will dominate in the low T limit described above. More specifically, if we scale δΩ → 3 δΩ and δω → 3 δω along with the scalings of T, ω * and k 2 * above, then in the limit → 0 the denominator is The right hand side vanishes for a suitable δΩ, which means that G εε is dominated by a pole at the corresponding location in the complex ω plane. More specifically, for a given fixed wavenumber k = k * (i.e. δ(k 2 ) = 0) we can replace δΩ by δω using equation (A22) and then solve for δω(ω * ). Recalling that ω * and k * are related by Ω(ω * , k * ) = 0, we obtain the pole location ω(k) = −i 3 2 quoted in equation (11) in the main text.
The result (A37) contains information about the low T limit of the hydrodynamic dispersion relation. For example, by comparing (A37) to (A6) we would obtain the correct low temperature expansion of the diffusivity. But the result (A37) is not simply the hydrodynamic expansion of the dispersion relation: the terms with higher powers of k 2 are not being ignored but are formally subleading in the expansion we have described. Furthermore, although we have assumed that Im(ω) > −2πT until now, we will shortly see that the result (A37) continues to hold beyond this. In fact, even for values of k outside the radius of convergence of the hydrodynamic expansion of the dispersion relation, we will see that there is still a pole at the location indicated by equation (A37).

c. Extension to general Im(ω)
In our derivation of the result (A37) we assumed that Im(ω) > −2πT in order that we could rewrite an integral appearing in the perturbative solution for δψ as (A31). When this inequality is not obeyed, we must take more care as the integral formally diverges near the horizon. To deal with this, we will manually separate out the divergent part of the integral before imposing ingoing boundary conditions in an analogous manner to the steps above.

36
The upshot of this is that G εε continues to be dominated by a pole with dispersion relation (A37) independently of the value of Im(ω) (within the low T limit under consideration).
The only exceptions to this are near the points (A28) that we previously highlighted, which we address in the following section.
To extend the results to lower in the complex ω plane, we write the inner integral appearing in the solution (A30) for δψ as where We can formally expand G(r) in a Taylor series near the horizon with coefficients G n and then write the integral as where the integral on the right hand side is finite provided that Im (ω * ) > −2πT (2 + N ).
With this expression we now evaluate the contribution of the δψ term in the perturbative expansion of the Greens function denominator (A14), anticipating that it will again be anomalously large in the relevant low T limit, and find The terms on the first line will vanish provided that Im (ω * ) < −2πT (1 + N ). Assuming this, along with the previous condition, we find a contribution to the Green's function denominator given by There are an extra set of terms in comparison to the corresponding term in the Green's function denominator (A33) for Im(ω * ) > −2πT .
By setting N = 0, 1, 2, . . . and evaluating this integral, we can extend our previous results to successively lower regions of the complex ω plane. To isolate the divergent contribution to the integral arising from the near-horizon region we follow the same procedure as before: take the low T , near-horizon scaling limit of the integrand and integrate the leading term up to a UV cutoff Λ T of the near-horizon region. After doing this explicitly for N = 0, 1, 2, 3 we find the same result (A35) as before, and we expect that this will continue to be the case for higher N . As a consequence, the results for the Green's function derived previously (including the dispersion relation (A37) for the pole) are also valid in the region Im(ω) < −2πT .

Analytic description of pole collisions
We now turn to the exceptional points (A28) for which the analysis in the previous section must be modified. Recall that for these points the leading order ingoing solution is not given by (A26), despite the fact that they lie on the line Ω = 0. In fact, at these points the ingoing solution is no longer uniquely defined and thus the perturbative expansion is more subtle.
These 'pole-skipping points' correspond to points where a line of poles of the field theory Green's function intersects with a line of zeroes [71-74]. Here, we will take a direct approach and just describe in practice how to obtain the appropriate perturbative solutions.
The main result of this calculation is that in the low T limit, there are in fact two poles of G εε in the vicinity of each exceptional point with n ≥ 0. One passes directly through the point while the other is separated from it by a distance ∼ T 2 in ω space. These two poles collide for specific complex values of k, that become real as T → 0. Assuming that the collision occurring at the smallest value of |k| characterizes the breakdown of diffusive hydrodynamics (which we will confirm numerically in the next section), this allows us to quantitatively extract ω eq and k eq (or equivalently τ eq and v eq ) in the low T limit.

a. Leading order ingoing solution
Our starting point is the result (A25) for the general solution for ψ, at a point (ω * , k * ) lying exactly on the line Ω(ω * , k * ) = 0. At the frequencies ω n given in (A28), this solution obeys ingoing boundary conditions at the horizon for any values of the constants. The corresponding values of k n are found by solving Ω(ω n , k n ) = 0. This equation has two distinct solutions for k 2 n and ultimately we will be interested in the one with k 2 n ∼ T at low T . Using equation (A3) to trade T and r 0 , explicitly this is where on the second line we have written the small T expansion.
To fix the boundary conditions at these points we must consider points in Fourier space infinitesimally away from (ω n , k n ), where the ingoing solution is unique. Specifically, we take the equation of motion (A12) at frequencies ω = ω n + δω and wavenumbers k 2 = k 2 n + δ(k 2 ) and construct a Taylor series solution for ψ near the horizon. This solution is ingoing by definition, and is unique (up to an overall normalisation). We subsequently expand the Taylor series solution at small δω ∼ δ(k 2 ), finding the leading order result ψ(r, ω n + δω, k n + δk) → 1 + c n δΩ δω (r − r 0 ) 2+n + . . . , where the coefficients c n are independent of both δω and δ(k 2 ). The important aspect of this solution is that, although it is unique, it depends on the ratio δΩ/δω that parameterizes exactly how we have moved away from the point (ω n , k n ) in Fourier space. By expanding our general leading order solution (A25) near the horizon and matching it to the ingoing solution (A44), we obtain the values of the constants C 0 and C 1 that correspond to an ingoing solution at the point (ω n + δω, k n + δk n ). Explicitly, the ingoing solution at leading order is then ψ(r, ω n , k n ) = 1 + i2πT α (1) where we have defined the constants and the integrals Correction to the leading order solution Having established the leading order solution (A45) near these special points in Fourier space, we will now repeat the same procedure at one higher order in the expansion at small δω ∼ δ(k 2 ) to obtain the first correction to this result. After substituting the leading order result (A45) into the equation (A24), we can trivially integrate the resulting equation to obtain the expression δψ(r, ω n , k n ) = dr e 2iωnr * (k 2 for the correction. There are two arbitrary integration constants arising from the two indefinite integrals in the first term, and these need to be fixed in accordance with the ingoing boundary conditions.
Before doing this, note that for all n ≥ 0 the inner integral of the first term diverges near the horizon and so it will be convenient to separate out the divergent terms by hand to leave a finite integral. 13 Specifically, we re-write it as where β n is the integration constant. Substituting this into (A48), we obtain the solution 13 For the special case n = −1, the summation terms in (A49) and subsequent equations can be set to zero.
in the form δψ(r, ω n , k n ) = δΩ r r 0 dr e 2iωnr * (k 2 n + r 3 f ) 2 In this expression we fixed one integration constant such that ψ(r 0 , ω, k) = 1. The remaining integration constant β n has been split up into parts proportional to δΩ, δΩ 2 /δω and δΩδ(k 2 )/δω for later convenience.
The integration constants β (i) n are fixed by imposing ingoing boundary conditions. To determine these conditions explicitly we again use the unique Taylor series solution for ψ(r, ω n + δω, k n + δk) near the horizon, and find the correction to equation (A44) at the next order in the small δω ∼ δ(k 2 ) expansion. At this order, the series has the form ψ(r, ω n + δω, k n + δk) = 0 + . . . + (r − r 0 ) 2+n γ (0) n δΩ + γ (1) The coefficients γ (i) n become increasingly complicated for higher values of n and so we will only explicitly present the n = 0 expressions, which are those relevant for the breakdown of To impose ingoing boundary conditions on the solution (A50), we expand it near the horizon and match it to the ingoing solution (A51) to fix the values of the three integration constants n . Explicitly, these integration constants are then and (A55) In summary, the first correction to the ingoing solution in the small δω ∼ δ(k 2 ) expansion is given by (A50) and (A53).

c. Low temperature limit of the Green's function
The ingoing solution near (ω n , k n ), including the first correction in the small δω ∼ δ(k 2 ) expansion, is given by the sum of (A44) and (A50). It is convenient to rewrite this as ψ(r, ω n + δω, k n + δk) = 1 + Φ (0) n (r) where the ellipsis represents higher order terms in the small δω ∼ δΩ expansion, and Substituting into (A14), the denominator of the energy density Green's function near (ω n , k n ) is (−r 2 ψ + iωψ) r→∞ ∝ lim r→∞ δΩ −r 2 Φ (0) n (r) + iω n Φ (0) n (r) + 2πT (2 + n)δω + iδω 2 42 The first important result can be seen from this formal expression: it vanishes as δω, δΩ → 0 and therefore there is always a pole of the Green's function passing through the point (ω n , k n ). The dispersion relation of the pole in the immediate vicinity of this point can be found by taking the O(δω) and O(δΩ) terms in (A58) and solving for δω as a function of δ(k 2 ).
As the O(δ 2 ) terms are formally small, one might naively conclude that the Green's function is dominated by this single pole in the vicinity of (ω n , k n ). However the second important result is that, in the low T limit, the coefficient of the δΩδω term becomes anomalously large for all n ≥ 0. This indicates that there is an additional pole of the Green's function near (ω n , k n ) in these cases. By carefully evaluating the expression (A58) for the Green's function at low T we will further be able to show that for each n ≥ 0 these two poles undergo a collision in complex Fourier space near (ω n , k n ). The collision closest to the origin of k space (n = 0) involves the hydrodynamic diffusion mode and thus characterizes the breakdown of hydrodynamics.
To take the low temperature limit, we formally rescale T ∼ and then evaluate the scaling of each term in (A58) in the limit → 0. 14 For the n = −1 case, the coefficient of However, for the cases n ≥ 0, the scaling of the coefficients in the low T limit is different: the coefficient of the δΩ term is ∼ T 0 while the coefficients of the δωδΩ and δΩ 2 terms are ∼ T −2 and ∼ T −1 respectively. Therefore the quadratic corrections are important at frequencies δω ∼ T 2 , indicating that the single-pole approximation to the Green's function breaks down here due to the presence of a second pole. In the limit T → 0, these two poles are parametrically closer to one another in Fourier space than the typical separation ∼ T between poles.
To make this more quantitative, we can formally rescale δΩ ∼ δω ∼ 2 and then compute (A58) in the low → 0 limit for n ≥ 0. If we keep only the leading order terms in this limit (O( 2 )), the Green's function denominator near (ω n , k n ) has the simple factorised form where we replaced δΩ using (A22), D n = 3/2m −1 , and τ n is given in equation (15) in the main text. 15 The expression (A59) clearly shows that the Green's function near the points (ω n , k n ) with n ≥ 0 is dominated by two poles. Within this approximation the poles will coincide at the exactly imaginary frequency δω = −iτ −1 n and real wavenumber δ(k 2 ) = τ −1 n D −1 n . But this is not yet really a pole collision 16 : the poles pass through each other undisturbed, and the dispersion relations δω(δ(k 2 )) extracted from (A59) have an infinite radius of convergence.
To observe the finite radius of convergence characteristic of a pole collision, we must go beyond the result (A59) and include subleading corrections in the expansion. One source of these are finite T corrections to τ n and D n . Although these corrections are the largest in magnitude (e.g. the leading correction to τ n is suppressed by T log T ), we will ignore them as they will not alter the functional form of (A59) and so will not produce a pole collision.
The other source of corrections are the δω, δΩ 2 and δω 2 terms present in the full expression (A58). The first two corrections enter at O( 3 ) while the latter is O( 4 ). Of the first two, we will only need to consider that proportional to δω as it is clear from inspecting the first line of (A59) that adding a term ∝ δΩ 2 will still result in two nicely factorised poles whose dispersion relations δω(δ(k 2 )) each have an infinite radius of convergence. Upon adding the leading δω term correction to (A59), we obtain the expression for the energy density Green's function denominator presented in the main text.
The expression (A60) no longer factorizes in a simple manner and solving explicitly for the dispersion relations δω(δ(k 2 )) of the two poles results in expressions with a nonzero discriminant term. The poles will collide at the wavenumber where the discriminant vanishes, which is when Expanding out these expressions for δω c and δ(k 2 c ) at low T , and recalling that ω = ω n + δω and k 2 = k 2 n + δ(k 2 ), gives the collision location for every n ≥ 0. We have shown the leading order real and imaginary terms in each expression, and the ellipses denote small T corrections to these. Note however, that there are corrections to the leading real (imaginary) term in k 2 c (ω c ) that are larger than the leading order imaginary (real) term (for example, coming from T log T suppressed corrections to τ n ).
There are two collision points in the complex ω plane with opposite signs of Re(ω) and these collision points asymptotically approach the imaginary (real) axis as T → 0. From (A62), the phase of the collision wavenumber and frequency are where we define these quantities using the collision in the upper right quadrant of the complex k c plane.
In summary, we have derived an analytic description of a series of pole collisions happening in the vicinity of the points (ω n , k n ) for all n ≥ 0. The n = 0 collision characterizes the breakdown of diffusive hydrodynamics (as can be seen from the numerical results in the following section) and from (A62) the local equilibration scales are therefore given by as quoted in equation (16) in the main text.

Comparison with numerical results
In this section we present exact results for the locations of poles of G εε obtained numerically. In addition to verifying our analytic results in the appropriate regions of Fourier space, these numerical results provide a global picture of the pole structure and therefore show how the different analytic results we have presented connect together.
To obtain the exact locations of poles of G εε one needs to solve for the relevant linear fluctuations around the black hole geometry (A2). As explained above one can determine G εε by solving for a suitably defined variableψ (see (A7)) that satisfies a decoupled equation of motion. For our numerical calculations we followed [73] and solved instead for which satisfies the equation of motion d dr We solved this equation numerically using a dimensionless radial variable z = r 0 /r. We integrated from the horizon at z = 1, where we imposed ingoing boundary conditions Ψ 1 = , to the boundary at z = 0 where we read off the leading and subleading contributions of Ψ 1 = ψ 0 +ψ 1 z +O(z 2 ). We used 8 terms in the near-horizon expansion, an IR cutoff of 1 − z = 10 −5 and working precision 50 in Mathematica. 17 Finally, the poles of G εε were identified as zeroes of the ratio ψ 0 /ψ 1 [73].
For small wavenumbers k, it is the hydrodynamic pole of G εε with dispersion relation (7) that lies closest to the origin of the complex ω plane. In order to characterize the breakdown of hydrodynamics we must first identify the non-hydrodynamic poles lying near to the origin with which this pole could collide. At small enough T /m (and for small real k), a family of poles with imaginary frequencies are the non-hydrodynamic poles closest to the origin. In Figure 8 we show the locations of the longest-lived such pole in ω space at small k, demonstrating agreement with (13) as T /m → 0. G εε of course exhibits other non-hydrodynamic poles besides these, but at low T /m these are far from the origin and so will not play any role in the rest of our discussion.
We will now examine more carefully how these infra-red poles move in the complex ω plane as real k is varied, and how they fit together with the diffusive hydrodynamic mode.
In the low T limit, we have established that there is generically a pole with the dispersion relation (A37) but that near each of the points (ω n≥0 , k n≥0 ) there is an additional pole.
Given that the infra-red modes lie at ω n≥0 as k → 0, the simplest explanation would be that the infra-red modes have a very weak k-dependence and are themselves the additional  poles. This explanation is correct, as is shown in Figure 4. Note that while this Figure   shows multiple instances of poles coming very close to one another, the collisions themselves cannot be seen as they occur at the complex values of k (A62).
The final part of our numerical analysis is to examine more closely the dynamics of the poles near (ω 0 , k 0 ) where we expect a collision between the hydrodynamic diffusion mode and the longest-lived infra-red mode to occur, characterising the breakdown of hydrodynamics.
Based on our calculations in the previous section, we expect the collision to occur for the complex value of the wavenumber (A62). In Figure 3 we show pole trajectories in the complex ω plane as |k| is varied at fixed T and phase of k, illustrating that two poles collide for an appropriate |k|.
In Figure 9 we quantitatively compare features of this pole collision to the analytic results obtained in the previous section, showing excellent agreement as T → 0. The upper two panels show exact results for the local equilibration scales ω eq and k eq , while the lower panels show the ratio Dk 2 eq /ω eq and the phase of the collision wavenumber φ k .
In summary, the low T numerical results show that as real k is increased, the hydrodynamic diffusion mode moves down the imaginary ω axis before coming into the vicinity of the longest-lived infra-red mode near ω = −i4πT . The breakdown of hydrodynamics is characterized by a collision between these two modes near this point when k has a small imaginary part. The low T collision is captured quantitatively by the analytic expressions derived in the previous section.
We finish this section with some observations on what happens after hydrodynamics breaks down. After the collision signaling the breakdown, one may expect the character of the collective modes of the system to change. But this does not happen. After the collision occurs, there continues to be a diffusive-type mode with dispersion relation (A37) as is apparent from Figure 4 in the main text. In fact, this mode then undergoes a second collision near ω = −i6πT before quickly re-appearing and so on through a whole series of collisions near the frequencies ω n≥0 . So while hydrodynamics formally breaks down, diffusion-like transport continues to occur. momentum with diffusivity D Π . We focus on the latter two modes, whose diffusivities can be computed from Einstein relations [14] as .

(B6)
Each hydrodynamic mode corresponds to a quasinormal mode of one of the four decoupled bulk variables. The sound modes arise as quasinormal modes of the longitudinal variable ψ L − , diffusion of energy and charge arises as a quasinormal mode of the longitudinal variable ψ L + , and transverse diffusion of momentum arises as a quasinormal mode of the transverse variable ψ T − . 18 From these observations, we obtain ∆ ε (0) = 2 and ∆ Π (0) = 1 for the k → 0 scaling dimensions of the variables that exhibit the respective hydrodynamic diffusion modes.

Numerical results
The linearised perturbation equations for the AdS 4 -RN black brane solution are significantly more complicated than for the neutral translation-breaking model of Appendix A.
We will therefore establish the results by solving these equations numerically to extract the relevant information about the poles of G εε and G ΠΠ . The poles of the Green's functions have previously been studied in [76,77,79,[92][93][94][95] and more recently the breakdown of hydrodynamics was looked at in [34,42]. The breakdown of hydrodynamics in the AdS 5 case was studied in [41].
To obtain our numerical results for the poles of G εε we used the coupled variables and equations described in [94]. 19 Although in these variables it is not manifest that the poles can be separated into two decoupled sectors, on occasion it will be useful for us to take advantage of this fact to simplify the results that we present. As for the perturbations of the neutral translation-breaking model we work with a dimensionless radial variable z = r 0 /r and integrated from the horizon at z = 1 (where we imposed ingoing boundary conditions) to the boundary at z = 0. To determine the poles of G εε we followed the determinant method of [96]. We used 12 terms in the near-horizon expansion, an IR cutoff of 1 − z = 10 −5 and working precision 60 in Mathematica. For the lowest temperatures and largest wavenumbers we worked with 1 − z = 10 −6 and working precision 80. 18 Our definition of the variables ψ L ± is the opposite of that in [79], and instead agrees with that in [75]. 19 Denoting the chemical potential in [94] byμ, µ = 2μ due to a different normalisation of the action.
Our numerical results for the G ΠΠ were obtained by solving the perturbation equation for the decoupled variable ψ T − , as in [76]. The poles arising from ψ T + will not be relevant to the results we present. In this case, we used 9 terms in the near-horizon expansion, an IR cutoff of 1 − z = 10 −5 , and working precision 50 in Mathematica. As before, for the lowest temperatures and largest wavenumbers we decreased the cutoff to 1 − z = 10 −6 and increased the working precision to 70.
For this state it is well-established that at low temperatures the longest-lived nonhydrodynamic poles of each Green's functions are a family of poles with imaginary frequencies ∼ T [76, 79]. In Figure 10 we show the motion of the diffusive pole of each Green's function, along with that of the longest-lived non-hydrodynamic poles, as function of real k for a fixed low temperature. For both Green's functions the diffusive mode moves down the imaginary ω axis and is extremely well approximated by the quadratic approximation to diffusive hydrodynamics, while the infra-red modes are approximately k-independent and lie very close to the locations (13) of the poles of the infra-red Green's functions. Note that G εε can be separated exactly into the sum of two independent pieces, each associated to one of the variables ψ L ± . In Figure 10, we show only poles of G εε that are associated to the variable ψ L + . The reason is that generically a pole of ψ L − cannot have any effect on the hydrodynamic diffusion pole (or its radius of convergence) even if these poles coincide, as they are solutions to independent differential equations. 20 In the vicinity of the frequencies ω n≥0 , the diffusive mode and an infra-red mode become very close, but a detailed inspection ( Figure 11) shows that there are no collisions near these points for any real k. To observe the pole collisions that characterize the breakdown of hydrodynamic diffusion, we must consider the motion of the poles for complex values of k.
For G ΠΠ , the first pole collision happens close to ω = −i2πT and plots of ω eq , Dk 2 eq /ω eq and φ k for this collision are shown in Figures 5 and 6 in the main text. In [42] it was shown that k eq (T → 0) → 0 (both by identifying the location of the pole collision we have described, and independently by explicitly computing the radius of convergence of the hydrodynamic series), which is consistent with the more quantitative results we have presented. At sufficiently high T , the nature of the collision characterising the breakdown 20 An exception to this arises at particular values of k where the equations for the two decoupled variables become equivalent. Such collisions can be relevant to the breakdown of hydrodynamics [34,41,42], but not in the low T regime which we study. of hydrodynamics undergoes a qualitative change [34,42].
For G εε the breakdown of hydrodynamics is characterized by a collision near ω = −i4πT between the diffusive mode and the first infra-red mode of ψ L + . In Figures 5 and 6 in the main text we show corresponding plots of ω eq , Dk 2 eq /ω eq and φ k . In Figure 12 we show the motion of the two poles near the collision point, thereby confirming that this is a branch point of the hydrodynamic dispersion relation. Our results for k eq are consistent with those in [42], which observed that k eq → 0 as T → 0. As for G ΠΠ , the collision that characterizes   the breakdown of diffusive hydrodynamics in G εε changes qualitatively at sufficiently high T [41,42].
Finally, we mention that it is manifest from Figure 10 that even far outside the formal regime of applicability of hydrodynamics (i.e. far past the first collision of the hydrodynamic pole), the system still supports modes whose dispersion relations are extremely well approximated by the quadratic approximation to diffusive hydrodynamics, with diffusivities