Longitudinal Sound and Diffusion in Holographic Massive Gravity

We consider a simple class of holographic massive gravity models for which the dual field theories break translational invariance spontaneously. We study, in detail, the longitudinal sector of the quasi-normal modes at zero charge density. We identify three hydrodynamic modes in this sector: a pair of sound modes and one diffusion mode. We numerically compute the dispersion relations of the hydrodynamic modes. The obtained speed and the attenuation of the sound modes are in agreement with the hydrodynamic predictions. On the contrary, we surprisingly find disagreement in the case of the diffusive mode; its diffusion constant extracted from the quasi-normal mode data does not agree with the expectations from hydrodynamics. We confirm our numerical results using analytic tools in the decoupling limit and we comment on some possible reasons behind the disagreement. Finally, we extend the analysis of the collective longitudinal modes beyond the hydrodynamic limit by displaying the dynamics of the higher quasi-normal modes at large frequencies and momenta.


Heraclitus
Hydrodynamics and effective field theories (EFTs) are powerful tools which can be used to describe the macroscopic and low energy dynamics of physical systems. EFTs have been widely and successfully applied in a broad variety of research areas, from cosmology to particle physics, and finally to condensed matter. Symmetries are the fundamental pillars upon which such effective frameworks are built; furthermore, the symmetries of a system strongly constrain the nature of its collective excitations and their low energy dynamics. However, in realistic settings, and especially in the context of condensed matter systems, part of the fundamental symmetries appear to be broken. Paradigmatic examples are the spontaneous breaking of the U (1) symmetry in superconducting materials, as well as the spontaneous breaking of translational (and rotational) invariance in crystals with long range order, i.e. systems modelled by periodic lattices. Nevertheless, even in cases where some of the symmetries are broken, effective field theory frameworks still provide an important guidance for theoretical descriptions [1]. In particular, the various symmetry breaking patterns can be used to classify the corresponding phases of matter [2]. Historically, such methods of classification have been of great importance for the understanding of the dynamics of fluids, liquid crystals, nematic crystals and ordered crystals [3,4].
Nevertheless, for systems of the types mentioned above, it is challenging to set up the hydrodynamic theory, i.e. finding a consistent and complete gradient expansion. Over the last decade, holographic dualities (also known as AdS/CFT or gauge/gravity dualities) have provided a powerful and effective tool set which can be used to tackle strongly coupled condensed matter systems [5][6][7]. In particular, holography has provided novel insights into hydrodynamics concerning new transport coefficients; bounds of transport coefficients; and the applicability of hydrodynamics in the presence of large gradients. Moreover, holography goes beyond hydrodynamics in the sense that it allows for direct computations of transport coefficients, which in the EFT description appear only as undetermined parameters.
Historically, holography has been focused on the description of strongly coupled fluids [8,9]. More recently, the attention switched to the description of strongly coupled solids or, in general a more general sense, systems with elastic properties. The longterm scope of this programme is to bring the theoretical framework closer to realistic systems in order to make concrete predictions which are testable in the lab.
The propagation of sound in compressible materials is a direct macroscopic manifestation of the spontaneous symmetry breaking (SSB) of translational invariance. Sound propagation admits the effective quasiparticle description by acoustic phonons, which are the Goldstone bosons corresponding to the spontaneous breaking of translational symmetry [10]. Sound can propagate in two different forms, transverse and longitudinal sound, which correspond to transverse and longitudinal phononic vibrational modes with dispersion relation ω = ±c T,L k − i D T,L k 2 , (1.1) where c T,L denotes the sound speed and D T,L denotes the diffusion constant. Such modes are tightly correlated to the mechanical properties of the material, in particular to its shear and bulk elastic moduli G, K, respectively [11,12]. The speed of transverse sound is determined by the simple expression where G is the shear elastic modulus and χ ππ the momentum susceptibility. The property of supporting propagating shear elastic waves is usually considered the main distinction between solids and fluids. 1 Longitudinal sound appears both in solids and in fluids and, in two spatial directions, propagates with speed where K is the bulk elastic modulus, which is the coefficient determining the response of the system under a compressive strain. The diffusion constants appearing in expression (1.1) are determined by the dissipative mechanisms of the system; in a conformal system, where the bulk viscosity vanishes, they are both given in terms of the shear viscosity η. In particular, for a relativistic CFT with two spatial directions, we have D T = η χ ππ , D L = 1 2 η χ ππ , with χ ππ = ε + p = s T, (1.4) which can be derived using both relativistic hydrodynamics [19] and holographic techniques [8]. The presence of a diffusive term in the dispersion relation of the phonons (1.1) indicates that the system has also a viscous component, i.e. that it is a viscoelastic material. In an unrealistic perfect crystal, with no dissipation nor effective viscosity, sound would propagate indefinitely, which is not realistic.
In a system with spontaneously broken translational symmetry the longitudinal sector contains an additional diffusive mode with dispersion relation with diffusion constant D Φ . This diffusive mode is intimately connected with the presence of Goldstone degrees of freedom, in fact its hydrodynamic nature is simply 1 This statement is imprecise and true only at low frequency and low momentum ω/T, k/T 1. It is well known that fluids at large frequency, i.e. beyond the so-called Frenkel frequency, support propagating shear waves [13]. Simultaneously, it is suggested by recent experiments [14,15] and theoretical developments [16][17][18] that the same might happen at large momentum, beyond the usually called k-gap. related to the conservation of the phase of the Goldstone bosons, i.e. the phonons. The presence of such a mode can be predicted by performing an effective hydrodynamic description for crystals [3,4,20]. Notably, the underlying physics is very similar to that of superfluids/superconductors [21], where a Goldstone diffusion mode appears as well.
All the features just mentioned can be formally derived using hydrodynamic techniques [3,4,20]. The analysis of the transverse modes within the holographic framework and its successful matching with the hydrodynamic predictions have already been presented in [22]. The scope of this paper is to extend the analysis to the longitudinal sector, in particular to the study of the diffusion constants of the sound and diffusive modes, in a simple holographic system with spontaneously broken translational invariance. Concretely, we use the holographic massive gravity theory introduced in [23,24] to implement the SSB breaking of translations as explained in [22].
The hydrodynamic degrees of freedom present in the longitudinal sector have previously been investigated in a more complicated model in [25], with more emphasis on the electric transport properties of the dual field theory. In particular the longitudinal sound mode and the diffusive mode have been identified numerically but without any systematic study or theoretical background. Related results have recently appeared in [26] in the presence of electromagnetic interactions and in [27] for holographic fluid models. Other recent works have analysed the origin of this extra diffusive mode in similar holographic models with a global U (1) symmetry [28,29]. In this work, thanks to the simplicity of the model at hand, we improve the analysis. In particular: 1. We study in the dynamics of the longitudinal collective modes dialing the SSB scale.
2. We compare numerical results obtained from the holographic model with the hydrodynamics expectations.
4. We numerically describe the higher modes, extending the analysis beyond the hydrodynamic limit towards large frequencies and momenta, ω/T, k/T 1.
From the numerical study we obtain the expected hydrodynamic modes, namely the longitudinal sound and the diffusive mode. Both modes display the correct dispersion relation. Moreover, the propagation speed and the attenuation constant of the sound mode are in agreement with the hydrodynamic predictions [3,4,20] at every value of the SSB scale. Instead, we find a disagreement between the diffusion constant of the diffusive mode and the value predicted by hydrodynamics. We will comment on this puzzle below. Our results indicate that a complete understanding of these holographic models, in terms of an hydrodynamic effective description, may still be lacking.

