A Holographic Superfluid Symphony

We study the hydrodynamic excitations of backreacted holographic superfluids by computing the full set of quasinormal modes (QNMs) at finite momentum and matching them to the existing hydrodynamic theory of superfluids. Additionally, we analyze the behavior of the low-energy excitations in real frequency and complex momentum, going beyond the standard QNM picture. Finally, we carry out a novel type of study of the model by computing the support of the hydrodynamic modes across the phase diagram. We achieve this by determining the support of the corresponding QNMs on the different operators in the dual theory, both in complex frequency and complex momentum space. From the support, we are able to reconstruct the hydrodynamic dispersion relations using the hydrodynamic constitutive relations. Our analysis rules out a role-reversal phenomenon between first and second sound in this model, contrary to results obtained in a weakly coupled field theory framework.


Introduction
The holographic superfluid model [1][2][3][4], often mislabelled as holographic superconductor, 1 is one of the most simple and widely studied setups in the context of applied holography [6][7][8][9]. Its potential application in the understanding of high-T c superconductors [10] has made this system the object of many thorough analyses in the years since its inception (see [11,12] for a review). Thus, the lack of a complete study of the quasinormal mode excitations corresponding to hydrodynamic modes of the holographic superfluid [13] is surprising. Until now, the study of the hydrodynamic modes has been limited to the so-called probe limit, where all the gravitational degrees of freedom (and in particular energy and momentum) are kept frozen [14][15][16][17].
In this work we fill this gap in the literature by carrying out the first complete study of the quasinormal modes of the holographic s-wave superfluid with backreaction [13]. 2 In particular, we numerically compute the spectrum of the lowest lying quasinormal modes for the original setup [13], and also for the modified system [30][31][32][33] where momentum dissipation is introduced using the so-called holographic axion model [34,35]. We successfully compare our results to the predictions of superfluid hydrodynamics [36] extending the results of [37,38]. This study of the QNMs will allow us to carry out a novel analysis of the support of the hydrodynamic modes of the holographic superfluid. Superfluids display a characteristic coexistence of two different sound modes [39] which can be easily captured within the two-fluid Tisza-Landau model [40,41]. In addition to the standard density fluctuations, giving rise to the common first sound excitation, superfluids support also a second sound mode which is driven by temperature fluctuations. First sound can be visualized as the in-phase motion of the normal fluid component and the superfluid one; on the contrary, second sound is displayed as the out-of-phase collective dynamics of the two. In [42], a field theoretical study of the nature of first and second sound as density or entropy waves or a mixture of the two was performed. Even though in superfluid 4 He, the nature of first and second sound is preserved for all temperatures, this seems not to be the case for ultra-cold atoms and weaklyinteracting bosons [42]. Here, we carry out a similar analysis within the holographic superfluid which stands as a putative toy model for a strongly coupled superfluid system: we determine the support of the hydrodynamic modes of the system on excitations of the different operators present in the dual theory.
In order to analyze the nature of the hydrodynamic excitations of the holographic superfluid we define a quantitative measure of the support that a quasinormal mode excitation has on the different dual operators. We then apply this procedure to the sound and diffusive modes present in our gravity dual. Our results rule out any role reversal between first and second sound and place the holographic model in the same phenomenological class of standard superfluid 4 He. This is in contrast to the weakly coupled bosonic field theory results of [42] described above. Naturally, an obvious difference between both setups is given by the inherently strongly coupled nature of the holographic description. Moreover, [42] is neglecting dissipating effects in their description.
In a final part of this work we study the response of the system from a different perspective. One is typically interested in the poles of (the Fourier-transformed) Green's functions in the complex frequency and real momentum space. These are the QNMs. Nevertheless one can also consider the poles in the space of real frequency and complex momentum. These modes carry information about the response of the system to a perturbation of the form e −iωt+ikx where ω now is a real number and k a complex one. Although this approach is not usual in holography, it is a standard practice in the study of propagation of waves, in particular in relation to absorption and reflection properties of media (e.g. EM waves in a medium [43]). Once momentum is decomposed as k = Re(k) + i Im(k), its real part defines the propagation wavelength in real space while its imaginary part the so-called penetration length λ ≡ 1/Im(k). In other words, this second setup corresponds to considering the decay of a wave in space rather than in time.
The curves ω(k) described by the QNMs in the complex plane feature a complicated structure (e.g. the presence of different sheets, and singularities). Hence, a priori it is not clear that the behavior of the modes at real frequency and complex momentum can be obtained straightforwardly from that of the standard QNMs. Indeed, notice that these two types of modes encode responses to perturbations with different boundary conditions. Consider for instance an elastic rod immersed in a viscous liquid. Considering complex frequency and real momentum corresponds to keeping the extreme of the rod fixed and creating an oscillating wave on it. The wave in that case will not propagate in space but just oscillate (and decay) in time. On the contrary, one could imagine exciting a wave at one of the extremes of the rod and studying its propagation along the rod. That is equivalent to considering real frequency and complex momentum.
In the context of holography, the modes at real frequency and complex momentum have been overlooked and only few works considered this second setup [44][45][46][47][48]. On a different note, recently, several studies of the modes in the complex frequency/complex momentum plane have been performed in relation to the so-called pole-skipping phenomenon [49] and the convergence radius of the linearized hydrodynamics expansion [50]. In the last part of our work, we study the excitations of a holographic superfluid at real frequency and complex momentum and discuss their physical interpretation in detail.
Structure of the paper.
In section 2, we first summarize the main results of the hydrodynamic description of a relativistic superfluid. We then introduce the holographic model and set the framework of the analysis we carry out in the remainder sections. In section 3, we compute the quasinormal modes of the backreacted holographic superfluid model numerically and use the results to match holographic and hydrodynamic predictions for the system. In section 4, we perform a new analysis of the support of the various hydrodynamic modes in terms of the dual field theory operators. Finally, in section 5, we compute the excitations of the system and their main properties by considering real frequency and complex momentum. We conclude in section 6 with a summary of our results and some comments for the future.

