Dispersion relations in non-relativistic two-dimensional materials from quasinormal modes in Hǒrava Gravity

We compute dispersion relations of non-hydrodynamic and hydrodynamic modes in a non-relativistic strongly coupled two-dimensional quantum field theory. This is achieved by numerically computing quasinormal modes (QNMs) of a particular analytically known black brane solution to 3+1-dimensional Hǒrava Gravity. Hǒrava Gravity is distinguished from Einstein Gravity by the presence of a scalar field, termed the khronon, defining a preferred time-foliation. Surprisingly, for this black brane solution, the khronon fluctuation numerically decouples from all others, having its own set of purely imaginary eigenfrequencies, for which we provide an analytic expression. All other Hǒrava Gravity QNMs are expressed analytically in terms of QNMs of Einstein Gravity, in units involving the khronon coupling constants and various horizons. Our numerical computation reproduces the analytically known momentum diffusion mode, and extends the analytic expression for the sound modes to a wide range of khronon coupling values. In the eikonal limit (large momentum limit), the analytically known dispersion of QNM frequencies with the momentum is reproduced by our numerics. We provide a parametrization of all QNM frequencies to fourth order in the momentum. We demonstrate perturbative stability in a wide range of coupling constants and momenta.


Introduction
Right after the discovery of the gauge/gravity correspondence [1], it was applied in parallel with hydrodynamics to describe strongly coupled fluids [2][3][4], which are relativistic. Most experiments, however, are conducted in non-relativistic materials, some of which may be described best as a non-relativistic strongly coupled fluid.
A holographic correspondence has been conjectured between Hořava Gravity [5][6][7][8] and non-relativistic quantum field theories [9,10]. An analytical Hořava black brane solution was found [11], and its fluctuations were studied analytically in the hydrodynamic limit [12]. A hydrodynamic momentum diffusion mode was found via a field redefinition in what we will refer to as axial sector of the theory. That field redefinition maps the axial sector of Hořava Gravity to the corresponding sector of Einstein Gravity. However, in the other sector, the polar sector, this map could only be performed in the special case of setting one

JHEP10(2019)087
Hořava coupling constant to zero, λ = 0. Two hydrodynamic sound modes were found as expected. But this failure of the map at λ = 0 motivates a closer study of the polar sector.
For our purposes, Hořava Gravity is Einstein Gravity with the addition of a scalar field, the khronon. This scalar field defines a preferred time-foliation and thus breaks Lorentz invariance. Generally, non-relativistic theories exhibit Lifshitz scaling, defined as different scaling in time, t, and spatial, x I coordinates, i.e. t → κ z t, x I → κx I , with constant κ and dynamical exponent z. The analytic Hořava black brane solution is peculiar in that it has z = 1. While this solution still breaks Lorentz boost invariance, and the theory is non-relativistic, see discussion in section 2, its Lifshitz scaling coefficient z = 1 is that of a relativistic theory. Furthermore, in the case of field theories with Lifshitz scaling z and dimension d it has been claimed that these contain only overdamped eigenmodes, i.e. quasinormal modes on the imaginary frequency axis, in their spectrum if d ≤ z + 1 at vanishing momentum [13,14]. The analytic Hořava black hole solution has z = 1, d = 3 and should thus not be overdamped, should it follow the behavior described in [13,14]. These references [13,14] consider a different system though, namely a probe scalar in a background with Lifshitz scaling. Therefore it is not clear if also quasinormal modes in Hořava Gravity should show this behavior. The peculiar value of z = 1, imposed by the one and only analytic Hořava black brane background available to us from the solutions described in [15], 1 motivates us to investigate the quasinormal modes of this solution in order to see if there are overdamped modes and to investigate the dispersion relations of this theory. It is a priori not clear, if they should be relativistic or non-relativistic.
In Hořava Gravity, fields can travel at distinct speeds that can even be infinitely large [11]. These speeds are set by the Hořava coupling constants. The existence of distinct speeds implies that fields experience different horizons, i.e. last points of return, in the presence of a black brane. In [12] it was conjectured that the relevant horizon for each set of fields is the sound horizon, which is identified as a regular singular point of the set of equations of motion of those fields. For the axial sector the relevant horizon was shown to be the horizon of the spin-2 graviton. In this work we confirm that in the polar sector the relevant horizon is the spin-0 graviton horizon, which here is identical to the universal horizon r h . The universal horizon determines the temperature in the dual field theory [11].
Hence, the main goal of this work is to calculate all quasinormal modes in the axial and polar sectors for a wide range of Hořava coupling values, mode numbers and momenta. From this, we extract the dispersion relation of each mode. The questions mentioned above will be answered in passing. In the polar sector, see section 3, the fluctuation equations for Hořava fields are very complicated. Hence, we use the equivalence between Hořava Gravity and Einstein-AEther Theory when the scalar field, the khronon, is required to be hypersurface-orthogonal [16,17]. In the axial sector we show very good numerical agreement of quasinormal modes from both theories. Compared to Einstein Gravity, in Hořava Gravity, an additional set of quasinormal modes is to be expected, namely those contributed by the khronon. These khronon modes turn out to be very special and we discuss them in detail in section 3.5. 1 These black brane solutions are smoothly connected to solutions with all possible two-derivative Hořava couplings non-zero [15].

JHEP10(2019)087
We study long-lived modes determining the behavior of the system at late times. Two types of modes are long-lived because their damping is small: hydrodynamic modes with small frequency and momentum, and non-hydrodynamic modes with large momentum. Analytic results in the large momentum (eikonal) limit [18][19][20], are a useful tool to check numerical accuracy and the structure of the quasinormal mode spectrum, see section 3.3. Our results are discussed in section 2.5 and 3.2, but visually summarized in figure 1 and figure 2. The parameter space of Hořava coupling constants λ and β is shown. The third coupling is set to zero in this work, α = 0. A summary of results and our conclusions are found in section 4. Equations of motion, and QNM data is collected in the supplementary material attached to this paper.
The two-dimensional materials we have in mind for application of our results could be thin films (e.g. NbSi [22]) or layered structures such as those found in high critical temperature cuprates, see e.g. [22,23] and more generally [24,25]. We assume that such materials -considered in an appropriate regime of their defining quantities -can be described by non-relativistic hydrodynamics [22]. In this regime, conserved quantities such as energy, momentum, and electric charge are conserved and lead to long-lived modes, for example sound modes. Outside of this hydrodynamic regime, such materials may contain non-hydrodynamic modes. Some of those may be long-lived as well. The quasinormal mode frequencies as a function of spatial momentum, which we compute in this work, holographically correspond to dispersion relations for the long-lived and short-lived modes which may occur in such materials. However, the system we study in this work is only a step towards non-relativistic descriptions of such materials. This is because we are working in a background which imposes that time scales just like spatial directions. But we will see that our system still singles out the time direction, allows for arbitrarily high velocities, and has further non-relativistic properties discussed below. So we are studying a putative two-dimensional material, sharing symmetries and hydrodynamic description with real materials. We are not attempting to describe the microscopic details of such real world materials.

