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

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.


Introduction
In 2009, Petr Hořava proposed a power-counting renormalizable ultraviolet (UV) completion of General Relativity (GR) [1] 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 → t ′ (t), x i → x ′i (t, x), (1.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.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 [2,3].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, 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 [4].A far from exhaustive list includes: understanding the infra-red (IR) dynamics of the theory [5,6]; the structure of constraints [7,8]; its perturbative UV dynamics [9][10][11][12][13][14][15]; the search for black-hole [16][17][18][19], regular [20,21], and cosmological solutions [22]; the interaction with matter fields [23,24], and low-energy signatures in observations [25][26][27][28][29][30][31][32][33][34], among others.However, despite incredible advances, the problem of classically solving and dynamically evolving equations of motion exhibiting the scaling (1.2) is still mostly unexplored, except in the simplest situations -such as for static and spherically symmetric space-times, or in the perturbative limit [35], 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 µ p µ = ω 2 − k 2 =0, with p µ = (ω, k i ) the four-momentum of the field1 ; the presence of higher spatial derivatives generically modifies it to where we have assumed the addition of an operator with 2z spatial derivatives, controlled by a dimensionful coupling Θ.Standard explicit numerical methods tend to develop instabilities when evolving equations of motion exhibiting this behaviour [36].For hyperbolic systems, a Von Neumann analysis of explicit methods -such as the standard Runge-Kutta algorithmsshows that the corresponding evolution is stable as long as the dimensionless Courant number satisfies C ≤ C max , where ∆t and ∆x are the time and space steps used for discretization, and v m is the maximum speed of propagation allowed by the equation.The value of C max generically depends on the method, but for explicit schemes it usually takes the value C max = 1.
In the dispersive case instead, a similar von Neumann analysis yields a more restrictive condition for stability, of the form ∆t (∆x) p ≤ α, (1.5) where p > 1 and α > 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 (1.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 2 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 3, while an application to the Lifshitz scalar field around a Lorentz-violating black hole is discussed in detail in section 4. Finally, we show our results in section 5, drawing conclusions and future research directions in section 6.

The Lifshitz scalar field
Field theories exhibiting the anisotropic scaling (1.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 [37].Its action in a 3 + 1 flat space-time with z = 3 reads where ∂ 2 = ∂ i ∂ i , Λ 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 (2.1) contains all possible parity-preserving spatial derivative operators with z ≤ 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 [1].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 (2.1) reads which leads to a modified dispersion relation of the form (1.3), that is This becomes explicitly Lorentz-violating for momentum scales k ≳ Λ. Stability of the degrees of freedom demands ω 2 > 0, which in particular requires the parabola y(x) = 1 + c 2 x + c 3 x2 to be concave, leading to c 3 > 0. Regarding the value of c 2 , ω 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 ≤ 4c 3 .Putting both conditions together, we get In order to extend this structure to curved space-times, we couple the scalar field to Hořava gravity [1], which implements the anisotropic scaling (1.2) through a foliation in co-dimension one space-like hypersurfaces.This is described by supplementing the space-time curved metric g µν with a hypersurface-orthogonal unit-norm and timelike vector U µ , called the aether [38].Time derivatives are identified with derivatives along U µ , while the spatial Laplacian ∂ 2 is replaced by ∆ γ = γ µν ∇ µ ∇ ν , where γ µν is the orthogonal projector γ µν = −g µν + U µ U ν .The action thus takes the form [39] where L U is the Lie derivative along U µ .The corresponding equation of motion is where acceleration of the aether, and K = γ µν K µν the trace of the extrinsic curvature of the foliation leafs K µν = (L U γ µν )/2.This equation can be then thought of as a generalization of (2.2) 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 (1.2), in particular regarding the integration of the equations of motion with numerical methods.Equation (2.5) 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 ≥ 0 and c 2 ≥ 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 [40].Notice that these conditions include the ones required before for stability of the degrees of freedom.
A first attempt at evolving (2.5) on top of a black hole space-time was performed in [35], 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 [36].Assuming that the coefficients of the equations are locally constant (both in space and time), the eigenmodes of the discretized equations take the form ϕ n j = ξ n e ikj∆x , where ϕ n j is the numerical approximation of the solution at time step t n = n∆t and grid point x j := j∆x, with ∆t and ∆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 ξ 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 ξ n , the discretized equations will be stable if and only if |ξ n | ≤ 1 for all k.For hyperbolic equations, this analysis implies the well known Courant-Friedrichs-Lewy (CFL) condition, which relates ∆t and ∆x with the maximum characteristic speed in the equation.For dispersive equations such as (2.5) instead, the anisotropic scaling between space and time leads to a stronger upperbound for the time step ∆t, which scales in the form given in (1.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.the specifics of Eq. (2.5), and consider instead a generic partial differential equation (PDE) in 1 + 1 dimensions for a scalar field ψ(t, x), of the form where F is a sufficiently smooth function, linear in ψ and in all its derivatives, and ψ 0 , ψ 1 are given initial conditions.The indices j and k span the ranges 1 ≤ j ≤ J, and 0 ≤ k ≤ K, with J ≥ 3, and K ≥ 0. With these conditions, Eq. (3.1) 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 ω 2 = p(k), where p(k) is a polynomial of order greater than or equal to three.Notice that the general problem (3.1) does include the Lifshitz field equation (2.5).
We aim to numerically evolve Eq. (3.1), for which we consider a first-order reduction in time, with dynamical variables ψ and Π := ∂ t ψ; namely, We consider a uniform spatial grid of N points x i = x 0 + (i − 1)∆x, i = 1, . . ., N , with step ∆x = L/(N − 1), being L the length of the spatial domain.We also discretize time as t n = n∆t, n = 0, 1, 2, . .., with time step ∆t, and define the grid functions as ψ n i := ψ(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 ψ at the grid point i, centered at a time step n + 1 2 , is computed as (Dψ) For the second derivative, one gets and similarly for higher-order derivatives.The same average is also performed for evaluating the field ψ For the time derivative, instead, we just keep the forward finite difference 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: 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 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 , ℓ j } and { dj , ūij , lj } are grid functions that also depend on ∆x and ∆t, but not on the dynamical fields.The system (3.10) can be solved with standard numerical methods for band-diagonal linear systems.These usually require O(N ) iterations at each time step, unlike Gaussian elimination, which requires O(N 3 ) iterations (see [36] for details on the corresponding algorithm).Finally, notice that the generic scheme (3.8)-(3.9) is also valid for the case in which F is a non-linear function of ψ, which would clearly lead to a non-linear system of coupled equations for the grid functions ψ n i and Π n i .Nevertheless, such a system can be solved e.g. by means of standard iterative numerical methods.

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 (2.5) on top of a spherically symmetric, static, and asymptotically flat black hole space-time [19], solution to the equations of motion of khronometric gravity in vacuum [5,41].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.

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 aether in Schwarzschild coordinates {τ, r, θ, φ} read with ) where µ is the mass of the black hole, dΩ 2 is the S 2 line element, and contains the only free parameter in the action, c 13 .Note that since the solution is static, there exists a Killing vector χ µ = (1, 0, 0, 0), with norm χ 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. (2.5).In this chart of coordinates, the Laplacian in the orthogonal leafs, ∆ γ , also contains higher time derivatives, so the corresponding equation of motion obtained from Eq. (2.5) 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 µ , 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) so that the metric and aether now read ) Notice that after this transformation, the metric takes the Arnowitt-Deser-Misner (ADM) form [42] where N, N i and h ij are the lapse, shift and induced metric in the foliation leafs, respectively, and are given by This chart of coordinates, however, has a pathology whenever N ≡ 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 [16,19,43].From the relation (4.6), we see that this point lies at finite r, but it is mapped to t → +∞, signalling that the foliation cannot be globally extended smoothly beyond this point, although observers moving inwards along U µ can still cross it in finite proper time (because dτ ≡ U µ dx µ 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 χ µ is space-like in the vicinity of the point r = r U , while the product χ • U = H(r)/(2A(r)) precisely vanishes at this point, and becomes negative in the inner region 2 .This is enough to characterize this surface as a universal trapping surface, hence named universal horizon [16,43].The region r < r U always lies behind the Killing horizon -otherwise χ • U = 0 would not be possible, as U µ is time-like everywhere by definition.In our particular case, one has For a more detailed discussion on causality within space-times endowed with universal horizons, see [44].
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 [45].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.

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 3.
We solve Eq. (2.5) in spherical coordinates {t, r * , θ, ϕ}, where the radial tortoise-like coordinate r * is chosen to push the universal horizon to infinity, corresponding to r * → −∞.The relation between the areal radius r and the new coordinate r * is given by 3 (4.12) 2 Note 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.1), under which N → N (dt ′ /dt). 3In principle, any transformation of the form p with integer exponent p ≥ 1 equally pushes the universal horizon to r * → −∞.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.
Although complicated in general, one can see that it behaves as r * ∼ r for large r due to asymptotic flatness of the background solution, which implies f (r → ∞) → 1 and A(r → ∞) → 1. Close to the universal horizon, we have H(r)/(2A(r)) ∝ (r −r U ) and hence r * ∝ −(r −r U ) −1 , indeed placing r U at r * → −∞.The condition (4.12) is an ordinary differential equation, which can be solved numerically using the standard fourth-order Runge-Kutta method.We also define the parameters κ 2 = c 2 /Λ 2 and κ 3 = c 3 /Λ 4 for computational convenience.Thus, to ensure the stability and well-posedness of the problem, we require κ 3 ≥ 0 and κ 2 > −2 √ κ 3 .After this change, Eq. (2.5) reads Due to the spherical symmetry of the background, we consider the following ansatz for the scalar field: where Y ℓm (θ, φ) are the spherical harmonics.Plugging this into Eq.(4.13), we get an effective 1 + 1 dimensional differential equation for every mode ψ ℓm (t, r * ), which is independent of m, and reads The coefficients ζ ij and the effective potential V eff are functions of r * through the coordinate transformation (4.12), the metric functions A(r), H(r), and their derivatives up to fifth order.They also depend on κ 2 , κ 3 and ℓ.The explicit formulae for ζ ij are given in Appendix (B), while we show V eff later in Eq. (5.1).Let us note that V eff vanishes for ℓ = 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 3, Eq.(4.15) admits the implicit discretization where We implemented centered finite difference operators with second-order accuracy, with a stencil of seven grid points.In this case, the system (3.10) has five upper diagonals and seven lower ones.Table 1 shows the explicit form of the coefficients in the matrices.

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 where a 0 , r * c and σ 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 This corresponds to an exact ingoing wave-packet when κ 2 = κ 3 = 0, satisfying ∂ψ/∂t| t=0 = ∂ψ/∂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 ω is small.Hereinafter, we set ω = 1.

Boundary conditions
Setting boundary conditions for equations endowed with an anisotropic scaling (1.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, 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.(4.13).This can be seen by simply plugging the ansatz (4.20) into (2.2).Only when -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 [46,47] 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 ψ n i at each time step by ψ n i → L i ψ n i .We choose in particular the following function for the ADL: where s controls the slope of the function in the extremes of the numerical domain, while x L and x R correspond to the left and right midpoints of the decaying regions, respectively.An illustration of this function is given in Figure 1.

Results
We report here the results of our simulations, performed with different values of the coupling parameters κ 2 and κ 3 controlling the strength of the Lorentz-breaking terms in Eq. (4.13).We also investigate the behaviour of wave modes with angular number ℓ > 0 for fixed values of the couplings, as well as the effect of the centrifugal barrier induced by the effective potential V 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.

Evolution of the Lifshitz field
We start off by studying the dynamical evolution of the ℓ = 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µ, as shown in figure 2. We set κ 2 = 0.1 and vary κ 3 , taking the values κ 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 4 , 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 κ 3 .A comparison with the standard wave equation evolution (that is, when κ 2 = κ 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 ∼ 80µ, whose magnitude increases faster for larger values of κ 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 κ 2 , and κ 3 , as shown in figure 3.
Furthermore, in figures 4 and 5 we display the evolution of the solution for fixed κ 3 , varying instead κ 2 .In particular, we set κ 3 = 0.01 and consider κ 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 κ 2 , at least at early times.At later times (t ∼ 60µ), and once the solutions have reached the Killing horizon, the order of magnitude of the wave profiles with κ 2 = 0.5 and κ 2 = 1.0 remain similar, until a bump around the horizon forms and starts growing faster as we increase κ 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 5 .In particular, let us note that Hořava gravity propagates a scalar degree of freedom together with the usual transverse traceless graviton perturbation [5].When expanded around the background solution considered here, the dynamics of the scalar mode must exhibit, for consistency, a Lifshitz scaling of the form (1.2), and hence the equation of motion for scalar perturbations must unavoidably take the form (2.5), 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 (2.5), it will also develop in the gravitational case, therefore signaling a linear instability of the universal horizon [48].

Effect of the potential barrier
The coefficient multiplying the term linear in ψ ℓm in (4.14) acts as an effective centrifugal barrier for modes with ℓ ̸ = 0. Its explicit expression reads .1) 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 eff has a single maximum for fixed ℓ, in the vicinity of the Killing horizon.Its value grows quickly with ℓ, 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 κ 2 = 0.5 and κ 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.

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.(3.1) 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 (4.16)-(4.17),with second-order accurate finite-difference operators, is given by where α 1 , α 2 , and α 3 are functions on the grid, and O(∆ 3 ) denote all terms that are cubic in time and spatial steps, like (∆t)(∆r * ) 2 .Thus, by rescaling ∆r * → λ∆r * and ∆t → λ∆t for some λ ∈ R, one gets which indicates that the scheme is second-order accurate.
In order to confirm that Eq. (5.3) approximately holds during the numerical evolution, we perform three different runs with λ = 1 (low resolution), λ = 0.5 (medium resolution) and λ = 0.25 (high resolution).Then, we define the ratio [36] which behaves as Q(t) ∼ 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 κ 3 (fixing κ 2 = 0.1).We can see that the method is approximately second-order accurate, as expected, with convergence improving for smaller values of κ 3 .Finally, and as a last consistency check of the numerical code, an independent residual evaluator test was also done, whose results are reported in Appendix C. 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 κ 2 = 0.1, and κ 3 = 0.01.We test three different positions of the left wall of the layer, namely x L = {−300µ, −400µ, −500µ}, while keeping the right wall fixed far from the horizon, at x R = 300µ.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 ∼ 100µ), this analysis suggests that choosing x L ∼ −300µ 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.

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 (1.2).Our approach is based on a generalization of the Crank-Nicolson 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 highresolution 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, κ 2 and κ 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.
Spherically symmetric and static solutions to (A.1) can be searched for with the ansatz [16,19] where H(r) = 1 + f (r)A(r) 2 , and the form of U µ is chosen for convenience, automatically satisfying U µ U µ = 1.Accidentally, spherical symmetry and staticity automatically impose hypersurface orthogonality, and hence all solutions to Einstein-AEther 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-AEther gravity are not hypersurface orthogonal [50,51].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 [16], there exist two corners where analytic solutions can be attained [19].These correspond to c 14 = 0 and c 123 = 0, where we are using the notation c ij...k = c i + c j + • • • + c k .The latter case, although leading to simpler functions, is incompatible with observational bounds constraining the parameter space of the theory [34].The former case, instead, is perfectly compatible with current bounds and leads to The parameter r ae is in principle arbitrary.However, the solution displays a singularity at the universal horizon, unless 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 (A.1) by higher derivative spatial terms [5].The contribution of these operators -which are critical for achieving power-counting renormalizability -, and in particular their backreaction onto space-times of the form (4.1), 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.

B Coefficients of the evolution equation
In this appendix, we explicitly give the form of the coefficients in Eq. (4.15).
, where the prime denotes derivative with respect to r, and the coefficients are functions of r * through the transformation r(r * ) given in (4.12).

C Further details of the Lifshitz dynamics
In this appendix we show in detail some features of the dynamics of the Lifshitz field.These include the formation of the first oscillations at early times, as well as the case κ 3 = 1, compared with the choices κ 3 = 0.01 and κ 3 = 0.1 reported in Section 5. We also report the results of an independent residual evaluator test, which serves as a consistency check of the numerical method, and complements the convergence studies reported in sections 5 and 6.

C.1 Early-time dynamics
In figure 9 we show the early dynamics of the ℓ = 0 mode of the scalar field.There, we can observe the production of modes with different amplitudes and propagation speeds, as a consequence of the dispersive nature of the evolution equation (2.5).Due to the presence of higher spatial derivatives, generic propagating modes exhibit velocities with depend on their wavenumber.Indeed, we observe the coexistence of modes with different speeds, which combine to produce a "cascade" effect.We also note that the slowest mode present coincides approximately with a solution to the wave equation (i.e. the case κ 2 = κ 3 = 0), which propagates with constant speed.In figure 10 we display the profile of the scalar field at t = 100µ, for κ 2 = 0.1 and κ 3 = {0.01,0.1, 1}.As can be seen, the "bump" between the universal and Killing horizons increases together with κ 3 .One can also observe that the ringing far from the Killing horizon is robust against the choice of this parameter.This rules out the possibility that these features are an artifact of the numerical method presented in Section 3, and supports their presence as genuine physical features of the scalar field evolution.

C.3 Independent residual evaluator
As a final consistency check of our implicit scheme, we performed an independent residual evaluator test, by forcing the evolution equation to admit an exact solution u exact (t, r * ) at the cost of adding an extra source.Numerically evolving the new "sourced equation" with initial data u exact (0, r * ) should reproduce the original solution.
For a given spatial resolution ∆r We verified that the above condition actually holds for our scheme, which is of order p = 2, by choosing the exact solution  The obtained result is shown in Figure 11.As can be easily seen, the numerical approximation converges to the exact solution with the expected convergence rate, ensuring a correct behavior of the numerical scheme, even for small values of κ 3 .

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 L and x = x R , corresponding to the midpoints of the damping zones.

Figure 2 .
Figure 2. Dynamics of the Lifshitz field (ID type I).Snapshots of the evolution of the ℓ = m = 0 mode, for κ 2 = 0.1 and different values of κ 3 : κ 3 = 0.01 (light blue), κ 3 = 0.05 (blue) and κ 3 = 0.1 (dark blue).The initial data is the ID type I given by Eq. (4.18) (a static pulse far from the horizons), with parameters a 0 = 1, r * c = 150µ and σ = 2.As κ 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 κ 2 = κ 3 = 0), represented by the dotted black profile.

