Damping of hard excitations in strongly coupled $\mathcal N\,{=}\,4$ plasma

The damping of high momentum excitations in strongly coupled maximally supersymmetric Yang-Mills plasma is studied. Previous calculations of the asymptotic behavior of the quasinormal mode spectrum are extended and clarified. We confirm that subleading corrections to the lightlike dispersion relation $\omega({\bf q}) = |{\bf q}|$ have a universal $|{\bf q}|^{-1/3}$ form. Sufficiently narrow, weak planar shocks may be viewed as coherent superpositions of short wavelength quasinormal modes. The attenuation and evolution in profile of narrow planar shocks are examined as an application of our results.


I. INTRODUCTION
In a strongly coupled Yang-Mills plasma, such as that of maximally supersymmetric Yang-Mills (N = 4 SYM) theory, the typical time scale for relaxation of non-hydrodynamic perturbations is set by the inverse temperature T −1 . In a dual holographic description, this scale may be interpreted as the characteristic gravitational infall time for perturbations falling through the horizon of black brane geometries which describe near-equilibrium states [1,2].
However, even in the strong coupling limit, sufficiently high momentum excitations are only weakly damped. This may, for example, be seen in the large wavenumber asymptotics of the quasinormal mode (QNM) spectrum. At zero temperature in N = 4 SYM, Fourier transformed two-point correlation functions, viewed as functions of frequency ω at fixed wavenumber q, have branch cuts starting at the lightcone, ω = ±|q|. 1 At non-zero temperature, and in the N → ∞ limit, this branch cut breaks up into a closely spaced series of poles at locations ω = {ω ± n (q)} known as quasinormal mode frequencies [3][4][5]. Festuccia and Liu [6] studied the large-q asymptotics of the quasinormal mode spectrum for scalar perturbations (or helicity ±2 stress-energy perturbations) and found that as |q| → ∞, ω ± n (q)/|q| ∼ ± 1 + c n e ∓iπ/3 (πT /|q|) 4/3 , with real "spectral deviation" coefficients {c n } which are discussed below. The small imaginary part (relative to the real part), Im ω ± n (q)/ Re ω ± n (q) = O(T /|q|) 4/3 , demonstrates the weak damping for |q| T , and shows that these high energy, short wavelength excitations may, in some respects, be regarded as quasiparticles, i.e., excitations whose mean free path is much longer than their de Broglie wavelength. However, because |q| T , these are highly athermal excitations which are exponentially rare in equilibrium. Moreover, because the spacing in energy between successive quasinormal modes is comparable to their width, |ω ± n+1 (q) − ω ± n (q)| ∼ | Im ω ± n (q)|, the spectral densities of correlation functions, at large q, do not have distinct narrow peaks in frequency associated with each quasinormal mode; instead the contributions of multiple QNMs merge to produce a slowly varying spectral density [7].
The weak damping of high q excitations may also be seen in the behavior of planar shocks. 2 At zero temperature, planar shocks propagate at the speed of light without dispersion or attenuation. At non-zero but low temperatures (small compared to the inverse width of the shock), the shock experiences weak thermal drag [8]. This slowly attenuates the amplitude of the shock and introduces dispersion, but this weak damping vanishes as T → 0.
In this paper, we study the damping of high q excitations in N = 4 SYM theory in greater detail. In section II we perform our own WKB analysis of the large-q asymptotics of helicity ±2 quasinormal mode frequencies. We confirm the relative |q| −4/3 form (1) of the leading corrections to a lightlike dispersion relation, with a universal ∓π/3 phase. However, we find values for the coefficients {c n } of these corrections which disagree in two respects with the result given in ref. [6], which was c n = K FL (2n + 1) 4/3 , n = 1, 2, · · · , with K FL = [ √ π Γ 7 4 /Γ 1 4 ] 4/3 = 0.344127 · · · . The (2n + 1) 4/3 dependence on mode number is asymptotically correct for high-lying modes, but is not accurate for low order modes. Moreover, the coefficient K FL differs by a factor of 2 5/3 from the correct value in the large order asymptotic form, with The need for a 2 5/3 correction factor in the value of the coefficient K FL for AdS 5 black holes was noted earlier in ref. [9], 3 but the inaccuracy of the estimate (3) for low order modes seems not to have been previously appreciated.
Complementary numerical results confirming the WKB analysis, examining the approach to the asymptotic regime, and studying helicity 0 and ±1 modes in addition to helicity ±2, are presented in section III. We calculate accurate results for the lowest fifteen quasinormal modes in each helicity channel for wavevectors up to |q|/(πT ) = 160. This extends previous results given in ref. [5]. For helicity ±2 perturbations, comparison of the numerical results with the WKB asymptotics clearly confirms the validity of the asymptotic analysis and shows that for low order modes the large-q asymptotic form (1) becomes a good approximation starting at modest wavenumbers of a few times πT . For helicity ±1 and 0 perturbations (which satisfy significantly more complicated equations) we have not performed a full WKB asymptotic analysis. However, our numerical results for these helicities very clearly support the assertion that the asymptotic form (1) is equally valid for these perturbations. Moreover, our extrapolated numerical values for the first fifteen spectral deviation coefficients c n strongly suggest that in these helicity channels the large order asymptotic form is c n ∼ K (2n) 4/3 , helicity ±1; K (2n − 1) 4/3 , helicity 0, (4) with exactly the same prefactor K = 1 2 [ √ π Γ 7 4 /Γ 5 4 ] 4/3 as for helicity ±2. As an application of our results, we discuss the propagation of narrow planar shocks in section IV. A sufficiently weak shock may be viewed as a coherent superposition of quasinormal modes. As noted above, as the shock moves through the dispersive N = 4 SYM plasma at temperature T it experiences friction; the maximum amplitude will decrease and the longitudinal energy density profile will evolve. We specifically study narrow planar shocks whose quasinormal mode spectrum is dominated by wavevectors large compared to T and discuss characteristic features of the resulting evolution. The final section V contains a few concluding remarks. Appendix A presents tabular data for QNM frequencies. Three subsequent appendices provide details on the numerical analysis, transformation to infalling coordinates, and the large wavevector, large order asymptotic analysis.

