Dispersion relations in non-relativistic two-dimensional materials from quasinormal modes in Ho\v{r}ava Gravity

TL;DR

This paper computes and analyzes the dispersion relations of modes in a non-relativistic strongly coupled 2D quantum field theory using quasinormal modes of a black brane in Hořava Gravity, revealing new analytic and numerical insights.

Contribution

It introduces a method to compute QNMs in Hořava Gravity, including an analytic expression for the decoupled khronon mode and an extension of sound mode dispersion relations.

Findings

Reproduces known momentum diffusion mode numerically.

Extends analytic expression for sound modes across khronon couplings.

Demonstrates perturbative stability over various parameters.

Abstract

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 Ho\v{r}ava Gravity. Ho\v{r}ava 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 Ho\v{r}ava 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.

Tables

Table 1. Table 1: Sample comparison of axial QNM frequencies ν 𝜈 \nu computed from Hořava Gravity using the shooting method and QNMs computed from Einstein-Æther Theory using spectral methods from Sec.  3 . The deviation d 𝑑 d is defined in Eq. ( 13 ), the samples are the hydrodynamic mode at momenta q = 0.1 ,  0.5 ,  1.0 𝑞 0.1 0.5 1.0 q=0.1,\,0.5,\,1.0 , and the lowest non-hydrodynamic mode at q = 0.1 ,  0.5 𝑞 0.1 0.5 q=0.1,\,0.5 .

q	Hořava $\nu$	deviation $d$ from Einstein-Æther [in %]
0.1	$2.040790625396911\times {10}^{-12}-0.003336148565969461i$	$5.92777\times {10}^{-7}+8.05655\times {10}^{-8}i$
0.1	$1.850804328436749-2.6634620817084143i$	$-1.04789\times {10}^{-7}+9.46626\times {10}^{-8}i$
0.5	$6.311321466437316\times {10}^{-13}-0.08515846254765713i$	$2.56115\times {10}^{-10}+7.41127\times {10}^{-10}i$
0.5	$1.8825874031004601-2.654086452911871i$	$-1.04405\times {10}^{-7}+6.75017\times {10}^{-8}i$
1.0	$5.842004644938402\times {10}^{-13}-0.3665129458025562i$	$3.33525\times {10}^{-10}+1.59394\times {10}^{-10}i$

Table 2. Table 2: Shown here are the expansion coefficients for small momentum axial modes, where ν ​ ( q ) ≈ ∑ m = 0 4 q m ​ ν m 𝜈 𝑞 superscript subscript 𝑚 0 4 superscript 𝑞 𝑚 subscript 𝜈 𝑚 \nu(q)\approx\sum_{m=0}^{4}q^{m}\nu_{m} .

${\nu }_0$	${\nu }_1$	${\nu }_2$	${\nu }_3$	${\nu }_4$
$0.$	$0.$	$0.-0.333333i$	$0.$	$0.-0.028111i$
$\pm 1.84942-2.66385i$	$0.$	$\pm 0.138538+0.039062i$	$0.$	$\mp 0.023605-0.000358i$
$\pm 3.16126-4.91642i$	$0.$	$\pm 0.105299+0.027848i$	$0.$	$\mp 0.01297-0.000165i$
$\pm 4.46435-7.16754i$	$0.$	$\pm 0.08816+0.023133i$	$0.$	$\mp 0.008838-0.000102i$
$\pm 5.76525-9.41808i$	$0.$	$\pm 0.077319+0.020279i$	$0.$	$\mp 0.006816-0.000059i$