Hořava Gravity
In this section, we will summarize those aspects of Hořava Gravity of relevance to our calculation of QNMs. Our main interest are all the QNMs of fluctuations around a particular analytically known Hořava black brane solution found by Janiszewski [11]. We close this section with a review of the (numerical) shooting method with which then the Hořava QNMs of the axial sector are computed. In the hydrodynamic limit, axial and polar 2 transport quantities have been found analytically in [12], 3 which will serve as a check of our numerical solutions. 2 That previous analysis specialized to the case of vanishing coupling λ = 0 for polar QNMs. In this work, we lift this severe restriction and analyze QNMs for nonzero λ.

Action, coupling constants & field content
We will specialize to a particular low energy solution of classical Hořava Gravity [11]. The degrees of freedom are represented by G IJ , N I , and N , which are constituents of the ADM decomposition of the spacetime metric, g XY ; 4 where we will be using the mostly positive metric convention. 5 In terms of spacetime coordinates x X = {t, r, x, y} the line element then is given by In (3+1) dimensions the analytically known metric solution [11] satisfies the equations of motion generated by the variation of the following action of Hořava Gravity: K IJ , G IJ , N I and N are the extrinsic curvature tensor, spatial metric, shift vector, and lapse function, respectively. K is the trace of K IJ . G is the determinant of G IJ . The lowering and raising of spatial indices are carried out by contraction with G IJ and G IJ , respectively, while ∇ I is the covariant derivative defined with respect to the spatial metric G IJ . In order for the propagation speeds of the graviton to be strictly positive the coupling constants,(α, β, λ), must obey the following inequalities [16,21]: These allowed regions of parameter space are represented by white surfaces plus the blue line along the β-axis in the (β,λ)-parameter plots figure 1 and figure 2; forbidden regions are colored red.

Hořava black brane background solution
Taking α = 0 and Λ = −3/L 2 (as usual, the AdS radius L can be set to unity, L = 1, by scaling symmetries of the equations of motion), we here review the aforementioned asymptotically AdS 4 black brane solution to Hořava gravity found in [11]. G IJ , N I , and N are known analytically and given by 6

JHEP10(2019)087
which notably is independent from the coupling λ. The AdS radius is L = 1, the time coordinate is t ≡ x 0 , spatial coordinates x ≡ x 2 , y ≡ x 3 , and the radial coordinate is r ≡ x 1 . The AdS boundary lies at r = 0 and the universal horizon at r = r h . This horizon is a trapping surface from which none of the Hořava Gravity fields can escape. In particular, the universal horizon is also the sound horizon 7 for the spin-0 graviton, which travels with infinite speed in this particular solution. Since Hořava Gravity is non-relativistic, there exists another horizon: the spin-2 sound horizon at r = r h /2 1 3 , which is a trapping surface for the spin-2 graviton. There is also the Killing horizon located at r = r k [12]. The temperature is determined by the universal horizon [11] and is given by In general, solutions to Hořava Gravity can exhibit Lifshitz scaling symmetry under the scaling with a constant κ: t → κ z t, x I → κx I , r → κr with the dynamical exponent z, leading to the asymptotic scaling ds 2 ∼ dt 2 r 2z + dx I 2 +dr 2 r 2 [21]. The metric solution in eq. (2.5) is an example with z = 1, which means that time and space coordinates scale the same way. However, time and spatial coordinates are still distinct because there exists a time-like vector (or gradient of the khronon) specifying a preferred time-foliation. To see this symmetry explicitly it is helpful to write Hořava Gravity as Einstein-AEther Theory with a particular constraint (hypersurface orthogonality) [12,16], as we will see below in section 3. Hence, the solution (2.5), like any typical solution to Hořava Gravity, is only invariant under those diffeomorphisms preserving the time-foliation.

Hořava black brane perturbations
One may perturb the metric (2.1) around a background value of (2.5) to linear order by a metric perturbation h XY g p where 1. Requiring a vanishing variation of the action with respect to h XY generates 10 coupled linear equations of motion for the 10 independent components of h XY . Since our action (2.2) and background metric (2.5) have a translational symmetry in x, ydirections (2.5), we make the standard plane wave ansatz of Fourier expanding into momentum space modesh XY (r; ω, k). 8 Since action and metric also are invariant under spatial rotations in the x, y-plane, without loss of generality we choose the momentum vector to point into the y-direction h XY (t, r, x I , r) = dωdk e i(ky−ωt)h XY (r; ω, k) r 2 . (2.8) Once (2.8) is substituted into the equations of motion we find that the 10 equations of motion decouple into a set of 7 equations for metric components which are odd under 7 We define the sound horizon of a field as the location along the radial AdS-coordinate (excluding the AdS-boundary), at which the fluctuation equation for that field contains singular coefficients when the leading derivative is normalized to have coefficient 1. For example, φ (r) + b(r)/(r − rc)φ (r) + b(r)/(r − rc) 2 φ(r) = 0 has a sound horizon at r = rc, if b(r) and c(r) can each be expanded in a Taylor series. 8 Recall that x ≡ x 2 and y ≡ x 3 .

JHEP10(2019)087
parity (axial ) and a set of 3 equations for metric components which are even under parity (polar ) [12]. For the QNMs in the Hořava case we only concern ourselves with the set of 3 equations of motion for the axial fieldsh xy ,h xt , andh xr , which are odd under parity. After a radial diffeomorphism,h xr (r; ω, k) can be set to vanish. A linear combination of these fields turns out to be a gauge invariant master field, ψ(r; ω, k) = (ωh xy + kh tx ), and obeys the following single equation of motion 9 found in [12] where the following variables have been rescaled to be dimensionless Eq. (2.9) is the master equation of motion for axial perturbations which we are going to solve in the remainder of this section to extract axial QNMs.