Superfluids in hydrodynamics and holography
In this section, we first review the main features of the hydrodynamic description of a conformal superfluid. Next, we introduce the holographic dual of a superfluid and set up the tools for the analysis performed in the remaining of this work.
Let us start by considering a relativistic charged fluid with conserved momentum [64,65]. The longitudinal spectrum contains two sets of hydrodynamic modes. First, charge and energy conservation imply the presence of a thermoelectric diffusive mode ω = −iD q k 2 whose diffusion constant is explicitly given later on in eq. (2.5). Additionally, we have a pair of propagating first sound modes ω = ±c 1 k − i Γ 1 k 2 which result from the interplay of longitudinal momentum and energy. In a conformal system, their propagation speed is given by where and p are the energy density and the thermodynamic pressure, respectively. Additionally, the attenuation constant reads with η being the shear viscosity and χ ππ = + p the momentum susceptibility. Below the critical temperature, a superfluid may be viewed in the two-fluid picture as composed of two components: a normal component as the (relativistic) fluid just described; and a condensed or superfluid component that flows without dissipation. Each component has its own velocity and density. In the normal phase, the interplay of charge, energy and longitudinal momentum gives rise to a propagating mode: first sound, and a diffusive mode. In the superfluid phase T < T c , first sound remains in the spectrum, the diffusive mode gets a purely imaginary gap, and the Goldstone fluctuations, coupled to those of charge and energy give rise to a new propagating mode: second sound ω = ±c 2 k − i Γ 2 k 2 , known as second sound [39]. Its speed and attenuation constant are given by Here we have defined: the superfluid density ρ s , the normal density ρ n = ρ − ρ s , the reduced entropy ξ = s/ρ, the diffusivity κ and the superfluid diffusivity ζ. Notice the only difference with the results presented in [36], and valid for d = 4, is the numerical constant in front of the term proportional to the shear viscosity η.
In the normal phase, ρ s = 0, the first sound remains in the spectrum and the second sound becomes the thermoelectric diffusive mode [66,67] Once momentum conservation is broken, the first sound in the normal phase gets substituted by the energy diffusion mode ω = −iD e k 2 whose diffusion constant obeys the standard Einstein relation. At the same time, in the superfluid phase, the second sound mode is replaced by the so-called fourth sound first predicted by [68,69] where we keep the normal component stationary. The dispersion relation of fourth sound reads ω = ±c 4 k − i Γ 4 k 2 , where the velocity and attenuation constant are given by For completeness, we list here all the Kubo formulas used to extract the various transport and thermodynamic parameters entering in the dispersion relations above together with the expression for the electric conductivity For slow momentum dissipation this two-point function is modified to In this discussion, φ refers to the superfluid Goldstone mode which is dual to the phase of the complex bulk scalar (see next section for the details). Also, G AB is the retarded Greens function of the operators A and B. These Kubo formulas will be used in section 3 when we analyze the transport properties of holographic superfluids.

The holographic model
We consider the standard bulk action for a s-wave holographic superfluid [1,13] where ψ is a complex scalar, F ≡ dA, and the covariant derivative is defined as We consider the 4-dimensional bulk case d = 3, and fix Λ = −3 and 2κ 2 4 = 1. The AdS radius L has also been set to one for convenience. The resulting equations of motion are We consider the following ansatz for the background in Eddington-Finkelstein coordinates {t, z, x, y} with the AdS boundary located at z = 0 and the black brane horizon at z h = 1. Assuming this ansatz, the temperature and chemical potential of the dual field theory are given as For numerical convenience, we impose A t (1) = 0. Furthermore, we choose the value of the background condensate to be real by setting ψ 1 . (2.17) The analysis in sections 3 and 4 will rely on the study of the fluctuations about the background discussed above. In particular, we switch on the full set of longitudinal perturbations g mn where we choose the momentum to be aligned with the x-direction without loss of generality. From now on, we will add to all symbols referring to background quantities a hatÂ. Moreover, we assume the radial gauge 4 We decompose all our fluctuations as The quasinormal modes and all the transport coefficients can be obtained following standard methods. To compute the thermodynamic derivatives accurately, for example the one appearing in eq. (2.3), we re-write the derivative in terms of dimensionless quantities; we then compute the background on a Chebychev grid in that variable (for example T /µ) and use the spectral differentiation matrices to compute the thermodynamic derivative (see appendix B.4 in [70] for more details).
In order to introduce momentum dissipation and study the properties of fourth sound, we supplement the action (2.12) by a linear axion term [34] where the index I runs only along the spatial directions I = x, y, and the background profile for the scalars isφ In the dual field theory, this corresponds to switching on an operator which breaks translational invariance explicitly. For more details about this model we refer to the recent review [35]. In presence of momentum dissipation, we add the fluctuation of the longitudinal axion field δφ x to the fluctuations in eq. (2.18).
Finally, throughout all the manuscript we set q = 1 and M 2 = −2.
4 Note that since we work in Eddington-Finkelstein coordinates, the bulk solution for the scalar field is complex even though we impose the ψ (2) 2 to be zero.