Holographic Set-up
We consider generic viscoelastic holographic massive gravity models [23,24] defined by the following four-dimensional gravitational action in the bulk with X ≡ 1 2 g µν ∂ µ φ I ∂ ν φ I and m is a coupling which relates to the graviton mass [24]. Here, φ I are Stückelberg scalars which admit a radially constant profile, φ I = x I with I = 1, 2. In the dual field theory the Stückelberg fields represent scalar operators breaking the translational invariance of the system. Whether the translational symmetry breaking is explicit or spontaneous depends crucially on the bulk potential of the scalar fields. In particular, depending on the choice of the potential V (X), the bulk solution φ I = x I may either correspond to the non-normalisable mode and hence it acts as a source term for the dual operators Φ I ; or it corresponds to the normalisable mode and hence it gives rise to a non-zero expectation value Φ I . In the following we will concentrate on potentials of the type which in a simple way realise the spontaneous breaking of translational invariance [22]. The thermal state is dual to an asymptotically AdS black hole geometry which in Eddington-Finkelstein (EF) coordinates is described by the metric where u ∈ [0, u h ] is the radial holographic coordinate. The conformal boundary is located at u = 0 and the horizon u h is defined by f (u h ) = 0. Hence the emblackening factor f (u) takes the simple form The corresponding temperature of the dual QFT reads while the entropy density is s = 2π/u 2 h . For simplicity, and without loss of generality, we will fix u h = 1 in the rest of the manuscript.
The heat capacity is c V = ∂ε/∂T , with ε being the energy density and T being the temperature, and it has been studied in [17,30]. The dual field theory has been shown to possess viscoelastic features and was studied in detail in [22,[30][31][32][33][34][35]. 2 It is important to stress that our model only has one dimensionless scale controlling the strength of the spontaneous translational symmetry breaking. More specifically, we have SSB ∼ m/T , (2.6) where the limiting case m = 0 corresponds to the translationally invariant system. In the following we consider fluctuations around the thermal equilibrium state and compute the quasi-normal modes. The transverse sector of the fluctuations has been studied in detail in [22] and it revealed the presence of damped transverse phonons whose speed are in perfect agreement with the theoretical predictions from hydrodynamics. Within this manuscript, we continue this analysis by considering instead the quasi-normal modes within the longitudinal sector.

Hydrodynamics and Longitudinal Quasi-Normal Modes
In this section we will study the hydrodynamic modes appearing in the longitudinal sector of the quasi-normal mode spectrum. 3 Furthermore, we compare the dispersion relations of the quasi-normal modes to the corresponding, analytic, hydrodynamic expectations.
We defer the relevant equations of motion for the fluctuations and the technicality concerning the holographic computation to appendices A and B. Moreover, in appendix C we review and (partially) extend the hydrodynamic effective description of systems with broken translations following the lines of [20]. 4

Hydrodynamic Regime of Quasi-Normal Modes
First we analyse the quasi-normal modes in the low-frequency and long-wavelength regime. In particular the spectrum exhibits:

A pair of longitudinal damped sound modes with dispersion relation
with the speed c L and the diffusion constant D p . 5 The ellipsis stands for higher momenta corrections. The mode (3.1) is exactly a longitudinal damped phonon mode, which is expected both in fluids and solids.

A diffusive mode with dispersion relation
where D Φ is the diffusion constant. The ellipsis stands for higher momenta corrections. The presence of such a diffusive mode is typical of systems which break translational invariance spontaneously. This hydrodynamic mode can be viewed as the diffusive mode of the Goldstone parallel to the momentum.
The two modes presented above are shown in figure 2, for the specific potential N = 3 and for different values of the dimensionless parameter m/T . We obtain similar results for other values of N with N > 5/2 (see figure 1 for a schematic representation of the hydrodynamic modes of our system).

Correlators of the Goldstone Modes
Before comparing the numerically obtained dispersion relations of the low-lying quasinormal modes to the expected dispersion relations of hydrodynamics, let us discuss the retarded Green's functions which will be relevant for the future analysis. We are particularly interested in correlators appearing due to the SSB of translations, i.e. those containing the Goldstone fields Φ I which are dual to the Stückelberg fields φ I in the bulk.
To be more precise, we discuss the retarded Green's functions G π I Φ J (ω, k) and G Φ I Φ J (ω, k), where Φ I and π I (with I = 1, 2) are the two Goldstone modes and the momentum operators, respectively. Using the hydrodynamic description of phases with spontaneously broken translational symmetry, found in [3,20] and repeated in our Appendix C, we can derive the form of the previous Green's functions at low energy/frequency where χ ππ is the momentum susceptibility and ξ relates to the Goldstone diffusion.
The holographic correlators can be derived from the equations for the fluctuations of the bulk fields using the standard holographic dictionary [37]; see Appendix B for further details, in particular regarding the mapping between the scalar Stückelberg fields φ I and the Goldstone operators Φ I .
We have conducted numerical studies of the behaviour of the Green's functions (3.3) within the context of our holographic model and for various potentials V (X) = X N with N > 3. One example of the results is shown in figure 3. We find perfect agreement between the Green's functions obtained from holography and from hydrodynamics. The transport coefficient ξ can be read off from the imaginary part of the Green's function where, in this limit, it is no longer necessary to to distinguish between Φ (1) or Φ (2) . 6 A plot of the ξ parameter for various N is shown in figure 4. Furthermore, the numerical results for ξ are in perfect agreement with the horizon formula given in [38] and derived in [39], which here is used at zero chemical potential µ = 0. 7 6 At zero momentum, k = 0, and for isotropic backgrounds, there is no distinction between the various directions. 7 Notice the different notations. The map between the conventions is

Matching to Hydrodynamics
Let us compare the dispersion relations obtained from the quasi-normal mode analysis presented above with the theoretical expectations from hydrodynamics. For details regarding the hydrodynamic analysis see Appendix C.

Sound Mode
First we analyse the dispersion relation of the propagating longitudinal mode (3.1), starting with its speed c L . The hydrodynamic analysis performed in Appendix C predicts the speed c L to be given by where ε is the energy density; p is the thermodynamic pressure; and χ ππ is the momentum susceptibility. Moreover, κ and G are the bulk and shear elastic modulus, respectively. Both κ and G vanish in the absence of translational symmetry breaking and the formula (3.7) reduces to the well-known speed formula for sound Reminiscent of solids, we can re-express the speed of the longitudinal sound modes as where K is the thermodynamic bulk modulus obtained from the inverse of the compressibility where V is the volume of the system and P is the mechanical pressure defined by P ≡ T ii (no sum over i implied). Let us notice that, in our model which means that the thermodynamic and mechanical pressures coincide only at zero SSB, m/T = 0. In order to establish an equivalence between (3.7) and (3.8) we must demand In other words, the first term in the above formula is a contribution which is finite even in the absence of SSB, and it is indeed what determines the speed of sound in a pure relativistic CFT. The second term captures additional effects due to spontaneous symmetry breaking.
Let us now determine the coefficients G, K and κ within our model. The elastic shear modulus is given in terms of the Kubo formula 8 Re G TxyTxy (ω, k) , (3.12) and has been discussed in [22,24,33,34,40] within the context of the holographic model considered here. The bulk modulus κ can be determined in various ways. Using the underlying conformal invariance of the dual QFT, as well as the structure of the energy-momentum tensor, the bulk modulus κ reads [40] In fact, using χ ππ = 3ε/2 as well as 9 which is valid for two-dimensional conformal solids, we can check that the consistency relation (3.11) is indeed satisfied. Finally, we may also compute κ by a Kubo formula; for more details see appendix C.4. Collecting all the above formulae we end up re-discovering a relation between the transverse and longitudinal speeds of sound put forward in [42,43], i.e. that we have Depending on the choice of the potential, the speed of sound might become superluminal and the system can exhibit causality pathologies (see [41]). However, this is not the case for the choice V (X) = X N , provided N ≥ 3. In the right panel we compare the numerical data extracted from the QNMs and the theoretical formula written down in equation (3.7). Within the hydrodynamic limit k/T 1 the agreement is excellent.
Let us now move on and discuss the dissipative part of the phonon dispersion relation. Because of finite temperature effects, or in other words due to the presence of an event horizon and its associated shear viscosity η, the longitudinal phonons obtain a finite diffusive damping D p . The same phenomenon occurs also in the transverse sector, the effects of which were studied in [22]. The interplay of elasticity (a propagating term) The comparison between the numerical data (solid line) and the theoretical expectations for m/T 1.
and viscosity (damping) qualify the system at hand to represent the gravity dual of a viscoelastic material. Hydrodynamics provides us with the formula for the diffusion constant, which reads (3.16) where we define the specific heat c V ≡ ∂ε/∂T and the parameters η, ξ,κ 0 and γ 2 as where Φ is the Goldstone operator (parallel to k) of the dual field theory (not to be confused with the scalar bulk field φ).
The equation for the diffusion constant can be considered as a sum of leading and sub-leading terms. The leading contribution, at m/T 1, is where the dots signify the sub-leading terms. The above expression reduces to the known result for sound damping, D p = η/(2sT ), in the absence of SSB [9,44]. In figure 6 we compare the behaviour of the dimensionless damping constant D p T with the numerical results computed for the potential with N = 3, as functions of the symmetry breaking scale m/T . As expected, the leading order result (3.21) represents a good approximation only for small values of m/T . Note that at large m/T → ∞ the damping coefficient goes to zero and the phonons become purely propagating modes. In other words, at zero temperature no dissipation is present. The results are analogous for various N > 3, see figure 7. The tendency is a smooth decrease when going to small temperature, which is consistent with the hypothesis that no dissipation can be at work at T = 0; thus, in this limit our system can be considered a "perfect elastic solid." We verified numerically, using the Kubo formulas (3.19) and (3.20), that both coefficientsκ 0 and γ 2 are zero in our holographic model, which is in agreement with [45]. 10 Thus, the full expression for the diffusion constant of the sound modes reduces to In figure 8 we plot this improved formula next to the numerical data, for N = 5. It is evident that the complete formula, including the sub-leading terms in equation (3.22), is now in very good agreement with the numerical data.

Diffusive Mode
As mentioned above, spontaneous breaking of translational invariance gives rise to a non-propagating diffusive mode in the longitudinal sector. 11 For an effective hydrodynamic description of this mode we refer the reader to [3,4,20] or to our appendix C. This additional diffusive mode was first observed numerically in the context of holography [25] and later discussed in [26]. In this subsection we will study this mode, in particular its diffusive constant D Φ . The numerical results for the dimensionless form of the diffusion constant are shown in figure 9, for a range of values of the power N . We notice that the value of D Φ is finite even at m/T = 0, suggesting that its nature can be understood entirely in the limit when the Goldstone modes decouple. 12 Moreover, D Φ decreases monotonically when increasing the SSB parameter m/T .
In Appendix C we derive the following analytic expression for the diffusion constant In the case ofκ 0 = γ 2 = 0, equation (3.23) displays a leading contribution

24)
11 Notice that such a mode does not appear when translations are broken explicitly as in [46]. 12 Similar observations appeared in [38]. where the dots represent corrections. Surprisingly, we find that the complete hydrodynamic formula in expression (3.23) does not agree with the numerical quasi-normal mode analysis, which is evident in the plot in the right panel of figure 10. Such a discrepancy between the numerical data and the analytic expression for D Φ points towards a disagreement between the hydrodynamic predictions and the holographic results. Before continuing the analysis of this issue, let us first comment on the various parameters entering in the formula (3.23): We are able to gain some insight in to the puzzle by doing an analysis in the bulk. As stated in several works [26-29, 38, 39], the diffusive Goldstone mode originates in the scalar sector and couples to the gravitational sector via the graviton mass m. In the limit of very small graviton mass, i.e. at values of m/T 1, one can approximate that the scalar and gravitational sectors decouple, hence the scalar fields can be treated as probes on a fixed gravitational background. In this decoupling limit the equation of motion for the Goldstone mode in the longitudnal sector reads with the Schwarzschild emblackening factor given by As expected, the quasi-normal modes calculated in this limit show the presence of a diffusive mode with dispersion relation ω = −iD φ k 2 + . . . , (3.27) which is shown in figure 11. Given the simplicity of the decoupled regime, we are able to obtain the Green's function of the scalar operator analytically, using perturbative techniques. Implementing standard methods, we find the following expression for the diffusion constantD 28) where N is the power in the potential V (X) which defines the model. Continuing, in the limit of small graviton mass the constituents of the leading contribution to the hydrodynamic formula for D Φ , given by equation (3.24), can be written (3.29) Using the above formulae we find the leading order of D Φ , at m/T 1, to be given by the formula Comparing (3.28) and (3.30), it is evident that the leading contribution to the hydrodynamic formula for the diffusion constant of the diffusive mode does not match the corresponding quantity extracted from holography, even at m/T 1. Moreover we can see thatD where we in the second equality use the identity (κ + G)/G = 3/2; and D ⊥ = Gξ is the diffusion constant of the transverse diffusive mode, which has been successfully tested in [35,38]. The above relation can also be seen directly by comparing the equation of motion (3.25) to its transverse counterpart The plot shown in the left panel of figure 10 indicates that the formula (3.31) is in good agreement with the numerical data in the limit of soft SSB, roughly m/T 0.5. It is important to note that a missing factor of N is not able to fix the disagreement between hydrodynamics and the holographic results at arbitrary values of m/T . This is clear from the left panel of figure 10, where the discrepancy is still present at values of m/T 1. Taking the difference between (3.30) and (3.24) we find which means we at leading order in contribution have the relatioñ In terms of dimensionless quantities the above relation becomes Surprisingly, we have found that this discrepancy is independent of the choice of potential, at least for small m/T . The predictions of hydrodynamics combined with such a shift reproduces the numerical values of the diffusion constant very well for potentials of the form V (X) = X N , even until large values of m/T , see the right panel of figure  10.
We are currently unable to find a precise hydrodynamic resolution for the mismatch described above. Assuming the correctness of our computations, we cannot discard the possibility that the hydrodynamic description is missing some effect encoded in a novel transport coefficient. At the same time, it is not guaranteed that the holographic models considered in this work are the exact duals of a system with spontaneously broken translations, as those described by hydrodynamics. One possible reason for the discrepancy is that the susceptibility matrix defined in eq.(C.6) is not as simple and in particular it might contain off diagonal terms. A preliminary analysis shows that such off-diagonal terms will not appear in our model and that they will not affect the results at m = 0. A more detailed investigation in this direction is needed.
A last possibility is related to the fact that our vacuum does not minimize the free energy. In other words, our solution is not thermodynamically favourable and hence it can been seen as an excited state of the dual field theory. The situation is nevertheless more subtle, because our theory can not be constructed around the wouldbe thermodynamic favourable vacuum φ I = 0. More precisely, the theory we consider is strongly coupled around such a background and the results are therefore not trustable. In principle this should not affect the quantities we are discussing and indeed for all of them beside the Goldstone diffusion constant we find very good agreement. In summary, our results are surprising given the large amount of verifications in this direction.