II. QUASINORMAL MODE FREQUENCIES: LARGE-q ASYMPTOTICS
We wish to study the dynamics of linearized perturbations on the background geometry of an AdS 5 black brane, which is dual to the thermal equilibrium state (at vanishing chemical potentials) of N = 4 SYM theory. We find it convenient to use infalling coordinates (t, x, u) with x ≡ (x 1 , x 2 , x 3 ) denoting ordinary spatial coordinates and u an (inverted) bulk radial coordinate. Choosing to set the AdS curvature scale L equal to unity, the metric reads The conformal boundary is at u = 0 and the horizon lies at u h ≡ (πT ) −1 , with The metric is translationally invariant in the Minkowski directions (t, x). Hence, it is natural to Fourier decompose the dependence of perturbations on these directions and, for non-zero wavevectors, to classify according to the helicity of the perturbation under the SO(2) little group [10]. In this section we focus, for simplicity, on helicity ±2 perturbations. Choosing the wavevector q to lie along the x 3 -direction (with magnitude q), we consider a metric perturbation of the form with h an undetermined function of u. Factoring out two powers of u, as shown, is convenient as the appropriate boundary condition for h at u = 0 then becomes just regularity. Similarly, because our infalling coordinates are non-singular across the future horizon, ingoing boundary conditions at the horizon correspond to h also remaining regular at u = u h . With this choice of perturbation, the only non-trivial part of the linearized Einstein's equations is the xy component, and the resulting equation reads Henceforth, it is convenient to choose units such that m = 1 (or equivalently, to measure ω and q in units of πT ), so that the helicity ±2 perturbation equation becomes 4 It will also prove convenient to denote the frequency to wavevector ratio by This ratio will be complex and wavenumber dependent [i.e., s = s(q)], although this dependence will not always be indicated explicitly. From the quasinormal mode equation (9) it is apparent that if h(u) is a solution with frequency ω then h(u) * is also a solution with frequency −ω * , showing that quasinormal mode frequencies (which are not pure imaginary) come in pairs with opposite real parts. Hence, it is sufficient to focus on solutions with Re s ≥ 0. One may eliminate first derivative terms in the helicity ±2 equation (9) by suitably redefining the radial function. Let with and z ≡ u 2 . Thenĥ(z) satisfies a zero-energy Schrödinger equation, with Boundary conditions onĥ (corresponding to regularity of h at horizon and boundary) arê The six singular points of eq. (9) (at u = 0, u = ∞, and u 4 = 1) are reduced to four in eq. (13): z = 0, z = ∞, and z = ±1. The resulting equation (13) is thus of the Heun type.

A. Leading behavior
As mentioned earlier, in the large q (or low temperature) limit, where the spatial wavevector is arbitrarily large compared to πT , quasinormal mode frequencies should approach the zerotemperature branch points at ω 2 = q 2 . To demonstrate that this is indeed the case, we insert an ansatz for the asymptotic behavior of the ratio s = ω/q, with exponent α > 0 and the "dispersive correction" s α (q) a smooth function of q which approaches a finite non-zero limit, with corrections vanishing as an inverse power of q. First, to show that s 0 must equal ±1, we make a proof by contradiction: assume that s 2 0 = 1 and demonstrate that there are no solutions. Eq. (13) becomes where An appropriate ansatz for a WKB approximation to the solution iŝ Subsequent terms in the exponent involve higher integer powers of q −1 and q −α . The ordering of the terms will be explained a-posteriori, when we find that α is non-integer and 1 < α < 2.
Inserting the expansion (20) into the radial equation (18) and collecting like powers of q produces the conditions: Solving for T 0 , T 1 , and T α yields two solutions (due to the sign ambiguity in √ Q 0 ). One choice gives where we define Q 0 (z) as the branch which approaches +is The resulting WKB approximations for two linearly independent solutions, which we denote byĥ ± WKB , have the formĥ The most general solution is an arbitrary linear combination ofĥ ± WKB . Subleading terms in these WKB approximations are negligible provided that |1 − s 2 0 − z 2 | q −α and |1 − s 2 0 − z 2 | z q −2 . The first condition ensures that the Q 0 (z) term in eq. (18) is large compared to Q α (z), while the second condition ensures that Q 0 (z) also dominates the Q 2 (z) term.
Near the horizon, 1−z 1, we have Only the behavior ofĥ + WKB matches the required near-horizon condition (15b), so this is the solution of interest.
Near the boundary, z 1, we have (1−s 2 0 )/z, with 1−s 2 0 defined to be positive just above the branch cut running from −1 to 1.
The form (25) cannot, however, be directly compared with the required boundary condition (15a), as the WKB approximation (23) is not valid all the way to z = 0; as noted above, the WKB approximation is limited to z q −2 /|1−s 2 0 |. Therefore, we must match the WKB solution to a suitable near-boundary approximation. 5 Provided s 2 = 1, the Schrödinger equation (13) forĥ simplifies near the boundary, z 1, tô with solutions given by regular or irregular modified Bessel functions, These forms are valid for z 1, regardless of the size of q 2 z, up to relative corrections of order z 2 . In the overlap region 1 z q −2 /|1−s 2 0 |, both WKB and near-boundary approximations are valid. Within this region, the arguments of the Bessel functions in the near-boundary approximations (27) are large and these solutions behave as 6 , then the WKB approximation is valid for all z ∈ [ , 1] for any q −2 /|1−s 2 0 |. But if s 2 0 is real and lies inside the interval [0, 1] then there is a quadratic turning point at z * = 1−s 2 0 . The WKB approximation (23) is not accurate in a neighborhood of this turning point. Nevertheless, this does not invalidate the following argument matching WKB and near-boundary approximations, as one may deform the contour in z along which one works from the real interval [0, 1] to a complex contour which runs from 0 to 1 but avoids the turning point at z * . This contour deformation argument does not apply when s 2 0 → 1, as the endpoints of the contour in z are necessarily fixed at 0 and 1. 6 These asymptotic forms, and the following argument, are valid provided (1−s 2 )z has positive real part. As the phase of √ 1−s 2 varies away from zero, it is convenient to perform the matching on the ray arg z = − arg √ 1−s 2 , along which the arguments of the modified Bessel functions remain real.
Comparing these forms to the WKB behavior (25), one sees thatĥ + WKB is proportional to the near-boundary solutionĥ irr , not toĥ reg . However, only the regular near-boundary solutionĥ reg satisfies the boundary condition (15a) requiring O(z 3/2 ) behavior as z → 0. The irregular solution h irr diverges as O(z −1/2 ) as z → 0, violating the required regularity condition. Consequently, the assumption that s 2 0 = 1 is inconsistent with the boundary conditions (15); solutions which satisfy the boundary condition at one end of our interval in z fail to satisfy the required boundary condition at the other end. Therefore, the only solutions which satisfy both boundary conditions must have s 2 0 = 1, implying that quasinormal mode frequencies approach ±q as q → ∞.

