Black hole horizons at the extremal limit in Lorentz-violating gravity

Lorentz-violating gravity theories with a preferred foliation can have instantaneous propagation. Nonetheless, it has been shown that black holes can still exist in such theories and the relevant notion of an event horizon has been dubbed `universal horizon'. In stationary spacetimes the universal horizon has to reside in a region of spacetime where the Killing vector associated with stationarity is spacelike. This raises the question of what happens to the universal horizon in the extremal limit, where no such region exists anymore. We use a decoupling limit approximation to study this problem. Our results suggest that at the extremal limit, the extremal Killing horizon appears to play the role of a degenerate universal horizon, despite being a null and not a spacelike surface, and hence not a leaf of the preferred foliation.

Lorentz-violating gravity theories with a preferred foliation can have instantaneous propagation. Nonetheless, it has been shown that black holes can still exist in such theories and the relevant notion of an event horizon has been dubbed "universal horizon". In stationary spacetimes the universal horizon has to reside in a region of spacetime where the Killing vector associated with stationarity is spacelike. This raises the question of what happens to the universal horizon in the extremal limit, where no such region exists anymore. We use a decoupling limit approximation to study this problem. Our results suggest that at the extremal limit, the extremal Killing horizon appears to play the role of a degenerate universal horizon, despite being a null and not a spacelike surface, and hence not a leaf of the preferred foliation.

I. INTRODUCTION
From a perspective of causality, there seem to be at least two distinct ways one can violate boost invariance -an essential part of Lorentz symmetry. The first way is to introduce a preferred frame and have different types of excitations propagate at different speeds in that preferred frame, perhaps some of them superluminal and some of them subluminal. The second way is to introduce a preferred foliation. This allows for the following: (i) a more general modification of the dispersion relation that can include higher-order terms of the type ω 2 ∝ k 2 + α k 4 /M 2 + . . ., where ω is the frequency, k is the wave number, M is a characteristic mass scale associated with Lorentz symmetry breaking, and α is a dimensionless parameter; (ii) degrees of freedom that carry no time derivatives at all [1][2][3]. Either way, one can have infinitely fast propagation (either for very large momenta or for all momenta respectively).
Superluminal propagation at finite or infinite speed requires one to rethink the concept of a black hole. The existence of black holes in the Universe might appear to be the ultimate vindication of the causal structure of general relativity and of Lorentz symmetry. However, before making such a statement, one actually needs to understand whether and how the concept of black holes fits into Lorentz-violating gravity theories. In fact, this is in its own right a very strong motive for studying Lorentzviolating theories in the first place.
In a Lorentz-violating theory in which excitations propagate with different but finite speeds in a preferred frame, it is rather straightforward to suitably generalise the definition of a black hole. Such excitations can be thought of as propagating along null cones of different effective metrics [4]. A stationary black hole spacetime will possess a succession of horizons, one for each excitation that travels at a unique speed. All of these horizons will be null surfaces of the corresponding metric. They will all cloak the singularity and the innermost one, the one corresponding to the fastest mode, will act as the boundary of the region that is causally disconnected from the exterior [4,5]. This scenario is well studied in the context of Einstein-aether theory (ae-theory), which is the most general theory of a unit timelike vector (aether) coupled to Einstein gravity, containing only two derivatives of the unit vector u µ [6,7]. Because of the unit constraint g µν u µ u ν = 1, the aether field u µ never vanishes but instead always defines a preferred frame. The theory propagates spin-2, spin-1 and spin-0 modes that propagate along null trajectories of the effective metrics g The issue of the existence and definition of black holes is more subtle in the case of a theory with a preferred foliation, since infinitely fast propagation seems particularly hard to reconcile with our conventional understanding of horizons. Additionally, notions like null trajectories, null surfaces, and null infinity become irrelevant for causality, which is now defined by the leaves of the foliation that correspond to constant preferred time surfaces (see Ref. [8] for a more detailed discussion). Remarkably, a suitable notion of black hole does exist in such theories and it is related with the behaviour of the foliation instead of the metric. In a generic asymptotically flat spacetime, the leaves of the foliation extend to spacelike infinity. If a leaf fails to satisfy this requirement it actually bounds a region that cannot communicate with infinity, since any signal can only penetrate a leaf in a single direction while traveling into the future. Hence, such a leaf defines a black hole (see Ref. [8] for a rigorous definition) and it is dubbed universal horizon [2,5].
The concept of a universal horizon was first uncovered in the context of static, spherically symmetric solutions in Einstein-aether theory and the infrared limit of Hořava gravity [2,5]. Hořava gravity [9][10][11] is actually a quantum gravity candidate. Its infrared limit, on which we will focus here, bears similarities with Einstein-aether theory, as will be discussed in more detail below. A key difference though is that unlike Einstein-aether theory, Hořava gravity has a preferred foliation. That is, the equations are second order in time derivatives only in a specific foliation, which is thus causally preferred [1,8,12]. Moreover, the theory exhibits instantaneous propagation even in the infrared limit [2,3]. Hence, Hořava gravity is the ideal test bed for a quantitative study of the concept of a universal horizon.
A foliation by spacelike hypersurfaces corresponds to the level surfaces of a scalar field whose gradient is everywhere timelike. Let us call this field T . The unit one-form allows one to describe the foliation in a coordinateinvariant fashion and without the need to label the surfaces, as u µ is invariant under reparametrizations of the type T → f (T ). The reason that we chose u µ for both the normalised gradient of T here and the aether in Einsteinaether theory will become clear shortly. In stationary spacetimes it has been shown that one can characterise the universal horizon by the local condition u µ χ µ = 0 provided that a µ χ µ = 0, where χ is the Killing vector associated with stationarity and a µ ≡ u ν ∇ ν u µ is the acceleration of u µ [8]. The condition u µ χ µ = 0 bears a strong similarity with the local characterisation of event horizons as Killing horizons in general relativity and the condition χ 2 = 0. The a µ χ µ = 0 requirement came up as a technical restriction in the proof given in Ref. [8], but it appears to have a physical interpretation as a nondegeneracy condition for the universal horizon. In particular, it can be shown that a µ χ µ is constant along the universal horizon [8] and it has been argued to play the role of surface gravity [13,14]. Since u µ is timelike per definition, the condition u µ χ µ = 0 can only be satisfied in a region of spacetime where χ µ is spacelike. Indeed, in all of the known stationary solutions the universal horizon lies in such a region [2,5,[15][16][17][18][19][20][21]. This raises the question of what happens to the universal horizon in the limit where a rotating or a charged black hole becomes extremal. In GR, in the extremal limit the Killing vector is everywhere timelike except on the (degenerate) Killing horizon, where it is null. Hence, it appears that no such solution can have a universal horizon. If instead one considers a nonextremal solution with a universal horizon, then the latter is expected to lie between two Killing horizons [8]. As one moves closer to extremality, the two Killing horizons should get closer to each other and to the universal horizon and there is no reason to expect that any singularity will appear even for an arbitrarily small deviation from extremality.
Motivated by the above, we study here the behaviour of the foliation in extremal and nearly extremal black holes in the infrared limit of Hořava gravity. To be able to restrict ourselves to spherical symmetry that simplifies the analysis, we consider charged instead of rotating black holes. As we will discuss below, the assumption of spherical symmetry has the added advantage that our analysis will apply to Einstein-aether theory as well [12,16]. An additional restriction we impose is that we use the decoupling limit: i.e. we consider only the behaviour of the foliation in a fixed black hole background that solves Einstein's equations and ignore the backreaction of the foliation-defining field on the geometry. Though this might seem drastic, we believe that it provides technical simplification without compromising the qualitative understanding of the problem at hand.

II. EINSTEIN-AETHER THEORY AND INFRARED LIMIT OF HOŘAVA GRAVITY
Below we will use the infrared limit of Hořava gravity. For the purposes of our analysis it is instructive to introduce the corresponding action and field equations through their relation to ae-theory. The action of aetheory reads where R is the Ricci scalar and with c i being dimensionless coupling constants whereas G ae is related to Newton constant G N as measured in a Cavendish experiment by G N = G ae /(1 − c 14 /2) and we are following the convention c ij = c i + c j . The aether satisfies the unit timelike constraint (1) that can be imposed explicitly in (3) as a Lagrange multiplier term ζ(g µν u µ u ν − 1). The equations of motion for the metric and aether fields can be obtained by varying the action with respect to g µν and u µ and using the unit constraint to eliminate the Lagrange multiplier. This process yields where G µν ≡ R µν − 1 2 Rg µν is the Einstein tensor, is the aether stress energy tensor, As discussed in the Introduction, the unit constraint (1) implies that the aether cannot vanish in any spacetime, including the Minkowski spacetime. Instead, the aether defines a preferred threading by timelike trajectories and hence it introduces a preferred frame.
One could further restrict the aether to be hypersurface orthogonal [12]. Taking into account the unit constraint as well, hypersurface orthogonality can be imposed by the local condition (2). The fact that u µ is timelike implies that T will always have a timelike gradient throughout spacetime. That is, T defines a foliation by spacelike hypersurfaces and u µ is the normal to the leaves of this foliation.
For the rest of this paper we will focus on the theory described by action (3) together with the condition (2). This condition is imposed before the variation and T is now the fundamental field. One can show (see Ref. [17] for a detailed discussion) that the equation of motion for T is The field equation that one obtains when varying with respect to the metric is given by Eq. (5) with u µ satisfying Eq. (2). If one decides to use T as a time coordinate and write the theory in the foliation defined by T , then the action coincides with the action of Hořava gravity, truncated to two derivatives in that foliation [12]. Hence, we refer to the theory as the infrared limit of Horava gravity and the action (3) supplemented by the condition (2) is its manifestly covariant formulation. In its full glory, Hořava gravity includes terms that contain higher order spatial derivatives. These terms lead to higher-order dispersion relations and infinitely fast propagation, as discussed in the Introduction. At the same time they improve the UV behaviour of the theory [9]. In principle one needs to include in the action all such terms that have up to six spatial derivatives and are compatible with the symmetries of the theory, i.e. diffeomorphisms that respect the foliation, t →t(t) and There is a large number of such terms that makes the theory less tractable and raises concerns about predictability. As a result, various restricted versions that consistently reduce the number of UV terms have appeared [22][23][24][25][26] but this reduction usually comes at the cost of infrared viability [27][28][29][30][31][32][33][34].
Here we will not be concerned with the UV completion of Hořava gravity. Hence, we will ignore the higher order terms and focus on the infrared limit as defined above. Already in that limit, the theory has the characteristics that we want: a preferred foliation and instantaneous propagation. This last statement might not seem obvious, because the infrared limit does away with the higher order terms in the dispersion relations and renders them linear. Some intuition can be gained by inspecting Eq. (9). This equation is fourth order in derivatives of T , as AE µ contain two derivatives of u µ , u µ contains a further derivative of T and there is a extra explicit derivative. However, if one chooses T as a time coordinate then the only nonvanishing component of u µ is u T = N , where N is the lapse of the foliation defined by T = const surfaces. Moreover, AE ν u ν = 0 identically: i.e., AE ν is orthogonal to u ν and it lies entirely on T = const surfaces.
Hence, the derivative that appears explicitly in Eq. (9) becomes a spatial divergence on a leaf. These suggest that the theory is second order in time derivatives only in the preferred foliation and that there is some elliptic equation that needs to be solved on each slice. Indeed, more detailed analysis has shown that both statements are true [1][2][3].

A. The Reissner-Nordström metric
In general relativity, the Reissner-Nordström (RN) metric describes a spherically symmetric object with electromagnetic charge. The parameters needed to describe this solution are its mass M and electric charge Q. The metric reads: where f (r) = 1 − rs r + r 2 q r 2 , r s = 2G N M and r 2 q = Q 2 G N /(4πε 0 ). The RN metric admits a Killing vector that is timelike at infinity and it is given in our coordinate as χ µ = (1, 0, 0, 0). The condition for the Killing vector being null determines the position of the Killing horizons that are located at The cosmic censorship conjecture [35] states that in order to avoid naked singularities r q ≤ r s /2. Thus, the extremality condition occurs as r q = r s /2 for which the two radii of Eq. (11) coincide at r = r ex = r s /2. We introduce the parameter ε, useful to express the deviation from extremality

B. Determining the foliation at decoupling
As stated in the Introduction, instead of considering the full field equations (5) and (9), we will work in the decoupling limit. That is, we will neglect entirely the backreaction that T has on the metric and we will solve Eq. (9) on the RN black hole background that was reviewed above. It is worth mentioning that Eqs. (9) and (6) share all static, spherically symmetric solutions with flat asymptotics [16]. Hence, our solutions will also determine fully the configuration of the aether in ae-theory in the decoupling limit.
Assuming that u µ respects the Killing symmetries of the metric, the most generic form for the scalar field T is with C(r) a generic function of r. The four-velocity vector will be and then, imposing condition (1) yields Recall that a µ ≡u µ . With this ansatz one has (15) As already discussed above, for static, spherically symmetric solutions with flat asymptotics, Eqs. (9) and (6) yield the same solutions [16] and in our setup they can be reduced to the following equation

IV. FOLIATING EXTREMAL AND NEARLY EXTREMAL BLACK HOLES
A. Extremal case (ε = 0) As we have already discussed in the Introduction, it appears to be impossible to satisfy the condition u·χ = 0 in an extremal RN spacetime. The Killing vector χ µ is everywhere timelike or null, but never becomes spacelike. u ν is always timelike and the inner product of a timelike vector with a null or timelike vector χ µ is never zero.
Recall that u · χ = 0 serves as a local characterisation of a universal horizon r UH , provided that a · χ = 0, where a µ is the acceleration of u µ . The condition a · χ = 0 has been suggested as a non-degeneracy condition, similar to requiring non-vanishing surface gravity for Killing horizons. Extremal horizons might well be degenerate (as in GR) and hence the arguments above regarding u · χ are not enough to conclude that there will not be a universal horizon in the extremal limit. In fact, there is no known local characterisation of an extremal universal horizon, so it is not possible to check without having the full solution and the global causal structure.
There is actually a rather trivial solution to Eqs. (6) and (9) when ε = 0, T (t, r) = t. It is very tempting to dismiss this solution, due to the behaviour of the T (t, r) = t foliation on the Killing horizon. At r = r ex = r s /2 a t = const surface is null, so it cannot be a leaf of a regular foliation by spacelike hypersurfaces. Said otherwise, the normal vector to these surfaces is parallel to the Killing vector χ, and (if normalized) it satisfies aether's equation of motion (6), but at the same time it needs to be null at r = r ex = r s /2. This implies that the foliation is not well defined there and cannot be smoothly continued past that surface.
Consider, however, t = const surfaces in the exterior of the extremal RN spacetime, which corresponds to the analytic solution of Eq. (16), The foliation leaves are everywhere spacelike and they have the right asymptotic behaviour N 0 (0) = 1. Indeed, for r → ∞, the aether u µ aligns with the killing vector χ µ and this corresponds to the condition N (0) = 1. Moreover, N 0 becomes zero at ξ = 1 (corresponding to r = r ex ). This last features implies that all leaves asymptote and pile up on the r = r ex = r s /2, which is exactly the behaviour they are expected to have as they approach a universal horizon. Hence, this seemingly unphysical, trivial solution is actually compatible with the naive expectation that the universal horizon and the Killing horizon coincide in the extremal limit. In the next section, we will demonstrate numerically that the exterior of the universal horizon in nearly extremal RN black holes does indeed approach this analytic solution as ε → 0.

B. Nearly extremal case
We now move on to nonextremal black holes. The plan is to generate the foliation numerically by integrating Eq. (16). First it is important to understand how many independent parameters our solutions will have. Clearly, given that Eq. (16) is a second-order ODE, generically there will be a 2-parameter family of solutions for a black hole of given mass and electric charge. However, these solutions are expected to be singular whenever the denominator of the right-hand side of Eq. (16) vanishes. The locations of these singularities are determined by the roots of the following equation which we collectively denote as ξ c . As discussed earlier, the spin-0 mode propagates along null geodesics of the effective metric, g µν ≡ g µν − (s 2 0 − 1)u µ u ν . This can be explicitly seen by performing perturbative analysis around an arbitrary background for Eq. (5) and either Eq. (6) or Eq. (9). The square of the speed of the spin-0 mode that is present in both ae-theory and Hořava gravity is .
At the decoupling approximation that we are using here, backreaction is neglected and one solves Eq. (6) or Eq. (9) in a fixed background. One can now determine the effective metric just by linearising Eq. (6) or Eq. (9). This yields g (0) µν ≡ g µν −(s 2 −1)u µ u ν , where s is as in Eq. (17). That is, s is the speed of the spin-0 mode at decoupling. Indeed, in the limit where all c i → 0, one has that T ae µν → 0 and s 0 → s, as expected. Killing horizons of g (0) µν are referred to as spin-0 horizons and act as causal boundaries for the spin-0 mode. Interestingly, their locations are also determined by roots of Eq. (19), ξ c , in our setup. In fact the denominator of the right-hand side of Eq. (16) is the norm of the Killing vector that generates the spin-0 horizon. This behaviour is expected, as it has been seen on all previous works (e.g. [2,4,5]).
This suggests that the most general solution will be singular on one or more spin-0 horizons. * There is always a 1-parameter subfamily of solutions that is regular on at least one spin-0 horizon of choice. To see this, recall that a solution to Eq. (16) can be generated starting from any given radius and selecting two pieces of "initial data", the values of N and N there (that correspond to the two parameters of the general family). One can then choose to start from a spin-0 horizon located at ξ c and denote N c ≡ N (ξ c ). Assuming that N c = 0 and imposing that N (ξ c ) ≡ N c satisfies the bond makes the numerator of the right-hand side of Eq. (16) vanish and one can integrate outwards and inwards and generate a solution that is regular on the spin-0 horizon and its vicinity. Given that the bond above relates the two pieces of initial data, there will be a 1-parameter family of such solutions.
For what comes next, we will apply this prescription to the outermost spin-0 horizon. The solutions of the 1parameter family one obtains when integrating outwards toward infinity (ξ → 0) do not generically satisfy the asymptotic condition N (0) = 1. Imposing this condition requires tuning the remaining piece of initial data.
In practice, we impose the asymptotic condition via shooting. We impose that N − 1 vanishes to a desired tolerance at a large but finite radius r max ≡ r s /(2ξ min ). We set the tolerance to be O(ξ min ), as the general expectation is that N will have a 1/r decay asymptotically * Note that here we are interested in solutions with a universal horizon. The latter will always be cloaked by a spin-0 horizon [8]. In static, spherically symmetric solutions this follows easily from the fact that u · χ = 0, which has to hold on the universal horizon, requires χ µ to be spacelike with respect to the spin-0 metric (or any spin-i metric as defined above). Since χ µ is timelike at infinity, there needs to be a spin-0 horizon outside the universal horizon.
(ξ min is chosen in the range [10 −4 , 10 −6 ], depending on the solution). Imposing the regularity condition (21) is not trivial when integrating numerically, because one needs to start the integration from a point where both the numerator and the denominator of the right-hand side of Eq. (16) vanishes. To avoid this problem, we expand N (ξ) in a Taylor series in δ = ξ − ξ c (δ can be positive or negative depending on the direction we are solving the equation) and we integrate analytically from ξ c to ξ c + δ, up to fifth order, with δ = 10 −5 . We then start the numerical integration from this point. The error of the perturbative solution at ξ c + δ is then O(δ 5 ) ∼ 10 −25 , comparable to the machine accuracy, i.e. 24 significant digits for our real numbers. In general, it is conceivable that a boundary value problem gives rise to multiple solutions with the same boundary conditions. According to Eq. (19), N c can only take values in the interval [0, ε/ (s 2 − 1)(1 − ε 2 )] for s ≥ 1. We have explored this interval numerically (∼ 100 points) and we did not find another solution that satisfies the asymptotic condition N (0) = 1. Though this is not a rigorous proof that our solution is unique, considering numerical limitations the outcome is rather reassuring.
Starting with the same values of N and N that gave a solution with the right asymptotic behaviour and regular on the outermost spin-0 horizon, one can then integrate inwards and generate the interior down to a desired radius. This process is quite straightforward up to (and including) the universal horizon, which is the location where N vanishes. In the next section we present and discuss these solutions, i.e. the exterior of the universal horizon. If one attempts to continue this integration well past the universal horizon, a complication arises: eventually one will encounter a second spin-0 horizon. † The solution does not need to be smooth there and no further regularity condition can be imposed. So, they are essentially two options: (a) the solution is "accidentally" smooth; (b) one will encounter a singularity. In Sec. IV B 2 we will present strong evidence that supports the latter case, but we will also argue that this does not affect the solution for the exterior of the universal horizon, even when ε → 0.

Exterior of the universal horizon
There is degeneracy between changing the mass of the RN black hole and rescaling the radial coordinate in our equations. This implies that one can find a solution for a fixed mass and then obtain all other solutions by rescal- † In most cases that can be referred to as the inner spin-0 horizon, but in some cases there might be more than two spin-0 horizons in total.
ing r. We have generated exterior solutions ‡ for different values of s but there is little qualitative difference between them. Hence, we will present here only solutions for s = 1.5 and for different values of ε. The T = const surfaces (surfaces orthogonal to aether) correspond to leaves of the preferred foliation. T is defined from Eq. (12) as T = t + C(r). In outgoing Eddington-Finkelstein coordinates (v, r * ), where v = t + r * and r * is the tortoise coordinate defined by dr * = dr/f (r), one has In Fig. 1 we show sets of foliation leaves for different values of ε in a (r, v) plane. The leaves have a different behaviour depending on the value of ε. As expected, as extremality is approached, the Killing horizon (red vertical dashed line) moves closer and closer to the universal horizon (green vertical solid line). In all non-extremal cases the leaves of the foliation asymptote to the universal horizon in the same fashion qualitatively. However, as one moves closer to extremality the foliation becomes denser, i.e. one would need to cross more leaves to move radially from a fixed r to another, smaller fixed r.
Another illuminating way to compare the foliations is to define a new spatial coordinate such that the locations of the Killing horizon of the metric and universal horizon r UH do not change for different values of ε. Indeed, one can define such a coordinate as where w(r UH ) = 0 and w(r + ) = 1.
In Fig. 2 we show one leaf from three different foliations (ε = 0.1, 0.02, 0.001) in a (w, v) plane. The leaves cross at a point (w 0 , v 0 ) but we remark that, due to the ε-dependent definition of w, (w 0 , v 0 ) does not correspond to a single point in (r, v) coordinates. One can reproduce the entire foliation simply by shifting up and down the leaves we are showing in this plot. The range of the plot has been chosen rather freely, approximately twice the distance between the Killing horizon and the universal horizon. This plot clearly exhibits the different ways leaves approach the universal horizon as one moves closer to extremality.
The bottom panel of Fig. 1, corresponding to the extremal case and the analytic solution (18) , shows remarkable similarity with the panel right above it that corresponds to ε = 0.001. A remarkable difference is that the analytic solution breaks down on the extremal horizon and, hence, there is no foliation leaf that coincides with (or crosses) this location. In what follows we attempt to ‡ From here onwards we will be referring to solutions that describe the exterior of the universal horizon simply as "exterior solutions". quantify the agreement between the numerical solutions and the analytic one as ε → 0. Fig. 3a shows the difference between the analytic, extremal solution N 0 (ξ) and the numerical solutions N ε (ξ). ξ s and ξ ex are the corresponding values to r s and r ex , respectively, i.e., ξ s = 0.5 and ξ ex = 1. For all values of ε, N 0 (ξ) serves as a good approximation at large radii, as expected. There is deviation as one approaches the black hole horizon, but it clearly decreases as ε decreases. In Fig. 3b one can clearly read off the scaling of these deviations. Indeed, this plot suggests that the numerical solutions differ from N 0 (ξ) by an ε 2 correction.
We close this section with a comment about backreaction, which we have ignored above, as we have used the decoupling limit. We can employ the approximation N (ξ) ≈ 1 − ξ − ε 2 2 ξ 2 (1 + ξ) we derived earlier, which should work well for small ε, in order to evaluate the stress-energy tensor for u µ or T field (7). Up to second order in ε The resulting stress-energy tensor is manifestly parametrically small which justifies the earlier assumption of ignoring backreaction.

Interior of the universal horizon
We now want to describe the behavior of the aether inside the universal horizon. As argued before, when integrating from the outer spin-0 horizon inwards, upon reaching another spin-0 horizon, one can encounter two possibilities: (a) the solution is "accidentally" smooth; (b) a singularity arises on the inner spin-0 horizon. It is not trivial to check numerically which of the above is true. If the solution were accidentally smooth, both the numerator and the denominator in the right-hand side of (16) would vanish at the spin-0 horizon. This would introduce an uncontrollable error in the numerical integration which would be hard to distinguish from an actual singularity.
To overcome this difficulty, we proceed as follows. We generate a new solution, starting from the second spin-0 horizon. We select a value for N and Eq. (19) is used to determine the value of ξ there. Regularity uniquely determines N and then we integrate outwards, in precisely the same fashion as we did for the outermost horizon. This generates a 1-parameter family of solutions. This parameter can be thought of as the value of N on the second spin-0 horizon. By changing N , we then attempt to match this solution to the solution we used for the exterior at some radius. We chose to do the match at the radius of the extreme horizon (or ξ ex ), since this location always lies between the two spin-0 horizons. Since Eq. (16) is a second order differential equation, N and N have to match at the matching point to have a smooth solution. However, we have only one parameter to tune to achieve a matching, the value on N on the second spin-0 horizon. If a solution exists that is accidentally smooth on the second spin-0 horizon, then when one matches N at the matching point, N should show no discontinuity as well. This is the hypothesis we have checked and it turns out not to be valid. Fig. 6 shows the difference between the derivatives of two solutions on the matching point. It is important to mention that we have attempted to change the accuracy of the matching in N or the overall accuracy of our solutions and recalculated the mismatch in N . ∆N has shown no sensitivity to such changes, so it seems sensible to conclude that there is a genuine discontinuity.
The above implies that asymptotically flat solutions with a regular outermost spin-0 horizon will exhibit a singularity on the second spin-0 horizon. Note that the decoupling limit that we are using here cannot be trusted anymore near the singularity. Within our approach this appears to be a singularity in the foliation only, because backreaction is neglected. However, one expects that, beyond decoupling, backreaction will render the spin-0 metric singular on the second spin-0 horizon.
Since the second spin-0 horizon is located inside both the outermost spin-0 and the universal horizon, this singularity poses no threat to causality -it is not naked in any sense. However, one can be rightfully worried about the role of this singularity as ε → 0, as the second spin-0 horizon approaches the universal horizon and the outermost spin-0 horizon. In fact, this is precisely the limit we have studied in the previous section. Fig. 6 clearly shows that ∆N actually decreases as ε → 0. For larger value of ε, ∆N is well above numerical accuracy, but it scales as ε 3 . The distance between the two spin-0 horizon scales as ε. Remarkably, for any given accuracy for numerical solutions, one can always find some sufficiently small ε, such that ∆N would appear to be continuous within numerical tolerance. This strongly suggests that the singularity on the second spin-0 horizon does not affect the numerical calculation that has been used to generate the exterior solution as ε → 0.

V. DISCUSSION
We have studied black holes in the extremal limit in Einstein-aether theory and the low-energy limit of Hořava gravity. We have used a decoupling limit approximation which reduces the problem to determining the preferred foliation of a Reissner-Nordström spacetime. For a black hole of given mass and charge and for flat asymptotics there exists a foliation that has a universal horizon and is regular on it and everywhere in its exterior. Our numerical results strongly suggest that this foliation is unique. However, it also appears to be singular at a specific location in the interior of the universal horizon. This location corresponds to an inner spin-0 horizon: a Killing horizon of the effective metric that defines the propagation cones of spin-0 excitations. This singularity, and the universal horizon, are always cloaked by an outermost spin-0 horizon.
The decoupling approximation that we have used here clearly breaks in the vicinity of the singular inner spin-0 horizon. However, our results strongly indicate that the decoupling approximation is valid for the exterior of the universal horizon and that our exterior solution is unaffected by the existence of this singularity, even at the extremal limit.
As expected, as extremality is approached, the universal horizon, the Killing horizons of the metric, and the outer and inner spin-0 horizons, all move toward the location of the extremal Killing horizon. The exterior of the universal horizon smoothly approaches an analytic solution in which preferred time slices are constant t slices, where t is the standard time coordinate in the Reissner-Nordström line element. Clearly, such a solution is not well defined on the extremal Killing horizon and on the extremal universal horizon, yet it offers a well defined description of the exterior.
Our results provide important hints for the appropriate definition of an extremal universal horizon. In every nearly-extremal solution, u · χ vanishes at constant radius surface, whose radius is larger than the radius of the Killing horizon in the extremal solution. This surface can be identified as the universal horizon. It is a spacelike surface and a leaf of the preferred foliation that fails to reach spatial infinity, well in line with the definitions and results of Ref. [8]. Moreover, a·χ, which can be thought of as analogous to surface gravity [8,13,14], does not vanish on that surface for any non-extremal solution.
It is important to examine the behaviour of the foliation in the extremal limit near the extremal Killing horizon: as Q 2 /M 2 → 1, a · χ → 0 everywhere in the exterior and u · χ → 0 as r → r ex . That is, all leaves that span the exterior appear to asymptote and "pile up" on the Killing horizon at the extremal limit. Remarkably, the extremal Killing horizon appears to satisfy the requirements of the global definition of a universal horizon, as given in Ref. [8]. The fact that a · χ → 0 in the extremal limit suggests that such a universal horizon is degenerate in the appropriate sense, as expected. However, most of the rest of the definitions and proofs of Ref. [8] are inapplicable to the extremal solution. A key assumption of Ref. [8] is that every point of the manifold belongs to a leaf of the foliation and our analysis indicates that points that lie on the extremal Killing horizon cannot satisfy this property. Indeed, a Killing horizon cannot be a leaf of the foliation, as it is a null surface. Moreover, when a · χ = 0, the proof of Ref. [8] that a universal horizon must be a leaf (Proposition 2 in Sec. 4) does not apply. § § Note that Theorem 4 in the same section relies heavily of Proposition 2.
It would be interesting to study nearly extremal and extremal black holes with universal horizons without using a decoupling approximation. It is also important to consider spinning as opposed to electrically charged black holes. It is also important to understand the role of extremal universal horizons in the context of universal horizon thermodynamics (see e.g. [13-15, 36, 37].