Quasi-Normal Modes Beyond Hydrodynamics
In the previous section we investigated hydrodynamic modes, i.e. poles of the retarded Green's function located at ω/T 1, k/T 1. Hence, these modes determine the late time and long distance behaviour of the system.
Using holographic techniques we can study further, so-called non-hydrodynamic, poles of the Green's functions. These poles correspond to higher quasi-normal modes which we determine numerically. In particular we are interested in the behaviour of the higher quasi-normal modes as functions of the SSB scale m/T , and momentum k/T .
The results for the potential V (X) = X 3 are shown in figure 12. The left panel of the figure displays the tower of quasi-normal modes at zero momentum while dialing the dimensionless parameter m/T . For any given ratio m/T there is a clear gap in the system. Moreover, as the mass m/T is increased all the higher non-hydrodynamic modes move away from the origin of the complex plane and the imaginary part of the The dispersion relation of the longitudinal sound mode beyond the hydrodynamic limit. In the hydrodynamic limit, k/T 1, the behaviour is described by (3.1). At large momenta, i.e. for k/T 1, the real part of the frequency asymptotes to the UV dispersion relation Re(ω) = ±k indicated by the dashed black line. The imaginary part of the frequency is small compared to the real part even for large momenta. frequency grows. Hence the late time and long distance behaviour is determined by the three hydrodynamic modes which we investigated in the previous section. The other poles do not give rise to quasi-particles due to their short lifetime.
In the right panel of figure 12 we show the quasi-normal mode spectrum as a function of the momentum k. In addition, the more refined results for the first seven modes are shown in figure 13. Note that the dynamics of the higher quasi-normal modes as a function of the momentum appears to be quite complicated. Especially noticeable is the large imaginary part of the QNM depicted in blue in figure 13 (corresponding to the diffusive mode in the hydrodynamic regime) upon increasing k/T . Finally, we analyse the large momentum behaviour of the quasi-normal modes corresponding to the two hydrodynamic sound modes. In particular, we consider the case k/T 1 which clearly exceeds the hydrodynamic regime. As discussed above, at small momenta k/T 1 the expected dispersion relation is given by (3.1). On the contrary at very large momenta we obtain the dispersion relation of a sound mode within the relativistic UV AdS fixed point. In this limit, the real part of the dispersion relation is given by Re(ω) = ± k. In contrast to the sound mode in the transverse sector (see [35]), the real part of the frequency of the sound mode never becomes zero; in fact, as shown in 13 the real part of the dispersion relation of the sound mode interpolates smoothly between Re(ω) = ±c L k in the hydrodynamic limit and Re(ω) = ± k in the large momentum limit.