B. Subleading behavior
Specializing (without loss of generality) to the case of s 0 = +1, the integrals appearing in the WKB functions (22) may be explicitly evaluated and give: Hence, the relevant WKB solution has the form As discussed above, neglected higher order terms are negligible provided z 2 q α 1 and z q 2 1. Once again, this solution will need to be matched, within a suitable overlap region, to an appropriate near-boundary solution. For z 1, T 0 (z) ∼ i 3 z 3/2 and (with no assumption on the size of z compared to inverse powers of q), the WKB solution (30) behaves aŝ We now turn to the near-boundary region. Non-uniformity between the small z and s 2 0 → 1 limits cause the near-boundary behavior for s 2 0 = 1 to be qualitatively different from the previously discussed s 2 0 = 1 case. So we must redo the analysis starting from eq. (13) and specializing to s 0 = 1. Assuming z 1 and q 1 (but making no assumptions about products of the form z a q), the Schrödinger equation (13) simplifies tô It is helpful to introduce a rescaled coordinate, so thath(y) ≡ĥ(z(y)) satisfies In terms of this rescaled coordinate, the small-z form (31) of the WKB solution becomeŝ and is valid for . All digits shown are accurate.
then we have a consistent description for asymptotically large q: the WKB solution has a universal small-z form,ĥ WKB (z(y)) ∼ y −1/4 exp i 3 y 3/2 − is ∞ α y −1/2 , valid for y 1, which can smoothly match onto a solutionh(y) of the q-independent near-boundary equation, To determine allowed values for the constant s ∞ α , one must find solutions to eq. (37) which are O(y 3/2 ) as y → 0 and, up to an irrelevant overall constant, approach y −1/4 exp i 3 y 3/2 − is ∞ α y −1/2 when y 1. Although it may seem most natural to work on the ray with arg y = 0 (corresponding to the original physical domain of z ∈ [0, 1]) when performing this matching, this is not required. For reasons which will momentarily become apparent, it is more convenient to work along the rotated ray arg y = π/3. So we define with w real and positive. Then h(w) ≡h(y(w)) satisfies In other words, by rotating the contour, our desired solution now vanishes exponentially for large argument. Moreover, with these boundary conditions eq. (39) is a self-adjoint eigenvalue problem. Specifically, λ is an eigenvalue of the self-adjoint positive operator From the form of the effective potential appearing in this operator, it is clear that it has a pure point spectrum. So the eigenvalues {λ n } must form a discrete set of real, positive values. Consequently, the subleading asymptotic coefficient s ∞ α must have the form with a real, positive sequence of values {c 1 , c 2 , · · · } equal to twice the eigenvalues {λ n }. The Schrödinger equation (39) has an irregular singular point at w = ∞ along with a regular singular point at w = 0. An analytic solution does not appear to be possible, but solving this equation numerically is relatively straightforward. We describe our numerical techniques in appendix B and present the resulting values for the first 25 spectral deviation coefficients {c n } in table I. The values of c n rapidly increase with increasing mode number n. For n 1, one may use a further WKB approximation to find the large n asymptotics of these coefficients. When the eigenvalue λ is large, a simple WKB approximation for high order eigenfunctions is valid in regions where the potential 1 4 w−λw −1 + 3 4 w −2 is sufficiently slowly varying. One must appropriately match to a near-boundary approximation (given by a Bessel function) for small w, and also match across the linear turning point at w ≈ 2 √ λ. Details of this exercise are presented in appendix D. One finds that solutions satisfying the required boundary conditions exist when with Figure 1 shows a comparison of this asymptotic form with the numerical results in table I. For the lowest n = 1 mode, the deviation from the asymptotic scaling (41) is approximately 6% (far larger than the precision of the results in table I). But by n = 5 the asymptotic form is accurate to about half a percent. The rapid approach to the asymptotic form (41) could have been anticipated from the fact that already for modest values of n the coefficients c n become quite large compared to unity. Examination of the rate of convergence confirms the expected 1/n 2 scaling of the deviation.
To summarize, we have shown that helicity ±2 quasinormal mode frequencies, for large wavenumbers, have the form (16) with α = 4/3 and s ∞ α having phase −π/3. Continuing the WKB analysis, it is straightforward to show that the next term in the large-q expansion is O(q −1 ). Therefore, for large wavenumbers, helicity ±2 quasinormal mode frequencies are given by plus reflected frequencies −ω n (q) * , with the real coefficients {c n } shown in table I. These coefficients have the large order asymptotic form (41). Restoring factors of πT gives the result (1) quoted in the introduction.