Shooting method
We intend to find QNMs, that is, we search for those solutions to fluctuation equations which satisfy two conditions: (i) modes are ingoing at the sound horizon, and (ii) vanish at the boundary, i.e. satisfy a Dirichlet boundary condition. With the rescaled-z coordinate (2.10), the sound horizon is located at z = 1, where the equation of motion (2.9) has a regular singular point; in fact, this regular singular point is the reason for us to call this location a sound horizon. We make the following near sound horizon ansatz separating the regular part F and irregular part of ψ(z; ν, q) from each other. Solving the equation (2.9) with the ansatz (2.11), one finds two α's that satisfy eq. (2.9) at the first nonzero order in the near-horizon expansion: α = −1 or α = − 2iν 3 . We choose the latter which corresponds to the ingoing mode. At the other expansion orders one can recursively solve for the f n (ν, q) coefficients. As usual, the series in eq. (2.11) is asymptotic, though one can assume validity for a small region around the sound horizon. Since there are singularities at z = 1 and z = 0 we numerically solve eq. (2.9) with Mathematica's NDSolve function on the restricted computational domain of 10 shown that for exceedingly small values of dr h the (2.11) ansatz fluctuates rapidly and can create large numerical errors [26], which guides our choice of dr h here. The boundary conditions, ψ(1 − dr h ) and ψ (1 − dr h ), are then provided by the horizon expansion (2.11).
For an arbitrary value of ν and q, this solution is not a quasinormal mode, i.e. ψ(dr b ) = 0. We use Mathematica's FindRoot to find ν such that ψ(dr b ; q) = 0. When compared to the pseudospectral method's modes (figure 3 and figure 6), the shooting method used here was found to be numerically less stable, especially for QNMs with larger momentum. Hence, in the next section we will switch to pseudospectral methods.

Hořava Gravity axial QNMs
It turns out that the axial QNMs ν found with the shooting method are numerically equal the QNMs one would find for an asymptotically AdS 4 Schwarzschild black brane within Einstein Gravity, i.e. ν Hořava axial QNM = ν Einstein axial QNM . However, the QNM frequencies ν are scaled by a factor of √ 1 + β and β vanishes in GR. So, up to numerical errors we empirically find the relation ω Hořava axial QNM = 1 + β ω Einstein axial QNM . (2.12) This particular β factor is in fact the speed of the spin-2 graviton [12]. At β = 0 the spin-2 graviton travels at unit speed, and with respect to just the axial perturbations, the theory returns to being relativistic. At small momentum our lowest lying mode agrees with the diffusion mode found in a hydrodynamic approximation given in eq. (3.31) through (3.34) of [12]. It must be mentioned that we attempted to find polar QNMs, however Mathematica's NDSolve used in the shooting method failed to find a converging solution to the fluctuation equations of motion. 11 It is difficult to determine indicial exponents in general for this coupled system, and to then separate the singular from the regular part of the fluctuations.
To circumvent these problems, we decided to calculate axial and polar QNMs in an equivalent theory, Einstein-AEther Gravity, using a different technique, namely pseudospectral methods. The equivalence of these two theories holds under the constraint of hypersurface orthogonality, discussed in section 3. Hence, QNMs found in the two theories are expected to be identical. 12 A comparison between axial QNMs computed in both theories is displayed in table 1. For the shooting 9 orders in the horizon expansion have been taken into account, the horizon and boundary cutoffs were chosen as r = 1 − 10 −3 and r = 10 −3 , respectively. For the computation of the Einstein-AEther QNMs a grid of size N grid = 80 was compared to one with N grid = 100 in order to determine convergent quasinormal modes, as described in appendix B. The percentage of deviation d of the Hořava Gravity QNM frequencies from the Einstein-AEther Theory QNMs, for that calculation see section 3, is given by which will be discussed in the next section.

Einstein-AEther theory
Einstein-AEther Theory differs from General Relativity by the addition of a scalar matter field φ, which is termed AEther field. We need a particular Einstein-AEther Theory, namely one that is equivalent to Hořava gravity. This requires the AEther field to define a timelike vector field which is hypersurface orthogonal. That goal is achieved by the following definition of the timelike vector which is a timelike hypersurface orthonormal unit vector matter field [4,17]. This vector field also breaks diffeomorphism invariance down to foliation preserving diffeomorphisms. In this context, the scalar φ is often referred to as khronon, determining the time foliation. Lower case Greek indices µ, ν denote spacetime indices. As coordinates, we choose x µ with x 0 = v, x 1 = r, x 2 = x, x 3 = y, and use the mostly positive metric convention, (-,+,+,+). A (3+1)-dimensional Einstein-AEther action has been constructed [17], which includes four quadratic derivative terms of u µ , 13

JHEP10(2019)087
where g is the determinant of the metric g µν , and R is the Ricci Scalar of the metric g µν . c 1 is redundant once we constrain u to be hypersurface orthogonal by construction. This allows us to rearrange (3.2) whose coupling constants then appear as linear combinations of c 1 , c 2 , c 3 , and c 4 as c 13 = c 1 + c 3 , c 2 , and c 14 = c 1 + c 4 [16], and we set c 1 = 0 without loss of generality. This equates Einstein-AEther theory to low energy Hořava Gravity after the identification −N = δ µ 0 u µ , and with the following coupling constants for the respective theories [27] Analogously to our perturbative treatment in the previous section around eq. (2.1), also here we will investigate linear perturbations around the black brane solution with α = 0 and Λ = −3 [11]. Our coordinates are similar to Eddington-Finkelstein coordinates [11], as seen in the metric where r h is the radius of the universal horizon. The the sign choice on f (r) and in u µ corresponds to the choice of infalling (lower signs) or outgoing (upper signs) Eddington-Finkelstein-like coordinates [11]. For this paper we choose the negative signs (f (r) = −1/r 2 and u r = −a(r)). This leads to infalling modes which are regular at their respective sound horizons. Note that eq. (3.4) reduces to the Schwarzschild AdS 4 metric with the Schwarzschild horizon located at r = r Schwarzschild = r h /2 1/3 when all remaining AEther couplings are set to zero, c 3 = 0 = c 2 , and we choose the vector field to vanish u µ = 0. As a side note, remarkably, one can express φ explicitly by integrating eq. (3.1) to obtain (3.5)