Table 3. Table 3: Shown are the expansion coefficients for small momentum q 𝑞 q polar modes where ν ​ ( q ) ≈ ∑ m = 0 4 q m ​ ν m 𝜈 𝑞 superscript subscript 𝑚 0 4 superscript 𝑞 𝑚 subscript 𝜈 𝑚 \nu(q)\approx\sum_{m=0}^{4}q^{m}\nu_{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 𝑞 q , e.g. integer multiples of “ ν = − 2.381 ​ i 𝜈 2.381 𝑖 \nu=-2.381i ” at q = 0 𝑞 0 q=0 .

${\nu }_0$	${\nu }_1$	${\nu }_2$	${\nu }_3$	${\nu }_4$
$0.+0.00006i$	$\mp 0.70163$	$\mp 0.07719-0.16993i$	$\pm 0.14775$	$\mp 0.10496+0.02857i$
$0.-2.381i$	$0.$	$0.$	$0.$	$0.$
$1.8495-2.66401i$	$0.$	$0.1374+0.04113i$	$0.$	$0.02127-0.01031i$
$-1.84948-2.66401i$	$0.$	$-0.13769+0.04113i$	$0.$	$-0.021-0.01031i$
$0.-4.7619i$	$0.$	$0.-0.20987i$	$0.$	$0.+0.00993i$
$3.16133-4.91639i$	$0.$	$0.10448+0.02638i$	$0.$	$0.01176+0.0062i$
$-3.16137-4.91639i$	$0.$	$-0.10406+0.02638i$	$0.$	$-0.0122+0.0062i$
$4.46443-7.16784i$	$0.$	$0.08739+0.02837i$	$0.$	$0.00831-0.01216i$
$-4.46445-7.16784i$	$0.$	$-0.08719+0.02837i$	$0.$	$-0.00852-0.01216i$
$5.76501-9.41813i$	$0.$	$0.07778+0.01646i$	$0.$	$0.00564+0.01621i$
$-5.76504-9.41822i$	$0.$	$-0.0776+0.02054i$	$0.$	$-0.00542-0.00365i$
$7.06506-11.6683i$	$0.$	$0.06959+0.01718i$	$0.$	$0.00521+0.0028i$
$-7.06509-11.6683i$	$0.$	$-0.06951+0.01718i$	$0.$	$-0.00486+0.0028i$

Table 4. Table 4: A Sample of the Einstein-Æther/Hořava Axial Quasinormal-Modes, found using pseudospectral methods.

Momentum ($q$)	Frequency ($\nu$)	$n^{th}Mode$
$0.1$	$-0.0033361486i$	0
$2.5$	$-6.5109332704i$	0
$0.1$	$\pm 1.8508043279-2.6634620863i$	1
$2.5$	$\pm 2.5648819977-1.7693681441i$	1
$20.0$	$\pm 20.2726721388-0.9160839493i$	1
$0.1$	$\pm 3.1623087258-4.9161394118i$	2
$2.5$	$\pm 3.3930208613-4.2840047943i$	2
$20.0$	$\pm 20.794701895-2.4963744671i$	2

Table 5. Table 5: A Sample of the Einstein-Æther Polar Quasinormal-Modes, found using pseudospectral methods.

Momentum ($q$)	Frequency ($\nu$)	$n^{th}Mode$
0.25198	$\pm 4.227043-0.006645i$	0
2.51984	$\pm 2.144945-0.454586i$	0
7.55953	$\pm 7.570777-7.55953i$	0
0.25198	$\pm 2.331225-3.355934i$	1
2.51984	$\pm 2.941333-4.227044i$	1
7.55953	$\pm 8.473500-3.621822i$	1
0.25198	$\pm 5.022401-7.803214i$	2
2.51984	$\pm 5.535820-7.707342i$	2
7.55953	$\pm 9.780225-7.023156i$	2
0.25198	$-3.788146i$	1 (Purely Dissipative )
4.78770	$-6.099892i$	1 (Purely Dissipative )

Equations

$$d{s}^2=g_{XY}d{x}^{X}d{x}^Y=-N^2{dt}^2+G_{IJ}({dx}^I+N^{I}dt)({dx}^J+N^{J}dt).$$

$$S_{\text{Horˇava}}=\int \frac{d^{4}x\sqrt{|G|}}{16\pi G_H}\left(K_{IJ}{K}^{IJ}-(1+\lambda )K^2+(1+\beta )(R-2\Lambda )+\alpha \frac{({\nabla }_{I}N)({\nabla }^{I}N)}{N^2}\right),$$

$$K_{IJ}\equiv \frac{1}{2N}({\partial }_t{G}_{IJ}-{\nabla }_I{N}_J-{\nabla }_J{N}_I).$$

$$\beta >-1\text{and}0\le \alpha \le 2(1+\beta )\text{and}\lambda \ge 0\text{or}\lambda \le -\frac{2}{3}.$$

$$G_{IJ}=\left(\begin{array}{ccc}{\left(\frac{r_h^3}{r(r_h^3-r^3)}\right)}^2& 0& 0\\ 0& \frac{1}{r^2}& 0\\ 0& 0& \frac{1}{r^2}\end{array}\right)N=\frac{1}{r}\left(1-\frac{r^3}{r_h^3}\right)N_I=\left(\frac{r\sqrt{1+\beta }}{r_h^3-r^3}\right),$$

$$T=\frac{3\sqrt{1+\beta }}{4\pi r_h}.$$

$$g_{XY}^p=g_{XY}+ϵh_{XY}(t,r,x^I),$$

$$h_{XY}(t,r,x^I,r)=\int d\omega dk{e}^{\text{boldmath}i(ky-\omega t)}\frac{{\tilde{h}}_{XY}(r;\omega ,k)}{r^2}.$$

$$\bigg{[}q^{4}z^{2}(-2+z^{3})^{2}(-1+z^{3})+2{\nu}^{2}\big{[}-2{\nu}^{2}z^{2}-i\nu z^{4}(-5+z^{3})+(-2+z^{3})^{2}(1+2z^{3})\big{]}\\ +q^{2}\big{[}{\nu}^{2}z^{2}(8-8z^{3}+z^{6})-2(2-3z^{3}+z^{6})^{2}+i\nu z^{4}(-10+6z^{3}+z^{6})\big{]}\bigg{]}\psi(z)\\ +z(-2+z^{3})\big{[}2q^{2}(-1+z^{3})(2+z^{3}(-3-i\nu z+z^{3}))\\ +{\nu}^{2}(4+z^{3}(-12-2i\nu z+z^{3}))\big{]}\psi^{\prime}(z)\\ z^{2}(-2+z^{3})^{2}(-1+z^{3})\big{[}{\nu}^{2}+q^{2}(-1+z^{3})\big{]}\psi^{\prime\prime}(z)=0\,,$$

$$z\equiv \frac{2^{1/3}r}{r_h},q\equiv \frac{r_h}{2^{1/3}}k,\nu \equiv \frac{r_h}{2^{1/3}\sqrt{1+\beta }}\omega .$$

$$\psi (z;\nu ,q)\equiv (1-z{)}^{\alpha }F(z;\nu ,q)=(1-z{)}^{\alpha }\sum _{n=0}^{\infty }{f}_n(\nu ,q)(1-z{)}^n,$$

$${\omega }_{\text{axial QNM}}^{\text{Horˇava}}=\sqrt{1+\beta }{\omega }_{\text{axial QNM}}^{\text{Einstein}}.$$

$$d=2\frac{{\nu }_H-{\nu }_{\text{AE}}}{{\nu }_H+{\nu }_{\text{AE}}}\times 100\%.$$

$${\nu }_{\text{QNM}}^{\text{Horˇava}}={\nu }_{\text{QNM}}^{\text{hypersurface-orthogonal Einstein-Æther}},$$

$$u_{\mu }=\frac{{\partial }_{\mu }\varphi }{\sqrt{-({\partial }_{\nu }\varphi )({\partial }^{\nu }\varphi )}},$$

$$S_{\text{\AE ther}}=\frac{1}{4\pi G_{ae}}\int d^{4}x\sqrt{-g}\bigg{(}R-2\Lambda+c_{4}(u^{\mu}{\nabla}_{\mu}u^{\nu})(u^{\sigma}{\nabla}_{\sigma}u_{\nu})-c_{3}({\nabla}_{\mu}u^{\nu})({\nabla}_{\nu}u^{\mu})\\ -c_{2}({\nabla}_{\mu}u^{\mu})^{2}-c_{1}({\nabla}^{\mu}u^{\nu})({\nabla}_{\mu}u_{\nu})\bigg{)}\,,$$

$$\frac{G_H}{G_{ae}}=1+\beta =\frac{1}{1-c_3},1+\lambda =\frac{1+c_2}{1-c_3},\alpha =\frac{c_4}{1-c_3}.$$

$$g_{\mu \nu }$$

$$u_{\mu }$$

$$e(r)$$

$$\phi(r,\nu)=\Bigg{(}\int^{r}\frac{2a(\rho)^{2}f(\rho)}{a(\rho)^{2}e(\rho)+f(\rho)^{2}}\,d\rho\Bigg{)}-\nu\,.$$

$$g_{\mu \nu }^p$$

$${\varphi }^p$$

$$h_{\mu \nu }(x^{\sigma })$$

$$\chi (x^{\sigma })$$

$$\nu =\left(\frac{3}{4\pi 2^{1/3}}\right)\frac{\omega }{T},$$

$${\omega }_{\text{Horˇava}}^{\text{axial}}=\sqrt{1+\beta }{\omega }_{\text{Einstein}}^{\text{axial}},$$

$$\underset{q\to 0^{+}}{\lim }{\nu }_h(q)=0.$$

$${\nu }_h(q)=-iD{q}^2+\mathcal{O}(q^4),$$

$${\omega }_{\text{Einstein-Æther}}^{\text{polar}}=\sqrt{1+\beta }{\omega }_{\text{Einstein}}^{\text{polar}}(\lambda =0).$$

$$\underset{q\to 0^{+}}{\lim }{\nu }_{hs}(q)=\pm v_{s}q-i\Gamma q^2=\pm \frac{1}{\sqrt{2}}q-i\frac{1}{6}{q}^2+\mathcal{O}(q^3),$$

$$\underset{k\to 0^{+}}{\lim }{\omega }_{hs}(q)=\pm \frac{1}{\sqrt{2}}\sqrt{1+\beta }k-i\frac{1}{6}\frac{r_h}{2^{1/3}}\sqrt{1+\beta }{k}^2+\mathcal{O}(k^3).$$

Full text

1 Introduction

Right after the discovery of the gauge/gravity correspondence Maldacena:1997re , it was applied in parallel with hydrodynamics to describe strongly coupled fluids Policastro2002 ; Policastro:2001yc ; Bhattacharyya:2008jc , 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 Horava:2009uw ; Horava:2009vy ; Anderson:2011bj ; horava2009 and non-relativistic quantum field theories Janiszewski:2012nf ; Janiszewski2013 . An analytical Hořava black brane solution was found Janiszewski2015 , and its fluctuations were studied analytically in the hydrodynamic limit Davison2016 . 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 Hořava coupling constant to zero, $\lambda=0$. Two hydrodynamic sound modes were found as expected. But this failure of the map at $\lambda\neq 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\to\kappa^{z}t,\,x^{I}\to\kappa x^{I}$, with constant $\kappa$ 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 Sec. 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\leq z+1$ at vanishing momentum Sybesma2015 ; Gursoy:2016tgf . 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 Sybesma2015 ; Gursoy:2016tgf . These references Sybesma2015 ; Gursoy:2016tgf 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 Janiszewski:2015ura 111These black brane solutions are smoothly connected to solutions with all possible two-derivative Hořava couplings non-zero Janiszewski:2015ura ., 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 Janiszewski2015 . 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 Davison2016 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 Janiszewski2015 .

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 Sec. 3, the fluctuation equations for Hořava fields are very complicated. Hence, we use the equivalence between Hořava Gravity and Einstein-Æther Theory when the scalar field, the khronon, is required to be hypersurface-orthogonal BHATTACHARYYA2013 ; Jacobson2000 . 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 Sec. 3.5.

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 Festuccia:2008zx ; Morgan:2009vg ; Fuini2016 , are a useful tool to check numerical accuracy and the structure of the quasinormal mode spectrum, see Sec. 3.3. Our results are discussed in Sec. 2.5 and 3.2, but visually summarized in Fig. 2 and Fig. 2. The parameter space of Hořava coupling constants $\lambda$ and $\beta$ is shown. The third coupling is set to zero in this work, $\alpha=0$. A summary of results and our conclusions are found in Sec. 4. Equations of motion, and QNM data is collected in four ancillary files.

The two-dimensional materials we have in mind for application of our results could be thin films (e.g. NbSi Hartnoll:2007ih ) or layered structures such as those found in high critical temperature cuprates, see e.g. Hartnoll:2007ih ; Blauvelt:2017koq and more generally Hartnoll:2009sz ; McGreevy:2009xe . We assume that such materials —considered in an appropriate regime of their defining quantities— can be described by non-relativistic hydrodynamics Hartnoll:2007ih . 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.

2 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 Janiszewski2015 . 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 polar222That previous analysis specialized to the case of vanishing coupling $\lambda=0$ for polar QNMs. In this work, we lift this severe restriction and analyze QNMs for nonzero $\lambda$. transport quantities have been found analytically in Davison2016 333In Davison2016 , our axial modes are there called vector modes, and our polar modes are there called scalar modes., which will serve as a check of our numerical solutions.

2.1 Action, coupling constants & field content

We will specialize to a particular low energy solution of classical Hořava Gravity Janiszewski2015 . 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}$444Capitalized roman letter indices $I,J,K,\ldots$ are non-temporal indices (ie not $t$). Capitalized roman letter indices $X,Y,Z$ are bulk spacetime indices.; where we will be using the mostly positive metric convention555The mostly positive convention (-,+,+,+).. In terms of spacetime coordinates $x^{X}=\{t,r,x,y\}$ the line element then is given by