Conclusions
In this paper we studied the longitudinal quasi-normal mode spectrum of a simple holographic system with spontaneously broken translational symmetry, at zero charge density. We successfully identify the three hydrodynamic modes: the pair of longitudinal sound modes and a diffusion mode. We compared the numerical results obtained from the holographic model with the predictions of hydrodynamics given in [20] and in our appendix C.
The numerical results for the propagation speed and the attenuation constant of the sound modes are in agreement with the hydrodynamic formulas at any strength of the SSB, at arbitrary graviton mass. Surprisingly, we find disagreement between the diffusion constant of the diffusion mode obtained numerically and the formula given by hydrodynamics. The discrepancy is corroborated with analytic methods in the decoupling limit and it is present even at zero graviton mass. At this stage we are not able to identify the reason behind this puzzle.
Some comments are in order. Mistakes present in our computations are of course a possibility. A second, more interesting, option is that the hydrodynamic framework considered is missing some effect. A lacking transport coefficient, for example, could potentially explain the discrepancy we observe. A last scenario is that the holographic model discussed in this work is not the gravity dual of the type of systems described by the hydrodynamic picture, i.e. a dissipative and conformal system with spontaneously broken translational invariance. Needless to say, finding the reason for this discrepancy is certainly the most important open question that our manuscript poses. It would be interesting to perform the same check in similar holographic models [25,38,40,47].
As a final remark, it would be beneficial to gain a better understanding regarding the propagation of sound and the elastic properties of quantum critical materials, both from EFT methods [43] and holographic methods [41]. There are several recent hints suggesting that phonons, viscoelasticity and glassy behaviours may play an important role in quantum critical systems [48] and High-Tc superconductors [49]. Moreover, electron-phonon coupling might play an important role in the high-Tc puzzle [50].
In conclusion, we have in this paper continued the study of simple holographic massive gravity models with spontaneously broken translational symmetry. We have tested several predictions obtained from the hydrodynamic theory of such systems. Finally, we provided food for thought in the form of a discrepancy between the hydrodynamic and holographic theories concerning the dispersion relation of the diffusive mode. In this case, we would be very happy to be proven wrong, and even more happy to find a novel and interesting physical reason behind the mismatch we observed. discussions; Oriol Pujolas and Victor Cancer Castillo for collaboration on closely related topics and for sharing with us unpublished results. Finally, we would like to sincerely thank Blaise Gouteraux for several discussions and suggestions.

