Gravitational quasinormal modes for Lifshitz black branes

We study the scalar and vector channels of gravitational quasinormal modes for Lifshitz black branes emerging in Einstein-Maxwell-Dilaton and Einstein-Proca theories in four and five dimensions, finding significant differences between the two models. In particular, rather surprisingly, in the Einstein-Maxwell-Dilaton model the dispersion relations for the shear and sound modes are given by $\omega_{shear} \sim -i\,k^4$ and $\omega_{sound}\sim-i \,k^2$, while in the Einstein-Proca model they take the more conventional form $\omega_{shear} \sim -i\,k^2$ and $\omega_{sound}\sim k$, the proportionality constants depend on the dynamical exponent and the appropriate factors of temperature. Through the holographic duality, this calculation provides information about the relaxation of the momentum and energy flux operators in a putative dual Lifshitz field theory. Comparing with the dispersion relations obtained directly by considering Lifshitz hydrodynamics suggest that the mass density of the equilibrium state in the Einstein-Maxwell-Dilaton model is infinite.


Introduction
Black holes, unlike many idealised physical systems, are intrinsically dissipative due to the presence of an event horizon. Thus, when considering the characteristic oscillations of these systems, instead of carrying out a standard normal-mode analysis, one reverts to the computation of quasinormal modes (QNMs). The latter have in principle complex frequencies, with the real part representing the actual frequency of the oscillation and the imaginary part representing the damping.
QNMs have been studied in various contexts. On the one hand, in asymptotically flat spacetimes they have been used in the context of black hole spectroscopy to infer the mass and angular momentum of the final black hole created after a binary merger as well as for testing no-hair theorems. On the other hand, in asymptotically AdS spacetimes, QNMs of black branes have been used for studying the near-equilibrium behaviour of strongly coupled plasmas with a dual gravity description revealing intriguing connections between the dynamics of horizons and relativistic hydrodynamics. For a recent review on QNMs for asymptotically flat and asymptotically AdS black holes/branes see [1].
In this paper we study QNMs in black brane geometries that asymptote to the so-called Lifshitz geometry 1 described by the line element (for a review see [3]).
where z is the critical exponent and i = 1, . . . , D − 2. These geometries manifestly realise the Lifshitz symmetry Lif D−2 (z) (rather than relativistic invariance) which comprises of temporal (H) and spatial (P i ) translations, spatial rotations (L ij ) as well as a scaling symmetry D z H : t → t = t + a P i : When z = 1 the metric is AdS and has full relativistic symmetry, but for z = 1 the system obtains anisotropic scaling between space and time and is thus non-relativistic. For z ≥ 1, Lifshitz geometries satisfy the strong energy condition R mn u m u n ≥ 0 for any future directed timelike vector u m , as well as the null energy condition G mn k m k n ≥ 0 for any future directed null vector k m . Therefore there are no obstruction to supporting the Lifshitz geometry with physically reasonable matter for z ≥ 1. For z < 1, some pathologies emerge caused by the violation of the null energy condition [4]. Lifshitz geometries have an anisotropic curvature tensor and solve the Einstein equations with a non-trivial stress energy tensor. They were first constructed in [5] and have been mainly realised in Einstein-Maxwell-Dilaton (EMD) [2,6], Einstein-Proca (EP) [2] and higher derivative gravity theories [2,[7][8][9]. Lifshitz black holes have been constructed in various dimensions and for different values of z, with some analytic and many numerical examples, see e.g. [10][11][12].
According to the holographic duality, these geometries are dual to strongly-coupled non-relativistic field theories with Lifshitz symmetry. However, the corresponding holographic dictionary is less well-understood than in the case of AdS and in fact a real-time formulation is still not fully developed. Holographic renormalisation, which is what allows one to determine the independent sources and their corresponding expectation values, is more elaborate [13][14][15]. These references find that the operator content of the dual theory corresponds to a non-relativistic stress tensor complex, consisting of an energy density E, energy flux E i , momentum density P i and spatial stress-tensor π ij . In the absence of sources, these satisfy the usual Ward identities In addition, a common feature present in holographic realizations of Lifshitz is the appearance of a scalar operator with irrational scaling dimensions. This can be associated to the longitudinal polarization of the massive vector in EP, or the dilaton in EMD theories. Recent studies of a specific z = 2 Lifshitz geometry arising as reduction of higherdimensional theories [16] suggested that Lifshitz field theories couple generically to Newton-Cartan geometries on the boundary -note that z = 2 is special as the symmetry group can be augmented by Galilean boosts. This study has lead to some resurgence of interest in nonrelativistic holography, including studies of intrinsically non-relativistic gravity theories such as Horava-Lifshitz -this is to be contrasted with what has been stated above, where one considers relativistic gravity theories that admit non-relativistic solutions. Scalar field QNMs in Lifshitz black brane backgrounds were considered previously in the literature. In particular, the cases that have been studied include D = 3, z = 3 in NewMassiveGravity [17], D ≥ 4, z = 2 in R 2 gravity [18] and D ≥ 2, z = 2 in a R 3 gravity [19]. Additionally, the case of D = 4, z = 2 has been studied in the Einstein-Proca-Scalar theory [20], in the EMD setup [21] and in a topological black hole in a Einstein-Maxwell-Proca background [22][23][24][25]. All these quasinormal modes were found to be purely imaginary. This is in contrast with the results of [26] for the EMD model that showed that, at zero momentum, these modes have purely imaginary frequencies for z ≥ D − 2, while for z < D − 2 they pick up a real part. In all case, all quasinormal modes are situated in the lower half-plane of complex frequencies, indicating stability. Furthermore, in the context of holographic superconductors in a Lifshitz background, [27] found hydrodynamic modes in the QNM spectrum of the charged scalar field close to the critical temperature.
In this work we study gravitational QNMs, focusing on the EP and the EMD model for D = 4, 5 and general z. These QNMs are more interesting than scalar ones because the corresponding fluctuations couple to conserved currents in the dual field theory. Currently, these have only been calculated for the EP model for D = 4, z = 2 at zero temperature [28] and negative imaginary frequencies were found. Our results for the leading behaviour of the dispersion relations of the gapless (hydrodynamic) modes are summarised in the table below where ω, k 1 are the frequency and the momentum of the modes respectively. Note that all the z dependence is hidden in the constant of proportionality that we have calculated numerically. We see substantial differences in the relaxation of the two theories, and in the case of the EMD theory, we also see significant deviation from the dispersion relations of the hydrodynamic modes in AdS D , which take the form where u s is the speed of sound, Γ is the attenuation and ν is the diffusion constant. In addition to the above hydrodynamic modes, we also find a tower of non-hydrodynamic modes the real parts of which go to zero at z = D − 2, just like for the scalar QNMs [26]. All the quasinormal modes found are located in the lower half plane and we thus conclude that the system is stable, which is in agreement with the non-linear time evolution of [29] for the EMD model. Through the holographic duality, QNMs correspond to poles of the (retarded) thermal correlators of dual (D − 1)-dimensional strongly interacting quantum field theories. The lowest QNM frequencies of black branes have a direct interpretation as dispersion relations of hydrodynamic excitations in the dual field theory, which in our case enjoys Lifshitz symmetry. Lifshitz hydrodynamics have been developed in [30,31] and more recently in [32]. This has been carried out in two competing approaches: both start with relativistic hydrodynamics but break the Lorentz symmetry in different ways along the hydrodynamic expansion. In particular, the difference lies on whether Galilean boosts are broken at the perfect fluid level or at the first dissipative order. In both cases new transport coefficients were identified and the bulk viscosity is found to vanish [30,32]. Dispersion relations for the hydrodynamic sound, shear and diffusion modes have been studied in [32] and a new expression for the speed of sound has been obtained [32,33]. The sound mode dispersion relation in a higher-derivative gravity theory for z = 3, D = 3 at finite temperature in the hydrodynamic limit was also studied in [34].
The remaining on this paper is structured as follows. In section 2 we give a brief overview of the computation of gravitational QNMs and in sections 3 and 4 we discuss, respectively, EP and the EMD theories and the corresponding results that we obtained. Then in section 5 we carry out a comparison with the Lifshitz hydrodynamics developed in [32] and finally in section 6 we conclude with some discussion and future directions.

Review of quasinormal mode computation
A detailed analysis of the quasinormal spectra for AdS-Schwarzschild and AdS-Reissner-Nordstrom black branes was discussed previously in the literature [35]. It is well known that the electromagnetic and gravitational perturbations split into the tensor (for D ≥ 5), vector and scalar sectors depending on their transformation properties. The scalar sector contains the sound and charge diffusion fluctuations, the vector contains the shear and transverse gauge field fluctuations and the tensor modes are decoupled scalar equations. Out of these fluctuations only the shear, sound and charge diffusion contain hydrodynamic modes, meaning modes that obey dispersion relations such that the frequency approaches zero as the momentum is decreased. These hydrodynamic modes correspond precisely to the shear and sound modes in the hydrodynamic limit of the dual CFT at finite T we are interested in.
In the case of black branes with Lifshitz asymptotics, the electromagnetic and gravitational fluctuations once again split into the three sectors described above. In this work we focus only on the scalar and vector fluctuation, which contain gapless (hydrodynamic) modes, meaning modes that obey dispersion relations such that the frequency approaches zero as the momentum is decreased. Note that in our case, even though we have a vector field in the bulk, there is no U(1) current in the dual theory, so only sound and shear modes associated to the stress tensor are expected. In particular, we do not expect to find charge diffusion hydrodynamic mode in the set-ups we are considering.
In particular, we consider linearised fluctuations around the Lifshitz backgrounds of the form where is a small expansion parameter. Note that here we have chosen the momentum k to point in the x direction without lose of generality. The precise form of the Lifshitz background, denoted by g (Lif z) , A (Lif z) and φ (Lif z) is described in detail in the sections below.
In the vector sector, the non-trivial metric and gauge field fluctuations considered are where i denotes the spatial directions transverse to x, in which we retain isotropic e.g. δg xy = δg xz in D = 5.
In the scalar sector, the non-trivial metric and gauge field fluctuations considered are The equations of motion for these fluctuations carry a lot of redundant information due to gauge invariance. In particular, under an infinitesimal coordinate transformation x µ → x µ + ξ µ , where ξ µ is an arbitrary function of r, the fluctuations transform as Typically this gauge freedom is dealt with by arranging the fluctuations into gauge invariant combinations [35]; this is the approach we follow in the majority of our calculations. For the scalar channel of the Einstein-Proca model we instead follow a method first discussed in [36] that mimics the DeTurck trick 2 . Specifically, we add a gauge fixing term ∇ (ν τ µ) to Einstein's equations, where and a posteriori check that it vanishes. Note that such term is only added for the metric and not for the gauge field. Since the gauge field is massive, the equation of motion completely determines A, which means that one can not do any gauge transformations. Comparing the two approaches, considering gauge invariant combinations is preferred as it boils down to a smaller number of equations to be solved numerically.
The final equations are then solved numerically subject to boundary conditions, namely ingoing boundary conditions at the horizon and fast enough fall-off close to the UV boundary compatible with the absence of sources.

The model
In this section we consider the Einstein-Proca model, described by the following bulk action where F = dA. The corresponding equations of motion are given by and it is easy to show that for they admit a solution which corresponds to a dual field theory at zero temperature with Lifshitz symmetry.
In order to put this system at finite temperature, one needs to consider configurations with regular horizons that approach the above the solution close to the boundary. This was achieved in [6], using the ansatz which gives rise to two second order equations for F, a t and an algebraic equation for R. These equations can be solved numerically using a shooting method subject to appropriate UV boundary conditions and regularity at the horizon (located at r = r h ). In particular, close to the boundary (r → ∞), the fields fall off according to where c 1 , c ± are undetermined constants, c 2 , c 3 , c 4 are functions of z (which we omit for simplicity) and In our UV boundary conditions we require that c − = 0, following the analysis of [13]. The temperature of the solution is finite and given by In Fig 1 we plot the temperature dependence with z for the numerical solutions we have constructed.  Figure 1: Temperature as a function of z for Lifshitz black branes in the Einstein-Proca model for D = 4 (blue) and D = 5 (red) . Note that the curves approach the AdS values T D=4 = 3r h /(4π), T D=5 = r h /π at z = 1.

Vector channel
Given the vector channel fluctuations, the equations of motion reduce to a set of two second order ODEs for the gauge invariant quantities δA y , δH xy = r 4 ∂ r (r −2 δg xy ) , which we omit for simplicity. The independent terms in the boundary asymptotics for the gauge invariants can be written as It is easy to see that setting the sources to zero requires a (0) = H (0) = 0 [37]. At the horizon we require ingoing boundary conditions, which implies that the fields behave as where T is the temperature and a admit regular power series expansions in the near horizon region. We determine the spectrum of QNMs by discretising the differential equations and solving the corresponding matrix problem in Mathematica. We find one gapless mode ω shear = −iν(z)k 2 + . . .

Scalar channel
We now consider the scalar channel fluctuations around the Lifshitz background (3.5). As explained in the previous section, we follow the gauge fixing approach of [36], which requires introducing an extra term in the equations of motion for the metric -we have verified that this term vanishes within numerical precision for the modes of interest, so that these fluctuations are indeed solutions of the linearised Einstein's equations.
The equations of motion for these fluctuations are solved numerically by discretisation subject to ingoing boundary conditions at the horizon and fast enough fall-off close to the boundary compatible with the absence of sources; see appendix A for more details. We find one gapless mode characterised by the standard dispersion relation

The model
In this section we consider the theory described by the bulk action which gives rise to the following equations of motions As it was shown in [2] (for generalisations see [39]), this theory admits an analytic solution black hole configuration given by .

(4.4)
It is easy to see that close to the boundary r → ∞ this solution approaches the Lifshitz metric with critical exponent z, while the black hole horizon is located at r = r h corresponding to temperature where r h is the radius of the horizon and in our units r h = 1. It is worth mentioning that these configurations are not smoothly connected to AdS branes in the limit z → 1. This is easily seen by noting the divergence in the coupling λ which appears explicitly in the equation of motion for the dilaton.

Vector channel
We consider the vector channel fluctuations around the background (4.3) in the gauge where δg ri = 0 and we focus our attention on the gauge invariant quantities δA y , δH ty = a 0 r D−z+1 h ∂ r (r −2 δg ty ) , (4.6) which obey second order ODEs. The independent terms in the boundary asymptotics for the gauge invariants can be written as The leading terms parametrise the field theory sources so we set them to zero. Given these boundary conditions in the UV and ingoing boundary conditions at the horizon, we proceed to solve the ODEs numerically to find the spectrum of QNMs. We find one gapless mode with dispersion In Fig 5(a) show the dispersion relation for D = 4, z = 3/2, which behaves line ∼ k 4 in the hydrodynamic limit -this was also confirmed by studying the logarithmic derivative of the dispersion kIm[ω] /Im[ω]. We find an analogous behaviour for D = 5. The dispersion constant,ν, depends on the value of z in a way which we depict in Fig 5(b).

Scalar channel
In this section we consider the scalar sector fluctuations around the background (4.3). In particular, in the gauge where δg tr , δg rr , δg rx , δa r = 0, these fluctuations are combined into the following gauge invariant combinations The equations of motion for these master fields take the form of three second order ODEs of degree 6 in the frequency. At the horizon we impose ingoing boundary conditions where Z (reg) i admit regular power series expansions. A mode analysis close to the boundary of the form reveals modes with scaling dimensions More concretely, the boundary conditions for the master fields take the form where c a , c b , c c , c d , c e , c f , c g are known functions of z -the UV asymptotics of the field fluctuations are recorded in appendix B. Given the above expansion, by settingc 1 , c 2 , c 3 = 0 we are demanding that the boundary sources vanish. Imposing these boundary conditions, we solve the corresponding eigenvalue problem to find the following dispersion relation for the sound mode ω sound = −iΓ(z)k 2 + . . . , (4.14) indicating that for the Einstein-Maxwell-Dilaton case the speed of sound is zero. In Fig.  6(a), we plot the dispersion relation for z = 3/2, D = 4 and we see a clear quadratic scaling -this was also confirmed by studying the logarithmic derivative of the dispersion kIm[ω] /Im[ω]. We find analogous behaviour for other values of z and D. We depict the behaviour ofΓ(z) in Fig. 6(b). We confirmed these results by repeating the calculation using a gauge fixing term instead of master fields.

Analytic results for z = D − 2 in the vector channel
In this section we carry a perturbative analysis in an attempt to get an analytic handle on the numerical results we obtained for the shear channel. Due to the non-analytic behaviour of the metric functions, the hydrodynamic expansion can only be pushed analytically for z = D − 2. In particular, we consider a perturbative expansion in the momentum k of the form where u = r h /r, β = z + D − 2 and the prefactors were chosen in such a way to simplify the IR boundary conditions. Note that, in units of r h = 1, we have 4πT = z + D − 2. At each order, i, we find two second order equations, which we can recast as a single higher order ODE for δa (i) y (u). Solving these equations, we find that at zeroth order while at order 1 (and in fact for all odd orders) Finally, we were able to solve the corresponding ODE at quadratic order and determinē ω (2) = 0 through imposing boundary conditions. Focusing on D = 4, in the IR the solution looks like (4.19) and in the UV δa (2) y (u) = c 0 π 2 96k 2 + c 2 + c 1 u 2 72 (9 − 14k 2 + 12k 2 log u) + . . . , (4.20) δH (2) ty (u) = − c 0 36 u 2 (9 − 14k 2 + 12k 2 log u) + . . . . (4.21) Similar results apply for D = 5. Determining analytic solutions to the higher order equations was not possible. To summarise, from this perturbative solution we have been able to conclude that the quadratic piece vanishes,ω (2) = 0, which is consistent with our numerical results. Unfortunately we were not able to solve the equations to a sufficiently high order in perturbation theory to determine the value of the first non-trivial term, namelyω (4) .

Comparing with Lifshitz hydrodynamics
The hydrodynamic modes of homogeneous and isotropic non-relativistic fluids with generic dynamical exponent z was discussed in [30][31][32], with the main difference between the two approaches being that in [30,31] boosts are broken at the first dissipative level while in [32] they are already absent at the ideal level. In particular, [32] derived the following hydrodynamic dispersion relations for sound and shear where ρ is the mass density, η is the charge density and the zero index denotes equilibrium. They also expressed the speed of sound u s and the attenuation Γ in terms of equilibrium quantities as follows where d = D − 1,˜ is the internal energy, p is the pressure andπ is the dissipative part of the thermal conductivity.
In the case of the Einstein-Proca model the dispersion relations we found numerically is consistent with the behaviour above for non-vanishing values of all parameters, in a way which smoothly connects to the AdS values as z → 1. However, for the EMD model we numerically find that the speed of sound as well as the leading term in the shear dispersion relation vanish. The latter was also confirmed analytically for a particular choice of z. To reconcile the two results, ρ 0 → ∞ together withπ 0 ∼ ρ 2 0 while the rest of the coefficients are order one. Checking this explicitly is left for future work.

Discussion and Future Direction
In this work we have considered linearised perturbations around Lifshitz black branes, with a focus on electromagnetic and gravitational fluctuations in the scalar and vector channels.
We found that in the case of the Einstein-Proca model, the dispersion relations of the hydrodynamic shear and sound mode have the same structure as in asymptotically AdS black branes, namely they have quadratic and linear dispersions, respectively. On the other hand, we found that the dispersions of the shear and sound modes in Einstein-Maxwell-Dilaton model are quartic and quadratic, respectively. This is a significant difference between the two models, signalling that the two models approach equilibrium differently. It would be interesting to understand this better from the field theory perspective. As mentioned in Sec. 4, the z → 1 limit of the EMD system does not smoothly approach the relativistic Einstein-Maxwell theory, at least taken while keeping all other quantities regular in (z − 1). This could explain the root of the qualitative change we observe in the dispersion relations of the hydrodynamic modes.
Preliminary analysis of the near-horizon expansion of the perturbation indicates the phenomenon of pole-skipping does occur in asymptotically Lifshitz geometries. Specifically, we find that the Matsubara frequencies exist in the lower half plane at the exact same locations as in the relativistic case. It would be interesting to check if the Lifshitz hydrodynamic sound mode, when driven to instability by a choice of a specific value of imaginary momentum that is well outside the standard regime of validity of hydrodynamics, exhibits connections with chaos through the phenomenon of pole-skipping [40]. In addition, in the spirit of [41,42], it would be interesting to further investigate the radius of convergence of the non-relativistic hydrodynamic expansion and, through resummation, test whether it is possible to extract non-hydrodynamic QNMs from the hydrodynamic ones. This will be shed light on the thermalisation properties of strongly coupled nonrelativistic fluids -for relativistic fluid thermalisation, it is known that the hydrodynamic expansion becomes applicable very early during dynamical processes.
An exciting recent technical development with the potential of leading to rapid progress in both gravitational physics and holography involves the limit of large number of dimensions (D) of general relativity [43,44]. This tool has so far been applied to relativistic theories of gravity, realised on geometries that are either asymptotically flat or asymptotically AdS, providing a number of very interesting results. By treating D as a free parameter, one can use this approach to perform a perturbative expansion in 1/D, which leads to a drastic simplification of the theory. The non-trivial black hole dynamics are localised within a distance 1/D from the horizon and it is thus possible to capture them with an effective theory given by a set of constraints that depend solely on the directions parallel to the horizon. It would be interesting to explore the applicability of this tool in spacetimes with Lifshitz asymptotics. The first step in this direction would be to see if the tower of gravitational QNMs splits into two subgroups, one controlling the dynamics in the near horizon region and the other in the far region, when D → ∞.
form [37] δg tt =r 2z F (r) g tr r −z−D+1 + . . . , δg tx =r 2z g tx r −z−D+2 + . . . where γ ± = − Note that (A.2) is meant to capture the general form of the expansion and not all coefficients are independent. Furthermore, just like in the DeTurk trick, the non-analytic terms, γ ± , are an artefact of the way we fix the gauge and the free constants that multiply it are in fact zero on the actual solution. In this expansion we identify the leading coefficients g tx r −z−D+2 + . . .
tx r z−D + g tx (c) r −z−D+2 + . . . , δg xx =r 2 g (0) xx + g (v) xx r 2−D−z + . . . , δg =r 2 g (0) yy + g (v) yy r 2−D−z + . . . , which, at least for the majority of the metric fields, looks similar to (A.2). Note that (B.1) is meant to capture the general form of the expansion and not all coefficients are independent. Close to the UV boundary we impose boundary conditions that kill the leading modes g x , φ (0) as they correspond to boundary sources.