Transport & hydrodynamic modes
In this section, we explicitly verify that the hydrodynamic framework described in the previous section (see [36] for the original reference, albeit in one dimension more) is indeed the correct hydrodynamic description of the holographic superfluid model [1]. Surprisingly, even though partial checks may be found in the literature [36][37][38], a complete matching of the two pictures considering the full backreaction [2] in the gravitational model has not appeared yet.  We start by plotting the normal and superfluid densities as function of T /T c in the superfluid phase. The normal and superfluid density, respectively, may be extracted numerically from eq. (2.10) (translationally invariant case) and eq. (2.11) (broken translations). While in the latter case, the superfluid density is simply given by the i/ω pole in the conductivity, we have to use that ρ n = ρ − ρ s and solve eq. (2.10) for ρ s in the translationally invariant case. The results are shown in fig. 1 for both the translational invariant model and the model where translations are broken explicitly using axion fields [35,71]. As it should, the superfluid density vanishes at the critical temperature for both models. In the translational invariant case, we also observe that the normal density vanishes at zero temperature. This last feature holds in superfluid Helium and BCS superconductors, while a residual normal density at zero temperature has been reported in some high-Tc superconductors [72]. Interestingly, [73] has shown that in quantum critical superfluids the nature of the IR fixed point of the theory is crucial for the fate of the normal density: it vanishes for fixed points with an emergent Lorentz symmetry as in our translation invariant case [74] , while some quantum critical points of the hyperscaling-Lifshitz type can support a nonvanishing normal density [75]. The model breaking translations (right panel of fig. 1) shows a slower decay of the normal density for T → 0 and a fit of the numerical data shows a scaling of ρ n /ρ ∼ 0.05265 + 3.14843(T /µ) 1.2583 for T /µ ∈ [5 · 10 −7 , 5 · 10 −6 ] indicating a non-zero normal density in the zero temperature limit. It would be interesting to characterize the IR geometry of our model and analyze it in the light of the results of [73]. We leave a more accurate analysis of this feature for the future.
We then move to the study of the sound modes. On the gravity side these correspond to quasinormal modes of the black hole geometry. As we review in section 4.1, with further technical details in appendix D, we obtain them by numerically solving the linearized equations of motion for the fluctuations (2.18). This will allow us to compute the speed and attenuation of the different sound modes and compare them against the hydrodynamic predictions of section 2.1. The propagation speed and attenuation constant of the first sound mode are shown in fig. 2. The propagation speed is fixed by conformal symmetry to the expected value c 2 1 = 1/2 and it is continuous across the phase transition at T = T c . Its dimensionless attenuation constant Γ 1 T is shown in the right panel of fig. 2 and it is in perfect agreement with the hydrodynamic formula in eq. (2.2) in which all the parameters are independently computed using thermodynamic relations and transport Kubo formulas. One should notice that the dimensionless attenuation constant vanishes at zero temperature, indicating that at T = 0 first sound is an infinite living excitation. This is tantamount to saying that the shear viscosity of the system vanishes at T = 0. Finally, we observe that the attenuation constant of first sound is continuous at T = T c even though its first derivative is not. T/μ  2) and the colored circles the numerical data extracted from the QNMs. Blue color indicates the normal phase. The plot for the speed is not extended into the normal phase since the speed is trivially constant in the full range of T due to the conformal symmetry.
We perform an analogous analysis for second sound in fig. 3. We remind the reader that the second sound mode exists only in the superfluid phase below the critical temperature T < T c . Interestingly, we observe a very non-monotonic behavior for the second sound speed which vanishes both at the critical temperature T = T c and at zero temperature T = 0, with a maximum around T = 0.6 T c . Apparently, this behavior is not universal and it rather depends on the charge and the conformal dimension of the dual scalar operator [37]. Again the hydrodynamic formula (2.3) fits the numerical data very well.
The attenuation constant of second sound is more interesting. First, it vanishes at zero temperature. Second, across the phase transition it is connected to the imaginary part of the gapped scalar modes present in the normal phase [76] and responsible for the superfluid instability at T = T c . Nevertheless, the data display a clear and sharp jump around T = T c where our hydrodynamic description is no longer applicable as is evident in the right panel of fig. 3. Moreover, the attenuation constant of second sound is not continuously connected to the diffusion constant which we present in the inset in the same figure 3. Our numerical observation was confirmed by the authors of [77] in a simpler setup (without chemical potential) using analytic arguments. T/μ Finally, when momentum is dissipated and the normal component is prevented from flowing, a new excitation called fourth sound appears in the superfluid [68,69]. We obtained the speed of propagation and the attenuation constant numerically in fig. 4 and compared them with the hydrodynamic expressions in eq. (2.6). As expected the speed of fourth sound interpolates between the speed of second sound at the critical temperature T c to the speed of first sound for T → 0 [69,78]. Additionally, notice how the attenuation constant of fourth sound is much larger than those of first and second sounds. Finally, we observe that the speed and attenuation constant of fourth sound have the same qualitative behavior as those of the sound mode obtained in the probe limit in [15]. We do expect the two to coincide in the limit of very fast momentum dissipation, α/T 1.
To conclude this first part, let us iterate that the hydrodynamic description presented in the previous section perfectly matches the numerical results obtained from the holographic superfluid model both in the normal and superfluid phases with and without momentum dissipation.

Tomography of superfluid sound modes
After verifying the hydrodynamic theory in our holographic model explicitly, we now move on to investigate the so-called support of the hydrodynamic modes.

Theory
A quasinormal mode is a solution to the linearized equations of motion in a black hole background with in-going boundary conditions on the horizon and normalizable ones at the boundary of AdS [79,80]. In general, such a solution will consist of a collection of fields which we denote by Φ I . For a single QNM all these fields oscillate with the same frequency and decay exponentially with the same decay rate. At late times and close to equilibrium, we may write a generic QNM solution as a superposition of modes The bulk wave functions of a QNM A I n ( k, z) do not contain the leading (nonnormalizable) terms of the asymptotic boundary condition. We can therefore multiply the field φ I with the appropriate power of z such that A I n (0, k) is the coefficient of the non-normalizable mode. According to the holographic dictionary the boundary value of Φ I n encodes the expectation value of an operator O I (t, x) (see appendix A). Since the solution (4.1) evaluated at the boundary represents a propagating and attenuated sound wave in the dual field theory, it is instructive to think about the quasi-normal mode as corresponding to a hydrodynamic sound wave. One should notice that sound waves consist of more than one perturbation, for example the perturbations of the energy density and the pressure p (as is clear from the thermodynamic formula for the speed of sound c 2 s = ∂p/∂ ). Hence, we expect the QNM solution dual to a sound mode to consist of more than one field Φ I and taking the form of eq. (4.1).
To simplify the analysis further, we assume that rotational symmetry is not broken. Then the quasinormal frequencies are function of | k| 2 only. From now on we write k for | k|. Since the collection Φ I includes vector fields and metric components the dependence of the amplitudes may be more complicated in general.
Demanding that the fields Φ I be real results in the following condition for the quasi-normal frequencies and amplitudes Combining this constraint with (4.1) leads to the condition that the quasinormal fequencies come in pairs (n, m) with or are purely damped with vanishing Ω n and A I * n (z, − k) = A I n (z, k). In view of this structure, the minimal way of constructing a real QNM solution for Φ I is to excite the mode n at wave number k and the mirror mode m with Ω m = −Ω n at wave number − k.
If we evaluate this at the boundary, we get the space-time evolution of the set of operators O I as where we take O I (0, k) to be positive since the phase α I n accounts for the sign of O I (t, x) . We need to take into account, however, that (4.4) represents a solution to the linearized equations of motion. This means that we can multiply it with an arbitrary complex number and obtain another valid solution. Therefore only relative amplitudes and phases contain physical information. Moreover, since a generic solution will involve operators O I with different dimensions, we should only compare dimensionless expectation values that can be easily obtained by normalizing with the corresponding power of the temperature, namely We also denote the dimensionless absolute values of the amplitudes as R I . We can then write our normalized operators as A measure of the relative support on the different operators O I of a given QNM is thus given by the ratio of their amplitudes R I n . Namely, r JI n = R J n /R I n measures the relative support of a mode on the operator O I with respect to the operator O J . In addition, the relative phases φ IJ = α I − α J encode the phase delay between the excitations of operators O I and O J . 5 A useful way of thinking about these ratios and relative phases of the expectation values, would be to consider a (thought) experiment where a device couples directly to the operator O I and only to that one. This device would act as a source of O I alone among all other operators. But the operator mixing inherent to the system would result in a state where operators other than O I are also excited. Then, the ratios r JI n and phases φ JI will precisely characterize the response of the operator O J when a source is switched on for the operator O I .
A (thorough) study of relative amplitudes and phases of QNMs is lacking in the literature, although closely related analyses of the residues of holographic Green's functions were carried out in [45,81,82]. We will do this here for a particularly interesting example, namely the first and second sound modes of a holographic superconductor. Notice that a similar study of relative phases and amplitudes for sound modes of a weakly coupled superfluid has been carried out in [83], and we will 5 We could also consider taking the logarithm of the ratios r IJ n and definingr IJ n = log(R J n ) − log(R I n ). This has the advantage that ther JI n are antisymmetric in JI as are the relative phases φ IJ . eventually compare our results to it.
Finally, we shall comment on the difference between the ratios of amplitudes that we will determine here, and the residues of a certain mode ω(k). The residue of a mode quantifies the contribution of that mode to a given response: i.e. the weight of its corresponding pole in the spectral representation of a given Green's function. Instead, the amplitudes we discuss here represent the weight of a given operator in the collective excitation corresponding to a specific mode. In a sense, these observables provide complementary information about the system.

