# Well-posed evolution of field theories with anisotropic scaling:
the Lifshitz scalar field in a black hole space-time

M. E. Rubio<EMAIL_ADDRESS>SISSA, Via Bonomea 265, 34136 Trieste, Italy and
INFN Sezione di Trieste  IFPU - Institute for Fundamental Physics of the
Universe, Via Beirut 2, 34014 Trieste, Italy  Á. D. Kovács<EMAIL_ADDRESS>SISSA, Via Bonomea 265, 34136 Trieste, Italy and INFN Sezione di Trieste  IFPU
- Institute for Fundamental Physics of the Universe, Via Beirut 2, 34014
Trieste, Italy  M. Herrero-Valea<EMAIL_ADDRESS>Institut de Fisica d’Altes
Energies (IFAE), The Barcelona Institute of Science and Technology, Campus
UAB, 08193 Bellaterra (Barcelona) Spain  M. Bezares
<EMAIL_ADDRESS>SISSA, Via Bonomea 265, 34136 Trieste,
Italy and INFN Sezione di Trieste  IFPU - Institute for Fundamental Physics of
the Universe, Via Beirut 2, 34014 Trieste, Italy  Nottingham Centre of
Gravity, University of Nottingham, University Park, Nottingham, NG7 2RD, UK
School of Mathematical Sciences, University of Nottingham, University Park,
Nottingham, NG7 2RD, UK  E. Barausse<EMAIL_ADDRESS>SISSA, Via Bonomea
265, 34136 Trieste, Italy and INFN Sezione di Trieste  IFPU - Institute for
Fundamental Physics of the Universe, Via Beirut 2, 34014 Trieste, Italy

###### Abstract

Partial differential equations exhibiting an anisotropic scaling between space
and time – such as those of Hořava-Lifshitz gravity – have a dispersive
nature. They contain higher-order spatial derivatives, but remain second order
in time. This is inconvenient for performing long-time numerical evolutions,
as standard explicit schemes fail to maintain convergence unless the time step
is chosen to be very small. In this work, we develop an implicit evolution
scheme that does not suffer from this drawback, and which is stable and
second-order accurate. As a proof of concept, we study the numerical evolution
of a Lifshitz scalar field on top of a spherically symmetric black hole space-
time. We explore the evolution of a static pulse and an (approximately)
ingoing wave-packet for different strengths of the Lorentz-breaking terms,
accounting also for the effect of the angular momentum eigenvalue and the
resulting effective centrifugal barrier. Our results indicate that the
dispersive terms produce a cascade of modes that accumulate in the region in
between the Killing and universal horizons, indicating a possible instability
of the latter.

##  I Introduction

In 2009, Petr Hořava proposed a power-counting renormalizable ultraviolet (UV)
completion of General Relativity (GR) Horava (2009) by endowing space-time
with a preferred foliation in space-like hypersurfaces. This comes at the cost
of reducing the local symmetry of the theory down to foliation preserving
diffeomorphisms (FDiff), namely

$t\rightarrow t^{\prime}(t),\quad x^{i}\rightarrow x^{\prime i}(t,x),$ (1)

where $t$ is the preferred time direction along the foliation, and $x^{i}$ is
a chart parametrizing the orthogonal leafs.

Action functionals invariant under (1) allow for time and spatial derivative
operators to be included independently, thus breaking local Lorentz invariance
(LLI), but also attaining a faster convergence of UV integrals in Quantum
Field Theory Fujimori  _et al._ (2016); Anselmi and Halat (2007). Ostrogradsky
ghosts are avoided simply by fixing the number of time derivatives to two,
with the action nonetheless still admitting higher spatial derivative
operators. Power-counting renormalizability is achieved in four-dimensional
gravity in this way by adding up to six spatial derivatives in a FDiff
invariant way.

A remarkable consequence of this construction is that, due to the different
order of derivatives along distinct directions, the UV dynamics of the theory
becomes invariant under an anisotropic Lifshitz scaling with critical exponent
$z$,

$\displaystyle t\rightarrow b^{z}t,\quad x^{i}\rightarrow bx^{i},\quad
b\in\mathbb{R},$ (2)

where the number of spatial derivatives in the UV is $2z$. This is the scaling
symmetry that allows for the modified power counting of UV divergences leading
to renormalizability – with $z=3$ in the case of gravity in four space-time
dimensions.

Since its formulation, countless works have explored the consequences of
Hořava’s proposal. A far from exhaustive list includes: understanding the
infra-red (IR) dynamics of the theory Blas  _et al._ (2011, 2010); the
structure of constraints Donnelly and Jacobson (2011); Bellorin  _et al._
(2022); its perturbative UV dynamics Barvinsky  _et al._ (2016, 2017a, 2019,
2022); Radkovski and Sibiryakov (2023); Griffin  _et al._ (2017); Benedetti
and Guarnieri (2014); the search for black-hole Barausse  _et al._ (2011);
Barausse and Sotiriou (2013a, b); Berglund  _et al._ (2012), regular Lara  _et
al._ (2021); Mazza and Liberati (2023), and cosmological solutions Mukohyama
(2010); the interaction with matter fields Pospelov and Shang (2012); Blas and
Lim (2015), and low-energy signatures in observations Audren  _et al._ (2015);
Barausse (2019); Bettoni  _et al._ (2017); Bonetti and Barausse (2015);
Cornish  _et al._ (2017); Emir Gümrükçüoğlu  _et al._ (2018); Yagi  _et al._
(2014a, b); Franchini  _et al._ (2021); Gupta  _et al._ (2021), among others.
However, despite incredible advances, the problem of classically solving and
dynamically evolving equations of motion exhibiting the scaling (2) is still
mostly unexplored, except in the simplest situations – such as for static and
spherically symmetric space-times, or in the perturbative limit Oshita  _et
al._ (2021), where the coefficients accompanying the higher derivative terms
are small enough.

Due to the presence of higher-order spatial derivatives, equations exhibiting
an anisotropic scaling have a dispersive nature. While in GR the dispersion
relation of massless fields takes a universal form
$p^{\mu}p_{\mu}=\omega^{2}-k^{2}$=0, with $p^{\mu}=(\omega,k^{i})$ the four-
momentum of the field111Hereinafter we use a mostly minus signature for the
metric. Greek indices denote the full set of space-time components, while
Latin ones are restricted to spatial directions only.; the presence of higher
spatial derivatives generically modifies it to

$\displaystyle\omega^{2}=k^{2}+\Theta k^{2z},$ (3)

where we have assumed the addition of an operator with $2z$ spatial
derivatives, controlled by a dimensionful coupling $\Theta$. Standard explicit
numerical methods tend to develop instabilities when evolving equations of
motion exhibiting this behaviour (Press  _et al._ , 2007). For hyperbolic
systems, a Von Neumann analysis of explicit methods – such as the standard
Runge-Kutta algorithms – shows that the corresponding evolution is stable as
long as the dimensionless Courant number

$C=\frac{v_{\text{m}}\Delta t}{\Delta x}$ (4)

satisfies $C\leq C_{\text{max}}$, where $\Delta t$ and $\Delta x$ are the time
and space steps used for discretization, and $v_{\text{m}}$ is the maximum
speed of propagation allowed by the equation. The value of $C_{\text{max}}$
generically depends on the method, but for explicit schemes it usually takes
the value $C_{\text{max}}=1$.

In the dispersive case instead, a similar von Neumann analysis yields a more
restrictive condition for stability, of the form

$\frac{\Delta t}{(\Delta x)^{p}}\leq\alpha,$ (5)

where $p>1$ and $\alpha>1$ depend both on the equation and on the scheme. This
is very inconvenient if long-time simulations with high spatial resolution are
sought, as the time step needed to keep the numerical error bounded has to be
chosen very small.

In this work, we circumvent this issue by providing a robust numerical scheme
that allows for evolving initial value problems exhibiting a UV Lifshitz
scaling (2) in non-trivial curved backgrounds. We do this by building an
implicit numerical method, where the usual finite difference discretization of
derivative operators is substituted by averages between two consecutive time
steps. This results in an integration scheme whose stability analysis leads to
no bounds for the time and spatial steps, allowing for long-time evolution
regardless of the desired spatial resolution. We show the robustness of our
methods by studying a situation of physical interest – the scattering of a
Lifshitz scalar field on top of a static, spherically symmetric, and
asymptotically flat black hole space-time, which we choose to be a solution to
the field equations of Hořava gravity at low energies.

This paper is organized as follows. In section II we introduce the formulation
of the Lifshitz scalar field theory both in flat and curved space-times,
obtaining the equations of motion and discussing their dispersive character.
An implicit method for the numerical evolution of equations of this kind is
introduced in section III, while an application to the Lifshitz scalar field
around a Lorentz-violating black hole is discussed in detail in section IV.
Finally, we show our results in section V, drawing conclusions and future
research directions in section VI.

##  II The Lifshitz scalar field

Field theories exhibiting the anisotropic scaling (2) have been known for
quite a long time. The simplest among them is the Lifshitz scalar field, used
in studies of quantum phase transitions in various strongly correlated systems
Ardonne  _et al._ (2004). Its action in a $3+1$ flat space-time with $z=3$
reads

$\displaystyle S=\frac{1}{2}\int
dtd^{3}x\left((\partial_{t}\phi)^{2}-\phi(-\partial^{2})\phi-\frac{c_{2}}{\Lambda^{2}}\phi(-\partial^{2})^{2}\phi-\frac{c_{3}}{\Lambda^{4}}\phi(-\partial^{2})^{3}\phi\right),$
(6)

where $\partial^{2}=\partial_{i}\partial^{i}$, $\Lambda$ is an energy scale,
all the $c_{i}$ are dimensionless, and we have fixed the speed of light to
one, which is always possible by rescaling the spatial coordinates. Note that
the action (6) contains all possible parity-preserving spatial derivative
operators with $z\leq 3$, while keeping two time derivatives to avoid
Ostrogradsky ghosts. This structure allows for a dynamical flow of $z$ with
energy, which changes from $z=1$ in the deep IR to $z=3$ in the UV, thus
mimicking the behavior expected in Hořava gravity Horava (2009). We have
omitted a possible potential for the scalar field, since our goal is to focus
on the effect of the higher derivative operators.

The equation of motion obtained from (6) reads

$\displaystyle\left[\partial_{t}^{2}-\partial^{2}+\frac{c_{2}}{\Lambda^{2}}(\partial^{2})^{2}-\frac{c_{3}}{\Lambda^{4}}(\partial^{2})^{3}\right]\phi=0,$
(7)

which leads to a modified dispersion relation of the form (3), that is

$\displaystyle\omega^{2}=k^{2}+\frac{c_{2}k^{4}}{\Lambda^{2}}+\frac{c_{3}k^{6}}{\Lambda^{4}}.$
(8)

This becomes explicitly Lorentz-violating for momentum scales
$k\gtrsim\Lambda$. Stability of the degrees of freedom demands $\omega^{2}>0$,
which in particular requires the parabola
$y(\hat{x})=1+c_{2}\hat{x}+c_{3}\hat{x}^{2}$ to be concave, leading to
$c_{3}>0$. Regarding the value of $c_{2}$, $\omega^{2}$ is positive for all
$c_{2}>0$, while for $c_{2}<0$ we must require the minimum of the parabola to
be above zero, corresponding to $c_{2}^{2}\leq 4c_{3}$. Putting both
conditions together, we get $c_{2}>-2\sqrt{c_{3}}$.

In order to extend this structure to curved space-times, we couple the scalar
field to Hořava gravity Horava (2009), which implements the anisotropic
scaling (2) through a foliation in co-dimension one space-like hypersurfaces.
This is described by supplementing the space-time curved metric $g_{\mu\nu}$
with a hypersurface-orthogonal unit-norm and timelike vector $U^{\mu}$, called
the _æther_ Jacobson and Mattingly (2001). Time derivatives are identified
with derivatives along $U^{\mu}$, while the spatial Laplacian $\partial^{2}$
is replaced by $\Delta_{\gamma}=\gamma^{\mu\nu}\nabla_{\mu}\nabla_{\nu}$,
where $\gamma^{\mu\nu}$ is the orthogonal projector
$\gamma^{\mu\nu}=-g^{\mu\nu}+U^{\mu}U^{\nu}$. The action thus takes the form
Barvinsky  _et al._ (2017b)

$\displaystyle S_{\phi}=\frac{1}{2}\int dtd^{3}x\sqrt{|g|}\left(({\cal
L}_{U}\phi)^{2}-\phi(-\Delta_{\gamma})\phi-\frac{c_{2}}{\Lambda^{2}}\phi(-\Delta_{\gamma})^{2}\phi-\frac{c_{3}}{\Lambda^{4}}\phi(-\Delta_{\gamma})^{3}\phi\right),$
(9)