Einstein-AEther black brane perturbations
Similar to how we perturbed the Hořava black brane (2.7) we perturb the metric (3.4) by adding a "small" linear term where 1, we choose where χ(x σ ) is a scalar field. Since φ p is still a scalar field after the perturbation, replacing φ → φ p in eq. (3.1) ensures hypersurface orthogonality and normalization of the now JHEP10(2019)087 perturbed u vector. The h µν fields and χ obey eleven coupled linear equations. A Fourier transformation similar to (2.8), is applied to (3.6), yielding Using (3.7) in the eleven equations of motion, they decouple to two sets of three equations of motion and eight equations of motion which depend on fields (h xt ,h xr ,h xy ), and (h xx ,h yy ,h yt ,h yr ,h rr ,h tt ,h rt ,χ), respectively. The reason for this decoupling is that the fields (h xt ,h xr ,h xy ) are odd under parity, hence we refer to them as axial. However, the remaining eight fields are even under the parity transformation x → −x, and we refer to them as polar. Since the perturbation equations are linear in perturbations, they can not couple fields of different symmetry properties [28,29]. The coupled fluctuation equations are lengthy, hence we include them in in the supplementary Mathematica [30] files attached to this paper.
In the rest of this section, we obtain all QNM results using pseudospectral methods. This method turns out to be more efficient in finding QNMs compared to the shooting method described in section 2. We apply the general techniques described well in [31]. More specifically, the recent Mathematica package for finding AdS quasinormal modes [32] has been very useful for generating and checking our code.
Pertaining to the polar modes, there are eleven coupled equations of motion, while there are three for axial modes. We convert each set of equations to a linear algebra statement using a Gauss-Lobatto grid, with 80 to 100 grid points. More specifically, the linear algebra problem is a generalized eigenvalue problem where the complex eigenvalues are the quasinormal mode frequencies we seek to find [32]. A spectral matrix is constructed as a representation of the problem, and Mathematica's Eigenvalues[. . . ] is used to find the generalized eigenvalues, i.e. the QNMs. We note here that the determinant of the relevant matrix vanishes in general, which obstructs inverting that matrix. Hence the treatment as a generalized eigenvalue problem. The procedure for finding convergent quasinormal modes is outlined in appendix B.
In order to determine which horizon is relevant to each sector, we determine the regular singular points of the linear system of differential equations. A perturbative analysis of the coefficients in these equations reveals that the coefficients become singular at a certain radial coordinate value, which is the sound horizon relevant for this sector. Our domain of integration stretches from the AdS-boundary r = 0 to the relevant sound horizon, which we set to r = 1. For the axial modes, this is achieved by fixing r h = 2 1/3 , because the relevant horizon for axial modes is the spin-2 sound horizon at r = r h /2 1/3 . Note, that this fixes the temperature, eq. (2.6) to For the polar modes, the horizon is set to r = 1 by choosing r h = 1, because the relevant horizon for polar modes is the universal horizon, which is also the sound horizon for the spin-0 graviton, r = r h . This fixes the temperature T = 3/(4π √ 1 − c 3 ) = 3 √ 1 + β/(4π) for the polar modes, which is distinct from the axial case temperature by a factor of 2 1/3 .

JHEP10(2019)087
Since the temperature is the only scale we fix here (except for the AdS-radius L fixed using scale symmetries of the equations of motion), our QNM results are still general. 14 With these definitions, all our frequencies will be expressed collectively in units of temperature, as nicely realized by the frequency definition eq. (2.10): for axial and polar sector.

Quasinormal mode results
We find two sets of QNMs coming from the axial and polar sector, respectively, corresponding to the decoupled equations of motion found in the previous section. Plotted in the complex frequency plane, these QNM frequencies are symmetric about the imaginary axis and have negative imaginary values, which shows perturbative stability of the background (in the large parameter region of couplings, λ, β, and momenta k, which we computed). In Einstein Gravity the speeds and horizons of all polar and axial excitations are identical. However, here in Einstein-AEther or equivalently Hořava Gravity, the excitations travel at distinct speeds and hence "see" distinct horizons. Axial QNMs are characterized by speed, √ 1 + β, which is the speed of the spin-2 graviton [11,33] and its sound horizon is located at r s = r h 2 1/3 [11]. The polar sector contains the spin-0 graviton (or khronon) and the remaining components of the metric, all traveling at infinitely large speed as α → 0 [11,21,33], and the corresponding horizon is r h , the universal horizon [11]. In Einstein Gravity, the horizon radius of the black brane solution naturally sets a scale. However, with various horizons, we have a choice where to apply boundary conditions, or to which horizon we normalize other scales, such as our QNM frequencies. As stated above, the computational domain for the axial modes reaches from the boundary r = 0 to the spin-2 sound horizon r s . For polar modes, the relevant horizon is the spin-0 sound horizon r h .

Axial modes
For axial modes in this section the relevant causally connected radial domain stretches from the boundary r = 0 to the spin-2 sound horizon r s = r h /2 1/3 . We choose to work at a fixed temperature by r = r h /2 1/3 = 1. In figure 3, we show axial QNMs corresponding to the axial perturbationsh xt ,h xr , andh xy . The color indicates the magnitude of the dimensionless momentum q in the range 0 ≤ q ≤ 6. For all axial QNMs, we find that their values measured in the dimensionless frequency ν agree numerically with the values of QNMs of a Schwarzschild black brane in Einstein Gravity [28,34]. This empirical evidence implies that the dimensionful frequencies are related as follows ω axial Hořava = 1 + β ω axial Einstein , (3.9) 14 We have checked this explicitly by redefining the radial coordinate, showing that r h disappears from the equations of motion. Hence the QNM frequencies we report are independent of this choice of horizon location in the units we are using. because ν = r h ω/(2 1/3 √ 1 + β), where β = 0 for Einstein Gravity. Therefore, also the dispersion relations ω(k) of this theory are related to those of Einstein Gravity by eq. (3.9).