Results
For concreteness, we study the longitudinal spectrum of excitations, to which first and second sounds belong. This amounts to solving for the set of fluctuations in eq. (2.18), as we describe at the end of section 2.2. Moreover, we will consider solutions corresponding to QNMs: these are dual to configurations in which all the sources of the dual operators vanish, while (some of) the expectation values do not. In appendix A, we work out the map between fluctuations and dual operators. As a result, in the dual theory we will be dealing with the following set of operators corresponding to the expectation values of the fluctuations δT tt : energy density, In eq. (A.6) we express these expectation values in terms of the asymptotics of the fluctuations (2.18). Notice that in the superfluid phase we are rewriting the operator dual to the fluctuation of the complex scalar in terms of its modulus and phase, which we denote by δH and ϕ respectively. Moreover, for the model breaking translations, we switch on an extra fluctuation, δφ x . The meaning of all of these quantities is well known apart from δφ x , corresponding to the longitudinal axion operator. The colored circles in (4.7) display the color scheme used in all the figures in this section.
In the remainder of this section we will present results for the amplitudes and phases of the expectation values (4.7) for the different QNMs of the system. These are all computed in Fourier space (see eq. (2.20)), at finite frequency and momentum. Finally, we restrict our analysis of the amplitudes and phases to the following hydrodynamic modes 1. first and second sounds in the superfluid phase (with preserved translations); 2. fourth sound and energy diffusion in the superfluid phase with broken translations; 3. first sound and a purely diffusive mode in the normal phase (with preserved translations); 4. thermo-electric diffusion modes (energy and charge diffusion in the decoupling limit) in the normal phase with broken translations.
Additionally, in the normal phase we shall also consider the gapped scalar mode that drives the superfluid phase transition at T = T c .
Finally, let us explain our normalization of the amplitudes. We may formulate the linearized equations for the fluctuations in terms of a generalized eigenvalue problem (see appendix D for more details) where A and B are differential operators and x is a vector consisting of the fluctuations. Solving this generalized eigenvalue problem numerically, we find a set of eigenfunctions x n (consisting of the fluctuations) for each eigenvalue ω n (the quasinormal frequency). We will refer to the set of eigenfunctions as eigenvector from now on. First, we make all expectation values of the eigenvectors consisting of the fluctuations (4.7) dimensionless by dividing by the appropriate power of the temperature (depending on the conformal dimension of the corresponding operator). Then, we normalize the sum of all absolute values of the dimensionless expectations values of the eigenvectors to one. This means that an amplitude of 0.5 implies for example that 50% of the support of such a mode is controlled by this type of fluctuation. With this normalization, the amplitudes are invariant under the scaling symmetry inherent to the linear analysis which acts on all fluctuations as δf i → λ δf i .
In order to make the reader familiar with the language, we start our analysis in the normal phase which corresponds to considering a standard relativistic charged fluid in the field theory dual. There, the two hydrodynamic modes are first sound and the purely diffusive mode. Nevertheless, for completeness we will consider also the gapped scalar modes which are the responsible for the superfluid instability at 0.01  T ≤ T c [76]. Notice that in the normal phase the background value of the scalar is zero and therefore the fluctuations of the phase are not a well defined object. Thus, in the normal phase, we will consider the expectation value of the fluctuations of the imaginary part of the scalar operator dual to the complex bulk field ψ.
The results for the amplitudes of the normal fluid are shown in fig. 5. As expected, we find that the gapped scalar mode is solely carried by the fluctuations of the scalar charge operator. The purely diffusive mode (left panel of fig. 5) is carried by the fluctuations of the charge density δJ t , as one can derive analytically from hydrodynamics (see appendix (C.3)). The hydrodynamic treatment shows that in first order hydrodynamics this mode consists solely of δJ t (and of the by a Ward identity related δJ t ). Moreover, the second order correction includes δT tt with a relative amplitude of δT tt /( δJ t T ) = k 2 /T ρ/(µρ + s T ). We shall point out that our result that the diffusive mode is solely carried by δJ t is not in contradiction with the findings of [66,84] where it was observed that the diffusive mode is carried by the 'incoherent charge' δQ diff ≡ δJ t − ρ/(µρ + sT ) δT tt 6 . As we detail in appendix C, our computation describes the situation where only the diffusive mode is excited in the system. The authors of [66,84] consider a generic perturbation in which all modes are excited instead. We show in that appendix that one can relate both situations via the hydrodynamic framework and match them to the holographic computation.
Finally, first sound is dominated by the fluctuations of energy density and is almost entirely supported on the fluctuations of energy density, pressure and longitudinal momentum. The ratio between pressure and energy fluctuations is approximately 0.5 as we expect from the Ward identities that constrain the expectation values given in eq. (A.10). Note that the pressure fluctuations contribute approximately equally via δT xx and δT yy to first sound as is obvious from the tracelessness of the energy-momentum tensor and the Ward identities (A.10). We furthermore observe  fig. 5 and in eq. (B.3). The simple relation eq. (B.3) is Fick's law for diffusion. This discussion shows that we can use the amplitudes and the Ward identities to compute the diffusion constant in the hydrodynamic dispersion relation and vice versa. In the case of first sound a similar analysis holds albeit it is a bit more complicated. We refer the interested reader to appendix B. In the normal phase, the scalar sector decouples and the fluctuations of the scalar operator do not enter into either of the sound modes in the normal phase.   Figure 6: The amplitudes of the various operators (4.7) for the thermo-electric diffusive modes in the normal phase with momentum dissipation. We fixed k/T = 0.01 and T /µ = 0.025. The color scheme is that illustrated in (4.7).
As already mentioned, when the axion fields are switched on and momentum is not a conserved quantity anymore, first sound is destroyed at large wavelength and the only two hydrodynamic modes in the normal phase are the thermo-electric diffusive modes [85]. We show the amplitudes for these two modes in fig. 6. At α = 0, the left panel corresponds to the diffusive mode shown in fig. 5. Indeed, for α/T 1, such a mode is totally dominated by charge fluctuations. Interestingly, around α/T ≈ 8, we observe a smooth crossover to a regime where the original diffusive mode is now dominated by energy fluctuations. Moreover, the ratio between pressure and energy fluctuations is approximately 0.5 for all α/T , as we already observed for the first sound mode which is a consequence of the tracelessness of the energy-momentum tensor (since δT tt = δT xx + δT yy ). Notice that in that large α/T 1 regime the contribution from charge fluctuations drops to zero indicating a drastic change in the nature of the diffusive mode.
As shown on the right panel of fig. 6, the situation is different for the other diffusive mode present in the normal phase of the model with momentum dissipation. For small values of momentum dissipation, α/T 2, the momentum chosen for the extraction of the amplitudes, k/T = 0.01, is too large and the mode is not diffusive anymore but already sound-like, i.e. propagating. There, we find the same features as for first sound in fig. 5. Around α/T ≈ 2, momentum is not a well-defined quantity anymore, and its contribution to the mode drops rapidly to zero. At very large value of momentum dissipation, the mode in the right panel of fig. 6 is totally decoupled from everything else and it is dominated by energy and pressure. Note that for α/T 200, we enter the regime where the scalar modes are unstable since the critical temperature (which is a function of the momentum dissipation strength) is now larger the temperature we chose. Similarly to the diffusive mode, we observe that the ratio between energy and pressure fluctuations is approximately 0.5 for all α/T .   Before moving on to the broken phase, let us investigate the normal phase with broken translations further; in particular, we are interested in the two thermoelectric diffusion modes. At intermediate values of T /µ and α/T , energy diffusion and charge diffusion are coupled and therefore their support is mixed between the various operators. This is no longer true in the so-called incoherent limit [66], at which the momentum dissipation rate becomes the dominant scale in the system, i.e. α/T, α/µ 1. In such a regime, we do expect the two diffusive modes to decouple again [86] and the corresponding diffusion constants acquire the simple values where κ is the thermal conductivity, c v the specific heat, σ the electric conductivity and χ ρρ the charge susceptibility. For the linear axion model considered in this work, the formulas for the diffusion constants above are known analytically and can be found in the literature (see for example [87]). In fig. 7, we repeat the analysis of fig. 6 but this time going into the aforementioned incoherent limit. As expected, in this case, the contributions to the support of the two diffusive modes from charge and energy fluctuations crosses around α/T ∼ 10, which corresponds exactly to the location at which the diffusion constants approach the values mentioned in eq. (4.9), i.e. the edge of the incoherent limit. Going towards α/T → ∞, the support of the two modes becomes solely dominated by either charge or charge and energy fluctuations, as expected from their decoupling. This result is explained in more detail in appendix C.2 where we show the relative coefficients in fig. 14 explicitly. Note that we extracted the diffusion constants from the amplitude data at one fixed k/T using eq. (B.3). Not only that, but the amplitudes of these two diffusive modes can be derived analytically at any value of the parameters from hydrodynamics (see appendix (C.3)), matching perfectly the numerical results presented here. We now move on to the analysis of the superfluid phase of our system. We study the amplitudes of the operators (4.7) for the two hydrodynamic modes in the longitudinal sector of this phase: first and second sound. The results are plotted in fig. 8. Notice that in this condensed phase we represent the fluctuations of the scalar operator in terms of its amplitude (the Higgs mode δH) and the time derivative of its phase which carries information about the norm of the Goldstone boson since in our background |A µ − i ∂ µ ϕ| ∼ µ − i ∂ t ϕ + O(ϕ 2 ). As one can see on the left panel of fig. 8, the behavior of the amplitudes of the first sound mode is rather featureless and dominated by energy, pressure, and longitudinal momentum fluctuations across the phase transition and down to the lowest temperatures. The fluctuations of the Goldstone mode are highly suppressed, and those of charge and current decrease as T is lowered. We show the results for second sound on the right panel of 8. This mode is dominated by the fluctuations of the superfluid condensate (the Higgs mode). In view of these results one can rule out a role reversal between first and second sound T/μ in which the nature of the fluctuations supporting those modes is interchanged as the temperature is lowered across the superfluid phase as in [83]. Finally, notice that the fluctuations of the scalar operator are not shown above the critical temperature in the left panel of fig. 8. They decouple from the other operators and therefore do not contribute to the support of the first sound mode above T c . T/μ