where ${\cal L}_{U}$ is the Lie derivative along $U^{\mu}$. The corresponding
equation of motion is

$\displaystyle\left[\square+\frac{c_{2}}{\Lambda^{2}}\Delta_{\gamma}^{2}-\frac{c_{3}}{\Lambda^{4}}\Delta_{\gamma}^{3}\right]\phi=0,$
(10)

where
$\Box=\mathcal{L}^{2}_{U}+a^{c}\nabla_{c}+K\mathcal{L}_{U}-\Delta_{\gamma}$ is
the standard four-dimensional diffeomorphism invariant D’Alembertian, with
$a_{\mu}=U^{\nu}\nabla_{\nu}U_{\mu}$ the acceleration of the æther, and
$K=\gamma^{\mu\nu}K_{\mu\nu}$ the trace of the extrinsic curvature of the
foliation leafs $K_{\mu\nu}=(\mathcal{L}_{U}\gamma_{\mu\nu})/2$. This equation
can be then thought of as a generalization of (7) to curved space-times, but
also as a useful toy model that captures many of the subtleties brought by the
presence of the anisotropic scaling (2), in particular regarding the
integration of the equations of motion with numerical methods.

Equation (10) admits a well-posed initial value problem; that is, for an
appropriate initial data set there exists (locally) a unique solution, which
depends continuously on the data. This condition holds when $c_{3}\geq 0$ and
$c_{2}\geq 0$ (or if $c_{2}<0$ and $|c_{2}|$ is sufficiently small). The proof
of this statement, as well as a more detailed study of the Cauchy problem in
Hořava-Lifshitz theories will be explored in more detail in a forthcoming work
Kovács . Notice that these conditions include the ones required before for
stability of the degrees of freedom.

A first attempt at evolving (10) on top of a black hole space-time was
performed in Oshita  _et al._ (2021), using a finite difference method.
Although their scheme is convergent, stability is only guaranteed as long as
the time step used for simulations is small enough, which can be problematic
when aiming for long run-times. This can be seen from a Von Neumann analysis
of the local stability of the numerical scheme Press  _et al._ (2007).
Assuming that the coefficients of the equations are locally constant (both in
space and time), the eigenmodes of the discretized equations take the form
$\phi^{n}_{j}=\xi^{n}e^{ikj\Delta x}$, where $\phi^{n}_{j}$ is the numerical
approximation of the solution at time step $t_{n}=n\Delta t$ and grid point
$x_{j}:=j\Delta x$, with $\Delta t$ and $\Delta x$ the time and spatial steps
used for numerical integration, and $k$ is the frequency of the mode (assuming
one spacial dimension for simplicity). The amplification factor $\xi^{n}$
generically depends on $k$ and controls the amplitude of the $k-$eigenmode.
Since time evolution of a single eigenmode scales as some power of $\xi^{n}$,
the discretized equations will be stable if and only if
$\left|\xi^{n}\right|\leq 1$ for all $k$. For hyperbolic equations, this
analysis implies the well known Courant-Friedrichs-Lewy (CFL) condition, which
relates $\Delta t$ and $\Delta x$ with the maximum characteristic speed in the
equation. For dispersive equations such as (10) instead, the anisotropic
scaling between space and time leads to a stronger upper-bound for the time
step $\Delta t$, which scales in the form given in (5), for some $p>1$. This
issue complicates the running of long simulations with high spatial
resolutions, as computations quickly become very costly. Achieving stable
simulations of dispersive equations of the Lifshitz kind thus requires to
consider a different approach.

##  III An implicit method for equations with anisotropic scaling

In order to achieve stable, high-resolution numerical evolution, we present
here a fully implicit scheme, following the spirit of the well-known Crank-
Nicholson method for simple diffusion equations (see (Press  _et al._ , 2007)
for references). For the sake of generality, we will diverge for the moment
from the specifics of Eq. (10), and consider instead a generic partial
differential equation (PDE) in $1+1$ dimensions for a scalar field
$\psi(t,x)$, of the form

$\displaystyle\partial^{2}_{t}\psi$ $\displaystyle=$ $\displaystyle
F\left(t,x,\psi,\partial_{t}\psi,\\{\partial^{j}_{x}\psi\\}_{j=1}^{J},\\{\partial^{k}_{x}\partial_{t}\psi\\}_{k=0}^{K}\right),$
$\displaystyle\psi\mathcal{j}_{t=0}$ $\displaystyle=$
$\displaystyle\psi_{0}(x),$ (11)
$\displaystyle\partial_{t}\psi\mathcal{j}_{t=0}$ $\displaystyle=$
$\displaystyle\psi_{1}(x),$

where $F$ is a sufficiently smooth function, linear in $\psi$ and in all its
derivatives, and $\psi_{0},\,\psi_{1}$ are given initial conditions. The
indices $j$ and $k$ span the ranges $1\leq j\leq J$, and $0\leq k\leq K$, with
$J\geq 3$, and $K\geq 0$. With these conditions, Eq. (III) always contains
more spatial derivatives than time derivatives. Hence, it corresponds to a PDE
with a dispersive nature, leading in fact to a dispersion relation of the form
$\omega^{2}=p(k)$, where $p(k)$ is a polynomial of order greater than or equal
to three. Notice that the general problem (III) does include the Lifshitz
field equation (10).

We aim to numerically evolve Eq. (III), for which we consider a first-order
reduction in time, with dynamical variables $\psi$ and
$\Pi:=\partial_{t}\psi$; namely,

$\displaystyle\dot{\psi}$ $\displaystyle=$ $\displaystyle\Pi$ (12)
$\displaystyle\dot{\Pi}$ $\displaystyle=$ $\displaystyle
F\left(t,x,\psi,\Pi,\\{\partial^{j}_{x}\psi\\}_{j=1}^{J},\\{\partial^{k}_{x}\Pi\\}_{k=0}^{K}\right).$
(13)

We consider a uniform spatial grid of $N$ points $x_{i}=x_{0}+(i-1)\Delta x$,
$i=1,\ldots,N$, with step $\Delta x=L/(N-1)$, being $L$ the length of the
spatial domain. We also discretize time as $t_{n}=n\Delta t$,
$n=0,1,2,\ldots$, with time step $\Delta t$, and define the grid functions as
$\psi^{n}_{i}:=\psi(t_{n},x_{i})$.

While for standard explicit methods, derivatives are discretized using finite
difference operators, here we replace them with the average of the derivatives
in two consecutive time steps. This procedure leads to a fully implicit
scheme, meaning that in order to evaluate the solution at a given time step,
we need to solve an algebraic equation containing the information of the
solution both at the current and previous steps. More specifically, a second-
order accurate approximation of the first derivative of $\psi$ at the grid
point $i$, centered at a time step $n+\frac{1}{2}$, is computed as

$(D\psi)^{n+\frac{1}{2}}_{i}\equiv\frac{1}{2}\left[\frac{\psi^{n}_{i+1}-\psi^{n}_{i-1}}{2\Delta
x}+\frac{\psi^{n+1}_{i+1}-\psi^{n+1}_{i-1}}{2\Delta x}\right].$ (14)

For the second derivative, one gets

$(D^{2}\psi)^{n+\frac{1}{2}}_{i}\equiv\frac{1}{2}\left[\frac{\psi^{n}_{i+1}-2\psi^{n}_{i}+\psi^{n}_{i-1}}{(\Delta
x)^{2}}+\frac{\psi^{n+1}_{i+1}-2\psi^{n+1}_{i}+\psi^{n+1}_{i-1}}{(\Delta
x)^{2}}\right];$ (15)

and similarly for higher-order derivatives. The same average is also performed
for evaluating the field $\psi^{n+\frac{1}{2}}_{i}$, namely

$\psi^{n+\frac{1}{2}}_{i}\equiv\frac{1}{2}\left(\psi^{n}_{i}+\psi^{n+1}_{i}\right).$
(16)

For the time derivative, instead, we just keep the forward finite difference

$(\partial_{t}\psi)^{n+1}_{i}\equiv\frac{\psi^{n+1}_{i}-\psi^{n}_{i}}{\Delta
t},$ (17)

ensuring thus second-order accuracy also in time.

Therefore, by replacing all the averaged fields and derivatives into the
previous system, we get the following scheme:

$\displaystyle\frac{\psi^{n+1}_{i}-\psi^{n}_{i}}{\Delta t}$ $\displaystyle=$
$\displaystyle\Pi^{n+\frac{1}{2}}_{i}$ (18)
$\displaystyle\frac{\Pi^{n+1}_{i}-\Pi^{n}_{i}}{\Delta t}$ $\displaystyle=$
$\displaystyle
F\left(t_{n+\frac{1}{2}},x_{i},\psi^{n+\frac{1}{2}}_{i},\Pi^{n+\frac{1}{2}}_{i},\\{(D^{j}\psi)^{n+\frac{1}{2}}_{i}\\}_{j=1}^{J},\\{(D^{k}\Pi)^{n+\frac{1}{2}}_{i}\\}_{k=0}^{K}\right),$
(19)

where $D^{j}$ denotes the finite difference average operator of the $j-$th
spatial derivative. Due to the linearity of $F$ with respect to derivatives,
the above system can be recast in matrix form, obtaining a band-diagonal
linear system for the variables evaluated at time step $n+1$, namely

$\begin{bmatrix}d_{1}&u_{11}&0&\cdots&\cdots&0\\\
\ell_{1}&d_{2}&u_{12}&u_{22}&\cdots&\vdots\\\
0&0&d_{1}&u_{11}&\ddots&\vdots\\\ \ell_{3}&\ell_{2}&\ell_{1}&d_{2}&\ddots&0\\\
0&0&0&0&d_{1}&\vdots\\\ \vdots&\ddots&\ddots&\ddots&\ddots&\vdots\\\
0&\cdots&\cdots&0&\cdots&d_{2N}\end{bmatrix}\begin{bmatrix}\psi^{n+1}_{2}\\\
\Pi^{n+1}_{2}\\\ \psi^{n+1}_{3}\\\ \Pi^{n+1}_{3}\\\ \vdots\\\ \vdots\\\
\psi^{n+1}_{N-1}\\\ \Pi^{n+1}_{N-1}\\\
\end{bmatrix}=\begin{bmatrix}\bar{d}_{1}&\bar{u}_{11}&0&\cdots&\cdots&0\\\
\bar{\ell}_{1}&\bar{d}_{2}&\bar{u}_{12}&\bar{u}_{22}&\cdots&\vdots\\\
0&0&\bar{d}_{1}&\bar{u}_{11}&\ddots&\vdots\\\
\bar{\ell}_{3}&\bar{\ell}_{2}&\bar{\ell}_{1}&\bar{d}_{2}&\ddots&0\\\
0&0&0&0&\bar{d}_{1}&\vdots\\\ \vdots&\ddots&\ddots&\ddots&\ddots&\vdots\\\
0&\cdots&\cdots&0&\cdots&\bar{d}_{1}\end{bmatrix}\begin{bmatrix}\psi^{n}_{2}\\\
\Pi^{n}_{2}\\\ \psi^{n}_{3}\\\ \Pi^{n}_{3}\\\ \vdots\\\ \vdots\\\
\psi^{n+1}_{N-1}\\\ \Pi^{n+1}_{N-1}\\\ \end{bmatrix}+\begin{bmatrix}b_{11}\\\
b_{21}\\\ 0\\\ 0\\\ \vdots\\\ \vdots\\\ b_{1N}\\\ b_{2N}\\\ \end{bmatrix},$
(20)

where the second term on the right-hand side contains the information on the
boundary conditions, which need to be specified. The number of upper/lower
diagonal bands in the matrices depends on the stencil used – that is, on the
number of neighboring points needed to approximate the derivatives at a given
grid point. The coefficients $\\{d_{j},u_{ij},\ell_{j}\\}$ and
$\\{\bar{d}_{j},\bar{u}_{ij},\bar{\ell}_{j}\\}$ are grid functions that also
depend on $\Delta x$ and $\Delta t$, but not on the dynamical fields. The
system (20) can be solved with standard numerical methods for band-diagonal
linear systems. These usually require $\mathcal{O}(N)$ iterations at each time
step, unlike Gaussian elimination, which requires $\mathcal{O}(N^{3})$
iterations (see (Press  _et al._ , 2007) for details on the corresponding
algorithm). Finally, notice that the generic scheme (18)-(19) is also valid
for the case in which $F$ is a non-linear function of $\psi$, which would
clearly lead to a non-linear system of coupled equations for the grid
functions $\psi^{n}_{i}$ and $\Pi^{n}_{i}$. Nevertheless, such a system can be
solved e.g. by means of standard iterative numerical methods.

##  IV The Lifshitz field in a black hole space-time