A Equations of Motion
In this work we determine the spectrum of quasi-normal modes in the longitudinal sector. Choosing the momentum along the y-direction, the relevant time-dependent fluctuations are δφ 2 , h tt , h ty , h yu , h xx , h uu , h yy and h tu . We choose to work in radial gauge i.e. h µu = 0, with µ ∈ {t, u, x, y, z}. Since we work in in-falling Eddington-Finkelstein coordinates, the fluctuations already satisfy in-going horizon conditions. The Fourier transformed equations of motions to first order in the fluctuations read where we introduced the short-hand notations h x,s and h x,a for the symmetric combination h x,s = 1/2 (h xx + h yy ) and the anti-symmetric combination h x,a = 1/2 (h xx − h yy ), respectively. In terms of these combinations it is straightforward to solve equation (A.8) by integrating twice, h x,s = c 1 + c 2 h x,s u. In order to determine the integration constants, we compare the solution to the near the boundary expansion at the conformal boundary given by In the context of the QNM computations we do not allow for sources of fluctuations c 1 and c 2 have to vanish leading to the trivial solution h x,s = 0. 13 In this case, the equations of motion simplify to explained in [51,52]. Note, that the sources of the fluctuations have to vanish which we implement by a suitable redefinition of the fields. We checked that our solutions fulfill the equations of motions and all constraint equations. To check the convergence of the numerical solution, we monitor the change of the solution for finer discretizations. As depicted in the l.h.s. of figure 14, the change of the quasi-normal mode frequency decays exponentially with a growing number of grid points; the same is valid for the corresponding eigenfunctions. Another check for the numerical method are the Chebychev-coefficients of the solution, displayed in the r.h.s. of figure 14; the coefficients decay exponentially, indicating exponential accuracy of our numerical method.