JHEP10(2019)087
Two distinct types of QNMs appear: hydrodynamic and non-hydrodynamic ones, ν h and ν nh , respectively. The hydrodynamic QNM obeys the defining relation that its frequency vanishes as the momentum vanishes lim q→0 + ν h (q) = 0 . (3.10) Large momentum for this hydrodynamic momentum diffusion mode leads to increasing imaginary frequency indicating increasing dissipation, as expected. For the nonhydrodynamic modes, large momenta are leading to large real parts of the frequencies.
All these axial QNMs are identical to the Hořava QNMs of section 2 up to numerical errors, as indicated by the examples in table 1. The hydrodynamic QNM frequency has vanishing real part and its imaginary part monotonically increases with momentum. At sufficiently small momentum q < 1, our numerically computed frequencies agree well with the analytically predicted momentum diffusion [12] ν h (q) = −iDq 2 + O(q 4 ) ,   brane, eq. (3.4), evaluates to 1 3 (1 + β) −1/2 r h /2 1/3 [12]. The relation (3.11) has been verified numerically, and is also visualized in the fit shown in figure 4. That figure shows the imaginary part of the hydrodynamic mode Figure 5 shows a comparison of the dispersions of the two lowest non-hydrodynamic modes with the hydrodynamic diffusion mode. It is interesting to note that around a momentum of q ≈ 2 the diffusion mode has an imaginary part rapidly growing with momentum, indicating that diffusion modes with larger momentum are rapidly damped. The non-hydrodynamic modes on the other hand display monotonically decreasing imaginary part with increasing momentum. This leads to a crossing between the lowest non-hydrodynamic mode and the diffusion mode at q cross ≈ 1.9; this JHEP10(2019)087 occurs at Im(ν h (q cross )) = Im(ν nh (q cross )). While the late time behavior of the system for momenta q < 1.9 is governed by the diffusion mode, the lowest non-hydrodynamic mode dominates the late time behavior for excitations with q < 1.9. The real part of the nonhydrodynamic quasinormal modes grows linearly for momenta outside the hydrodynamic regime, i.e. with q > 2. All these observations mirror the behavior of relativistic dispersion relations extracted from holography by virtue of the empirical relation eq. (3.9). In table 2 we collect the expansion coefficients parametrizing the dispersion relations of the five lowest quasinormal modes (and for their mirror images over the imaginary frequency axis). We allow the expansion coefficients to be complex valued. The hydrodynamic diffusion mode has only imaginary coefficients in agreement with the requirement that the corresponding transport coefficient, the diffusion constant D, be real valued.

Polar modes
Polar QNMs correspond to the perturbationsh xx ,h yy ,h yt ,h yr ,h rr ,h tt ,h rt , andχ. The eight 15 lowest polar QNMs are displayed in figure 6. Similar to the axial QNMs, the polar QNMs can be scaled and expressed in dimensionless units (2.10). As indicated in figure 1, the QNM spectrum looks different for the two cases λ = 0 and λ = 0. With λ = 0, the polar QNMs are equivalent to black brane QNMs, up to a factor of √ 1 + β, which can be seen by comparison to the Einstein Gravity AdS 4 QNMs computed in [28,34]. Just like the axial QNMs, also the polar QNMs numerically agree with the Einstein Gravity QNMs of an AdS 4 Schwarzschild black brane when measured in the dimensionless frequency ν, so When λ = 0, additional purely dissipative non-hydrodynamic modes are found along the imaginary frequency axis. We refer to these as khronon modes, because they are associated with the fluctuations of the scalar field, χ. This can be confirmed by artificially setting the metric fluctuations to zero and the remaining frequencies found with the spectral method are indeed the same non-hydrodynamic frequencies found when λ = 0 and h µν = 0. Remarkably, the location of the khronon modes and the QNMs associated with the metric appear to be independent of each other, independent of the value of the scalar coupling λ. This implies that eq. (3.12) is true for the QNMs associated with metric fluctuations even at λ = 0.

JHEP10(2019)087
Similar to the axial QNMs, there are both hydrodynamic and non-hydrodynamic modes present in the polar sector. There are two polar hydrodynamic QNMs, ν hs (q), which obey the following dispersion relation up to numerical errors: which can be rewritten in terms of physical frequency and momentum (3.14) Eq. (3.14) agrees exactly with the analytic sound dispersion found in [12] for λ = 0. Our numerical data demonstrates that eq. (3.14) is valid for all values of λ and β.
Unlike the axial hydrodynamic diffusion QNM, the polar hydrodynamic QNM frequencies have a non-zero real part, see figure 7. The mode is propagating with a speed v s = 1/ √ 2. This value is identical to the conformal speed of sound, 1/ √ d − 1, in a relativistic d-dimensional field theory with two spatial dimensions [35]. The imaginary part in eq. (3.13) contains the sound attenuation coefficient Γ = 1/6. This is also consistent JHEP10(2019)087  [35] when factors of the spin-2 sound velocity √ 1 + β are re-instated, as already mentioned in [12]. Up to a rescaling of frequency with the speed factor √ 1 + β, our system has the same dispersion relation as a relativistic theory, see eq. (3.12), regardless of the value for λ. Hence, as expected, our fluid, and in particular the sound attenuation Γ, do not receive the corrections computed in [36]. It is noteworthy, that the sound modes do not grow at large momenta, q > 3, see figures 7 and 8. The latter figure in particular shows that a potential cross-over between imaginary parts of the hydrodynamic sound mode (Polar Hydro Mode) and the lowest non-hydrodynamic mode (1st Non-Hydro Mode) would have to occur at large momentum, q 5. In table 3, we present a parametrization of the dispersion relations of the lowest 14 QNMs at small momentum q < 1. It is suspicious that the higher khronon modes are approximately integer multiples of the lowest khronon frequency, e.g. at vanishing momentum ν khronon ≈ −i 2.381 n ≈ −i 3 2 1/3 n , n = 0, 1, 2, 3, . . . . (3.15) This point is discussed in section 3.5.