As a proof of concept of the numerical method introduced in the previous
section, we evolve the Lifshitz equation (10) on top of a spherically
symmetric, static, and asymptotically flat black hole space-time Berglund  _et
al._ (2012), solution to the equations of motion of khronometric gravity in
vacuum Blas  _et al._ (2011); Jacobson (2014). This corresponds to the low-
energy limit of Hořava gravity, as discussed in appendix A. We are thus
implicitly assuming that the backreaction of the scalar field onto the
geometry is negligible, and that gravitational perturbations are suppressed.
It remains unknown whether this is a solid assumption, precisely because a
proper understanding of the gravitational dynamics in Hořava gravity at all
energies would require a more sophisticated version of the methods that we are
pioneering here.

###  IV.1 The background solution

As already commented, we describe the derivation of the space-time solution in
appendix A, but report it here for practical purposes. The metric and æther in
Schwarzschild coordinates $\\{\tau,r,\theta,\varphi\\}$ read

$\displaystyle{\rm d}s^{2}=f(r){\rm d}\tau^{2}-\frac{{\rm
d}r^{2}}{f(r)}-r^{2}{\rm d}\Omega^{2},\quad U_{\mu}{\rm
d}x^{\mu}=\frac{H(r)}{2A(r)}{\rm d}\tau+\frac{1-f(r)A(r)^{2}}{2A(r)f(r)}{\rm
d}r$ (21)

with

$\displaystyle f(r)$
$\displaystyle=1-\frac{2\mu}{r}-c_{13}\frac{r_{\text{\ae}}^{4}}{r^{4}},$ (22)
$\displaystyle A(r)$
$\displaystyle=\frac{1}{f(r)}\left(-\frac{r_{\text{\ae}}^{2}}{r^{2}}+\sqrt{f(r)+\frac{r_{\text{\ae}}^{4}}{r^{4}}}\right),$
(23) $\displaystyle H(r)$ $\displaystyle=1+f(r)A(r)^{2},$ (24)

where $\mu$ is the mass of the black hole, ${\rm d}\Omega^{2}$ is the
$\mathbb{S}^{2}$ line element, and

$r_{\text{\ae}}=\frac{\mu}{2}\left(\frac{27}{1-c_{13}}\right)^{1/4}$ (25)

contains the only free parameter in the action $c_{13}$. Note that since the
solution is static, there exists a Killing vector $\chi^{\mu}=(1,0,0,0)$, with
norm $\chi^{2}=f(r)$. The latter flips sign at the surface $r=r_{K}$, where
$r_{K}$ is the positive solution of $f(r_{K})=0$, signaling the position of a
Killing horizon.

Although this form of the metric is useful for solving the equations of
motion, it can be problematic for numerically evolving Eq. (10). In this chart
of coordinates, the Laplacian in the orthogonal leafs, $\Delta_{\gamma}$, also
contains higher time derivatives, so the corresponding equation of motion
obtained from Eq. (10) does not fit the cases discussed in the previous
section. Nevertheless, this issue can be solved by aligning the time direction
with the integral curves of $U^{\mu}$, describing Eulerian observers in the
preferred frame of Hořava gravity. This can be achieved simply by introducing
another time coordinate $t$ (the preferred time) satisfying

$\displaystyle{\rm d}\tau={\rm d}t-\frac{1-f(r)A(r)^{2}}{f(r)H(r)}{\rm d}r,$
(26)

so that the metric and æther now read

$\displaystyle{\rm d}s^{2}=f(r){\rm d}t^{2}-2\frac{1-f(r)A(r)^{2}}{H(r)}{\rm
d}t{\rm d}r-\frac{4A(r)^{2}}{H(r)^{2}}{\rm d}r^{2}-r^{2}{\rm d}\Omega^{2},$
(27) $\displaystyle U_{\mu}{\rm d}x^{\mu}=\frac{H(r)}{2A(r)}{\rm d}t.$ (28)

Notice that after this transformation, the metric takes the Arnowitt-Deser-
Misner (ADM) form Arnowitt  _et al._ (1959)

$\displaystyle{\rm d}s^{2}=(N^{2}-N_{i}N^{i}){\rm d}t^{2}-2N_{i}{\rm
d}x^{i}dt-h_{ij}{\rm d}x^{i}{\rm d}x^{j},$ (29)

where $N,N^{i}$ and $h^{ij}$ are the lapse, shift and induced metric in the
foliation leafs, respectively, and are given by

$\displaystyle N=\frac{H(r)}{2A(r)},\quad N_{i}{\rm
d}x^{i}=\frac{1-f(r)A(r)^{2}}{H(r)}{\rm d}r,\quad h_{ij}{\rm d}x^{i}{\rm
d}x^{j}=\left(\frac{2A(r)}{H(r)}\right)^{2}{\rm d}r^{2}+r^{2}{\rm
d}\Omega^{2}.$ (30)

This chart of coordinates, however, has a pathology whenever $N\equiv 0$,
which corresponds to $H(r_{U})=1+f(r_{U})A(r_{U})^{2}=0$ for some value
$r_{U}$ of the radial coordinate Blas and Sibiryakov (2011); Barausse  _et
al._ (2011); Berglund  _et al._ (2012). From the relation (26), we see that
this point lies at finite $r$, but it is mapped to $t\rightarrow+\infty$,
signalling that the foliation cannot be globally extended smoothly beyond this
point, although observers moving inwards along $U^{\mu}$ can still cross it in
finite proper time (because ${\rm d}\tau\equiv U_{\mu}{\rm d}x^{\mu}$ remains
finite). Remarkably, the surface ${r=r_{U}}$ represents a trapping surface for
all modes, regardless of their propagation speed. This can be seen by noting
that the Killing vector $\chi^{\mu}$ is space-like in the vicinity of the
point $r=r_{U}$, while the product $\chi\cdot U=H(r)/(2A(r))$ precisely
vanishes at this point, and becomes negative in the inner region222Note
importantly that this cannot be avoided by a change of coordinates, since once
in the preferred frame the symmetry group of the theory is restricted to FDiff
(1), under which $N\rightarrow N(dt^{\prime}/dt)$.. This is enough to
characterize this surface as a universal trapping surface, hence named
_universal horizon_ Blas and Sibiryakov (2011); Barausse  _et al._ (2011). The
region $r<r_{U}$ always lies behind the Killing horizon – otherwise $\chi\cdot
U=0$ would not be possible, as $U^{\mu}$ is time-like everywhere by
definition. In our particular case, one has

$r_{U}=\frac{3\mu}{2}.$ (31)

For a more detailed discussion on causality within space-times endowed with
universal horizons, see Bhattacharyya  _et al._ (2016).

For our purposes here, the universal horizon implies a limitation. Since our
dynamical equations will be evolved in preferred time, we can only cover the
region $r>r_{U}$, as the foliation and the time coordinate $t$ do not extend
into the inner region. However, from the practical point of view of an
observer sitting at a large radius (the “asymptotic infinity”) this is enough,
since they cannot observe anything coming from inside the universal horizon.
Note, however, that the same is not true for the Killing horizon. While the
surface sitting at $f(r_{K})=0$ is a trapping surface in GR, it is not the
case here anymore, since causal modes can move at speeds larger than unity
Cropp  _et al._ (2014). The universal horizon is the only true trapping
surface within this geometry. However, the region between the two horizons
still encodes important features of the dynamics of the system, due to the
character change of the Killing vector, which is associated with the only
notion of conserved energy in the system.

###  IV.2 Numerical Implementation

We now provide details about the specifics of our simulations. In particular,
we introduce the ansatz for the solution, as well as the initial data and
boundary conditions. We also discuss the numerical scheme, following the
general approach described in Section III.

We solve Eq. (10) in spherical coordinates $\\{t,r^{*},\theta,\phi\\}$, where
the radial tortoise-like coordinate $r^{*}$ is chosen to push the universal
horizon to infinity, corresponding to $r^{*}\rightarrow-\infty$. The relation
between the areal radius $r$ and the new coordinate $r^{*}$ is given by333In
principle, any transformation of the form
$\frac{\text{d}r}{\text{d}r^{*}}=\left(\frac{H(r)}{2A(r)}\right)^{p}$ with
integer exponent $p\geq 1$ equally pushes the universal horizon to
$r^{*}\to-\infty$. We choose here $p=2$ so that the transformation decays in a
smoother way when approaching the universal horizon, avoiding localized high-
frequency instabilities that would otherwise require the addition of
artificial dissipation in the equations.

$\frac{\text{d}r}{\text{d}r^{*}}=\left(\frac{H(r)}{2A(r)}\right)^{2}.$ (32)

Although complicated in general, one can see that it behaves as $r_{*}\sim r$
for large $r$ due to asymptotic flatness of the background solution, which
implies $f(r\rightarrow\infty)\rightarrow 1$ and
$A(r\rightarrow\infty)\rightarrow 1$. Close to the universal horizon, we have
$H(r)/(2A(r))\propto(r-r_{U})$ and hence $r_{*}\propto-(r-r_{U})^{-1}$, indeed
placing $r_{U}$ at $r_{*}\rightarrow-\infty$. The condition (32) is an
ordinary differential equation, which can be solved numerically using the
standard fourth-order Runge-Kutta method.

We also define the parameters $\kappa_{2}=c_{2}/\Lambda^{2}$ and
$\kappa_{3}=c_{3}/\Lambda^{4}$ for computational convenience. Thus, to ensure
the stability and well-posedness of the problem, we require $\kappa_{3}\geq 0$
and $\kappa_{2}>-2\sqrt{\kappa_{3}}$. After this change, Eq. (10) reads

$\displaystyle\left[\square+\kappa_{2}\Delta_{\gamma}^{2}-\kappa_{3}\Delta_{\gamma}^{3}\right]\phi=0.$
(33)

Due to the spherical symmetry of the background, we consider the following
ansatz for the scalar field:

$\phi(t,r^{*},\theta,\varphi)=\sum_{\ell=0}^{\infty}\sum_{m=-\ell}^{\ell}{\psi_{\ell
m}(t,r^{*})Y_{\ell m}(\theta,\varphi)},$ (34)

where $Y_{\ell m}(\theta,\varphi)$ are the spherical harmonics. Plugging this
into Eq. (33), we get an effective $1+1$ dimensional differential equation for
every mode $\psi_{\ell m}(t,r^{*})$, which is independent of $m$, and reads

$\left[\partial^{2}_{t}+\zeta_{10}\partial_{t}+\zeta_{11}\partial_{t}\partial_{r^{*}}+\sum_{j=1}^{6}{\zeta_{0j}\partial^{j}_{r^{*}}}+V_{\text{eff}}\right]\psi_{\ell
m}=0.$ (35)

The coefficients $\zeta_{ij}$ and the effective potential $V_{\rm eff}$ are
functions of $r^{*}$ through the coordinate transformation (32), the metric
functions $A(r),\,H(r)$, and their derivatives up to fifth order. They also
depend on $\kappa_{2}$, $\kappa_{3}$ and $\ell$. The explicit formulae for
$\zeta_{ij}$ are given in Appendix (B), while we show $V_{\rm eff}$ later in
Eq. (42). Let us note that $V_{\rm eff}$ vanishes for $\ell=0$, and
corresponds to an effective centrifugal barrier, which plays an important role
in the behavior of the solution close to the Killing horizon, as we will
discuss later.

Following the scheme introduced in Section III, Eq. (35) admits the implicit
discretization

$\displaystyle\frac{\psi^{n+1}_{i}-\psi^{n}_{i}}{\Delta t}$ $\displaystyle=$
$\displaystyle\Pi^{n+\frac{1}{2}}_{i}$ (36)
$\displaystyle\frac{\Pi^{n+1}_{i}-\Pi^{n}_{i}}{\Delta t}$ $\displaystyle=$
$\displaystyle-\zeta^{i}_{10}\Pi^{n+\frac{1}{2}}_{i}-\zeta^{i}_{11}(D\Pi)^{n+\frac{1}{2}}_{i}-\sum_{j=1}^{6}{\zeta^{i}_{0j}(D^{j}\psi)^{n+\frac{1}{2}}_{i}}-V^{i}_{\rm
eff}\psi^{n+\frac{1}{2}}_{i},$ (37)

where $\zeta^{i}_{jk}:=\zeta_{jk}(r^{*}_{i})$ and $V^{i}_{\rm eff}:=V_{\rm
eff}(r^{*}_{i})$.

We implemented centered finite difference operators with second-order
accuracy, with a stencil of seven grid points. In this case, the system (20)
has five upper diagonals and seven lower ones. Table 1 shows the explicit form
of the coefficients in the matrices.