Correlation Functions
In order to calculate the correlation functions, outlined in equation (3.3), we have to source the scalar field explicitly. Note, that for k = 0 it is sufficient to consider the δφ 2 − h ty sector since the h tt and h x,a decouple. The corresponding asymptotic expansions are of the form 14 where we set h ty,source = 0 and φ source = 1. Note, that the logarithmic terms are not present for N = 3 and N / ∈ Z/2. In order to ensure convergence in presence of logarithmic term (e.g. in the cases N = 4, 5, 6, ...), we use a suitable mapping for the radial coordinate u, for instance u → u 2 (for more details see [51,52]).

Green's Function
The retarded Green's function of operators O a and O b is given by The Fourier-transformed retarded Green's functionG OaO b (ω, k) reads 15 and it describes the linear response to δ O a (ω, k) for a given perturbation j b (ω, k) by 16 Using holographic techniques the retarded Green's functions are given by where h µν vev , h ρσ source , φ vev and φ source are the boundary values as given in (B.2). The prefactors of the Green's function G φφ may be explained as follows: the multiplicative factor 2∆ − d arising in holographic renormalisation reduces to 2N − 5 for the (2 + 1)dimensional QFT considered here. The factors N m 2 are due to the non-canonical normalisation of the kinetic term for the fluctuations of the scalar field. In particular, the prefactor of the kinetic term reads m 2 V (X) whereX denotes X evaluated on the background.
Next we determine the correct normalization of the Goldstone operator. Note that the Goldstone operator Φ I has to satisfy G TtyΦ (2) (ω, k = 0) = i/ω, while on the gravity side we numerically confirmed Hence we can establish the following map between the Goldstone operator Φ I , on the quantum field theory side, and the Stückelberg field φ I on the gravity side,