Large momentum dispersion relations (eikonal limit)
The QNMs at large momentum are rather difficult to find. The amount of computation time used by the pseudospectral method code increases as momentum increases. At large momentum, q 1, for both the axial and polar modes, the real part of the nonhydrodynamic frequencies are to leading order linear in the momentum. 16 This tendency is already seen in figure 5 and figure 8, and we confirm numerically Re(ν) ≈ q. It was shown in [18] that QNMs of a scalar field in AdS 4 Einstein Gravity at large momentum q at higher mode number n take the form  Table 3. Shown are the expansion coefficients for small momentum q polar modes where ν(q) ≈ 4 m=0 q m ν m . The two hydrodynamic sound modes are collected in the first entry. Empirically we find that the purely imaginary khronon modes are integer multiples of the lowest khronon frequency, regardless of the value of the momentum q, e.g. integer multiples of "ν = −2.381i" at q = 0. where s n is the q −1/5 coefficient for the nth mode, with a phase Arg(s n ) = ±π2/5. Eq. (3.16) was numerically shown to be also true for metric QNMs in AdS 4 [19], approximately even at lower mode numbers n = 1, 2. Analytically, eq. (3.16) was shown for AdS 5 metric QNMs in [20]. The linear behavior ν ∼ q indicates a light-like propagation and is not surprising in a relativistic theory, or in our case in a theory related to a relativistic one by mere scaling factors of the frequency and momentum, given in eq. (2.10). More interesting is the universal correction at order q −1/5 . For smaller values of q, we expect its coefficient to change with momentum, i.e. s n (q), but at large momentum, we expect it to asymptote to a constant independent of momentum lim q→∞ s n (q) → s n . Indeed, this is what we find, as seen in figure 9 for axial modes. All non-hydrodynamic axial QNMs approach each a value s n , labeled by the mode number n = 1, 2, 3, 4. Their phases approach a common value of Arg(s n ) ≈ −1.27, which is approximately −2π/5 as expected. For comparison, we also show the coefficient s diffusion of the hydrodynamic diffusion mode, which does not approach any constant value at large momentum. For the polar modes, figure 10, there are two types of behavior. First, there is the set of modes associated with the metric perturbations, which behave according to eq. (3.16), approaching constant values s n and a common phase Arg(s n ) ≈ −1.27. Second, there are the purely imaginary khronon modes associated with the scalar field, which appear JHEP10(2019)087 to be damped with a power different from q −1/5 , as indicated by the two trajectories not asymptoting to any constant at large momenta in figure 10. Their phase appears to approach −π, which indicates a negative sign for the subleading correction to the linear behavior. A WKB analysis similar to [18,20] may yield an analytic expression for the behavior of these khronon modes at large momentum. However, the equations of motion contain many terms and are coupled to each other, so analytically it is difficult to perform a WKB analysis. For now, we observe that the khronon modes do not display the large momentum behavior expected from QNMs of metric or scalar fields with higher mode numbers n. They rather resemble the behavior of the hydrodynamic diffusion mode shown in figure 9.

"Semi-AEther" field QNMs
As an interesting aside, while analyzing the polar and axial equations of motion, we considered two additional types of perturbations of the AEther field, which preserve its time-like unit normality but not its hypersurface orthogonal condition: One could claim that (3.18) is more "correct" than (3.17), because (3.18) has the correct number of degrees of freedom and preserves the hypersurface orthogonality. In addition to introducing (3.17) and (3.18), we also have to include the time-like unit constraint in the Einstein-AEther action eq. (3.2), we do so by adding it with a Lagrange multiplier λ AE (u 2 + 1). Then, we find λ AE must be set to λ AE = 3c 3 2c 3 r 6 (c 3 − 1) r 6 h + 1 , (3.19) in order to satisfy the unit constraint on the u p µ (x σ ) field. Conveniently (3.19) works for both perturbations, (3.17) and (3.18). We can derive a new set of equations of motion generated by these perturbations (3.17), (3.18). Applying a Fourier transformation and utilizing pseudospectral methods, we have two sets of axial and polar QNMs.
In both the polar and axial sectors, the semi-AEther QNMs found with the new AEther field perturbations (3.17) and (3.18) are numerically indistinguishable from those QNMs found using eq. (3.6). This suggests that the requirement of hypersurface orthogonality does not change the QNMs in our system at hand.
It should be noted that we find the metric fluctuation QNMs (coupled to the scalar khronon) to converge on ten significant figures at the grid size we work with, N grid = 80, 100, see figure 3 and figure 5). At the same grid size, the purely imaginary khronon modes were found to only converge on four significant figures in the case of the (3.17).

Khronon modes
In this section, we discuss the khronon modes. In particular we discuss the question if the khronon modes we find are fake modes or true QNMs. As discussed before, the JHEP10(2019)087 khronon modes are those modes in the polar sector of Einstein-AEther Theory, which have purely imaginary (quasi)eigenfrequencies and are non-hydrodynamic taking on a nonzero frequency value at vanishing momentum. The khronon fluctuation,χ couples to the other six fluctuations in the polar sector,h xx ,h yy ,h yt ,h yr ,h rr ,h tt , andh rt . At vanishing couplings λ = 0 and α = 0, it is possible to analytically map the polar fluctuations to the corresponding Schwarzschild-AdS 4 fluctuations in Einstein Gravity using a field redefinition [12]. At nonzero λ, we find no way of decoupling the system of linear differential equations analytically.
However, when solving the coupled system with pseudospectral methods as a generalized eigenvalue problem, we find that forcing the khronon fluctuation to vanish,χ = 0, does not affect the values of the other polar QNM frequencies (up to numerial errors). In turn, when forcing the metric fluctuations to vanishh xx =h yy =h yt =h yr =h rr =h tt = h rt = 0, we find the khronon eigenfrequencies unaffected (up to numerical errors). 17 Hence, we conclude that the khronon numerically decouples from the metric fluctuations in the polar sector.
We furthermore observe that the khronon eigenfrequencies, or khronon modes, assume purely imaginary values, which are integer multiples of the lowest khronon mode, ν khronon = −i 2.381 n with n = 0, 1, 2, 3, . . . , as stated in eq. (3.15). In fact, if we change our frequency normalization by a factor toν = ν2 1/3 , thenν khronon = −i 3 n. Such integer value solutions are known from various analytical solutions for quasinormal mode frequencies [38]. However, such behavior is also known from various fake modes [15]. The latter are normally revealed because they either do not converge to any value as accuracy is improved (grid size is increased), and such fake modes normally do not change their frequency with changing momentum. However, we find that the khronon modes converge to fixed frequencies, although not as quickly as the metric QNMs. The khronon modes converge to four significant figures while the metric QNMs converge to ten at a grid size of N grid = 100. Furthermore, the khronon modes do move with momentum, as illustrated in figure 10 and figure 8. The latter figure indicates a quadratic rise at small momentum q < 2 and a linear rise thereafter. Hence, the convergence and momentum dependence indicate that the khronon modes are not fake modes.
The khronon is a scalar field, but its equation of motion is not written in the standard Klein-Gordon form. This comes from the fact, that the khronon enters the Einstein AEther action (3.2) through dynamical terms for the vector u µ ∝ ∂ µ φ, which are quadratic in derivatives on u µ . This leads to a fourth order equation of motion for φ or its fluctuatioñ χ, see eq. (3.6). So it is worthwhile analyzing this fourth order equation separately by forcing the metric fluctuations to vanish, see footnote 17. Our near-horizon analysis reveals that this fourth order equation has a regular singular point at the universal horizon r = r h = 1. There are four indicial exponents for the khronon fluctuation near the horizoñ χ ≈ χ 0 (1 − r) α , which are all of the form