Element | Expression
---|---
$u_{11}$ | $\displaystyle\frac{1}{\Delta t}$
$u_{12}$ | $\displaystyle\frac{\zeta_{01}(\Delta x)^{5}+2\zeta_{02}(\Delta x)^{4}-2\zeta_{03}(\Delta x)^{3}-8\zeta_{04}(\Delta x)^{2}+5\zeta_{05}\Delta x+30\zeta_{06}}{4(\Delta x)^{6}}$
$u_{13}$ | $\displaystyle\frac{\zeta_{11}}{4\Delta x}$
$u_{14}$ | $\displaystyle\frac{\zeta_{03}(\Delta x)^{3}+2\zeta_{04}(\Delta x)^{2}-4\zeta_{05}\Delta x-12\zeta_{06}}{4(\Delta x)^{6}}$
$u_{15}$ | $\displaystyle\frac{\zeta_{05}\Delta x+2\zeta_{06}}{4(\Delta x)^{6}}$
$d_{1},\,\bar{d}_{1}$ | $\displaystyle\frac{1}{2}$
$d_{2}$ | $\displaystyle\frac{\zeta_{10}}{2}+\frac{1}{\Delta t}$
$\bar{d}_{2}$ | $\displaystyle-\frac{\zeta_{10}}{2}+\frac{1}{\Delta t}$
$\ell_{1}$ | $\displaystyle\frac{V_{\rm eff}}{2}-\frac{\zeta_{02}(\Delta x)^{4}-3\zeta_{04}(\Delta x)^{2}+10\zeta_{06}}{(\Delta x)^{6}}$
$\ell_{2}$ | $\displaystyle-u_{13}$
$\ell_{3}$ | $\displaystyle-u_{12}$
$\ell_{4}$ | $\displaystyle-u_{14}$
$\ell_{5}$ | $\displaystyle-u_{15}$
Table 1: Explicit form of the matrix elements for the scheme (36)-(37). Each
of them is a grid function, which should be evaluated at the grid point
$r^{*}_{i}$. Here we have omitted the corresponding indices for the sake of
notation clarity. For instance, $u_{11}$ corresponds to $u_{11}(r^{*}_{i})$.
Similarly for the coefficients $\zeta_{ij}$, which represent the grid
functions $\zeta_{ij}(r^{*}_{i})$.

####  IV.2.1 Initial data

As initial data for the Lifshitz field, we considered two different profiles.
The first one is a static Gaussian pulse, which we refer to as ID Type I,
given by

$\displaystyle\psi(0,r^{*})$ $\displaystyle=$ $\displaystyle
a_{0}\exp{\left[-\frac{(r^{*}-r^{*}_{c})^{2}}{\sigma^{2}}\right]}\,,$ (38)
$\displaystyle\Pi(0,r^{*})$ $\displaystyle=$ $\displaystyle 0\,,$

where $a_{0}$, $r^{*}_{c}$ and $\sigma$ are the Gaussian amplitude, mean and
variance, respectively.

The second profile that we consider, which we refer to as ID Type II, is given
by an (approximately) ingoing wave-packet

$\displaystyle\psi(0,r^{*})$ $\displaystyle=$
$\displaystyle\exp{\left[-\frac{(r^{*}-r^{*}_{c})^{2}}{\sigma^{2}}\right]}\cos{(\omega
r^{*})}\,,$ (39) $\displaystyle\Pi(0,r^{*})$ $\displaystyle=$
$\displaystyle-\left(\frac{2(r^{*}-r^{*}_{c})}{\sigma^{2}}\cos{(\omega
r^{*})}+\omega\sin{(\omega
r^{*})}\right)\exp{\left[-\frac{(r^{*}-r^{*}_{c})^{2}}{\sigma^{2}}\right]}\,.$

This corresponds to an exact ingoing wave-packet when
$\kappa_{2}=\kappa_{3}=0$, satisfying $\partial\psi/\partial
t|_{t=0}=\partial\psi/\partial r^{*}|_{t=0}$. In the Lifshitz case, it will
also contain outgoing modes, but we expect those to be negligible far enough
from the gravitational well, as long as the energy of the Killing energy of
the wave-packet $\omega$ is small. Hereinafter, we set $\omega=1$.

####  IV.2.2 Boundary conditions

Setting boundary conditions for equations endowed with an anisotropic scaling
(2) is highly non-trivial, due to the different order in derivatives along
distinct directions. While in the two-dimensional wave equation (regardless of
the boundary conditions imposed), the general solution can always be written
as the superposition of left-moving and right-moving waves,

$\displaystyle\Psi(t,r)=f_{L}(r-ct)+f_{R}(r+ct),$ (40)

with $f_{L},f_{R}$ arbitrary functions and $c$ the propagation speed of the
waves, this is no longer true for the case at hand, and in particular for Eq.
(33). This can be seen by simply plugging the ansatz (40) into (7). Only when
$f^{\prime}_{L}(x)\propto f_{L}(x)$, and $f^{\prime}_{R}(x)\propto f_{R}(x)$ –
where a prime denotes differentiation – solutions can be decoupled, as long as
$c$ takes the right value. This poses a problem for setting up a successful
evolution scheme. While in the case of the wave equation one can always impose
pure ingoing or outgoing boundary conditions, simply by selecting left or
right movers at the appropriate boundary surface, such procedure is not
possible here. This is of particular relevance at the universal horizon, which
is a semi-permeable surface that only allows for ingoing modes.

One possibility to face this issue that has been recently explored in other
contexts (Dima and Barausse, 2020; Dolan, 2013) is to add a perfectly matching
layer (PML) covering a small region close to the boundaries. This introduces
artificial dissipation suppressing spurious reflected waves, by modifying the
kinetic term in the equations of motion. Although this can be implemented
systematically when only second order derivatives are involved, we have not
found a way to extend it to the case with higher derivatives.

Nevertheless, following the spirit of the PML method, we implement instead an
artificial dissipative layer (ADL), controlled by a function $L(x)$, which
suppresses waves exiting the domain of interest during the numerical
evolution. The layer takes the shape of a function whose value is unity within
the physical domain, but decays smoothly to zero in the regions close to the
boundary. We implement it in our work by replacing the numerical solution
$\psi^{n}_{i}$ at each time step by $\psi^{n}_{i}\to L_{i}\psi^{n}_{i}$.

We choose in particular the following function for the ADL:

$L(x)=\left\\{\begin{array}[]{rcl}\displaystyle\frac{1-\tanh\left[s\left(x-x_{\rm
R}\right)\right]}{2},&\mbox{if}&x\geq\displaystyle\frac{x_{\rm L}+x_{\rm
R}}{2}\\\ \\\ \displaystyle\frac{1-\tanh\left[s\left(x_{\rm
L}-x\right)\right]}{2},&\mbox{if}&x<\displaystyle\frac{x_{\rm L}+x_{\rm
F}}{2}\end{array}\right.\,\,,$ (41)

where $s$ controls the slope of the function in the extremes of the numerical
domain, while $x_{\rm L}$ and $x_{\rm R}$ correspond to the left and right
midpoints of the decaying regions, respectively. An illustration of this
function is given in Figure 1.

Figure 1: Artificial dissipative layer. Example of $L(x)$ with parameters
$s=0.2$, $x_{L}=-75$, and $x_{R}=275$. The green-shaded rectangle represents
the physical region where the evolution is unaltered. The layer aims to mimic
boundary conditions describing waves escaping the numerical domain by
artificial damping at the extremes. The red vertical dotted lines are located
at $x=x_{\rm L}$ and $x=x_{\rm R}$, corresponding to the midpoints of the
damping zones.

##  V Results

We report here the results of our simulations, performed with different values
of the coupling parameters $\kappa_{2}$ and $\kappa_{3}$ controlling the
strength of the Lorentz-breaking terms in Eq. (33). We also investigate the
behaviour of wave modes with angular number $\ell>0$ for fixed values of the
couplings, as well as the effect of the centrifugal barrier induced by the
effective potential $V_{\rm eff}$. Finally, we probe the validity of the
numerical code by performing convergence tests, and show the robustness of our
results against a change in the position of the ADL zone.

###  V.1 Evolution of the Lifshitz field

We start off by studying the dynamical evolution of the $\ell=m=0$ mode for
the two different initial data introduced in the previous section.

We first evolve the initial data ID type I (static pulse) from $t=0$ to
$t=100\mu$, as shown in figure 2. We set $\kappa_{2}=0.1$ and vary
$\kappa_{3}$, taking the values $\kappa_{3}=\\{0.01,0.05,0.1\\}$ (from light
to dark blue in the plot panels). As it can be seen from the figure, we
observe the formation of a rapid cascade of modes nearing the black hole as
time increases. This is produced by wave modes travelling at different speeds
depending on their frequency, a behavior that is expected to occur due to the
dispersive character of the equation induced by the higher derivative terms.

From early evolution times, modes with faster speeds rapidly escape the
numerical domain, unlike slower ones, which stay longer within the physical
region. Moreover, we clearly observe that the magnitude of the modes within
the cascade grows faster for larger values of $\kappa_{3}$. A comparison with
the standard wave equation evolution (that is, when $\kappa_{2}=\kappa_{3}=0$)
is shown in the dotted black profile. As expected, the solution does not
exhibit a dispersive behavior in this case, as propagation speeds are bounded
and independent of the frequency. Notice also that due to the Lorentz
violating character of the equation, the solution can smoothly penetrate the
Killing horizon, travelling towards the universal horizon, which is located at
negative infinity. Finally, we observe a bump in the neighborhood of the
Killing horizon forming at later times, around $t\sim 80\mu$, whose magnitude
increases faster for larger values of $\kappa_{3}$. A similar behavior is
found in the case of the initial data ID Type II (ingoing pulse) for the same
choice of the parameters $\kappa_{2}$, and $\kappa_{3}$, as shown in figure 3.

Furthermore, in figures 4 and 5 we display the evolution of the solution for
fixed $\kappa_{3}$, varying instead $\kappa_{2}$. In particular, we set
$\kappa_{3}=0.01$ and consider $\kappa_{2}=\\{0.1,0.5,1.0\\}$. A cascade also
develops in this case, but we observe a decrease in the propagation speed of
the slowest modes when increasing $\kappa_{2}$, at least at early times. At
later times ($t\sim 60\mu$), and once the solutions have reached the Killing
horizon, the order of magnitude of the wave profiles with $\kappa_{2}=0.5$ and
$\kappa_{2}=1.0$ remain similar, until a bump around the horizon forms and
starts growing faster as we increase $\kappa_{2}$.

The appearance of such a bump in the region between the universal horizon and
the Killing horizon is interesting, since it might signal an instability of
the space-time background solution under certain assumptions. In particular,
let us note that Hořava gravity propagates a scalar degree of freedom together
with the usual transverse traceless graviton perturbation Blas  _et al._
(2011). When expanded around the background solution considered here, the
dynamics of the scalar mode must exhibit, for consistency, a Lifshitz scaling
of the form (2), and hence the equation of motion for scalar perturbations
must unavoidably take the form (10), with the parameters $c_{2}$ and $c_{3}$
somehow related to the couplings in the gravitational action. Hence, we can
conjecture that, if the bump found in our numerical experiments is a generic
feature of (10), it will also develop in the gravitational case, therefore
signaling a linear instability of the universal horizon (Herrero-Valea and
Rubio, ).

Figure 2: Dynamics of the Lifshitz field (ID type I). Snapshots of the
evolution of the $\ell=m=0$ mode, for $\kappa_{2}=0.1$ and different values of
$\kappa_{3}$: $\kappa_{3}=0.01$ (light blue), $\kappa_{3}=0.05$ (blue) and
$\kappa_{3}=0.1$ (dark blue). The initial data is the ID type I given by Eq.
(38) (a static pulse far from the horizons), with parameters $a_{0}=1$,
$r^{*}_{c}=150\mu$ and $\sigma=2$. As $\kappa_{3}$ increases, a faster cascade
of high frequency modes grows towards the horizon (orange dotted vertical
line). This behavior can be compared with the corresponding evolution of the
wave equation (i.e., by setting $\kappa_{2}=\kappa_{3}=0$), represented by the
dotted black profile.

Figure 3: Dynamics of the Lifshitz field (ID type II). Evolution of the
$\ell=m=0$ mode from the ID type II (ingoing pulse) given in Eq. (39), with
parameters $r^{*}_{c}=150\mu$, $\sigma=3$ and $\omega=1$. We fixed
$\kappa_{2}=0.1$, and set $\kappa_{3}=0.01$ (light red curve) and
$\kappa_{3}=0.1$ (brown curve). A cascade towards the Killing horizon
(represented by the gray dotted vertical line) behaves similarly to the
previous case shown in Figure 2. The evolution of the wave equation is also
shown in dotted black, for comparison.

Figure 4: Lifshitz evolution of a static pulse, for different values of
$\kappa_{2}$ and fixing $\kappa_{3}=0.01$. The green curve corresponds to
$\kappa_{2}=0.1$, the blue one to $\kappa_{2}=0.5$, while the red curve
corresponds to $\kappa_{2}=1$. At early times, the propagation speed of the
slowest modes decreases when increasing $\kappa_{2}$. At late times, instead,
the bump at the Killing horizon (orange dotted vertical line) grows when
increasing $\kappa_{2}$. The evolution of the wave equation from the same
static pulse is shown in dotted black. The parameters of the initial static
pulse are the same as those considered in Figure 2.

Figure 5: Lifshitz evolution of an (approximately) ingoing pulse, for fixed
$\kappa_{3}=0.01$, $\kappa_{2}=0.1$ (blue curve) and $\kappa_{2}=1$ (red
curve). The dotted black profile corresponds to the solution of the wave
equation, for comparison. The initial data parameters are the same as those
considered in Figure 3.

###  V.2 Effect of the potential barrier

The coefficient multiplying the term which linear in $\psi_{\ell m}$ in (34)
acts as an effective centrifugal barrier for modes with $\ell\neq 0$. Its
explicit expression reads

$\displaystyle V_{\text{eff}}(r){}$
$\displaystyle=\frac{H^{2}\ell(\ell+1)}{32A^{9}r^{6}}\left(8A^{7}\left(r^{4}+\kappa_{2}\ell(\ell+1)r^{2}+\kappa_{3}\ell^{2}(\ell+1)^{2}\right)+4A^{5}H\left(rH^{\prime}\left(\kappa_{2}r^{2}+3\kappa_{3}\ell(\ell+1)\right)\right.\right.$
(42)
$\displaystyle\left.\left.-H\left(\kappa_{2}r^{2}+7\kappa_{3}\ell(\ell+1)\right)\right)+A^{3}H\kappa_{3}\left(12H^{3}-r^{3}\left(H^{\prime}\right)^{3}+Hr^{2}H^{\prime}\left(11H^{\prime}-4rH^{\prime\prime}\right)-H^{2}r\left(r^{2}H^{{}^{\prime\prime\prime}}\right.\right.\right.$
$\displaystyle\left.\left.\left.-6rH^{\prime\prime}+22H^{\prime}\right)\right)+15H^{4}\kappa_{3}r^{3}\left(A^{\prime}\right)^{3}-4A^{4}H^{2}rA^{\prime}\left(\kappa_{2}r^{2}+3\kappa_{3}\ell(\ell+1)\right)+AH^{3}\kappa_{3}r^{2}A^{\prime}\left(A^{\prime}\left(23H\right.\right.\right.$
$\displaystyle\left.\left.\left.-25rH^{\prime}\right)-10HrA^{\prime\prime}\right)+A^{2}H^{2}\kappa_{3}r\left(11r^{2}A^{\prime}\left(H^{\prime}\right)^{2}+Hr\left(7rA^{\prime}H^{\prime\prime}+H^{\prime}\left(7rA^{\prime\prime}-34A^{\prime}\right)\right)+H^{2}\left(r^{2}A^{{}^{\prime\prime\prime}}\right.\right.\right.$
$\displaystyle\left.\left.\left.-6rA^{\prime\prime}+22A^{\prime}\right)\right)\right),$

where the prime denotes differentiation with respect to $r$, and we have
omitted the argument in the functions $A(r)$ and $H(r)$ for the sake of
simplicity. Despite this highly non-linear expression, $V_{\rm eff}$ has a
single maximum for fixed $\ell$, in the vicinity of the Killing horizon. Its
value grows quickly with $\ell$, as can be noticed from the left panel of
figure 6.

The effect of such a potential in the dynamics of the Lifshitz field is
stronger for larger values of the angular momentum eigenvalue. This behavior
has been verified numerically, setting the coupling constants to
$\kappa_{2}=0.5$ and $\kappa_{3}=0.01$. The results are reported in the right
panel of figure 6. As can be seen, higher harmonics are pushed away by the
barrier, which prevents them from penetrating the Killing horizon at late
times. This suggests that the physics of the field within the interior region
is captured by the first modes.

Figure 6: Effect of the centrifugal barrier. Left panel: Plot of the effective
potential $V_{\rm eff}$ (in log-scale) as a function of $r^{*}$ (in units of
the black hole’s mass, $\mu$), for angular numbers $\ell=1$ (black), $\ell=3$
(dotted golden), $\ell=7$ (orange) and $\ell=10$ (dotted yellow). As $\ell$
increases, the maximum of the potential grows by several orders of magnitude.
Right panel: Evolution of different wave modes of the Lifshitz field, with the
same angular numbers and colors as in the left panel. For higher harmonics,
the maximum of the effective potential is larger, acting as a stronger
centrifugal barrier, thus preventing the corresponding modes from penetrating
the horizon.

###  V.3 Convergence tests

In order to validate our numerical code, we performed several convergence
tests, aiming to assess the corresponding accuracy order of the simulations. A
simple analytic calculation for Eq. (III) in flat space — which matches the
characteristic structure of the full equation at high frequencies –, shows
that the difference between an exact solution and the numerical approximation
using the implicit scheme (36)-(37), with second-order accurate finite-
difference operators, is given by

$|\phi_{\text{exa}}-\phi_{\text{num}}|=\alpha_{1}(\Delta
r^{*})^{2}+\alpha_{2}(\Delta t)^{2}+\alpha_{3}(\Delta r^{*})(\Delta
t)+\mathcal{O}(\Delta^{3}),$ (43)

where $\alpha_{1},\alpha_{2}$, and $\alpha_{3}$ are functions on the grid, and
$\mathcal{O}(\Delta^{3})$ denote all terms that are cubic in time and spatial
steps, like $(\Delta t)(\Delta r^{*})^{2}$. Thus, by rescaling $\Delta
r^{*}\to\lambda\Delta r^{*}$ and $\Delta t\to\lambda\Delta t$ for some
$\lambda\in\mathbb{R}$, one gets

$|\phi_{\text{exa}}-\phi_{\text{num}}|=\mathcal{O}(\lambda^{2}),$ (44)

which indicates that the scheme is second-order accurate.

In order to confirm that Eq. (44) approximately holds during the numerical
evolution, we perform three different runs with $\lambda=1$ (low resolution),
$\lambda=0.5$ (medium resolution) and $\lambda=0.25$ (high resolution). Then,
we define the ratio Press  _et al._ (2007)

$Q(t)=\frac{\lVert\phi_{\text{low}}-\phi_{\text{med}}\rVert}{\lVert\phi_{\text{med}}-\phi_{\text{high}}\rVert},$
(45)

which behaves as $Q(t)\sim 2^{p}$, with $p$ the accuracy order of the desired
numerical scheme (i.e. $p=2$ in our case). The result of these tests is shown
in Figure 7, where the value of $p$ is computed as a function of time for
different values of $\kappa_{3}$ (fixing $\kappa_{2}=0.1$). We can see that
the method is approximately second-order accurate, as expected, with
convergence improving for smaller values of $\kappa_{3}$.

We also explore the dependence of our results on the position of the ADL. For
doing so, we perform three runs with the same initial parameters, using the
initial data ID Type I, and setting $\kappa_{2}=0.1$, and $\kappa_{3}=0.01$.
We test three different positions of the left wall of the layer, namely
$x_{\rm L}=\\{-300\mu,-400\mu,-500\mu\\}$, while keeping the right wall fixed
far from the horizon, at $x_{\rm R}=300\mu$. Results are shown in figure 8. We
notice that there is a minimum distance between the left wall of the ADL and
the Killing horizon, above which the dynamical features of the field remain
almost unchanged. If the layer is close enough to the boundary of the physical
region of interest, the dynamics at long times is altered by its presence,
contrary to what happens if the layer is placed far away. Such a minimum
position depends on the final evolution time, with the layer having to be
placed further and further away if longer evolution times are required.
Nevertheless, the computational cost of placing the layer far from the
physical region – i.e., the number of grid points needed in order to keep the
same spatial resolution – is kept reasonable thanks to the implicit character
of the method, as larger time steps are allowed, unlike in standard explicit
methods.

Focusing on the particular simulations considered here (where we evolve until
$t\sim 100\mu$), this analysis suggests that choosing $x_{\rm L}\sim-300\mu$
is optimal. Although there may be small effects from the layer – no matter how
far it is placed –, the general features of the solution remain robust. This
allows us to draw solid conclusions on the features of the Lifshitz field
independently of the layer, as discussed above.

Figure 7: Convergence tests. Accuracy order of the implicit scheme, with
coupling parameters $\kappa_{2}=0.1$ and $\kappa_{3}=0.01$ (green continuous
curve); $\kappa_{3}=0.05$ (red dotted curve) and $\kappa_{3}=0.1$ (blue dotted
curve), as a function of time (in units of the black hole mass, $\mu$). The
simulations are performed from the initial data ID type I, with initial steps
$\Delta t=0.1$ and $\Delta r^{*}=0.3$, taking $\lambda=1$ (low resolution),
$\lambda=0.5$ (medium resolution) and $\lambda=0.25$ (high resolution).

Figure 8: Effect of the dissipative layer. Evolution of the Lifshitz field for
different positions of the left side of the ADL: $x_{L}=-300\mu$ (black);
$x_{L}=-400\mu$ (maroon) and $x_{L}=-500\mu$ (coral). At early times, the
evolution remains almost unchanged, until $t\sim 60\mu$, when small
differences start showing up. Nevertheless, the global behavior of the field
is the same in all cases, even at late times. Notice in particular the
appearance of the bump close to the Killing horizon. The simulations were run
with $\kappa_{2}=0.1$, $\kappa_{3}=0.01$, $\Delta t=0.1$, and $\Delta
r^{*}=0.125$, using the initial data ID type I. The position of the right side
of the layer is fixed at $x_{R}=300\mu$.

##  VI Conclusions

In this paper we have introduced an implicit numerical scheme that allows us
to solve evolution equations with an anisotropic scaling between time and
space of the form (2). Our approach is based on a generalization of the Crank-
Nicholson method for diffusion equations, replacing the usual discretization
of spatial derivatives by finite-difference operators with averages of them
evaluated in two consecutive time steps. This allows for evading the stringent
stability bound implied by the CFL condition for standard explicit methods,
removing obstructions to high-resolution and long-time evolutions. Our
implicit scheme is free from the CFL constraint, and can therefore be evolved
for substantially long times with small grid sizes. Let us highlight that this
method is equally valid for linear and non-linear equations, as long as they
feature some form of anisotropic scaling.

As an application and proof of concept of our numerical scheme, we have
studied the case of a Lifshitz scalar field in four dimensions, propagating in
a spherically symmetric and static black-hole space-time, solution to the
equations of motion of Hořava gravity at low energies. In contrast to the case
of fields with a relativistic dispersion relation, the Lifshitz scalar field
can probe the region enclosed by the Killing horizon, freely escaping from it.
Instead, it is the universal horizon, sitting at a smaller radius, that
represents the inner semi-permeable boundary for the propagation of the field
within the background geometry.

We have performed simulations with varying values of the higher derivative
couplings, $\kappa_{2}$ and $\kappa_{3}$, with two classes of initial data: a
static Gaussian pulse and an approximately ingoing wave-packet. Being unable
to disentangle ingoing and outgoing modes exactly due to the dispersive
character of the equations, we have implemented boundary conditions by
introducing an artificial dissipative layer, essentially absorbing the modes
far from the physical region of interest.

Our results show a consistent picture, where UV modes of the field solution
develop a cascade near the Killing horizon, growing stronger with larger
values of the couplings accompanying higher derivatives. At late times, this
cascade accumulates in the region between the Killing and universal horizons,
producing a bump in the amplitude, which grows exponentially. This may have
important implications for the fate of the universal horizon, which has so far
been studied only in the low energy limit of Hořava gravity. The effect of
higher derivatives on its structure and stability is therefore unknown so far.
Provided that one could model the scalar mode contained in the dynamical
degrees of freedom of Hořava gravity as a Lifshitz field, our results seem to
strongly indicate a linear instability of the universal horizon. We have also
shown that a large centrifugal barrier with peak outside the Killing horizon
develops for higher harmonics, so that only the first few modes can
effectively penetrate the Killing horizon.

Finally, we have shown that our methods retain an approximately second-order
convergence along time evolution, and that the effect of the dissipative layer
is negligible within the physical region, which allows us to trust on the
robustness of our results.

The research presented here constitutes a first step within a larger program.
Although here we have focused on the dynamics of the Lifshitz scalar field,
which is linear, our implicit method is equally suitable for non-linear
equations. In particular, this seems the way to approach the non-linear and
non-perturbative evolution of the fully general equations of motion of Hořava
gravity, which also present a dispersive character with anisotropic scaling.
We hope that our method can provide stable simulations even in this more
challenging situation.

## Acknowledgements

We are grateful to Stefano Liberati, Shinji Mukohyama, Marc Schneider, and
Toby Wiseman for discussions throughout the realization of this work. M. H-V.
wants to thank the APP group at SISSA for their hospitality during the
completion of this work. M. R., A. K., M. B. and E. B. acknowledge support
from the European Union’s H2020 ERC Consolidator Grant “GRavity from
Astrophysical to Microscopic Scales” (Grant No. GRAMS-815673) and the EU
Horizon 2020 Research and Innovation Programme under the Marie Sklodowska-
Curie (Grant Agreement No. 101007855). The work of M. H-V has been supported
by the Spanish State Research Agency MCIN/AEI/10.13039/501100011033 and the EU
NextGenerationEU/PRTR funds, under grant IJC2020-045126-I; and by the
Departament de Recerca i Universitats de la Generalitat de Catalunya, Grant No
2021 SGR 00649. IFAE is partially funded by the CERCA program of the
Generalitat de Catalunya.

##  Appendix A Black holes in Einstein-Æther and Hořava Gravity

As discussed in the main text, the background space-time (21) is solution to
the equations of motion of Einstein-Æther gravity Jacobson and Mattingly
(2001), with action

$\displaystyle S_{\rm EA}=\frac{1}{16\pi
G}\int\mbox{d}^{4}x\sqrt{|g|}\left(-R+{\cal L}_{\rm
U}+\eta(U_{\mu}U^{\mu}-1)\right),$ (46)

where $R$ is the Ricci scalar, $\eta$ is a Lagrange multiplier implementing
the unit norm condition of the æther, and ${\cal L}_{\rm
U}=K^{\alpha\beta}_{\mu\nu}\nabla_{\alpha}U^{\mu}\nabla_{\beta}U^{\nu}$, with
$K^{\alpha\beta}_{\mu\nu}=c_{1}g^{\alpha\beta}g_{\mu\nu}+c_{2}\delta^{\alpha}_{\mu}\delta^{\beta}_{\nu}+c_{3}\delta^{\alpha}_{\nu}\delta^{\beta}_{\mu}+c_{4}U^{\alpha}U^{\beta}g_{\mu\nu}$
and couplings $c_{i}\in\mathbb{R}$.

At the two derivative level, this is the most general covariant action that
leads to Lorentz violations and modified dispersion relations for all degrees
of freedom in the gravitational action – the usual graviton, plus vector and
scalar modes contained in $U^{\mu}$ Jacobson and Mattingly (2004) –, as well
as for generic matter fields coupling to them – perhaps through
renormalization group flow in the latter case. As such, it includes specific
models as particular cases of the coefficients $c_{i}$. If one imposes the
æther to be hypersurface orthogonal at the level of the action, and identifies
$c_{1}=c_{3}$, $c_{2}=(1-\beta^{\prime})(\lambda-1)-\beta^{\prime}$,
$\beta^{\prime}\equiv c_{1}+c_{3}=1-\eta^{-1}$, and $c_{1}+c_{4}=\alpha$, (46)
reproduces the low energy limit of Hořava Gravity, also known as khronometric
gravity Blas  _et al._ (2011) – fully described by three couplings
$(\alpha,\lambda,\eta)$.

Spherically symmetric and static solutions to (46) can be searched for with
the ansatz Berglund  _et al._ (2012); Barausse  _et al._ (2011)

$\displaystyle{\rm d}s^{2}=f(r){\rm d}\tau^{2}-\frac{B(r)^{2}}{f(r)}{\rm
d}r^{2}-r^{2}{\rm d}\Omega^{2},\quad U_{\mu}{\rm
d}x^{\mu}=\frac{H(r)}{2A(r)}{\rm
d}t+\frac{B(r)(1-f(r)A(r)^{2})}{2A(r)f(r)}{\rm d}r,$ (47)

where $H(r)=1+f(r)A(r)^{2}$, and the form of $U^{\mu}$ is chosen for
convenience, automatically satisfying $U_{\mu}U^{\mu}=1$. Accidentally,
spherical symmetry and staticity automatically impose hypersurface
orthogonality, and hence all solutions to Einstein-Æther gravity with these
isometries are solutions to Hořava gravity as well (and vice versa). Hence,
both theories are generally studied together when discussing their features
and phenomenology in static spherically symmetric configurations. This is not
true when any of the conditions above are relaxed. In particular, axisymmetric
solutions in Einstein-Æther gravity are not hypersurface orthogonal Barausse
_et al._ (2016); Adam  _et al._ (2022). We also do not expect them to be
equivalent when higher derivatives are included in the gravitational action,
even in the spherically symmetric and static case.

The explicit form of the functions $f(r),B(r)$, and $A(r)$ must be found
solving the equations of motion explicitly. Although numerical solutions can
always be achieved at any generic point of the parameter space Barausse  _et
al._ (2011), there exist two corners where analytic solutions can be attained
Berglund  _et al._ (2012). These correspond to $c_{14}=0$ and $c_{123}=0$,
where we are using the notation $c_{ij\dots k}=c_{i}+c_{j}+\dots+c_{k}$. The
latter case, although leading to simpler functions, is incompatible with
observational bounds constraining the parameter space of the theory Gupta  _et
al._ (2021). The former case, instead, is perfectly compatible with current
bounds and leads to

$\displaystyle f(r)=1-\frac{2\mu}{r}-c_{13}\frac{r_{\text{\ae}}^{4}}{r^{4}},$
(48) $\displaystyle B(r)=1,$ (49) $\displaystyle
A(r)=\frac{1}{f(r)}\left(-\frac{r_{\text{\ae}}^{2}}{r^{2}}+\sqrt{f(r)+\frac{r_{\text{\ae}}^{4}}{r^{4}}}\right).$
(50)

The parameter $r_{\text{\ae}}$ is in principle arbitrary. However, the
solution displays a singularity at the universal horizon, unless

$\displaystyle
r_{\text{\ae}}=\frac{\mu}{2}\left(\frac{27}{1-c_{13}}\right)^{1/4}.$ (51)

This value is therefore chosen to ensure regularity of the solution everywhere
except for the central singularity, sitting at $r=0$.

At high energies, Hořava gravity departs from action (46) by higher derivative
spatial terms Blas  _et al._ (2011). The contribution of these operators –
which are critical for achieving power-counting renormalizability –, and in
particular their backreaction onto space-times of the form (21), are unknown.
Although here we assume that such effects are small and that gravitational
perturbations are negligible, we leave this as an open question for the
future.

##  Appendix B Coefficients of the evolution equation

In this appendix, we explicitly give the form of the coefficients in Eq. (35).

$\displaystyle\zeta_{01}={}$
$\displaystyle\frac{1}{64A^{11}Hr^{5}}\left\\{16\left(6\kappa_{3}\left(\ell^{2}+\ell\right)^{2}H^{3}+r\left(((6-H)H-8)r^{4}+2H^{2}\kappa_{2}\ell(\ell+1)r^{2}+3H^{2}\kappa_{3}\left(\ell^{2}+\ell\right)^{2}\right)H^{\prime}\right)A^{9}\right.$
$\displaystyle\left.-16Hr\left(4Hr^{4}-4r^{4}+H^{2}\left(-r^{4}+2\kappa_{2}\ell(\ell+1)r^{2}+3\kappa_{3}\left(\ell^{2}+\ell\right)^{2}\right)\right)A^{\prime}A^{8}-4H^{2}\left(12\kappa_{3}\ell(\ell+1)H^{3}\right.\right.$
$\displaystyle\left.\left.+r\left(\left(2\kappa_{2}r^{2}-8\kappa_{3}\ell(\ell+1)\right)H^{\prime}+r\left(2\left(\kappa_{2}r^{2}+\kappa_{3}\ell(\ell+1)\right)H^{\prime\prime}+r\left(\kappa_{2}r^{2}+3\kappa_{3}\ell(\ell+1)\right)H^{(3)}(r)\right)\right)H^{2}\right.\right.$
$\displaystyle\left.\left.-2r^{3}H^{\prime}\left(\kappa_{2}rH^{\prime}+\left(\kappa_{2}r^{2}+3\kappa_{3}\ell(\ell+1)\right)H^{\prime\prime}\right)H+r^{3}\left(\kappa_{2}r^{2}+3\kappa_{3}\ell(\ell+1)\right)\left(H^{\prime}\right)^{3}\right)A^{7}+4H^{3}r\left(\left(\left(2\kappa_{2}r^{2}\right.\right.\right.\right.$
$\displaystyle\left.\left.\left.\left.-8\kappa_{3}\ell(\ell+1)\right)A^{\prime}+r\left(2\left(\kappa_{2}r^{2}+\kappa_{3}\ell(\ell+1)\right)A^{\prime\prime}+r\left(\kappa_{2}r^{2}+3\kappa_{3}\ell(\ell+1)\right)A^{(3)}(r)\right)\right)H^{2}+r\left(r\left(\kappa_{2}r^{2}\right.\right.\right.\right.$
$\displaystyle\left.\left.\left.\left.+3\kappa_{3}\ell(\ell+1)\right)H^{\prime}A^{\prime\prime}+A^{\prime}\left(4\kappa_{3}\ell(\ell+1)H^{\prime}+r\left(\kappa_{2}r^{2}+3\kappa_{3}\ell(\ell+1)\right)H^{\prime\prime}\right)\right)H-r^{2}\left(\kappa_{2}r^{2}\right.\right.\right.$
$\displaystyle\left.\left.\left.+3\kappa_{3}\ell(\ell+1)\right)A^{\prime}\left(H^{\prime}\right)^{2}\right)A^{6}+H^{2}r\left(\kappa_{3}r^{4}\left(H^{\prime}\right)^{5}-2H\kappa_{3}r^{3}\left(H^{\prime}+2rH^{\prime\prime}\right)\left(H^{\prime}\right)^{3}+2H^{2}r^{2}\left(\kappa_{3}\left(2\left(H^{\prime}\right)^{2}\right.\right.\right.\right.$
$\displaystyle\left.\left.\left.\left.+r\left(2H^{\prime\prime}+rH^{(3)}(r)\right)H^{\prime}+2r^{2}\left(H^{\prime\prime}\right)^{2}\right)-2\left(\kappa_{2}r^{2}+3\kappa_{3}\ell(\ell+1)\right)\left(A^{\prime}\right)^{2}\right)H^{\prime}+H^{3}r\left(-8\left(\kappa_{2}r^{2}\right.\right.\right.\right.$
$\displaystyle\left.\left.\left.\left.+2\kappa_{3}\ell(\ell+1)\right)\left(A^{\prime}\right)^{2}-16r\left(\kappa_{2}r^{2}+3\kappa_{3}\ell(\ell+1)\right)A^{\prime\prime}A^{\prime}+\kappa_{3}\left(-H^{\prime\prime}\left(6H^{\prime\prime}+5rH^{(3)}(r)\right)r^{2}-8\left(H^{\prime}\right)^{2}\right.\right.\right.\right.$
$\displaystyle\left.\left.\left.\left.+H^{\prime}\left(H^{(4)}(r)r^{3}+6H^{\prime\prime}r\right)\right)\right)+H^{4}\kappa_{3}\left(4H^{\prime}+r\left(r\left(H^{(5)}(r)r^{2}+4H^{(4)}(r)r+2H^{(3)}(r)\right)-4H^{\prime\prime}\right)\right)\right)A^{5}\right.$
$\displaystyle\left.-H^{3}r\left(\kappa_{3}\left(4A^{\prime}+r\left(A^{(5)}(r)r^{3}+4A^{(4)}(r)r^{2}+2A^{(3)}(r)r-4A^{\prime\prime}\right)\right)H^{4}+\kappa_{3}r\left(r\left(r\left(H^{\prime\prime}\left(12A^{\prime\prime}+5rA^{(3)}(r)\right)\right.\right.\right.\right.\right.$
$\displaystyle\left.\left.\left.\left.\left.+5rA^{\prime\prime}H^{(3)}(r)\right)+2H^{\prime}\left(3A^{(4)}(r)r^{2}+8A^{(3)}(r)r+6A^{\prime\prime}\right)\right)+2A^{\prime}\left(r\left(3H^{(4)}(r)r^{2}+8H^{(3)}(r)r+6H^{\prime\prime}\right)\right.\right.\right.\right.$
$\displaystyle\left.\left.\left.\left.-12H^{\prime}\right)\right)H^{3}-r^{2}\left(12\left(\kappa_{2}r^{2}+3\kappa_{3}\ell(\ell+1)\right)\left(A^{\prime}\right)^{3}+\kappa_{3}\left(-24\left(H^{\prime}\right)^{2}+2r\left(8H^{\prime\prime}+rH^{(3)}(r)\right)H^{\prime}\right.\right.\right.\right.$
$\displaystyle\left.\left.\left.\left.+11r^{2}\left(H^{\prime\prime}\right)^{2}\right)A^{\prime}+\kappa_{3}rH^{\prime}\left(rA^{\prime\prime}H^{\prime\prime}-2H^{\prime}\left(2A^{\prime\prime}+3rA^{(3)}(r)\right)\right)\right)H^{2}+2\kappa_{3}r^{4}\left(H^{\prime}\right)^{2}\left(H^{\prime}A^{\prime\prime}+5A^{\prime}H^{\prime\prime}\right)H\right.\right.$
$\displaystyle\left.\left.-3\kappa_{3}r^{4}A^{\prime}\left(H^{\prime}\right)^{4}\right)A^{4}+H^{4}\kappa_{3}r^{2}\left(\left(3A^{\prime\prime}\left(6A^{\prime\prime}+5rA^{(3)}(r)\right)r^{2}+A^{\prime}\left(11A^{(4)}(r)r^{2}+32A^{(3)}(r)r+18A^{\prime\prime}\right)r\right.\right.\right.$
$\displaystyle\left.\left.\left.-16\left(A^{\prime}\right)^{2}\right)H^{3}+r\left(4\left(4H^{(3)}(r)r^{2}+7H^{\prime\prime}r+12H^{\prime}\right)\left(A^{\prime}\right)^{2}+r\left(A^{\prime\prime}\left(80H^{\prime}+29rH^{\prime\prime}\right)+46rH^{\prime}A^{(3)}(r)\right)A^{\prime}\right.\right.\right.$
$\displaystyle\left.\left.\left.+25r^{2}H^{\prime}\left(A^{\prime\prime}\right)^{2}\right)H^{2}-2r^{2}A^{\prime}H^{\prime}\left(A^{\prime}\left(6H^{\prime}+8rH^{\prime\prime}\right)-11rH^{\prime}A^{\prime\prime}\right)H+6r^{3}\left(A^{\prime}\right)^{2}\left(H^{\prime}\right)^{3}\right)A^{3}\right.$
$\displaystyle\left.-2H^{5}\kappa_{3}r^{3}A^{\prime}\left(\left(14\left(A^{\prime}\right)^{2}+2r\left(31A^{\prime\prime}+15rA^{(3)}(r)\right)A^{\prime}+35r^{2}\left(A^{\prime\prime}\right)^{2}\right)H^{2}+rA^{\prime}\left(85rH^{\prime}A^{\prime\prime}+A^{\prime}\left(32H^{\prime}\right.\right.\right.\right.$
$\displaystyle\left.\left.\left.\left.+15rH^{\prime\prime}\right)\right)H+5r^{2}\left(A^{\prime}\right)^{2}\left(H^{\prime}\right)^{2}\right)A^{2}+3H^{6}\kappa_{3}r^{4}\left(A^{\prime}\right)^{3}\left(A^{\prime}\left(26H+35rH^{\prime}\right)+70HrA^{\prime\prime}\right)A-105H^{7}\kappa_{3}r^{5}\left(A^{\prime}\right)^{5}\right\\}$

$\zeta_{10}=\frac{(H-2)HrA^{\prime}+A(2-H)rH^{\prime}}{4A^{3}r},\qquad\zeta_{11}=2-\frac{4}{H}$

$\displaystyle\zeta_{02}={}$
$\displaystyle\frac{1}{16A^{8}H^{2}r^{4}}\left\\{-16A^{8}\left(H^{2}\left(2\kappa_{2}\ell(\ell+1)r^{2}+3\kappa_{3}\left(\ell^{2}+\ell\right)^{2}-r^{4}\right)+4Hr^{4}-4r^{4}\right)+4A^{6}H^{2}\left(6H^{2}\kappa_{3}\ell(\ell+1)\right.\right.$
$\displaystyle\left.\left.+7r^{2}\left(H^{\prime}\right)^{2}\left(\kappa_{2}r^{2}+3\kappa_{3}\ell(\ell+1)\right)-2Hr\left(2rH^{\prime\prime}\left(\kappa_{2}r^{2}+3\kappa_{3}\ell(\ell+1)\right)+H^{\prime}\left(4\kappa_{2}r^{2}+\kappa_{3}\ell(\ell+1)\right)\right)\right)\right.$
$\displaystyle\left.+35H^{6}\kappa_{3}r^{4}\left(A^{\prime}\right)^{4}-A^{4}H^{2}r\left(H^{2}r\left(4\left(A^{\prime}\right)^{2}\left(\kappa_{2}r^{2}+3\kappa_{3}\ell(\ell+1)\right)+\kappa_{3}\left(22r^{2}\left(H^{\prime\prime}\right)^{2}+6\left(H^{\prime}\right)^{2}\right.\right.\right.\right.$
$\displaystyle\left.\left.\left.\left.+rH^{\prime}\left(21rH^{(3)}(r)+62H^{\prime\prime}\right)\right)\right)+31\kappa_{3}r^{3}\left(H^{\prime}\right)^{4}-4H\kappa_{3}r^{2}\left(H^{\prime}\right)^{2}\left(17rH^{\prime\prime}+10H^{\prime}\right)+2H^{3}\kappa_{3}\left(2H^{\prime}\right.\right.\right.$
$\displaystyle\left.\left.\left.-r\left(3r^{2}H^{(4)}(r)+11rH^{(3)}(r)+2H^{\prime\prime}\right)\right)\right)-2AH^{6}\kappa_{3}r^{3}\left(A^{\prime}\right)^{2}\left(35rA^{\prime\prime}+24A^{\prime}\right)+8A^{5}H^{3}r\left(H\left(2rA^{\prime\prime}\left(\kappa_{2}r^{2}\right.\right.\right.\right.$
$\displaystyle\left.\left.\left.\left.+3\kappa_{3}\ell(\ell+1)\right)+A^{\prime}\left(4\kappa_{2}r^{2}+\kappa_{3}\ell(\ell+1)\right)\right)-3rA^{\prime}H^{\prime}\left(\kappa_{2}r^{2}+3\kappa_{3}\ell(\ell+1)\right)\right)+A^{2}H^{4}\kappa_{3}r^{2}\left(8r^{2}\left(A^{\prime}\right)^{2}\left(H^{\prime}\right)^{2}\right.\right.$
$\displaystyle\left.\left.+HrA^{\prime}\left(3rA^{\prime\prime}H^{\prime}+A^{\prime}\left(4H^{\prime}-11rH^{\prime\prime}\right)\right)+H^{2}\left(14r^{2}\left(A^{\prime\prime}\right)^{2}+2\left(A^{\prime}\right)^{2}+rA^{\prime}\left(27rA^{(3)}(r)+70A^{\prime\prime}\right)\right)\right)\right.$
$\displaystyle\left.-A^{3}H^{3}\kappa_{3}r\left(12r^{3}A^{\prime}\left(H^{\prime}\right)^{3}+Hr^{2}H^{\prime}\left(5rA^{\prime\prime}H^{\prime}-A^{\prime}\left(15rH^{\prime\prime}+4H^{\prime}\right)\right)+H^{2}r\left(A^{\prime}\left(3r^{2}H^{(3)}(r)+4rH^{\prime\prime}\right.\right.\right.\right.$
$\displaystyle\left.\left.\left.\left.-4H^{\prime}\right)+r\left(3rA^{(3)}(r)H^{\prime}+4A^{\prime\prime}\left(H^{\prime}-2rH^{\prime\prime}\right)\right)\right)+H^{3}\left(6r^{3}A^{(4)}(r)+22r^{2}A^{(3)}(r)+4rA^{\prime\prime}-4A^{\prime}\right)\right)\right\\}$

$\displaystyle\zeta_{03}=\frac{1}{4A^{3}Hr^{2}}{}$
$\displaystyle\left\\{8A^{3}\left(2H\textit{k}_{2}r-3H^{\prime}\left(\textit{k}_{2}r^{2}+3\textit{k}_{3}\ell(\ell+1)\right)\right)+24A^{2}HA^{\prime}\left(\textit{k}_{2}r^{2}+3\textit{k}_{3}\ell(\ell+1)\right)\right.$
$\displaystyle\left.+H\textit{k}_{3}\left(-90r^{2}A^{\prime}\left(H^{\prime}\right)^{2}+Hr\left(45rA^{\prime}H^{\prime\prime}+H^{\prime}\left(45rA^{\prime\prime}+92A^{\prime}\right)\right)-\left(H^{2}\left(15r^{2}A^{(3)}+46rA^{\prime\prime}+6A^{\prime}\right)\right)\right)\right.$
$\displaystyle\left.+A\textit{k}_{3}\left(90r^{2}\left(H^{\prime}\right)^{3}-2HrH^{\prime}\left(45rH^{\prime\prime}+46H^{\prime}\right)+H^{2}\left(15r^{2}H^{(3)}+46rH^{\prime\prime}+6H^{\prime}\right)\right)\right\\},$

$\displaystyle\zeta_{04}=\frac{1}{A^{2}H^{2}r^{2}}{}$
$\displaystyle\left\\{4A^{4}\left(\textit{k}_{2}r^{2}+3\textit{k}_{3}\ell(\ell+1)\right)+A^{2}\textit{k}_{3}r\left(4H\left(5rH^{\prime\prime}+12H^{\prime}\right)-65r\left(H^{\prime}\right)^{2}\right)-25H^{2}\textit{k}_{3}r^{2}\left(A^{\prime}\right)^{2}\right.$
$\displaystyle\left.-2AH\textit{k}_{3}r\left(10HrA^{\prime\prime}+3A^{\prime}\left(8H-15rH^{\prime}\right)\right)\right\\},$

$\zeta_{05}=-\frac{12A\textit{k}_{3}\left(5HrA^{\prime}+A\left(2H-5rH^{\prime}\right)\right)}{H^{3}r},\qquad\zeta_{06}=-\frac{16\textit{k}_{3}A^{4}}{H^{4}},$

where the prime denotes derivative with respect to $r$, and the coefficients
are functions of $r^{*}$ through the transformation $r(r^{*})$ given in (32).

## References

  * Horava (2009) P. Horava, Phys. Rev. D 79, 084008 (2009), arXiv:0901.3775 [hep-th] .
  * Fujimori  _et al._ (2016) T. Fujimori, T. Inami, K. Izumi, and T. Kitamura, PTEP 2016, 013B08 (2016), arXiv:1510.07237 [hep-th] .
  * Anselmi and Halat (2007) D. Anselmi and M. Halat, Phys. Rev. D 76, 125011 (2007), arXiv:0707.2480 [hep-th] .
  * Blas  _et al._ (2011) D. Blas, O. Pujolas, and S. Sibiryakov, JHEP 04, 018 (2011), arXiv:1007.3503 [hep-th] .
  * Blas  _et al._ (2010) D. Blas, O. Pujolas, and S. Sibiryakov, Phys. Rev. Lett. 104, 181302 (2010), arXiv:0909.3525 [hep-th] .
  * Donnelly and Jacobson (2011) W. Donnelly and T. Jacobson, Phys. Rev. D 84, 104019 (2011), arXiv:1106.2131 [hep-th] .
  * Bellorin  _et al._ (2022) J. Bellorin, C. Borquez, and B. Droguett, Phys. Rev. D 106, 044055 (2022), arXiv:2207.08938 [hep-th] .
  * Barvinsky  _et al._ (2016) A. O. Barvinsky, D. Blas, M. Herrero-Valea, S. M. Sibiryakov, and C. F. Steinwachs, Phys. Rev. D 93, 064022 (2016), arXiv:1512.02250 [hep-th] .
  * Barvinsky  _et al._ (2017a) A. O. Barvinsky, D. Blas, M. Herrero-Valea, S. M. Sibiryakov, and C. F. Steinwachs, Phys. Rev. Lett. 119, 211301 (2017a), arXiv:1706.06809 [hep-th] .
  * Barvinsky  _et al._ (2019) A. O. Barvinsky, M. Herrero-Valea, and S. M. Sibiryakov, Phys. Rev. D 100, 026012 (2019), arXiv:1905.03798 [hep-th] .
  * Barvinsky  _et al._ (2022) A. O. Barvinsky, A. V. Kurov, and S. M. Sibiryakov, Phys. Rev. D 105, 044009 (2022), arXiv:2110.14688 [hep-th] .
  * Radkovski and Sibiryakov (2023) J. I. Radkovski and S. M. Sibiryakov, (2023), arXiv:2306.00102 [hep-th] .
  * Griffin  _et al._ (2017) T. Griffin, K. T. Grosvenor, C. M. Melby-Thompson, and Z. Yan, JHEP 06, 004 (2017), arXiv:1701.08173 [hep-th] .
  * Benedetti and Guarnieri (2014) D. Benedetti and F. Guarnieri, JHEP 03, 078 (2014), arXiv:1311.6253 [hep-th] .
  * Barausse  _et al._ (2011) E. Barausse, T. Jacobson, and T. P. Sotiriou, Phys. Rev. D 83, 124043 (2011), arXiv:1104.2889 [gr-qc] .
  * Barausse and Sotiriou (2013a) E. Barausse and T. P. Sotiriou, Phys. Rev. D 87, 087504 (2013a), arXiv:1212.1334 [gr-qc] .
  * Barausse and Sotiriou (2013b) E. Barausse and T. P. Sotiriou, Class. Quant. Grav. 30, 244010 (2013b), arXiv:1307.3359 [gr-qc] .
  * Berglund  _et al._ (2012) P. Berglund, J. Bhattacharyya, and D. Mattingly, Phys. Rev. D 85, 124019 (2012), arXiv:1202.4497 [hep-th] .
  * Lara  _et al._ (2021) G. Lara, M. Herrero-Valea, E. Barausse, and S. M. Sibiryakov, Phys. Rev. D 103, 104007 (2021), arXiv:2103.01975 [gr-qc] .
  * Mazza and Liberati (2023) J. Mazza and S. Liberati, JHEP 03, 199 (2023), arXiv:2301.04697 [gr-qc] .
  * Mukohyama (2010) S. Mukohyama, Class. Quant. Grav. 27, 223101 (2010), arXiv:1007.5199 [hep-th] .
  * Pospelov and Shang (2012) M. Pospelov and Y. Shang, Phys. Rev. D 85, 105001 (2012), arXiv:1010.5249 [hep-th] .
  * Blas and Lim (2015) D. Blas and E. Lim, Int. J. Mod. Phys. D 23, 1443009 (2015), arXiv:1412.4828 [gr-qc] .
  * Audren  _et al._ (2015) B. Audren, D. Blas, M. M. Ivanov, J. Lesgourgues, and S. Sibiryakov, JCAP 03, 016 (2015), arXiv:1410.6514 [astro-ph.CO] .
  * Barausse (2019) E. Barausse, Phys. Rev. D 100, 084053 (2019), [Erratum: Phys.Rev.D 104, 069903 (2021)], arXiv:1907.05958 [gr-qc] .
  * Bettoni  _et al._ (2017) D. Bettoni, A. Nusser, D. Blas, and S. Sibiryakov, JCAP 05, 024 (2017), arXiv:1702.07726 [astro-ph.CO] .
  * Bonetti and Barausse (2015) M. Bonetti and E. Barausse, Phys. Rev. D 91, 084053 (2015), [Erratum: Phys.Rev.D 93, 029901 (2016)], arXiv:1502.05554 [gr-qc] .
  * Cornish  _et al._ (2017) N. Cornish, D. Blas, and G. Nardini, Phys. Rev. Lett. 119, 161102 (2017), arXiv:1707.06101 [gr-qc] .
  * Emir Gümrükçüoğlu  _et al._ (2018) A. Emir Gümrükçüoğlu, M. Saravani, and T. P. Sotiriou, Phys. Rev. D 97, 024032 (2018), arXiv:1711.08845 [gr-qc] .
  * Yagi  _et al._ (2014a) K. Yagi, D. Blas, E. Barausse, and N. Yunes, Phys. Rev. D 89, 084067 (2014a), [Erratum: Phys.Rev.D 90, 069902 (2014), Erratum: Phys.Rev.D 90, 069901 (2014)], arXiv:1311.7144 [gr-qc] .
  * Yagi  _et al._ (2014b) K. Yagi, D. Blas, N. Yunes, and E. Barausse, Phys. Rev. Lett. 112, 161101 (2014b), arXiv:1307.6219 [gr-qc] .
  * Franchini  _et al._ (2021) N. Franchini, M. Herrero-Valea, and E. Barausse, Phys. Rev. D 103, 084012 (2021), arXiv:2103.00929 [gr-qc] .
  * Gupta  _et al._ (2021) T. Gupta, M. Herrero-Valea, D. Blas, E. Barausse, N. Cornish, K. Yagi, and N. Yunes, Class. Quant. Grav. 38, 195003 (2021), arXiv:2104.04596 [gr-qc] .
  * Oshita  _et al._ (2021) N. Oshita, N. Afshordi, and S. Mukohyama, JCAP 05, 005 (2021), arXiv:2102.01741 [gr-qc] .
  * Press  _et al._ (2007) W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery,  _Numerical Recipes 3rd Edition: The Art of Scientific Computing_, 3rd ed. (Cambridge University Press, 2007).
  * Ardonne  _et al._ (2004) E. Ardonne, P. Fendley, and E. Fradkin, Annals Phys. 310, 493 (2004), arXiv:cond-mat/0311466 .
  * Jacobson and Mattingly (2001) T. Jacobson and D. Mattingly, Phys. Rev. D 64, 024028 (2001), arXiv:gr-qc/0007031 .
  * Barvinsky  _et al._ (2017b) A. O. Barvinsky, D. Blas, M. Herrero-Valea, D. V. Nesterov, G. Pérez-Nadal, and C. F. Steinwachs, JHEP 06, 063 (2017b), arXiv:1703.04747 [hep-th] .
  * (39) A. D. Kovács, forthcoming .
  * Jacobson (2014) T. Jacobson, Phys. Rev. D 89, 081501 (2014), arXiv:1310.5115 [gr-qc] .
  * Arnowitt  _et al._ (1959) R. L. Arnowitt, S. Deser, and C. W. Misner, Phys. Rev. 116, 1322 (1959).
  * Blas and Sibiryakov (2011) D. Blas and S. Sibiryakov, Phys. Rev. D 84, 124043 (2011), arXiv:1110.2195 [hep-th] .
  * Bhattacharyya  _et al._ (2016) J. Bhattacharyya, M. Colombo, and T. P. Sotiriou, Class. Quant. Grav. 33, 235003 (2016), arXiv:1509.01558 [gr-qc] .
  * Cropp  _et al._ (2014) B. Cropp, S. Liberati, A. Mohd, and M. Visser, Phys. Rev. D 89, 064061 (2014), arXiv:1312.0405 [gr-qc] .
  * Dima and Barausse (2020) A. Dima and E. Barausse, Classical and Quantum Gravity 37, 175006 (2020), arXiv:1007.3503 .
  * Dolan (2013) S. R. Dolan, Phys. Rev. D 87, 124026 (2013), arXiv:1212.1477 .
  * (47) M. Herrero-Valea and M. E. Rubio, in preparation .
  * Jacobson and Mattingly (2004) T. Jacobson and D. Mattingly, Phys. Rev. D 70, 024003 (2004), arXiv:gr-qc/0402005 .
  * Barausse  _et al._ (2016) E. Barausse, T. P. Sotiriou, and I. Vega, Phys. Rev. D 93, 044044 (2016), arXiv:1512.05894 [gr-qc] .
  * Adam  _et al._ (2022) A. Adam, P. Figueras, T. Jacobson, and T. Wiseman, Class. Quant. Grav. 39, 125001 (2022), arXiv:2108.00005 [gr-qc] .


# Evolution of choices over time: The U.S. Presidential election 2012 and the
NY City Mayoral Election, 2013

Mukkai Krishnamoorthy, Wesley Miller
Rensselaer Polytechnic Institute, Troy, NY  Raju Krishnamoorthy
Columbia University, New York, NY

###### Abstract

We conducted surveys before and after the 2012 U.S. Presidential election and
prior to the NY City Mayoral election in 2013. The surveys were done using
Amazon Turk. This poster describes the results of our analysis of the surveys
and predicts the winner of the NY City Mayoral Election.

##  1 Introduction

Several authors have published analyses of poll and survey data of the U.S.
Presidential election 2012 [1][2], with vastly different predictions. Both
sampling differences in voter populations and differences in methods of
analysis may have accounted for the difference in the prediction. In the
current experiment, we perform a survey and an indirect analysis, using Amazon
Turk, as opposed to analysis of direct voter preferences using Amazon Turk
[4][3]. We have focused on voters’ two most pressing concerns, one national
and the other international, from a choice of 5 in each category.

A lot of data analysis of the presidential election in 2012 has been done. In
particular, please see [1]. Other data analysis [2] did not fare that well.
One of the reasons for such a difference could be the sampling of voter
population. Another reason could be due to data analysis. Please see these two
other papers for using Amazon Turk [4][3] in Election Prediction. As opposed
tp these two papers, we do an indirect analysis - not measuring direct voter
preferences.The perspective voters were asked to give on a pressing national
and international concern. Our approach is to do a seconday analysis,
concentrating more on the most important issue (one of a national concern and
the other of international concern among choices of 5 in each category) that
the surveyed people had in their minds.

These are the limitations on our data:

  * •

Did all those surveyed for the presidential election 2012, vote?

  * •

Were the same people surveyed before and after the election?

  * •

Will those surveyed for the NY City Mayoral Primary vote?

  * •

Sample size is extremely small, compared to the voting population (200 for the
US Presidential election and 100 for the NY City Mayoral)

Even with these limitations, we feel that it exhibits interesting patterns.

##  2 Data Collection

The data was collected through Amazon Mechanical Turk. Our pool of
participants was people of voting age in the United States, the voting
population. We selected 200 people before the election and asked each one of
them which domestic issue would affect his/her voting decision most when
considering the economy, healthcare, tax reform, education, and national
security. We then asked each one of them which international issue would have
the greatest effect on his/her voting decision when considering withdrawal
from wars, international security, global trade issues, more active
engagement, and international partners. Once the election was over and
President Obama had been re-elected, we posted the same survey for 200 more
people, asking this time which of the issues affected their decisions the most
when they voted. The third set of data was collected from a pool of 100
participants and asked the same set of questions (that we asked during the
presidential election).

##  3  Amazon Turk Experimental Results

Presidential Before
---
| National | International
---|---
| The Economy | 145
---|---
Healthcare | 30
Education | 11
Tax Reform | 12
National Security | 2
| Withdrawal from Wars | 109
---|---
International Security | 39
Global Trade Issues | 40
More Active Engagement | 6
International Partners | 6
Presidential After
| National | International
---|---
| The Economy | 147
---|---
Healthcare | 36
Education | 9
Tax Reform | 7
National Security | 1
| Withdrawal from Wars | 105
---|---
International Security | 43
Global Trade Issues | 42
More Active Engagement | 3
International Partners | 7
Mayoral Before
| City | National
---|---
| The Economy | 71
---|---
Healthcare | 15
Education | 7
Tax Reform | 7
City Security | 0
| Withdrawal from Wars | 39
---|---
National Security | 29
Global Trade Issues | 15
More Active Engagement | 11
State Partners | 6

##  4 Observations

The results of ordering all the surveys (shown in the previous section) remain
the same. This tells us about the general anxiety level of the people. Even
aftera gap of one year, rankings (of concerns) did not change. Even though our
sample size is small, we believe in the authenticity of data because the
survey after the election did not vary even slightly from those taken before
the election. What is more striking is the approximate percentage of people
prefering the choices (in almost all three cases). There is a small change of
percentages of healthcare and International Security (before and after
election) as both are hot button issues among the voting population.

Based on presidential outcome, we conjecture that Mr. de Blasio, the
democratic party’schoice, will win the November Mayoral election in New York
City (Nate Silver predicted Ms. Quinn, who lost in the democratic primaries
[5]).

##  5 Conclusion

There was evident variation in the sets of data; however, it is so small that
it can be attributed to the small sample size relative to the voting
population. There is no evidence of a change in voter opinion between the two
surveys and the mayoral election.

##  6 Acknowledgements

The Authors wish to thank Dr Janaki Krishnamoorthy and Prof Gurpur Prabhu with
their constructive suggestions and editing the manuscript. First author wishes
to thank Mr. Sean O’Sullivan for establishing Rensselaer Center for Open
Source Software where part of this work is carried out.

## References

  * [1] Nate Silver, _Five-Thirty-Eight Blogs from New York Times_. http://fivethirtyeight.blogs.nytimes.com/2012/11/05/in-ohio-polls-show-benefit-of-auto-rescue-to-obama/ November 5, 2012. 2012\.
  * [2] Dick Morris, _Prediction Romney 325, Obama 213_. http://www.dickmorris.com/prediction-romney-325-obama-213/ November 5, 2012.
  * [3] John Sides, _How Representative Are Amazon Mechanical Turk Workers?_ http://themonkeycage.org/2012/12/19/how-representative-are-amazon-mechanical-turk-workers/
  * [4] Paul Cuff, Sanjeev Kulkarni, Mark Wang and John Strum, _Voting Research- Voting Theory_ http://www.princeton.edu/~cuff/voting/theory.html
  * [5] Jacob Kornbluh, _Nate Silver: Quinn the most likely democrat, Jewish Press, Sept 6, 2013._ http://www.jewishpress.com/news/breaking-news/nate-silver-quinn-the-most-likely-democrat/2013/06/09/