T/μ
A norm T/T c Figure 9: The amplitudes of the various operators (4.7) for the diffusive mode (left) and 4 th sound (right) in the superfluid phase with momentum dissipation. We fixed α/T = 6 and k/T = 0.01. The color scheme is that illustrated in (4.7).
We can now perform the same analysis in the superfluid phase in presence of momentum dissipation. In this case, first sound is replaced by a diffusive mode and second sound is replaced by fourth. In fig. 9 we show the results for the amplitudes of the fluctuations for the diffusive mode (left) and fourth sound (right). Except at temperatures near the critical temperature, where the Higgs mode takes over, the diffusive mode is dominated by energy and pressure fluctuations. Fourth sound is always dominated by energy fluctuations. Interestingly, at low temperature the ratios between energy, longitudinal momentum, and pressure fluctuations asymptote to the same ratios observed for first sound since the speed of fourth sound tends to the speed of first sound as T → 0 there (while the attenuation constant tends to zero). Finally, we shall point out that we also observe a purely imaginary, gapped non-hydrodynamic mode which is solely supported by the scalar fluctuations. This higher quasi-normal mode might be the Higgs mode discussed in [88][89][90][91].

About the relative phases
As explained in detail in section 4.1, the value of the amplitudes of the various fluctuations (4.7) of the different (hydrodynamic) modes is not the only information we can extract from our analysis. We can also compute the phase differences between the various fluctuations present in a specific (hydrodynamic) mode. We used the phase of the charge density fluctuation Θ δJ t as reference phase and computed the phases of the other fluctuations with respect to Θ δJ t . Note that the phases are defined modulo 2π. The phase differences of the various fluctuations are listed in table 1. The phases differences remain almost constant throughout the phase diagram. However, it is important to notice that the "ideal" phase differences as listed in 1 are only reached in the k → 0 limit. As we show in eq. (B.4) and eq. (B.9), the relative phases between certain fluctuations are intimately related to damping, i.e. for vanishing phase differences we reach ideal hydrodynamics. In other words, the relative phases give us insights about the damping in the system and we may compute the diffusion and attenuation constants, respectively, from them.