Figure 3 .
Figure 3. Dynamics of the Lifshitz field (ID type II).Evolution of the ℓ = m = 0 mode from the ID type II (ingoing pulse) given in Eq. (4.19), with parameters r * c = 150µ, σ = 3 and ω = 1.We fixed κ 2 = 0.1, and set κ 3 = 0.01 (light red curve) and κ 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 .
Figure 4. Lifshitz evolution of a static pulse, for different values of κ 2 and fixing κ 3 = 0.01.The green curve corresponds to κ 2 = 0.1, the blue one to κ 2 = 0.5, while the red curve corresponds to κ 2 = 1.At early times, the propagation speed of the slowest modes decreases when increasing κ 2 .At late times, instead, the bump at the Killing horizon (orange dotted vertical line) grows when increasing κ 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 .
Figure 5. Lifshitz evolution of an (approximately) ingoing pulse, for fixed κ 3 = 0.01, κ 2 = 0.1 (blue curve) and κ 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.

Figure 6 .
Figure 6.Effect of the centrifugal barrier.Left panel: Plot of the effective potential V eff (in log-scale) as a function of r * (in units of the black hole's mass, µ), for angular numbers ℓ = 1 (black), ℓ = 3 (dotted golden), ℓ = 7 (orange) and ℓ = 10 (dotted yellow).As ℓ 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.

Figure 8 .
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µ (black); x L = −400µ (maroon) and x L = −500µ (coral).At early times, the evolution remains almost unchanged, until t ∼ 60µ, 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 κ 2 = 0.1, κ 3 = 0.01, ∆t = 0.1, and ∆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µ.

Figure 9 .
Figure 9.Early-time dynamics.Lifshitz evolution of the ℓ = 0 mode from t = 0 to t = 8.0µ, for the parameters κ 2 = 0.1 and κ 3 = 0.01.From t = 0.8µ, a "cascade-like" effect can be appreciated, as a consequence of the dispersive nature of the evolution equation.

Figure 10 .
Figure 10.Varying κ 3 .Snapshot of the evolution for the Lifshitz field at t = 100.0µ,for different values of the parameter κ 3 .The growth of the bump around the Killing horizon increases with κ 3 , and the ringing oscillations are still present at large distances.