III. QUASINORMAL MODE FREQUENCIES: NUMERICS
To validate the large-q asymptotics (43) and examine the accuracy of this approximation for intermediate ranges of wavenumber, we use pseudo-spectral methods [11] to solve numerically the quasinormal mode equations for a wide range of wavenumbers. 8 This extends previous work in ref. [5]. We consider first the helicity ±2 case, and then examine helicity ±1 and 0. A. Helicity ±2 As previously noted, frequencies for which the helicity ±2 quasinormal mode equation (9) has solutions satisfying the required regularity conditions at horizon and boundary come in pairs with opposite real parts (and identical imaginary parts): {ω n } and {−ω * n }. So it is sufficient to consider only the positive frequency spectrum, i.e., Re ω ≥ 0.
To apply spectral methods, it is convenient to return to the original form (9) of the helicity ±2 QNM equation. Representing h as a (truncated) series of Chebyshev polynomials, automatically satisfies the required regularity conditions at u = 0 and 1. Demanding that equation (9) [multiplied by u(1−u 2 )] be satisfied at each point u = u k on the collocation grid, for k = 0, · · ·, M , yields a finite set of linear equations of the form where, given an explicit choice of the wavevector q, A and B are numerical (M +1) × (M +1) matrices. The generalized eigenvalue problem (46) may be efficiently solved in O(M 3 ) time using standard methods. Results for the first fifteen helicity ±2 quasinormal mode frequencies ω n (q) [or rather, the deviation (ω n (q)−q)], for wavevectors q = 10, 20, 40, 80 and 160, are shown in table IV of appendix A. The real and imaginary parts of the first five helicity ±2 quasinormal modes are plotted in figure 2 for modest wavenumbers up to 4πT . We have verified that our quasinormal frequencies for q=2 agree with those given in app. B of ref. [5]. 9 To present results for larger wavenumbers in a manner which allows easy comparison with the asymptotic form (43), we use the definition (16) of the dispersive correction function (with α = 4/3), repeated here, The magnitude and phase of the dispersive correction s α,n (q) for the first 5 modes are shown in figure 3 for wavenumbers up to q/(πT ) = 50. One sees, as expected, that s α,n (q) approaches the 9 Note that Kovtun and Starinets [5] give results in units of 2πT , not πT . asymptotic value c n e −iπ/3 extracted from the WKB analysis. Lower modes converge faster than higher modes. The rapid rise of the magnitude |s α,n (q)| as q increases from zero is an artifact of definition (47) (since ω n (q) has a finite q → 0 limit). But the leveling off of the magnitude after this rise provides a clear visual indicator of the onset of the asymptotic regime. From the figure it might appear that the convergence of the magnitude of s α,n (q) toward its asymptotic value c n is considerably faster than the convergence of the phase to −π/3. This, however, is an illusion produced by the rather compressed range of the ordinate in the right hand plot (which was chosen to make the different phase curves visually distinct). One may parameterize the raw data in table IV of appendix A using the functional form and demanding that the result reproduce the values in table IV. This form is a truncation of the series which is generated by higher order asymptotic analysis. 10 The resulting values for the first coefficient A (1) n , when multiplied by e iπ/3 , provide independent estimates of the asymptotic coefficients {c n }. These estimates, based on what is effectively an extrapolation to q = ∞, are less accurate than the values listed in table I, but the agreement is quite good. The deviation is less than a part in 10 4 for the first few modes, but grows to about half a percent for n = 15. (This reflects the slower approach to the large-q asymptotic form of progressively higher modes.) Moreover, we have explicitly tested that using the parameterization (48) of the data in table IV, the resulting functions reproduce the directly calculated values of the quasinormal mode frequencies used to produce figure 3 (showing the range 10 ≤ q ≤ 50) to within a precision of two parts in 10 4 .