C The Hydrodynamic Theory
In this section we present the hydrodynamic calculations which lead to the dispersion relations used in this paper. We follow the reasoning of [20] and expand on their results, focusing on the case of zero density. First we present the effective hydrodynamic theory and then compute the dispersion relations in the transverse and longitudinal sector.

C.1 Set-up
We will consider the linearised hydrodynamics of small fluctuations around a 2 + 1 dimensional system at equilibrium. Without loss of generality we only allow excitations to depend on x and t, hence the momentum is in the x-direction. We also adopt the convention of denoting the x-component of a quantity by , and the y-component by ⊥, to indicate the orientation with respect to the momentum. The hydrodynamic variables which are of interest for our analysis are for, in order, energy density; momentum density parallel and transverse to the momentum; and λ (t, x), λ ⊥ (t, x) are the divergence and curl of the two-component Goldstone field Φ(t, x), respectively (not to be confused with the scalar bulk fields φ I ).
To be more precise, for the system at hand, the field Φ(t, x) is a vector of the form From the above Goldstone field we define x) in the case of yindependent fluctuations. The sources of the theory are which denote the temperature; the parallel and transverse components of v(t, x); and linear displacements parallel and transverse to the momentum, respectively. The ordering of the above sources indicates the corresponding variable when compared to the list (C.1).
The hydrodynamic variables, collectively denoted ϕ a , are related to the their corresponding sources, s b , through the thermodynamic susceptibilities χ ab via the following equation In the linear hydrodynamic regime the susceptibilities will be between small fluctuations of the variables and sources. Hence, for the system under consideration, the above relations can be expressed in matrix form as follows where we have assumed that the susceptibility between momentum and velocity is identical in the parallel and transverse sectors. In the above equation it is manifest that the translation from sources to variables is facilitated by the inverse of the expressed susceptibility matrix. Conservation of energy and momentum provides us with the equationṡ where a dot above a variable indicates time-derivative. The constitutive relation for τ i0 (t, x) reads, to first order in derivatives, whereκ 0 and γ 2 are transport coefficients. 17 Furthermore, the spatial component of the stress-tensor, τ ij (t, x), to first order of derivatives, is expressed by where p(t, x) is the thermodynamic pressure; κ and G are the bulk and shear elastic modulus, respectively; and is the dissipative term, with η being the shear viscosity. In addition to the conservation equations there are also two Josephson relations. For λ ⊥ (t, x) the Josephson relation is given bẏ and for λ (t, x) the relation iṡ where ξ ⊥ and ξ are transport coefficients. 18