Momentum dependence of the amplitudes
So far, all the fluctuations amplitudes have been computed for a single value of momentum, k/T = 0.01. In this section, we briefly discuss the momentum dependence of the amplitudes presented in section 4.2. The results are shown in fig. 10. For first sound, both the amplitudes and the relative phases are constant in a large range of momenta (we considered them until k/T ≈ 12). In the case of second sound, we observe that the biggest amplitude, the Higgs fluctuations, remains dominant until k/T ≈ 10 and the relative phases approximately constant. Away from the small momentum limit, we notice that the contribution from the fluctuations of the superfluid condensate becomes smaller and the mode is eventually dominated by the energy fluctuations which are the biggest contribution at large k. At the same time, at k/T ∼ 1, we observe a crossover between the fluctuations of the charge density (which is the first subleading contribution at small k, even thought it accounts for less than 10% of the total weight) and that of energy and longitudinal momentum without momentum dissipation π π π π 0 0 0 2 nd sound with Re(ω) > 0 (T < T c ) π/2 π/2 π/2 π/2 0 0 0 with momentum dissipation 4 nd sound with Re(ω) > 0 (T < T c ) π π π π 0 0 0 π/2 diffusion (T < T c ) π π/2 π π −π/2 π 0 π/2 Table 1: The phase differences modulo 2π with respect to δJ t for the various fluctuations for the normal fluid (T > T c ) and superfluid (T < T c ) phases for temperatures sufficiently far from the phase transition and in the ideal limit k → 0. The gray cells indicate that such a fluctuation is not present in the corresponding mode. In the case of the sound modes, the phases of the "partner" mode may be obtained by taking into account that the phases add up to 0 (for parity even operators) and π (for parity odd operators), respectively. Concretely, the computation of the phases for the momentum dissipation case has been performed for α/T = 6.
which become more relevant at large wave-vector. All in all, we can state that, in the hydrodynamic limit k/T 1, our results depend only minimally on the specific value of the momentum k.

Theory
Quasinormal modes can be understood as the response of the system to a perturbation that happened in the past. Thus, quasinormal modes describe how the system relaxes back to equilibrium at late times.
A different way of studying the response is by perturbing the system with a periodic source of frequency ω and looking for the corresponding wave with complex wave number q = k + iκ. If we were to excite only one particular mode then the response would be Φ I (z, t, x) =Ã I n (ω, z)e −iωt+iknx e −κn|x| .
The complex momentum is now a function of the real frequency k n = k n (ω) and κ n = κ n (ω). This was investigated in [45] where it was shown that one can isolate the modes that decay exponentially by choosing the integration in the complex momentum plane appropriately. Only a few other works considered this choice within the holographic framework [44,46,47].
We can now proceed with the same attitude as in the previous section and investigate the amplitudes and the relative phases. Whether these are different from the ones of the quasinormal modes is an open question and it is our main motivation for this last analysis. A naive expectation would be that the phases and amplitude ratios are the same as long as the mode in question has not gone through some crossing point with another mode in the complex plane. Generally, such a crossing point has complex frequency and complex momentum so its neither visible from the quasinormal mode (real momentum) nor from the complex momentum mode (real frequency).

Results
Before discussing the amplitudes and the phases of the various operators, let us briefly comment on the nature of the low-energy modes for complex momentum and real frequency. Let us start with the equation for a generic attenuated sound wave where as always v, Γ are the propagation speed and the attenuation constant, respectively. Let us now consider eq. (5.2) in terms of a real-valued frequency and a complex valued momentum. At low frequency, it is easy to verify that The same analysis can be repeated with a damped diffusive mode (such as the diffusive mode in the superfluid phase) The equations (5.3) and (5.5) represent the low-momentum behavior of the lowenergy modes in this new inverted picture. As expected, our numerical results shown in fig. 11 are in perfect agreement with this picture. We can now analyze the ampli- tudes of the fluctuations in (4.7) at constant frequency. We show them in fig. 12 for a fixed small value of (real) frequency ω/T = 0.01. The results for first sound are exactly identical to those obtained at complex frequency and real momentum (inset in the left panel of fig. 8). In the second sound case shown on the right panel of fig. 12 (to be compared with the right panel of fig. 8), the dominant contributions corresponding to the condensate and charge density fluctuations are identical while the subleading contributions are slightly different between the real and imaginary frequency modes. Still the qualitative behavior in the two frameworks is very similar. T/μ