B. Helicity ±1 and 0
To analyze perturbations with helicity ±1 and 0, it is convenient to use the gauge invariant linear combinations of metric perturbations introduced by Kovtun and Starinets [5]. With a 10 The powers of q in the A (1) n and A (2) n terms reflect the result (43) of sec. II B. When recast as an expansion of ω(q)/q, higher order terms involve products of positive integer powers of q −4/3 and q −2 arising from the decomposition (18) of the effective potential, and form a series in integer powers of q −2/3 . Fefferman-Graham form for the metric of the black brane geometry, helicity ±1 and 0 linear combinations are, respectively, Decoupled second order linear equations satisfied by these fluctuations were derived in ref. [5]. Converting to our preferred infalling coordinates leads to the following equations for these perturbations, with f (u) ≡ 1 − u 4 . Details of the transformation yielding these equations are given in appendix C. The required boundary conditions for the functions Z 1 (u) and Z 2 (u) are just regularity at both horizon (u=1) and boundary (u=0). Frequencies for which solutions satisfying these boundary conditions exist are either pure imaginary, or else come in pairs with opposite real parts, ω and −ω * . Therefore, without loss of generality, in the following discussion we consider Re ω ≥ 0. After multiplying the helicity ±1 equation (51) by its frequency-dependent denominator ω 2 − q 2 f , and likewise multiplying the helicity 0 equation (52) by 3ω 2 −(f +2) q 2 , both equations become cubic generalized eigenvalue problems of the form where each O i is a linear operator. By replicating the function space on which one works, this may be converted into a conventional generalized eigenvalue problem, A X = ω B X, where X ≡ (ω 2 Z, ω Z, Z) and 11 Applying pseudo-spectral methods to convert the linear radial differential operators O i into matrices and solving the resulting finite dimensional generalized eigenvalue problem proceeds in the same manner described previously. Results for the first fifteen helicity ±1 and 0 quasinormal mode frequencies ω n (q) [or rather, the deviation (ω n (q)−q)], for wavevectors q = 10, 20, 40, 80 and 160, are given in tables V and VI of appendix A. We first discuss helicity ±1 perturbations. The real and imaginary parts of the first five quasinormal frequencies are plotted in fig. 4 for q ≤ 4πT . 12 For modest wavenumbers, q 2.6 πT ,  the most weakly damped mode is a hydrodynamic shear mode whose frequency is pure imaginary and vanishes as q → 0. This frequency, which is shown by dashed lines in fig. 4, moves down the imaginary axis as q increases. As seen in the figure and noted in ref. [12], the frequency of this mode crosses the imaginary parts of other mode frequencies (having non-zero real parts) at various intermediate values of q. For q T , this mode becomes highly damped and is not among the minimally damped modes discussed below.
To examine larger values of q and the approach to the asymptotic regime, we plot in fig. 5 the magnitude and phase of the dispersive correction s α,n (q), defined via eq. (47), of the lowest five helicity ±1 modes (excluding the hydrodynamic shear mode) for q/(πT ) up to 20. Unlike the helicity ±2 case, one sees non-monotonic behavior in the lowest modes as q increases. Although we have not done an independent WKB calculation for helicity ±1 to determine asymptotic values directly, from the plots it certainly appears that the magnitudes |s α,n (q)| are approaching constant values while all phases are converging to a value near −π/3. Near-asymptotic behavior begins to be apparent for quite modest values of wavenumber, q/(πT ) ≈ 5.
One may parameterize the helicity ±1 data in  Although not a formal proof, the consistency and accuracy of the parameterization (48), when applied to our helicity ±1 data, strongly suggests that helicity ±1 quasinormal mode frequencies have the same large-q asymptotic form (43) as do helicity ±2 modes. The first coefficients {A (1) n } of the parameterization, when multiplied by e iπ/3 , directly give estimates of the asymptotic spectral deviation coefficients {c n } for helicity ±1 modes. Table II lists these estimates for the first fifteen modes. Within the accuracy of the parameterization (as determined by our tests with 10 < q < 20), the phases of the asymptotic coefficients c n are all compatible with zero.
We now turn to helicity 0 modes, whose behavior largely parallels that of the helicity ±1 modes just discussed. Fig. 6 plots the real and imaginary parts of the first five helicity 0 quasinormal modes for q/(πT ) ≤ 4. There is one hydrodynamic helicity 0 (sound) mode, whose frequency vanishes as q → 0 (with Re ω = O(q) and Im ω = O(q 2 )). As noted in ref. [13], the helicity 0 hydrodynamic sound mode smoothly evolves from small to large values of q and always remains the most weakly damped mode. Its damping, as measured by Im ω/ Re ω, rises linearly from q = 0, reaches a maximum at q/(πT ) ≈ 2.120, and then decreases as O(q −4/3 ) as q continues to grow. Fig. 7 plots the modulus and phase of the spectral deviation function s α,n (q) for these modes out to q/(πT ) = 20. From the figure one sees, once again, that the magnitudes |s α,n (q)| are nearly constant for q/(πT ) 5 and all phases approach a value close to −π/3.
As before, one may parameterize the helicity 0 data in   helicity 0 quasinormal mode frequencies for wavevectors throughout the range 10 ≤ q ≤ 20 to within a precision of a part in 10 4 . This consistency strongly suggests that helicity 0 quasinormal mode frequencies also have the same asymptotic form (43). Table III shows our resulting estimates, extracted from this simple parameterization, for the spectral deviation coefficients {c n } for the first fifteen helicity 0 modes. Within the accuracy of the parameterization, the phases of the asymptotic coefficients c n are, once again, all compatible with zero.
In summary, we have compelling evidence that, for all helicities of metric perturbations, quasinormal mode frequencies have the large q asymptotic form ω n (q) = q + c n e −iπ/3 q −1/3 + O(q −1 ), with O(q −4/3 ) relative corrections to a massless ω = ±q dispersion relation, and with real positive coefficients {c n } characterizing the dispersive correction. This large-q asymptotic form is reasonably accurate, at least for low order modes, down to q/(πT ) ≈ 5.

C. Large order asymptotics
Comparing the helicity ±2 spectral deviation coefficients listed in table I with our corresponding helicity ±1 or 0 values shown in tables II and III, it is clear by inspection that the helicity ±1 and 0 spectral deviation coefficients grow about as fast with increasing mode number as do the helicity ±2 coefficients. Given the known asymptotic behavior (41) of the helicity ±2 coefficients, it is natural to try fitting our estimates of helicity ±1 and 0 spectral deviation coefficients using a function of the form b (2n + a) 4/3 . It is also instructive, for comparison, to apply the same procedure to estimates of the helicity ±2 spectral deviation coefficients generated by the parameterization (48) applied to the data of table IV. In performing these fits, we set to zero the small residual phases in the c n estimates (which are all compatible to zero, within the accuracy of the parameterizations).
For all helicities, the resulting best fit value of the overall coefficient b coincides with the analytically known value (42) of the helicity ±2 coefficient K = 1 2 √ π Γ( 7 4 )/Γ( 5 4 ) 4/3 = 1.092 · · · to within a percent, and is quite insensitive to whether one fits all 15 modes or, for example, just modes 6 to 10. We take this as compelling evidence that our fitting function correctly describes the large order asymptotic behavior of spectral deviation coefficients for helicities ±1 and 0, as well as ±2, with the same overall coefficient K for all helicities. If one fixes the overall coefficient b = K, then the resulting best fit value for the shift a is very close to an integer, and is reasonably insensitive to the limits of the fitting range. For helicity ±2 the best fit value for the shift a equals the correct answer +1 to within four percent. For helicity ±1, the best fit value for the shift a is quite close to zero, somewhere in the interval −0.002 to −0.03 depending on the chosen limits of the fitting range. And for helicity 0, the best fit value for the shift a equals −1 to within a percent or two. Given that we are only fitting data up to n = 15, these results are nicely consistent with the known large order asymptotic form (41) for the helicity ±2 spectral deviation coefficients, and very clearly suggest corresponding large order asymptotic forms for helicity ±1 and 0 coefficients, as reported in the introduction. Explicitly, for helicity s, Figure 8 shows, for each helicity, a comparison of this asymptotic form with our numerical estimates for spectral deviation coefficients produced by using the functional form (48) to parameterize the data of tables IV, V and VI of appendix A. The uppermost curve shows helicity ±2, the middle curve is helicity ±1, and the lower curve shows helicity 0. Fast approach to the large order asymptotic form (55) is manifest. The helicity ±2 curve shown here agrees with the plot in fig. 8, which used the the highly accurate c n values of table I, up to a permille. Curiously, the approach to the asymptotic form is even faster for helicity ±1 and 0 compared to helicity ±2. For helicity ±1, the deviation is only 2% for the lowest n = 1 mode, and falls to half a percent or less for all higher modes. For helicity 0, the deviation is −7% for the lowest mode, but also falls to half a percent or less for all higher modes.