C.2 Dispersion Relations of Transverse Sector
Although the main analysis of this paper is concerned with the dispersion relations of the longitudinal sector of the hydrodynamics, it is nevertheless beneficial to establish the techniques used to compute the quantities under investigation by first considering the simpler case of the transverse sector.
In the transverse sector the hydrodynamic variables are π ⊥ (t, x) and λ ⊥ (t, x), with corresponding sources v ⊥ (t, x) and s ⊥ (t, x), respectively. To linear order in fluctuations the expansions are for the variables, and for the sources. The barred quantities represent the constant equilibrium values upon which we are perturbing. In the present sector the susceptibility relations (C.6) reduce to the two independent equations where G is the shear elastic modulus. The only conservation equation of relevance for the transverse sector is that of momentum conservation (C.8), which under the circumstances for the linearised system reads δπ ⊥ (t, x) + ∂ x τ xy (t, x) = 0, (C. 20) where, evidently, the only spatial derivative that contributes is with respect to x. Using the constitutive relation of the stress-tensor (C.10) we get where we have used ∂ x Φ (2) (t, x) = λ ⊥ (t, x); and the dissipative term σ xy (t, x) is given by equation (C.11), which yields At linear order the derivative of the stress-tensor is equal to Using the inverse of (C.18) the velocity fluctuation can be related to the momentum density, hence we arrive at The momentum conservation equation, up to linear order in fluctuations of the hydrodynamic variables, thus reads 25) In addition to the conservation equation we need to consider the Josephson relation (C.12). In the same way as for the stress-tensor, linearising the equation and expressing it in terms of the thermodynamic variables, yields The first step in the pursuit of the dispersion relations is to evaluate the derivatives in the equation of momentum conservation and the Josephson relation. This requires that we Fourier transform and identify coefficients. Doing so, we arrive at the system of equations where ω is the frequency and k is the momentum in the x-direction. The above system can be written in matrix form as which has non-trivial solutions only if the matrix on the left-hand side is non-invertible, i.e. if its determinant is zero. Setting this constraint means Solving the above equation for ω, taking the limit of small momentum, and comparing to the known form of a propagating mode with diffusion we identify the speed of sound of the transverse mode and diffusion constant This concludes the calculation of the dispersion relations in the transverse sector.

C.3 Dispersion Relations of Longitudinal Sector
The procedure of finding the dispersion relations in the longitudinal sector is in principle the same as for the transverse sector. The main differences stem from a larger set of hydrodynamic variables, which in turn provide us with more conservation equations and thus more hydrodynamic modes. The variables for the parallel sector are: energy density ε(t, x); momentum π (t, x); and field component λ (t, x). Respectively, the sources are the temperature T (t, x); velocity v (t, x); and displacement s (t, x). For the given set of variables and sources the conservation equations areε  where c V denotes the specific heat c V = ∂ε/∂T . We start by considering the equation of energy conservation (C.34). The temporal component of the stress-tensor is given by (C.9), hence at linear order we have Using the susceptibilities to write the stress-tensor in terms of hydrodynamic variables, Fourier transforming, and evaluating derivatives, gives the expression for the energy conservation equation −iω +κ 0 c V k 2 δε(ω, k) + ikδπ (ω, k) + γ 2 (κ + G)T k 2 δλ (ω, k) = 0. (C.45) Moving on, let us turn our attention to the equation of momentum conservation (C.35). The spatial stress-tensor components are given by equation (C.10), which at linear order gives ∂ x τ xx (t, x) = ∂ x δp(t, x) − (κ + G)∂ x δλ (t, x) − η∂ 2 x v (t, x) + O(δ 2 ), (C. 46) where we used ∂ x Φ (1) (t, x) = λ (t, x). In terms of the hydrodynamic variables the above equation is equal to Fourier-transforming and differentiating the momentum conservation equation thus yields −iω + η χ ππ k 2 δπ (ω, k) + v 2 u ikδε(ω, k) − ik(κ + G)δλ (ω, k) = 0. (C. 48) Finally, following all the above steps for re-expressing the Josephson relation, yields at linear order −iω + ξ (κ + G)k 2 δλ (ω, k) − ik χ ππ δπ (ω, k) + γ 2 c V k 2 δε(ω, k) = 0. (C. 49) In the same fashion as for the transverse sector, the equations (C.45), (C. 48) and (C.49) can be expressed as a matrix whose determinant must be zero in order to yield non-trivial solutions. This constraint on the determinant leads to a cubic equation in ω, hence the system we are considering will have three frequency modes. Two propagating modes (moving in opposite directions) and one diffusive mode. The generic dispersion relation of the longitudinal propagating mode is ω = ±c L k − iD p k 2 . (C.50) The speed of sound for this mode is given by Moreover, the dampening is D p = 1 2 η χ ππ + 1 2 c V (κ + G) 2 ξ − (κ + G)c V (∂p/∂ε)T γ 2 + (∂p/∂ε)κ 0 χ ππ − γ 2 (κ + G)χ ππ c V (κ + G + (∂p/∂ε)χ ππ ) .

C.4 Kubo Formulas
The susceptibilities and conservation equations presented in the above subsection can be used in the standard machinery for obtaining retarded Green's functions G R OaO b (ω, k). In turn, these Green's functions provide us with the following, especially relevant, Kubo formulas In the above formulas we have expressed the Green's functions in terms of Φ (1) (t, x) by relating the Goldstone field to λ (t, x), which is done by Fourier-transforming and evaluating the derivative of (C.3), together with the structure of (B.4). 20 Assuming that the speed of sound c 2 L is given by (C.51), the Kubo formula (C.61) provides us with a definition for κ via the expression κ = lim Further discussions regarding the results obtained from the above Kubo formulas can be found in the main text.