Conclusions
In this work we have revisited the holographic S-wave superfluid model [2] focusing on the hydrodynamic description of the low-energy modes and taking into account the backreaction on the gravity side. Despite the long history of this model, a complete analysis of this sort was missing. Through this approach we have explicitly verified the match between relativistic superfluid hydrodynamics and the holographic model, both in terms of the dispersion relations and, most importantly, in terms of all the first order transport coefficients. Our study has revealed that all the hydrodynamic transport coefficients are continuous at the phase transition. Moreover, differently from what [15] found in the probe limit, even the second sound attenuation constant is smoothly connected to the gapped scalar modes in the normal phase. We also analyzed the low energy modes for real frequency and complex momentum, finding agreement with the hydrodynamic description. Finally, we extended the analysis to a setup where momentum was not conserved and succesfully matched the peculiar fourth sound mode with its hydrodynamic counterpart. Interestingly, we find a nonzero normal density in the setup with broken translations, similarly to the results of [73,75] for hyperscaling-Lifschitz type geometries.
After demonstrating the validity of the hydrodynamic description of the holographic superfluid, we carried out a study of the support of the various hydrodynamic modes in terms of the dual field theory operators. In a nutshell, we have quantified "how much" a single boundary field theory operator contributes to a specific hydrodynamic mode. The analysis has been performed both in the superfluid and normal phases and also in presence of momentum dissipation. In contrast to the field theory results of [83], we do not observe a role-reversal phenomenon between first and second sound in terms of the support of those modes. From the amplitudes and phases, we can reconstruct the dispersion relations of the hydrodynamic modes. This allows us to compute the speed of sound and the attenuation/diffusion constants from our data of the amplitudes and relative phases at a fixed value of the momentum. We notice that the absolute values of the amplitudes give us insights about the speed of the sound modes, and the relative phases between the fluctuations are related to damping. Furthermore, the amplitudes and relative phases give us an (computational) easy recipe to compute the hydrodynamic dispersion relations since it is sufficient to determine them at only one point of the dispersion relation (within the radius of convergence). Beyond hydrodynamics, we observe a purely imaginary gapped mode in the superfluid phase which is solely supported by the scalar fluctuations; this could be the Higgs mode discussed in [88,90,91].
In the future, it would be interesting to enlarge the analysis in this work to scenarios like those in the following list • Consider more complex symmetry breaking patterns such as the pseudo-spontaneous breaking of the U (1) symmetry (as initiated in [77]). Even though the explicit breaking of a U (1) vector symmetry is definitely not a physical option, this could serve as a toy-model to improve our understanding of the pseudospontaneous breaking of translations in the homogeneous models [92][93][94][95][96][97][98][99] and the possible relation between phase relaxation and pseudo-spontaneous breaking [100][101][102][103][104][105][106]. Finally, if these two breaking mechanisms are analogous or even equivalent, one could try to tackle the phenomenology of superfluid current relaxation induced by vortices [107] within this effective construction.
• Compute the second order contributions to equilibrium transport, as initiated in [120,121], in superfluids and systems with second order phase transitions in general.
• Work out the hydrodynamics of holographic superfluids in presence of a background superfluid velocity [122][123][124]. The excitations in the probe limit have been already presented in [16] and the dynamical Landau instability in [125]. It would be interesting to extend them to the fully backreacted scenario and compare them to the hydrodynamic description.
• Holographic superfluid models in presence of disorder have been constructed in the past years [126,127]. It would be interesting to study the transport of those models and compare it with our results in the simpler axion model.
We leave these questions for the near future.
where γ mn denotes the induced metric at a z = hypersurface, while K mn is the extrinsic curvature at that hypersurface, K is its trace, and R (γ) mn the Ricci tensor. 7 In order to read the one-point functions from (A.2) we need the boundary expansions of the bulk fields. These take the generic form For the background solution (indicated with a hat symbol), we havê (A.5) Hence, we are taking the leading contribution of the scalar towards the boundary to be real. As usual this will correspond to picking a real value for the symmetrybreaking condensate.
Equations (A.2-A.5) result in the following expressions for the expectation values of the dual field theory operators up to first order in fluctuations i (x, t) , (i = x, y) , (A.6d) tt (x, t) , (A.6e) ty (x, t) , T xy = −3h (3) xy (x, t) , (A.6h) where the energy density is given by ε = 2 u 3 . The subindex (b) stands for the finite value of the different fields towards the boundary, hence once stripped of the powers of z (e.g. δγ mn(b) = δγ mn /z 2 ).
It is worth pointing out that the boundary expansions (A.4) fulfill the necessary relations for the relevant Ward Identities to be satisfied. First, they obey the constraint yy (x, t) = 0 , (A.7) 7 We write the induced metric and extrinsic curvature in the four-indices notation, see e.g. [128].
guaranteeing that T µ µ = 0, as required by conformal invariance. As for the Ward Identity ∇ µ T µν = 0, that results from translational invariance, it boils down to the following constraints: where we have assumed an harmonic spacetime dependence of the fluctuations as in (2.20). Finally, the Ward Identity ∇ µ J µ = 0, stemming from the global U (1) symmetry, is satisfied once k a (1) x + ω a (1) holds. By solving the equations (2.13) towards the boundary we have checked that indeed the constraints (A.8, A.9) are satisfied. From the Ward identities and the tracelessness of the energy-momentum tensor, we find the following relations between the expectation values of the fluctuations where the δ in the expectation values indicates that we are referring to the fluctuations.
Regarding the axion fields introducing momentum dissipation, the boundary expansion for the axion fluctuations reads we have (see [34] for the derivation) Hence, in addition to the expectation values (A.6), we have the expectation value of the axion fluctuation [34] φ I = 6 φ I (3) . (A. 13) for the operator dual to φ I . Since the axions break translational invariance, the Ward identities discussed above do no longer hold at the level of the fluctuations. For a discussion see [129]. Finally, by computing the one-point functions at sources on from (A.2), we can determine the form of the correlators relevant for the Kubo-formulas employed in the main text:

B Dispersion Relations from Amplitudes
From the amplitudes computed at a certain value of the momentum, we might extract the diffusion constants and the hydrodynamic dispersion relations. The dispersion relation of the purely diffusive mode in the normal phase is for example given by eq. (2.5) Using the Ward identity eq. (A.11), we find or in other words, the diffusion constant is given by which is nothing else than the very well-known Fick's law. The imaginary unit in eq. (B.3) simply implies that the amplitudes have a phase difference of π/2 where we required the diffusion constant to be real in the last step. Requiring the diffusion constant to be real and positive and thus the cosine to vanish fixes the relative amplitudes to π/2 mod 2π. From this equation it is obvious, that damping is intimately related to a phase difference between the fluctuations. We will see the same for sound waves. Note that the momentum dependence of the expectation values is more general, i.e. they also contain all higher orders in k.
For a sound mode with speed v s , attenuation constant Γ s and real momentum k, we notice with eq. (A.10) that In the case, where we know the speed of sound (first sound), we can simply take the absolute value of this equation and solve for the attenuation constant The easiest way to compute the speed is to use the information about the relative phases and decompose the ratio of the expectation values as and we thus find for the speed of sound according to eq. (B.5) and Euler's formula Note that in the case of first sound we observe a ratio of approximately 0.7 between | δT tx | and | δT tt | (see for example figure 5 or 8) which yields a conformal speed of 1/ √ 2. The attenuation constant is then simply given by Instead of the Ward identity between δT tt and δT tx , we can of course derive similar relations with δT xx and δT yy using the tracelessness of the energy-momentum tensor. This is in particular important in the case of 4th sound, where the Ward identity associated with energy momentum conservation no longer holds. Independently, we may also use eq. (A.11) to express the speed and attenuation constant in terms of the expectation values.

C Eigenvectors in Hydrodynamics
In this appendix, we discuss two simple examples within the hydrodynamic framework. They serve the purpose of putting our results for the amplitudes in section 4 into context and connecting them to the existing literature.