JHEP10(2019)087
with a momentum dependent real-valued function f(q), and we recall thatν = 2 1/3 ν. It is a novelty for the indicial exponent to depend on momentum as this is not the case for any QNM equation as far as we know. This momentum dependence can be traced back to the equation being fourth order. In general, f (q) is rather complicated. Let us consider first the case of vanishing momentum, q = 0. Then the four indicial exponents simplify to At this point, we recall that the equations of motion are written in terms of Eddington-Finkelstein coordinates, such that ingoing modes at the horizon are regular, others are singular. None of the modes in eq. (3.21) is regular at the horizon, except for special values of ν. Those special values areν = −i 3 n with n = 0, 1, 2, 3, . . . . Generalizing the momentum to q > 0, we find that regular modes appear at frequency valueŝ ν = −i 3 (n − f(q)) , n = 0, 1, 2, 3, . . . . (3.22) These are the khronon modes found by our pseudospectral method when solving the generalized eigenvalue problem. Eq. (3.22) explains the momentum dependence discussed above. At q = 0, eq. (3.22) also explains the observation that khronon mode frequencies are integer multiples of 3 when written in terms ofν. The numerical data indicates that this property also holds at q = 0, which implies that f(q) ∝ n. Now the question remains, if these khronon modes are to be regarded as true QNMs or as fake modes. 18 One defining property of a QNM is that it vanishes at the AdS-boundary. If the khronon modes are QNMs, then they have to vanish at the AdS-boundary. This can be checked by calculating the eigenvectors associated with the purely imaginary khronon frequencies in question. This analysis is technically very difficult in the full system, because the relevant matrix in the eigenvalue problem is not invertible. Since the khronon modes and the metric QNMs seem to numerically decouple, we hence restrict our eigenvector analysis to the case in which we force all metric fluctuations to vanish as above. In that case, the matrix is invertible, and we confirm that all khronon modes assume non-trivial values along the radial direction and all vanish at the AdS-boundary. 19 As a further test, we consider the large momentum q 1 limit, also called eikonal limit. Results for the two purely imaginary khronon modes were already discussed in section 3.3, and presented in figure 10. The observed phase is Arg(s n ) ≈ −3π/4 and the subleading correction is not of order q −1/5 . This behavior is neither that of a scalar nor that of a metric fluctuation. However, that may not be too surprising because the khronon does not satisfy a simple linearized scalar equation of motion of second order. It rather satisfies a fourth order equation and one should probably conduct a WKB analysis for the vector u µ and compare its numerical behavior with what is expected from that analysis. However, such a treatment is beyond the scope of this work. 18 It is a logical possibility that there exist ingoing solutions for the khronon, which are not regular in the coordinates and field definitions we have chosen here. However, it is not clear to us how to find such modes. 19 As an interesting aside, in this case we also observe that at larger momenta q > 3/2, the khronon spectrum contains both purely imaginary and also propagating modes with real and imaginary part to their QNM frequencies. This seems interesting in light of the observation that only one of these two types appeared [13] at vanishing momentum, depending on the dynamical scaling z and the number of dimensions. However, this behavior is not observed when the metric is allowed to fluctuate. Hence, this appears to be an artifact of the artificial decoupling.

JHEP10(2019)087
Our large momentum analysis ends up being inconclusive. However, the khronon modes (at least when decoupled from the metric fluctuations) satisfy the two defining relations of a QNM: they vanish at the AdS-boundary, and are ingoing at the horizon. Based on this, we decide to interpret the khronon modes as true QNMs which are part of the polar sector of the theory. We speculate that our limitation of α = 0 is forcing part of the khronon dynamics to vanish. This is plausible because the term in the AEther action (3.2) set to zero by α = 0 (equivalent to c 4 = 0) is essentially quadratic in a time derivative of the vector u µ . We speculate that α = 0 would allow the khronon modes to propagate. In that case, however, the analytic background solution is not valid anymore and one has to work with numerical background solutions [11], which is left for future work.

Summary & conclusions
In this paper we have calculated non-relativistic gravitational QNMs on an analytically known asymptotically AdS 4 black brane solution, eq. (2.5) and (3.4), with one of the three coupling constants vanishing, α = 0 [11]. The theory is comprised of two sectors, the parity even polar sector, and the parity odd axial sector. Each sector consists of gravitational fields which travel at a certain speed, either the spin-0 or the spin-2 speed. Correspondingly, the relevant horizon for the axial sector is the spin-2 sound horizon experienced by the spin-2 graviton. While the relevant horizon for the polar sector is the spin-0 sound horizon experienced by the spin-0 graviton. We presented QNMs up to mode number n = 5 for both sectors over a large range of Hořava couplings λ, β, and including large momentum up to q = 50. Our results are summarized in figure 1 and figure 2 for the polar and axial sector, respectively. Equations of motion, and QNM data is collected in the supplementary material attached to this paper.
In this work, we have shown numerically that all Einstein Gravity QNMs are contained in the set of QNMs of Hořava Gravity for any value of λ, β at α = 0, when expressed in appropriate units. At λ = 0, Hořava Gravity has an additional set of purely imaginary QNMs, the khronon modes. The khronon modes seem to numerically decouple from the metric modes and we conjecture an analytic dispersion relation, eq. (3.22), Furthermore, we conjecture an analytic relation between the QNMs of Einstein Gravity and all QNMs of Hořava Gravity at arbitrary λ, β, and at α = 0, except the khronon modes: where r Schwarzschild is the Schwarzschild horizon of a black brane in Einstein Gravity, and r sound is the sound horizon relevant for each sector of QNMs in the analytic Hořava Gravity black brane solution, eq. (3.4). That is the universal horizon r h in the polar sector, and the spin-2 sound horizon r h /2 1/3 in the axial sector. In other words, the QNM frequencies in Einstein Gravity and in the analytic Hořava black brane solution, measured in units of the respective horizon, are equal to each other except for a factor of √ 1 + β.