[TABLE]

In (3+1) dimensions the analytically known metric solution Janiszewski2015 satisfies the equations of motion generated by the variation of the following action of Hořava Gravity:

[TABLE]

with

[TABLE]

$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 $\nabla_{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,$(\alpha,\beta,\lambda)$, must obey the following inequalities BHATTACHARYYA2013 ; Griffin2013 :

[TABLE]

These allowed regions of parameter space are represented by white surfaces plus the blue line along the $\beta$-axis in the ($\beta$,$\lambda$)-parameter plots Fig. 2 and Fig. 2; forbidden regions are colored red.

2.2 Hořava black brane background solution

Taking $\alpha=0$ and $\Lambda=-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 Janiszewski2015 . $G_{IJ}$, $N^{I}$, and $N$ are known analytically and given by666The radial coordinate $r$, spatial momentum $k$, and frequency $\omega$ carry units of $length$, $length^{-1}$, and $length^{-1}$ respectively.

[TABLE]

which notably is independent from the coupling $\lambda$. The AdS radius is $L=1$, the time coordinate is $t\equiv x^{0}$, spatial coordinates $x\equiv x^{2}$, $y\equiv x^{3}$, and the radial coordinate is $r\equiv 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 horizon777We 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, $\phi^{\prime\prime}(r)+b(r)/(r-r_{c})\phi^{\prime}(r)+b(r)/(r-r_{c})^{2}\phi(r)=0$ has a sound horizon at $r=r_{c}$, if $b(r)$ and $c(r)$ can each be expanded in a Taylor series. 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^{\frac{1}{3}}$, which is a trapping surface for the spin-2 graviton. There is also the Killing horizon located at $r=r_{k}$ Davison2016 . The temperature is determined by the universal horizon Janiszewski2015 and is given by

[TABLE]

In general, solutions to Hořava Gravity can exhibit Lifshitz scaling symmetry under the scaling with a constant $\kappa$: $t\to\kappa^{z}t$,   $x^{I}\to\kappa x^{I},\,r\to\kappa r$ with the dynamical exponent $z$, leading to the asymptotic scaling $ds^{2}\sim\frac{dt^{2}}{r^{2z}}+\frac{{dx^{I}}^{2}+dr^{2}}{r^{2}}$ Griffin2013 . The metric solution in Eq. (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-Æther Theory with a particular constraint (hypersurface orthogonality) Davison2016 ; BHATTACHARYYA2013 , as we will see below in Sec. 3. Hence, the solution (5), like any typical solution to Hořava Gravity, is only invariant under those diffeomorphisms preserving the time-foliation.

2.3 Hořava black brane perturbations

One may perturb the metric (1) around a background value of (5) to linear order by a metric perturbation $h_{XY}$

[TABLE]

where $\epsilon\ll 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) and background metric (5) have a translational symmetry in $x,y$-directions (5), we make the standard plane wave ansatz of Fourier expanding into momentum space modes ${\tilde{h}}_{XY}(r;\omega,\vec{k})$.888Recall that $x\equiv x^{2}$ and $y\equiv x^{3}$. 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

[TABLE]

Once (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 parity (axial) and a set of 3 equations for metric components which are even under parity (polar) Davison2016 . For the QNMs in the Hořava case we only concern ourselves with the set of 3 equations of motion for the axial fields ${\tilde{h}}_{xy}$, ${\tilde{h}}_{xt}$, and ${\tilde{h}}_{xr}$, which are odd under parity. After a radial diffeomorphism, ${\tilde{h}}_{xr}(r;\omega,k)$ can be set to vanish. A linear combination of these fields turns out to be a gauge invariant master field, $\psi(r;\omega,k)=(\omega h_{xy}+kh_{tx})$, and obeys the following single equation of motion999Here $\nu$ and $k$ dependence in $\psi$ is suppressed. found in Davison2016

[TABLE]

where the following variables have been rescaled to be dimensionless

[TABLE]

Eq. (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.

2.4 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 (10), the sound horizon is located at $z=1$, where the equation of motion (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

[TABLE]

separating the regular part $F$ and irregular part of $\psi(z;\nu,q)$ from each other. Solving the equation (9) with the ansatz (11), one finds two $\alpha$’s that satisfy Eq. (9) at the first non-zero order in the near-horizon expansion: $\alpha=-1$ or $\alpha=-\frac{2i\nu}{3}$. We choose the latter which corresponds to the ingoing mode. At the other expansion orders one can recursively solve for the $f_{n}(\nu,q)$ coefficients. As usual, the series in Eq. (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 (9) with Mathematica’s NDSolve function on the restricted computational domain of101010z=1 is the location of the sound horizon. $z\in[{dr}_{b},1-{dr}_{h}]$ where ${dr}_{b},{dr}_{h}\ll 1$. It as been shown that for exceedingly small values of ${dr}_{h}$ the (11) ansatz fluctuates rapidly and can create large numerical errors Kaminski:2009ce , which guides our choice of $dr_{h}$ here. The boundary conditions, $\psi(1-{dr}_{h})$ and $\psi^{\prime}(1-{dr}_{h})$, are then provided by the horizon expansion (11). For an arbitrary value of $\nu$ and $q$, this solution is not a quasinormal mode, i.e. $\psi({dr}_{b})\neq 0$. We use Mathematica’s FindRoot to find $\nu$ such that $\psi({dr}_{b};q)=0$. When compared to the pseudospectral method’s modes (Fig. 3 and Fig. 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.

2.5 Hořava Gravity axial QNMs

It turns out that the axial QNMs $\nu$ found with the shooting method are numerically equal the QNMs one would find for an asymptotically AdS4 Schwarzschild black brane within Einstein Gravity, i.e. $\nu^{\text{Ho\v{r}ava}}_{\text{axial QNM}}=\nu^{\text{Einstein}}_{\text{axial QNM}}$. However, the QNM frequencies $\nu$ are scaled by a factor of $\sqrt{1+\beta}$ and $\beta$ vanishes in GR. So, up to numerical errors we empirically find the relation

[TABLE]

This particular $\beta$ factor is in fact the speed of the spin-2 graviton Davison2016 . At $\beta=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 Davison2016 .

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.111111This is due to an apparent pole at $r=\frac{r_{h}}{2^{1/3}}=r_{s}$. This still happened despite the polar mode equations of motion not having any factors of $(r-r_{s})$ in the perturbation coefficients, which would have indicated regular/irregular poles at $r_{s}$. 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-Æther Gravity, using a different technique, namely pseudospectral methods. The equivalence of these two theories holds under the constraint of hypersurface orthogonality, discussed in Sec. 3. Hence, QNMs found in the two theories are expected to be identical121212Only a smaller subset of the frequencies found with spectral methods were also found via the shooting method.. 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-Æther 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-Æther Theory QNMs, for that calculation see Sec. 3, is given by

[TABLE]

Up to numerical errors which are at worst on the order of $10^{-6}$ percent, we find agreement between QNMs of hypersurface-orthogonal Einstein-Æther Theory and Hořava Gravity in the five axial QNMs we have checked for many values of momentum $q$, and report only examples in table 1. This serves as a check of the equality of the QNM spectra

[TABLE]

which will be discussed in the next section.

3 Einstein-Æther Theory

Einstein-Æther Theory differs from General Relativity by the addition of a scalar matter field $\phi$, which is termed Æther field. We need a particular Einstein-Æther Theory, namely one that is equivalent to Hořava gravity. This requires the Æther field to define a timelike vector field which is hypersurface orthogonal. That goal is achieved by the following definition of the timelike vector

[TABLE]

which is a timelike hypersurface orthonormal unit vector matter field Jacobson2000 ; Bhattacharyya:2008jc . This vector field also breaks diffeomorphism invariance down to foliation preserving diffeomorphisms. In this context, the scalar $\phi$ is often referred to as khronon, determining the time foliation. Lower case Greek indices $\mu,\nu$ denote spacetime indices. As coordinates, we choose $x^{\mu}$ 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-Æther action has been constructed Jacobson2000 , which includes four quadratic derivative terms of $u_{\mu}$,131313We omit here the constraint term, $\lambda_{\text{\AE}}(u^{2}+1)$, since Eq. (15) incorporates this unit constraint by construction.

[TABLE]

where $g$ is the determinant of the metric $g_{\mu\nu}$, and $R$ is the Ricci Scalar of the metric $g_{\mu\nu}$. $c_{1}$ is redundant once we constrain $u$ to be hypersurface orthogonal by construction. This allows us to rearrange (16) 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}$ BHATTACHARYYA2013 , and we set $c_{1}=0$ without loss of generality. This equates Einstein-Æther theory to low energy Hořava Gravity after the identification $-N=\delta^{\mu}_{0}u_{\mu}$, and with the following coupling constants for the respective theories Barausse2011

[TABLE]

Analogously to our perturbative treatment in the previous section around Eq. (1), also here we will investigate linear perturbations around the black brane solution with $\alpha=0$ and $\Lambda=-3$ Janiszewski2015 . Our coordinates are similar to Eddington-Finkelstein coordinates Janiszewski2015 , as seen in the metric

[TABLE]

where $r_{h}$ is the radius of the universal horizon. The the sign choice on $f(r)$ and in $u_{\mu}$ corresponds to the choice of infalling (lower signs) or outgoing (upper signs) Eddington-Finkelstein-like coordinates Janiszewski2015 . 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. (18) reduces to the Schwarzschild $AdS_{4}$ metric with the Schwarzschild horizon located at $r=r_{\text{Schwarzschild}}=r_{h}/2^{1/3}$ when all remaining Æther couplings are set to zero, $c_{3}=0=c_{2}$, and we choose the vector field to vanish $u_{\mu}=0$. As a side note, remarkably, one can express $\phi$ explicitly by integrating Eq. (15) to obtain

[TABLE]

3.1 Einstein-Æther black brane perturbations

Similar to how we perturbed the Hořava black brane (7) we perturb the metric (18) by adding a “small” linear term where $\epsilon\ll 1$, we choose

[TABLE]

where $\chi(x^{\sigma})$ is a scalar field. Since $\phi^{p}$ is still a scalar field after the perturbation, replacing $\phi\to\phi^{p}$ in Eq.(15) ensures hypersurface orthogonality and normalization of the now perturbed $u$ vector. The $h_{\mu\nu}$ fields and $\chi$ obey eleven coupled linear equations. A Fourier transformation similar to (8), is applied to (20), yielding

[TABLE]

Using (21) 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 ($\tilde{h}_{xt}$, $\tilde{h}_{xr}$, $\tilde{h}_{xy}$), and ($\tilde{h}_{xx}$, $\tilde{h}_{yy}$, $\tilde{h}_{yt}$, $\tilde{h}_{yr}$, $\tilde{h}_{rr}$, $\tilde{h}_{tt},\tilde{h}_{rt}$, $\tilde{\chi}$), respectively. The reason for this decoupling is that the fields ($\tilde{h}_{xt}$, $\tilde{h}_{xr}$, $\tilde{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\to-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 Miranda2008 ; Kovtun:2005ev . The coupled fluctuation equations are lengthy, hence we include them in ancillary Mathematica Mathematica10 files with this submission.

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 Sec. 2. We apply the general techniques described well in Boyd2000 . More specifically, the recent Mathematica package for finding AdS quasinormal modes Jansen2017 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 Jansen2017 . A spectral matrix is constructed as a representation of the problem, and Mathematica’s $\mathsf{Eigenvalues[\dots]}$ 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. (6) to $T=3/(4\pi 2^{1/3}\sqrt{1-c_{3}})=3\sqrt{1+\beta}/(4\pi 2^{1/3})$. 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\pi\sqrt{1-c_{3}})=3\sqrt{1+\beta}/(4\pi)$ for the polar modes, which is distinct from the axial case temperature by a factor of $2^{1/3}$. 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.141414We 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. With these definitions, all our frequencies will be expressed collectively in units of temperature, as nicely realized by the frequency definition Eq. (10):

[TABLE]

for axial and polar sector.

3.2 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, $\lambda,\,\beta$, 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-Æther or equivalently Hořava Gravity, the excitations travel at distinct speeds and hence “see” distinct horizons. Axial QNMs are characterized by speed,$\sqrt{1+\beta}$, which is the speed of the spin-2 graviton Jacobson2004 ; Janiszewski2015 and its sound horizon is located at $r_{s}=\frac{r_{h}}{2^{1/3}}$ Janiszewski2015 . The polar sector contains the spin-0 graviton (or khronon) and the remaining components of the metric, all traveling at infinitely large speed as $\alpha\rightarrow 0$ Jacobson2004 ; Griffin2013 ; Janiszewski2015 , and the corresponding horizon is $r_{h}$, the universal horizon Janiszewski2015 . 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}$.

3.2.1 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 Fig. 3, we show axial QNMs corresponding to the axial perturbations $\tilde{h}_{xt}$, $\tilde{h}_{xr}$, and $\tilde{h}_{xy}$. The color indicates the magnitude of the dimensionless momentum $q$ in the range $0\leq q\leq 6$. For all axial QNMs, we find that their values measured in the dimensionless frequency $\nu$ agree numerically with the values of QNMs of a Schwarzschild black brane in Einstein Gravity Morgan:2009pn ; Miranda2008 . This empirical evidence implies that the dimensionful frequencies are related as follows

[TABLE]

because $\nu=r_{h}\omega/(2^{1/3}\sqrt{1+\beta})$, where $\beta=0$ for Einstein Gravity. Therefore, also the dispersion relations $\omega(k)$ of this theory are related to those of Einstein Gravity by Eq. (23).

Two distinct types of QNMs appear: hydrodynamic and non-hydrodynamic ones, $\nu_{h}$ and $\nu_{nh}$, respectively. The hydrodynamic QNM obeys the defining relation that its frequency vanishes as the momentum vanishes

[TABLE]

Large momentum for this hydrodynamic momentum diffusion mode leads to increasing imaginary frequency indicating increasing dissipation, as expected. For the non-hydrodynamic 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 Davison2016

[TABLE]

up to corrections of order $\mathcal{O}(q^{3})$, and with the diffusion coefficient $D=1/3$. That value is consistent with the relativistic value $D=\eta/(\epsilon+P)$ Herzog:2003ke , which for the Hořava black brane, Eq. (18), evaluates to $\frac{1}{3}(1+\beta)^{-1/2}r_{h}/2^{1/3}$ Davison2016 . The relation (25) has been verified numerically, and is also visualized in the fit shown in Fig. 4. That figure shows the imaginary part of the hydrodynamic mode

Fig. 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\approx 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_{\text{cross}}\approx 1.9$; this occurs at $\operatorname{Im}(\nu_{h}(q_{\text{cross}}))=\operatorname{Im}(\nu_{nh}(q_{\text{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 non-hydrodynamic 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. (23). 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.

3.2.2 Polar modes

Polar QNMs correspond to the perturbations $\tilde{h}_{xx}$, $\tilde{h}_{yy}$, $\tilde{h}_{yt}$, $\tilde{h}_{yr}$, $\tilde{h}_{rr}$, $\tilde{h}_{tt},\tilde{h}_{rt}$, and $\tilde{\chi}$. The eight 151515Counting includes mirror QNMs with opposite sign of the real part of the frequency. lowest polar QNMs are displayed in Fig. 6. Similar to the axial QNMs, the polar QNMs can be scaled and expressed in dimensionless units (10). As indicated in Fig. 2, the QNM spectrum looks different for the two cases $\lambda=0$ and $\lambda\neq 0$. With $\lambda=0$, the polar QNMs are equivalent to black brane QNMs, up to a factor of $\sqrt{1+\beta}$, which can be seen by comparison to the Einstein Gravity $AdS_{4}$ QNMs computed in Morgan:2009pn ; Miranda2008 . 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 $\nu$, so

[TABLE]

When $\lambda\neq 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, $\chi$. 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 $\lambda\neq 0$ and $h_{\mu\nu}=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 $\lambda$. This implies that Eq. (26) is true for the QNMs associated with metric fluctuations even at $\lambda\neq 0$.

Similar to the axial QNMs, there are both hydrodynamic and non-hydrodynamic modes present in the polar sector. There are two polar hydrodynamic QNMs, $\nu_{hs}(q)$, which obey the following dispersion relation up to numerical errors:

[TABLE]

which can be rewritten in terms of physical frequency and momentum

[TABLE]

Eq. (28) agrees exactly with the analytic sound dispersion found in Davison2016 for $\lambda=0$. Our numerical data demonstrates that Eq. (28) is valid for all values of $\lambda$ and $\beta$.

Unlike the axial hydrodynamic diffusion QNM, the polar hydrodynamic QNM frequencies have a non-zero real part, see Fig. 7. The mode is propagating with a speed $v_{s}=1/\sqrt{2}$. This value is identical to the conformal speed of sound, $1/\sqrt{d-1}$, in a relativistic $d$-dimensional field theory with two spatial dimensions Herzog:2003ke . The imaginary part in Eq. (27) contains the sound attenuation coefficient $\Gamma=1/6$. This is also consistent with the known relativistic formula $\Gamma=\frac{d-2}{d-1}\frac{\eta}{\epsilon+P}$ Herzog:2003ke when factors of the spin-2 sound velocity $\sqrt{1+\beta}$ are re-instated, as already mentioned in Davison2016 . Up to a rescaling of frequency with the speed factor $\sqrt{1+\beta}$, our system has the same dispersion relation as a relativistic theory, see Eq. (26), regardless of the value for $\lambda$. Hence, as expected, our fluid, and in particular the sound attenuation $\Gamma$, do not receive the corrections computed in deBoer:2017abi . It is noteworthy, that the sound modes do not grow at large momenta, $q>3$, see Fig. 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\gg 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

[TABLE]

This point is discussed in Sec. 3.5.

3.3 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\gg 1$, for both the axial and polar modes, the real part of the non-hydrodynamic frequencies are to leading order linear in the momentum.161616For the purely dissipative polar non-hydro modes, the real part obviously vanishes. This tendency is already seen in Fig. 5 and Fig. 8, and we confirm numerically $Re(\nu)\approx q$. It was shown in Festuccia:2008zx that QNMs of a scalar field in $AdS_{4}$ Einstein Gravity at large momentum $q$ at higher mode number $n$ take the form

[TABLE]

where $s_{n}$ is the $q^{-1/5}$ coefficient for the $n$th mode, with a phase $Arg(s_{n})=\pm\pi 2/5$. Eq. (30) was numerically shown to be also true for metric QNMs in $AdS_{4}$ Morgan:2009vg , approximately even at lower mode numbers $n=1,2$. Analytically, Eq. (30) was shown for $AdS_{5}$ metric QNMs in Fuini2016 . The linear behavior $\nu\sim 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. (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\limits_{q\to\infty}s_{n}(q)\to s_{n}$. Indeed, this is what we find, as seen in Fig. 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})\approx-1.27$, which is approximately $-2\pi/5$ as expected. For comparison, we also show the coefficient $s_{\text{diffusion}}$ of the hydrodynamic diffusion mode, which does not approach any constant value at large momentum.

For the polar modes, Fig. 10, there are two types of behavior. First, there is the set of modes associated with the metric perturbations, which behave according to Eq. (30), approaching constant values $s_{n}$ and a common phase $Arg(s_{n})\approx-1.27$. Second, there are the purely imaginary khronon modes associated with the scalar field, which appear 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 Fig. 10. Their phase appears to approach $-\pi$, which indicates a negative sign for the subleading correction to the linear behavior. A WKB analysis similar to Festuccia:2008zx ; Fuini2016 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 Fig. 9.

3.4 “Semi-Æther” field QNMs

As an interesting aside, while analyzing the polar and axial equations of motion, we considered two additional types of perturbations of the Æther field, which preserve its time-like unit normality but not its hypersurface orthogonal condition:

[TABLE]

One could claim that (32) is more “correct” than (31), because (32) has the correct number of degrees of freedom and preserves the hypersurface orthogonality. In addition to introducing (31) and (32), we also have to include the time-like unit constraint in the Einstein-Æther action Eq. (16), we do so by adding it with a Lagrange multiplier $\lambda_{\text{\AE}}(u^{2}+1)$. Then, we find $\lambda_{\text{\AE}}$ must be set to

[TABLE]

in order to satisfy the unit constraint on the $u^{p}_{\mu}(x^{\sigma})$ field. Conveniently (33) works for both perturbations, (31) and (32). We can derive a new set of equations of motion generated by these perturbations (31), (32). 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-Æther QNMs found with the new Æther field perturbations (31) and (32) are numerically indistinguishable from those QNMs found using Eq. (20). 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 Fig. 3 and Fig. 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 (31).

3.5 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 khronon modes are those modes in the polar sector of Einstein-Æther Theory, which have purely imaginary (quasi)eigenfrequencies and are non-hydrodynamic taking on a nonzero frequency value at vanishing momentum. The khronon fluctuation, $\tilde{\chi}$ couples to the other six fluctuations in the polar sector, $\tilde{h}_{xx}$, $\tilde{h}_{yy}$, $\tilde{h}_{yt}$, $\tilde{h}_{yr}$, $\tilde{h}_{rr}$, $\tilde{h}_{tt}$, and $\tilde{h}_{rt}$. At vanishing couplings $\lambda=0$ and $\alpha=0$, it is possible to analytically map the polar fluctuations to the corresponding Schwarzschild-$AdS_{4}$ fluctuations in Einstein Gravity using a field redefinition Davison2016 . At nonzero $\lambda$, 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, $\tilde{\chi}=0$, does not affect the values of the other polar QNM frequencies (up to numerial errors). In turn, when forcing the metric fluctuations to vanish $\tilde{h}_{xx}=\tilde{h}_{yy}=\tilde{h}_{yt}=\tilde{h}_{yr}=\tilde{h}_{rr}=\tilde{h}_{tt}=\tilde{h}_{rt}=0$, we find the khronon eigenfrequencies unaffected (up to numerical errors). 171717 Note that this choice of vanishing metric fluctuations is inconsistent. According to the Einstein equations, also the khronon has to vanish in this case. Potentially the khronon field can be shown to decouple after a field transformation or use of master fields Kodama:2000fa . 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, $\nu_{khronon}=-i\,2.381\,n$ with $n=0,1,2,3,\dots$, as stated in Eq. (29). In fact, if we change our frequency normalization by a factor to $\hat{\nu}=\nu 2^{1/3}$, then $\hat{\nu}_{khronon}=-i\,3\,n$. Such integer value solutions are known from various analytical solutions for quasinormal mode frequencies Starinets:2002br . However, such behavior is also known from various fake modes Janiszewski:2015ura . 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 Fig. 10 and Fig. 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 Æther action (16) through dynamical terms for the vector $u_{\mu}\propto\partial_{\mu}\phi$, which are quadratic in derivatives on $u_{\mu}$. This leads to a fourth order equation of motion for $\phi$ or its fluctuation $\tilde{\chi}$, see Eq. (20). 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 horizon $\tilde{\chi}\approx\chi_{0}(1-r)^{\alpha}$, which are all of the form

[TABLE]

with a momentum dependent real-valued function $\mathfrak{f}(q)$, and we recall that $\hat{\nu}=2^{1/3}\nu$. 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

[TABLE]

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. (35) is regular at the horizon, except for special values of $\nu$. Those special values are $\hat{\nu}=-i\,3\,n$ with $n=0,1,2,3,\dots$. Generalizing the momentum to $q>0$, we find that regular modes appear at frequency values

[TABLE]

These are the khronon modes found by our pseudospectral method when solving the generalized eigenvalue problem. Eq. (36) explains the momentum dependence discussed above. At $q=0$, Eq. (36) also explains the observation that khronon mode frequencies are integer multiples of 3 when written in terms of $\hat{\nu}$. The numerical data indicates that this property also holds at $q\neq 0$, which implies that $\mathfrak{f}(q)\propto n$. Now the question remains, if these khronon modes are to be regarded as true QNMs or as fake modes.181818It 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.

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.191919As 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 Sybesma2015 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.

As a further test, we consider the large momentum $q\gg 1$ limit, also called eikonal limit. Results for the two purely imaginary khronon modes were already discussed in Sec. 3.3, and presented in Fig. 10. The observed phase is $Arg(s_{n})\approx-3\pi/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_{\mu}$ and compare its numerical behavior with what is expected from that analysis. However, such a treatment is beyond the scope of this work.

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 $\alpha=0$ is forcing part of the khronon dynamics to vanish. This is plausible because the term in the Æther action (16) set to zero by $\alpha=0$ (equivalent to $c_{4}=0$) is essentially quadratic in a time derivative of the vector $u_{\mu}$. We speculate that $\alpha\neq 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 Janiszewski2015 , which is left for future work.

4 Summary & conclusions

In this paper we have calculated non-relativistic gravitational QNMs on an analytically known asymptotically $AdS_{4}$ black brane solution, Eq. (5) and (18), with one of the three coupling constants vanishing, $\alpha=0$ Janiszewski2015 . 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 $\lambda$, $\beta$, and including large momentum up to $q=50$. Our results are summarized in Fig. 2 and Fig. 2 for the polar and axial sector, respectively. Equations of motion, and QNM data is collected in four ancillary files.

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 $\lambda,\beta$ at $\alpha=0$, when expressed in appropriate units. At $\lambda\neq 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. (36),

[TABLE]

Furthermore, we conjecture an analytic relation between the QNMs of Einstein Gravity and all QNMs of Hořava Gravity at arbitrary $\lambda$, $\beta$, and at $\alpha=0$, except the khronon modes:

[TABLE]

where $r_{\text{Schwarzschild}}$ is the Schwarzschild horizon of a black brane in Einstein Gravity, and $r_{\text{sound}}$ is the sound horizon relevant for each sector of QNMs in the analytic Hořava Gravity black brane solution, Eq. (18). 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 $\sqrt{1+\beta}$.

In the axial sector, at any value of $\lambda$ and $\beta$, there is one hydrodynamic diffusion mode and a set of propagating (not overdamped) non-hydrodynamic QNMs, see Fig. 3. The absence of overdamped (purely imaginary modes) in this sector is in agreement with the claim from Sybesma2015 ; Gursoy:2016tgf . The hydrodynamic diffusion mode starts out having quadratic dispersion at small momentum in agreement with the analytic prediction, Eq. (25) and Fig. 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 Kaminski:2009ce . At large momentum, $q\gg 1$, the non-hydrodynamic modes dominate the system since their damping decreases and they are long lived as seen from Fig. 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, $\lambda=0$ and $\lambda\neq 0$, see Fig. 6. If $\lambda=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 Fig. 7, agree with the linear propagation and quadratic damping in Eq. (28), which was derived only at $\lambda=0$. Our analysis shows that this equation holds also at $\lambda\neq 0$. Again, at large momentum, $q\gg 1$, the system is likely dominated by the non-hydrodynamic modes, because again their damping decreases, as seen in Fig. 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 [math] and $-i$ for large momentum. In addition to that, in the other case, $\lambda\neq 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$ Sybesma2015 ; Gursoy:2016tgf . 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 Morgan:2009vg ; Festuccia:2008zx ; Fuini2016 , see Fig. 9 and 10. At large momentum, $q\approx 50$, all our metric (non-overdamped) QNMs follow the relation $\nu(q)\approx q+s_{n}q^{-1/5}$, asymptoting to a constant magnitude for $s_{n}$ and with a universal phase $Arg(s_{n})\approx-\pi 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 Fig. 9 and 10, their $s_{n}$ values do not asymptote to constants and their phase is not $\pm\pi 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 $\alpha$ 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) and (16). Moreover, also the metric modes in the polar sector will travel at a finite speed, which can be distinct from the khronon speed, depending on the Hořava couplings Jacobson2004 . 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, $\Gamma\propto\eta$, 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 $\alpha$.

In passing, we have verified the equivalence of axial QNMs derived from hypersurface orthogonal Einstein-Æther Theory and those derived from Hořava Gravity. Remarkably, the hypersurface orthogonality does not influence the QNMs in our system, as our semi-Æther results indicate, see Sec. 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 $\alpha$ coupling, or equivalently Einstein-Æther 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 $\alpha=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\leq 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. Jensen:2011xb , and then take a non-relativistic limit sending the speed of light to infinity Kaminski:2013gca , where the choice of the hydrodynamic frame is important Jensen:2014wha . A second way is to identify the non-relativistic data structures directly, as is done in the context of Newton-Cartan geometry Son:2013rqa ; Jensen:2014aia ; Jensen:2014ama . It has been shown that dynamical Newton-Cartan geometry gives rise to Hořava Gravity Hartong:2015zia . 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.

Acknowledgements.

MK thanks C. Uhlemann for helpful correspondence. We thank R. Davison for helpful comments on the manuscript. This work was supported, in part, by the U.S. Department of Energy grant DE-SC-0012447.

Appendix A Tabulated data

Exemplary QNM values for comparison can be found in tables 4 and 5. Much more extensive data has been included with this submission as ancillary files.

Appendix B Convergence

Following Jansen2017 , the pseudo-spectral method used in this paper is essentially an eigenvalue solver. Simply put, the expansion of fluctuations into polynomials of $n$th 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 $\Delta\omega_{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 Fig. 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 $\Delta\omega_{cutoff}=10^{-8}$ to test for convergence. The precision of our arithmetic is $\thicksim 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.

Appendix C Khronon scalar field

The Khronon scalar field takes the following form.

[TABLE]