C.1 The purely diffusive mode in normal fluids
Let us consider the longitudinal sector of a conformal charge fluid with conserved momentum. Following the same notations as [84], we may write the Fourier transformed hydrodynamic equations of motion in the form Note that δT tt and δT tx are not independent variables but related by the Ward identities (A.10). Instead of an independent hydrodynamic variable, δT tx should be viewed as an auxiliary variable to reduce the generalized eigenvalue problem to first order in ω as explained below eq. (D.1).
Solving the hydrodynamic equations as a generalized eigenvalue problem with respect to the frequency ω, we find a pair of sound modes and a purely diffusive mode in agreement with the results in the main text, i.e.
together with their corresponding eigenvectors Importantly, this means that the eigenvector corresponding to the purely diffusive mode, v 3 , consists only of δJ t to first order in k, exactly as we observe in the left plot of fig. 5. Furthermore, we see that the ratio of the contributions to first sound from δT tt and δJ t is given by µρ+sT ρ , exactly as we observe in the middle plot of fig. 5.
Let us now explain why our results are not in contradiction with those in [84] and how they relate to them. In [84] it was stated that the purely diffusive mode is carried by the incoherent current δQ diff ≡ δJ t − ρ sT +µρ δT tt . To get to this result, we have to perform a basis transformation and consider Notice that, in this basis, the diffusive mode corresponds toṽ 3 = δJ t − ρ sT +µρ δT tt . In the language of [84],ṽ 3 is the incoherent current. Note that the relative coefficient iñ v 3 is exactly the inverse of the relative coefficient between δT tt and δJ t in the sound mode. For simplicity, we define the coefficients c I appearing in the decomposition of the various eigenvectors in terms of the basis of operator fluctuations using: where n labels the eigenvectors and I the basis of operators. The coefficients c I are exactly the quantities discussed in [84] and displayed in figure 13, 14 and 15. The remaining vectors in (C.4) arẽ δT tx + . . .

(C.6)
and they correspond to the first sound mode. Using the Ward Identity relating δT tt and δT tx and assuming that the modes corresponding toṽ 1,2 (sound modes) are not excited, i.e.ṽ 1 ,ṽ 2 = 0, one immediately obtains the condition δT tt = 0. This implies that only charge can fluctuate and leads toṽ 3 = δJ t , which matches exactly our result presented in the main text shown in fig. 5. In a nutshell, the main physical difference between the approach of [84] and ours lies in considering the response of the system to different perturbations. While in section 4 we studied the support of the different hydrodynamic modes when only one of those modes is switched on, ref. [84] considered the question for a generic perturbation in which all modes are excited. Obviously, as shown in eq. (C.4), the information of the two approaches can be extracted from the same original hydrodynamic matrix (e.g. eq. (C.1)). As a proof of that, we can invert the matrix of the amplitudes of the eigenvectors, project it on the basis of the fluctuations and obtain exactly the ρ sT +µρ factor derived in [84]. The results of this procedure are shown in fig. 13.
(C.7) We have defined: σ the electric conductivity, χ ≡ ∂ρ/∂µ the charge susceptibility, κ the thermal conductivity,ᾱ the thermoelectric conductivity, ζ = ∂s/∂µ and the The relative amplitudes of δT tt and δJ t for the two thermoelectric diffusion modes at T /µ = 5 (in blue, corresponding to the data in fig. 7) and the hydrodynamic prediction for the ratio (dashed red line). Bottom: The normalized, dimensionless coefficients c δJ t T 2 and c δT tt T 3 in the basis of the eigenvectors from eq. (C.5). In the incoherent limit α µ, δJ + → δJ t and δJ − ∼ δT tt + γ − αδJ t . This implies that δJ − is purely δJ t and the other mode is still a mixture of δT tt and δJ t . The relative coefficient of the bottom right plot is presented in fig. 15. specific heat c µ = T ∂s/∂T . All these coefficients are known analytically for the linear axion model, see for example [87]. At this point, we can perform the same analysis as before and compute the eigenvectors J ± for the matrix eq. (C.7). The matrix is diagonal in the basis of the eigenvectors δJ ± = a ± (δT tt + γ ± α δJ t ), where a ± is just an overal normalization and γ ± = − 3 4α ρ 1 ± 1 + 16α 2 ρ 2 9 2 [86], with being the energy density and α the dimensionful strength of momentum dissipation. We present the amplitudes and coefficients c I in fig. 14. In fig. 15, we explicitly verify that the relative coefficient matches the analytic expression of [86].

D Numerical methods
In this section, we briefly describe the numerical methods we employed to solve the differential equations in this work following [130][131][132] (for a detailed introduction see [133]). Both, the background equilibrium solutions as well as the solutions to the linearized perturbations about the equilibrium state are obtained by means of pseudospectral methods. The idea of spectral methods is to approximate the numerical solution in terms of basis functions on a discretized grid. Throughout this work, we choose Chebychev polynomials as basis functions and a Chebychev-Lobatto grid to discretize the radial direction. Spectral methods solve the equations of motions globally which is a big advantage compared to shooting methods or finite-differences where we have to vary the initial conditions on one side of the domain until we find the desired boundary values at the other end of the interval. More importantly, spectral method are highly accurate and have a fast convergence rate and are thus perfectly suited for problems in numerical holography. Computing the static background and two-point functions in terms of pseudospectral methods is fairly standard and we thus only comment on how to obtain the quasi-normal modes. After obtaining the static background solution, we want to investigate time-and space dependent (linearized) fluctuations about that background to compute (a) the quasi-normal modes and (b) the retarded Green's functions for the transport coefficients. The coupled ordinary differential equations for the linearized fluctuations are generally of the form (C(k) ω 2 + A(k) ω − B(k)) x(ω, k) = 0, (F (ω) k 2 + D(ω) k − E(ω)) x(ω, k) = 0 (D.1) where A to F are differential operators and the vector x consists of all the fluctuations x = {h tt , h tx , h xx , h yy , a t , a x , δψ 1 , δψ 2 , δφ x }. We may recast both problems in terms of a generalized eigenvalue problem with respect to the quasi-normal frequency ω or the momentum k by introducing auxiliary functions of the formf = ω f andf = k f so that we find (Ã ω −B)x 1 = 0, (D k −Ẽ)x 2 = 0 (D. The tilde inx 1 andx 2 denote the auxiliary functions. By solving the generalized eigenvalue problem, we find for each eigenvalue ω n (or k n ) an eigenvectorx. Finally, to resolve the strong gradients at the horizon at low temperatures accurately we use up to 350 gridpoints in order to satisfy the equations of motion and constraints with accuracy better than 10 −12 . The Chebychev coefficients of the slowest convergent solution drop below 10 −13 in this case.