JHEP10(2019)087
In the axial sector, at any value of λ and β, there is one hydrodynamic diffusion mode and a set of propagating (not overdamped) non-hydrodynamic QNMs, see figure 3. The absence of overdamped (purely imaginary modes) in this sector is in agreement with the claim from [13,14]. The hydrodynamic diffusion mode starts out having quadratic dispersion at small momentum in agreement with the analytic prediction, eq. (3.11) and figure 4. However, then it increases faster in magnitude around q = 1. Around q = 2, it is damped more than the lowest non-hydrodynamic mode (with mode number n = 1). This has been observed before in relativistic theories [26]. At large momentum, q 1, the non-hydrodynamic modes dominate the system since their damping decreases and they are long lived as seen from figure 5. Dispersion relations for the lowest 9 axial QNMs are parametrized to fourth order in momentum in table 2.
In the polar sector, we distinguish two cases, λ = 0 and λ = 0, see figure 6. If λ = 0, then there are no modes associated with the khronon field, only those associated with the metric. Those are two hydrodynamic sound modes and a set of non-hydrodynamic QNMs. The sound modes at small momentum q < 1, see figure 7, agree with the linear propagation and quadratic damping in eq. (3.14), which was derived only at λ = 0. Our analysis shows that this equation holds also at λ = 0. Again, at large momentum, q 1, the system is likely dominated by the non-hydrodynamic modes, because again their damping decreases, as seen in figure 8. Although, this cross-over probably occurs at a much larger momentum than in the axial sector. This is because the polar hydro modes (sound modes) seem to asymptote to a constant value between 0 and −i for large momentum. In addition to that, in the other case, λ = 0, purely imaginary khronon modes are present. In that case, our QNM spectrum contains both overdamped and non-overdamped modes. The overdamped modes are associated with the scalar khronon field fluctuation, while the non-overdamped modes are associated with the metric fluctuations. This is interesting in the context of the claim that only one type, namely overdamped or non-overdamped modes should appear at a given combination of dynamical exponent z and number of dimensions d [13,14]. The latter works consider cases in which a massive scalar probe field does not couple to the metric fluctuations. Hence, it would be interesting in which form the claim needs to be generalized to coupled systems like the one we have studied here. Dispersion relations for the lowest 14 QNMs are parametrized to fourth order in momentum in table 3.
We have also performed a large momentum analysis (eikonal limit) and found to match the analytic expectation based on [18][19][20], see figures 9 and 10. At large momentum, q ≈ 50, all our metric (non-overdamped) QNMs follow the relation ν(q) ≈ q + s n q −1/5 , asymptoting to a constant magnitude for s n and with a universal phase Arg(s n ) ≈ −π2/5. Our overdamped modes, the khronon modes, however, do not show the large momentum dispersion expected from either a scalar or a metric perturbation QNM. As seen in figures 9 and 10, their s n values do not asymptote to constants and their phase is not ±π2/5.
It is interesting to speculate about why the khronon modes decouple from the other modes in the polar sector. The limit of infinite speed is likely the reason for this. A nonzero α leads to finite speed and a horizon for the khronon which will be different from the universal horizon. This case then allows for time derivatives of the khronon in the actions (2.2) and (3.2). Moreover, also the metric modes in the polar sector will travel at a finite speed,

JHEP10(2019)087
which can be distinct from the khronon speed, depending on the Hořava couplings [33]. It will be interesting to see the dynamics and interplay of fields in the polar sector in this case.
From a physical perspective, it is remarkable, that the relativistic relation between shear viscosity and sound attenuation, Γ ∝ η, still holds in this solution of Hořava Gravity. In the latter, gravitational fields (in the axial sector) giving rise to shear modes travel at a different speed than those (in the polar sector) giving rise to sound modes. So one could generally have expected that physical quantities from the polar sector have nothing to do with the axial sector. It would be interesting to check this relation at nonzero α.
In passing, we have verified the equivalence of axial QNMs derived from hypersurface orthogonal Einstein-AEther Theory and those derived from Hořava Gravity. Remarkably, the hypersurface orthogonality does not influence the QNMs in our system, as our semi-AEther results indicate, see section 3.4.
In summary, an obvious, though numerically challenging extension of this work would be the calculation of QNMs for Hořava Gravity with nonzero α coupling, or equivalently Einstein-AEther theory with nonzero c 4 . It will especially be interesting to see the full dynamics of the khronon field unfold. It is expected that QNMs will shift compared to the case α = 0, and that the polar sector should be truly coupled, with one common set of QNMs.
In that setting, one should check how the prediction of purely imaginary modes at d ≤ z +1 needs to be modified for a khronon fluctuations coupling to metric fluctuations. A technical improvement may simplify this computation: possibly gauge invariance could be used to define master fields, reducing the number of fields and field equations that need to be solved.
Relativistic hydrodynamics has been systematically constructed and restricted as an effective field theory over the past years. While Lorentz covariance serves as a fundamental construction principle in that case, it was less clear how to construct non-relativistic hydrodynamics systematically. One way is to start from a relativistic hydrodynamic description, e.g. [39], and then take a non-relativistic limit sending the speed of light to infinity [40], where the choice of the hydrodynamic frame is important [41]. A second way is to identify the non-relativistic data structures directly, as is done in the context of Newton-Cartan geometry [42][43][44]. It has been shown that dynamical Newton-Cartan geometry gives rise to Hořava Gravity [45]. Thus, it would be interesting to use Hořava Gravity as a framework for testing explicitly the proposals for non-relativistic hydrodynamics mentioned above. This may reveal inconsistencies or lead to the discovery of neglected effects.

B Convergence
Following [32], the pseudo-spectral method used in this paper is essentially an eigenvalue solver. Simply put, the expansion of fluctuations into polynomials of nth order, allows n eigenvalues. The number n increases with the size of the chosen grid N grid . Nevertheless, not all of the eigenvalues found are quasinormal frequencies. Some of them are fake modes. In order to determine which eigenvalues are quasinormal modes, we compare two sets of eigenvalues, found using two different grid sizes. The frequencies that do not move more than a set distance ∆ω cutoff in the complex plane are deemed to be convergent and are called a quasinormal frequency. For example, with k = 1 and c 3 = 0 the absolute values of the imaginary part and real part of eigenvalues are displayed in figure 11. For modes of physical interest (black dots in the plot), we can see not much change for increasing grid size.
In contrast to that, the blue dots representing fake mode eigenvalues, change appreciably. However, the number of grid points alone does not always determine convergence. For larger number of grid points used, we find that numerical errors from insufficient numerical precision in our calculation ruins convergence. So a higher numerical precision is required to find higher order modes. We use 80 and 100 grid points and a cutoff value of ∆ω cutoff = 10 −8 to test for convergence. The precision of our arithmetic is ∼ 10 −45 . For polar modes, we found some unphysical modes which did not change up to precision used with changes in momentum. Here, "did not change up to precision used" means that with the digits used, these modes seemed to not change as momentum varied. Physical modes are expected to change with momentum. We discard these non-varying modes due to this pathological behavior.

JHEP10(2019)087
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.