IV. PLANAR SHOCKS PROPAGATING IN N = 4 SYM PLASMA
The general solution for the time evolution of linearized perturbations to the metric of the AdS black brane geometry (with fixed boundary geometry and incoming conditions at the horizon), may be represented as a linear combination of quasinormal modes, where q ≡ |q|. The sum runs over all quasinormal modes (of metric perturbations) with the symmetries of interest, A n (q) is the amplitude of a given mode, and H n encodes the tensor structure of the mode, e.g., H n = dx 1 ⊗dx 2 +dx 2 ⊗dx 1 for helicity ±2 modes with the indicated polarization. Extracting a factor of u 2 , as shown, allows one to fix the normalization of the radial profiles h n (u; q) by requiring that they have a fixed boundary value, h n (0; q) = 1. The corresponding change in the expectation value of the energy-momentum tensor of the dual QFT is then [14] δT µν (t, where κ = m 4 L 3 /(4πG) with L the AdS curvature scale which, elsewhere, has been set to unity. In terms of field theory quantities, Similarly, if one considers a perturbation to the equilibrium state created by a time dependent source coupled to the QFT stress-energy tensor (i.e., a fluctuation in the spacetime geometry of the 4D field theory), then the induced response is given by a convolution with the retarded stress-energy correlator, Quasinormal mode frequencies correspond to the poles of the retarded Green's function [15], so evaluating the frequency integral using Cauchy's theorem yields where R n (q) denotes the residue of the retarded Green's function G R (ω, q) at ω n (q). Both expressions (57) and (59) represent the response as a sum of contributions from quasinormal modes, and make it obvious that at sufficiently late times the response will be dominated by those modes whose frequencies ω n (q) have the smallest (in magnitude) negative imaginary parts. Specifically, these are long wavelength hydrodynamic modes with q T , together with the short wavelength modes with q T discussed above. To examine qualitative features of the resulting evolution, it will be sufficient to use the asymptotic form (43) of quasinormal mode frequencies, repeated here for convenience, Consider a metric perturbation δg, represented in the form (56), which at time t = 0 has rapid spatial variation and a spatial Fourier transform concentrated near some wavevector q 0 with |q 0 | T . The asymptotic form (60) implies that the characteristic relaxation time of such an excitation will be of order τ (q 0 ) ≡ q 1/3 0 (πT ) −4/3 , with higher modes (larger n) damping faster than lower modes. At times comparable or larger than this relaxation time, dominant contributions will come from the n=1 mode with wavenumbers near q 0 . 13 Provided the perturbation is sufficiently small, so a linearized treatment is valid, there is no mode-mixing populating other ranges of wavevector. The resulting late-time energy-momentum tensor (57) is then well described by just the n = 1 contribution, with a damped amplitude and rapidly varying phase If one asks when a disturbance, initially localized near x = 0 at time 0, will reach a distant point x, the dominant contributions to the integral (61) come from the neighborhood of the stationary phase point where ∇φ(q) = 0. (provided A(q) is slowly varying on the scale of |x| −1 ). This yields the standard result that disturbances localized in wavenumber near q = q 0 travel at the group velocity, and hence arrive at position x = dq 0 at time t = d/v g (q 0 ). The damping (62) decreases monotonically with increasing wavenumber, while the group velocity (64) increases monotonically, asymptotically approaching the speed of light. Hence, shorter wavelength features attenuate slower, and propagate faster, than longer wavelength features. For disturbances with a significant range of wavenumbers, the overall evolution will reflect a combination of the wavenumber dependent damping (62) and the dispersive propagation (64).

B. Planar shocks at late times
The evolution of planar shocks in a strongly coupled N = 4 SYM plasma provides an interesting, concrete illustration of the above features. At zero temperature, planar shocks (composed of helicity 0 perturbations) move at the speed of light with no dispersion or damping. For a shock propagating in, say, the +x 3 direction with an arbitrary longitudinal energy density profile κ h(x 3 ), stress-energy components at time t are with all other components vanishing. Equivalently, with A(q 3 ) the Fourier transform of h(x 3 ). The dual geometry is an exact analytic solution of Einstein's equations [8]. Analogous planar stress-energy perturbations with helicity ±1 or ±2 tensor structures correspond to solutions of Einstein's equations linearized about AdS 5 , but analytic solutions at the non-linear level are not known. We are interested in the modification in the evolution of planar shocks induced by the presence of a background thermal plasma. 14 We assume that the amplitude of the shock is sufficiently small so that a linearized treatment is adequate. And, for simplicity, we assume that the perturbation has a single tensor structure corresponding to either helicity 0, ±1, or ±2.
Planarity of the shock implies that the expression (57) for the stress-energy (at times t ≥ 0) simplifies to a one-dimensional Fourier transform, with the coefficients {A n (q 3 )} characterizing the chosen shock profile at time t = 0. As discussed in the previous subsection, since large-q modes experience minimal damping, rapidly varying features in the longitudinal profile of the shock will outlive more slowly varying features. A particularly clear picture emerges for narrow shocks. A shock with narrow width w 1/T will have a broad Fourier spectrum extending from small wavenumbers at least up to |q 3 | ∼ 1/w. More precisely, the breadth of the Fourier spectrum reflects the (inverse) scale of variation of the sharpest spatial features. As a coherent superposition of many different wavevectors, small differences in the propagation of different wavenumbers will imprint themselves on the evolution of the shock profile. In particular, since the speed of propagation approaches the speed of light as |q 3 | → ∞, contributions from the highest-wavenumber modes will accumulate very near the light cone, forming an increasingly sharp leading edge. These sharp features will persist longer than contributions from lower wavenumbers (except for very small q hydrodynamic modes) which attenuate more quickly.
To illustrate this explicitly with simple examples, we consider perturbations which are dominated by the lowest quasinormal mode (of a given helicity), so that One may either regard this as fine-tuning the initial data, or the result of starting with a more general perturbation and waiting to sufficiently late times where higher modes are small compared to the lowest mode. For simplicity, we include only modes with Re(q 3 ω 1 (q 3 )) > 0 in the perturbation (68), with no corresponding contributions from the reflected modes with frequency −ω 1 (q 3 ) * ; this means that we are focusing on right-moving perturbations. We compare three different longitudinal profiles, 15 A Gauss Step" Step" [top to bottom]. Snapshots are taken at times t = 0, 1 2 , 1, 3 2 , 2 for helicity ±2, and at t = 0, 1, 2, 3, 4, 5 for helicity 0; the different time scales reflect the faster damping of helicity ±2 perturbations. For the helicity-2 perturbations one sees that the longest surviving features are associated with the steepest portions of the initial profile. This is especially apparent with the compactly-supported "blob" and "step" profiles. For helicity 0 there is, in addition, a visible slowly varying contribution from the hydrodynamic sound mode.
These are Fourier transforms of position space profiles h(x 3 ) = dq/(2π) e iqx 3 A 1 (q) which are, respectively, a Gaussian, a semicircular "blob," and a "top-hat" step function, each normalized to unit area: We choose σ = 1/10 for h Gauss , and σ = 1/5 for h Blob and h Step . For the dispersion relation ω 1 (q) we construct a cubic spline interpolating function from the numerical results of sec. III for low and intermediate wavenumbers, which smoothly connects to the large-q asymptotics of eq. (60). Since the large-q asymptotic form is already quite accurate for intermediate values of q, this procedure is straightforward. 16 Figure 9 shows plots of the resulting time evolution for perturbations with each of the above profiles, for the helicity ±2 tensor structure dx 1 ⊗ dx 2 as well as the helicity 0 structure (dx 3 −dt) ⊗ (dx 3 −dt). Comparing the plots, one clearly sees the stronger damping of helicity ±2 fluctuations relative to helicity 0, due to the larger values of the dispersive coefficient c 1 (cf. tables I and III). For the helicity ±2 perturbations, shown on the left side of the figure, the longest surviving features are associated with the steepest portions of the initial profile. This is especially apparent for the "blob" and "step" profiles which have compact support. At late times, one sees upward "spikes" which evolve from the portion of the initial profile with large negative gradient, and downward spikes evolving from the steeply rising part of the initial profile.
The helicity 0 evolution, shown on the right-hand side of the figure, shows similar sharp high-q features but with the addition of a slowly varying hydrodynamic (sound) contribution which moves slower than the leading features and gradually broadens. (The speed of sound in the conformal N = 4 plasma is 1/ √ 3 [5].) Hence, as time increases there is an increasingly clear separation between the attenuating high-q and low-q remnants.
The helicity 0 planar shocks have a conserved energy (and linear momentum). At late times, the energy and momentum of the shock is entirely carried by the q → 0 hydrodynamic contribution; for the profiles of fig. 9, the upward and downward high-q spikes cancel each other upon integration. More generally, the long-lived fine structure carries little net energy or momentum. This might seem odd, since high momentum quasiparticles (or short wavelength waves) in other contexts can transport energy and momentum. But stress-energy quasinormal modes in strongly coupled (and large N c ) SYM plasma are perturbations in which the energy density, momentum density, and/or stress of the fluid vary sinusoidally (for non-zero q) and hence unavoidably vanish upon integration. This is fundamentally different from, say, an electromagnetic wave in vacuum in which the EM field varies sinusoidally but the energy density is quadratic in the wave amplitude and always positive.

V. DISCUSSION
We have extended and clarified previous work on quasinormal mode frequencies for metric perturbations of AdS-Schwarzschild black branes, or equivalently stress-energy perturbations in strongly coupled N = 4 SYM plasma. The large wavenumber asymptotic behavior has the universal form (60), in all helicity channels, with mode-dependent spectral deviation coefficients shown in tables I, II and III. The relaxation rate of short wavelength quasinormal modes vanishes asymptotically as T 4/3 q −1/3 . We find that the large-q asymptotic form (60) of quasinormal mode frequencies agrees well with the exact values (for low order modes) already at rather moderate values of q/T , and thus provides a good approximation over a wide range of scales. The spectral deviation coefficients of high order modes approach the asymptotic form (55) but the coefficients of low order modes deviate from this simple expression.
In the strongly coupled SYM plasma it is noteworthy that hydrodynamic fluctuations are not the only arbitrarily long-lived excitations. The asymptotic vanishing of relaxation rates implies that there are two types of long-lived propagating excitations: long-wavelength sound waves, moving at v s = c/ √ 3, and short-wavelength fluctuations with group velocities asymptotically approaching c. As vividly seen in figure 9, illustrating examples of planar shock propagation, the decreasing attenuation plus increasing group velocity (as the wavevector q increases) combine to produce sharp, long-lived "spikes" which evolve from the most rapidly varying parts of an initial waveform. This same phenomena is evident in the study [13] by Chesler, Ho and Rajagopal of the "cyclotron radiation" produced by circular stirring of an SYM plasma (see fig. 3 of ref. [13]). Whether such long-lived "fine structure" could have observable phenomenological consequences in heavy ion collisions is an interesting open question.
There has been considerable discussion in the literature [7,[17][18][19][20][21][22] regarding "top-down" thermalization in strongly coupled SYM, as compared with "bottom-up" thermalization at weak coupling [23]. Evidence suggesting a top-down picture of thermalization at strong coupling comes from the decreasing relaxation times of quasinormal modes as the mode number increases at fixed wavenumber, and related observables probing similar physics. 17 However, interpreting this as topdown thermalization is, in our view, conflating the dephasing or decoherence of highly virtual off-shell excitations (|ω| 2 q 2 ), with relaxation of high momentum but near on-shell excitations (|ω − q| 2 q 2 ). It is these latter excitations, corresponding to large wavenumber but low order quasinormal modes, which should be considered when discussing top-down versus bottom-up thermalization. And these hard but low virtuality excitations thermalize slowly at both weak and strong coupling. 18 Finally, we limited our attention to metric (or stress-energy) perturbations. We expect quasinormal mode frequencies for perturbations in other supergravity fields to have the same universal asymptotic form (60), but it would be worthwhile to demonstrate (or disprove) this explicitly, especially for fermionic fluctuations. We leave such questions for future work.

Appendix B: Numerical techniques
For solving linear differential equations such as our quasinormal mode equations (9), (51) and (52), (pseudo)spectral methods are superior to traditional short-range discretization methods. Spectral methods converge faster, provide greater accuracy for a given number of discretization points, and allow one to easily enforce boundary conditions at either end of the computational domain without use of inefficient "shooting" techniques. 19 The basic approach, as sketched in section III A, is to represent the unknown function as a (truncated) expansion in a set of basis functions, and demand that the original differential equation be satisfied on a discrete set of points (the "collocation grid") within the computational interval. The optimal grid depends on the choice of basis functions. When using Chebyshev polynomials up to order M , the Chebyshev-Gauss-Lobatto grid points (45), consisting of the endpoints plus extrema of the highest order basis function, are an optimal grid. For functions which are analytic (in a neighborhood of the computational interval), the Chebyshev expansion converges exponentially rapidly with truncation size M .
To solve the helicity ±2 quasinormal mode equation (9), for a given numerical value of q, one may directly represent the unknown function h(u) as a Chebyshev sum (44), as the desired solution is regular at both u = 0 and u = 1. The radial equation (9)  The smallest eigenvalues (in absolute value) converge most rapidly as the truncation size M increases, with any given eigenvalue ω n (q; M ) showing exponential convergence for sufficiently large M . For any given value of M , the largest eigenvalues will always be sensitive to the truncation and hence dominated by discretization artifacts; at most some fraction of the smallest eigenvalues will be well converged. As the chosen value of the wavevector q increases, even the lowest quasinormal mode eigenfunction becomes highly oscillatory, and this necessitates the use of a truncation size M which grows linearly with q. For sufficiently large M , use of extended precision is also necessary to avoid excessive round-off error. For these reasons, a direct numerical solution of the quasinormal mode equation (9) becomes quite challenging for values of q beyond about 1000.
Such large-q computational difficulties are not present in the q-independent matching equation (37) which emerged from the WKB analysis of section II B. However, this equation needs to be solved on the positive halfline, and the equation has an irregular singular point at infinity plus a regular singular point at the origin. To find numerical solutions one may work either on the original halfline y ∈ R + , or on the rotated halfline (38) where y = e iπ/3 w with w ∈ R + . To be definite, we describe here our approach when working with the original form (37).
After writing the resulting equation in a form where all terms remain finite at u = 0 and 1, we arrive at Solutions of interest to eq. (B5) are now regular at both u = 0 and 1. Applying the same pseudospectral approximation scheme described above converts the differential equation to a generalized eigenvalue problem (with s ∞ α now the eigenvalue of interest). Before doing so, however, we make one final variable transformation, setting u = v 2 and using a Chebyshev-Gauss-Lobatto grid in v, as this was found to improve convergence of the spectral approximation. To obtain the results shown in table I, accurate to more than 12 digits, truncations up to M = 600 were used. 20 can match onto the decaying WKB-II solution at large v, so the solution within the turning point region is f TP (v) = Ai κ(v−2)/2 1/3 . (D10) For v−2 κ −1 , the asymptotic behavior of this Airy function coincides with the near turning point behavior (D9) of the WKB-II solution, as required. On the other side of the turning point, when 2−v κ −1 , the asymptotic behavior of the Airy function with large negative argument gives This agrees with the oscillatory behavior (D6) of the WKB-I solution near the turning point, up to a shift in the phase. For a consistent solution, these phase shifts must also agree modulo π (since a difference of π can be absorbed by a sign flip of an overall coefficient). Consequently, we require that for some integer n. Solving for the eigenvalue λ = κ 2 and inserting the value (D7) of I yields the next-to-leading approximation for large order eigenvalues, Inclusion of higher order terms in the WKB expansion will generate relative corrections to this result of order κ −3 ∼ n −2 . One may verify that the allowed region solution (D5) has n−1 nodes when λ = λ WKB n implying that, as written, n is the level number when counting starts from 1. Recalling [from eq. (39)] that the eigenvalue λ = 1 2 s ∞ α e iπ/3 = 1 2 c n one finds the result (41) quoted earlier for the large order behavior of the asymptotic coefficients {c n }.