Rotating black holes in Einstein-aether theory

While the aether has four parameters, treating the theory as a low energy effective description of gravity [22], these are constrained in order to have realistic phenomenology. Taking c_A very close to zero, the theory obviously behaves metrically as GR, and so passes phenomenological tests, provided the aether obeys basic constraints (for example, having a stable evolution, and ${s}_{0,1}^{2}\geqslant 1$ to avoid matter energy loss via Cherenkov radiation into aether modes), but does not behave in an interesting manner (unless one can observe the aether vector directly).

However, having imposed the LIGO constraint, so the spin-2 wavespeed is that of light, substantially changes the original observational constraints computed in [23]. There remain two phenomenologically viable parameter regions, one where the twist coupling may be large, so $\sim O(1)$, while the other couplings are very small in magnitude, and the other where both the twist and expansion couplings are large [7]. An important point is that these constraints do not involve strong field dynamics (apart from considering nucleosynthesis) and given the large couplings, this raises the fascinating possibility that observations of the strong field regime could differ from GR. For example, black holes in these parameter regions may strongly deviate from the Kerr black hole of GR, in the sense that geometric invariants of the black holes deviate from their Kerr counterparts by O(1) amounts for fixed mass and angular momentum. This in principle could allow the theory to be distinguished from GR on the basis of its strong field behaviour, such as the properties of black holes. Whether this indeed occurs is the key focus of this paper.

We now discuss these two allowed regions. The absence of decay of cosmic rays to Cherenkov radiation in the spin-0 and spin-1 modes implies that ${s}_{0,1}^{2}\geqslant 1$ to high precision. Further constraints come from agreeing with weak field tests of gravity [23]. The two PPN parameters α_1,2 are constrained to be small, with Solar System observations requiring |α₁| ≲ 10⁻⁴ and |α₂| ≲ 10⁻⁷. Written in terms of the couplings c_ω,θ,a together with the condition that c_σ = 0, ⁸

$$\begin{align}\hfill \begin{aligned}\hfill {\alpha }_{1}& =-4{c}_{a}\hfill \\ \hfill {\alpha }_{2}& \simeq -\frac{{c}_{a}}{2}\left(1-\frac{3{c}_{a}}{{c}_{\theta }}\right)\hfill \end{aligned}\end{align} \tag{ 2.8 }$$

where we have used that |α₁| ≲ 10⁻⁴ implies |c_a | ≪ 1. Nucleosynthesis constrains |c_θ | ≲ 0.3 (assuming standard matter content) ⁹ . Positivity of the kinetic terms implies,

$$\begin{equation}0\leqslant {c}_{a},{c}_{\theta }.\end{equation} \tag{ 2.9 }$$

Dissipative dynamics of binary (and ternary) pulsar systems has recently been found to improve the constraint on α₁ by a factor of 10, i.e. |α₁| ≲ 10⁻⁵ [8]. Another strong constraint $\vert {\hat{\alpha }}_{2}\vert < 2\times 1{0}^{-9}$ comes from spin precession of a millisecond pulsar [25, 26], where ${\hat{\alpha }}_{2}$ is a strong field analog of α₂. However the translation of the constraint on ${\hat{\alpha }}_{2}$ to one on the coupling parameters is complicated as it involves the neutron star sensitivities to higher order in velocity relative to the aether than they are currently known (see also [7] for a discussion of this point).

2.1.1. One large parameter: region I

Firstly we assume that c_ω ∼ O(1) and not small in magnitude—otherwise we will find ourselves in the trivial regime, |c_A | ≪ 1 where the aether decouples from the theory ¹⁰ . Then taking the remaining parameters as,

$$\begin{equation}{c}_{\omega }=O(1),\quad 1{0}^{-7}\ll {c}_{a}\lesssim 1{0}^{-5},\qquad \frac{{c}_{\theta }}{3}={c}_{a}+O(1{0}^{-7})\end{equation} \tag{ 2.10 }$$

satisfies the phenomenological constraints. Hence in this region c_ω is the only large O(1) parameter. All the other parameters are of order O(10⁻⁵) or smaller, so the aether action is dominated by the twist term. We restrict this region to refer to the parameter range where c_a ≫ 10⁻⁷ which implies c_θ ≃ 3c_a . We will call this parameter region I. We see the spin-0 wavespeed is close to that of light, ${s}_{(0)}^{2}\simeq 1$, while the spin-1 wavespeed diverges as ${s}_{(1)}^{2}\sim 1/{c}_{a}=O(1{0}^{5})$.

2.1.2. Two large parameters: region II

Taking both |α₁|, |α₂| ∼ 10⁻⁷, so that α₁ is smaller than necessary to satisfy current bounds, allows for two independent parameters to be large. Again taking c_ω ∼ O(1) and not small in magnitude, and further allowing c_θ to be an independent large O(1) parameter, these bounds are satisfied by,

$$\begin{equation}{c}_{\omega }=O(1),\quad 1{0}^{-7}\ll {c}_{\theta }\lesssim 0.3,\qquad {c}_{a}=O(1{0}^{-7}).\end{equation} \tag{ 2.11 }$$

The aether action is then dominated by the twist and expansion terms. We term this parameter region II. Now both the spin-0 and spin-1 speeds become very large as ${s}_{(0)}^{2},{s}_{(1)}^{2}\sim 1/{c}_{a}=O(1{0}^{7})$.

2.1.3. The transition between regions I and II

In figure 1 we depict these two regions. The transition region joining these two regions I and II in parameter space is given by,

$$\begin{equation}{c}_{\omega }=O(1),\qquad {c}_{a}=O(1{0}^{-7}),\qquad {c}_{\theta }=O(1{0}^{-7}).\end{equation} \tag{ 2.12 }$$

This joins to the region II taking c_θ very small rather than O(1). On the other hand this space joins with region I when c_a = O(10⁻⁷). A subtlety of this region is that if c_θ → 0 for fixed c_a then it is unclear the PPN analysis holds as the spin-0 wavespeed tends to zero. As this transition region is much smaller than the parameter regions I and II implied by the observational constraints we will not consider it further in this work.

Zoom In Zoom Out Reset image size

We begin by considering region II as it is simpler to analyse than region I. Here the aether action is dominated by the quadratic twist and expansion terms. When the shear and acceleration couplings are set to zero in the aether Lagrangian (2.2), leaving only the twist and expansion couplings, any solution to the usual GR Einstein equation admitting a maximal foliation is also an exact Einstein-aether solution, with the aether set equal to the unit normal to the foliation. The reason is that the twist and expansion of this aether vanish, and the action is quadratic in those quantities, so its variation away from such a configuration vanishes when the Lagrange multiplier for the unit norm constraint vanishes. If the acceleration coupling is instead small, rather than strictly zero, one may thus expect there to exist a solution that is a small perturbation of that in which it does vanish. We will now develop this more explicitly.

We begin by taking ≡ c_a = O(10⁻⁷) which we consider to be a small parameter. For simplicity let us consider a vacuum solution, and hence no matter. Let us attempt to write a consistent vacuum solution as an expansion in about a solution of vacuum GR and a twist free aether configuration, given in terms of a potential function f, as,

$$\begin{align}\hfill \begin{aligned}\hfill {g}_{\mu \nu }& ={\bar{g}}_{\mu \nu }+{\epsilon}{g}_{\mu \nu }^{(1)}+O({{\epsilon}}^{2})\hfill \\ \hfill {u}_{\mu }& ={\bar{u}}_{\mu }+{\epsilon}{u}_{\mu }^{(1)}+O({{\epsilon}}^{2}),\qquad {\bar{u}}_{\mu }=k{\partial }_{\mu }f\hfill \end{aligned}\end{align} \tag{ 3.1 }$$

where ${\bar{g}}_{\mu \nu }$ is a vacuum GR solution, so Ricci flat. Now at order O(⁰) the aether equation receives a contribution from the twist and expansion terms in the action. The twist term vanishes as our leading aether configuration is twist free, and since the term is quadratic, its contribution to the equations of motion and the aether stress tensor vanishes too. However from the aether equation (2.6) we see the expansion term leads to the condition;

$$\begin{equation}{c}_{\theta }{\bar{h}}^{\mu \nu }{\partial }_{\nu }\bar{\theta }=0\end{equation} \tag{ 3.2 }$$

where indices are raised/lowered with the leading metric ${\bar{g}}_{\mu \nu }$ and where $\bar{\theta }$ is the expansion of the leading order twist-free aether,

$$\begin{equation}\bar{\theta }=k{\bar{h}}^{\mu \nu }{\bar{\nabla }}_{\mu }{\partial }_{\nu }f\end{equation} \tag{ 3.3 }$$

where ${\bar{h}}_{\mu \nu }={\bar{g}}_{\mu \nu }+{\bar{u}}_{\mu }{\bar{u}}_{\nu }$ is the projection of the metric onto the constant f hypersurfaces. Thus the expansion $\bar{\theta }$ must be constant on the hypersurfaces of constant f.

Having taken a vacuum GR solution, then the Einstein tensor ${\bar{G}}_{\mu \nu }$ of the metric ${\bar{g}}_{\mu \nu }$ vanishes, and we must also check consistency with the Einstein equation. From (2.7) at leading O(⁰) order we find the expansion term backreacts as,

$$\begin{equation}0={\bar{G}}_{\mu \nu }=\frac{{c}_{\theta }}{3}\left(-2{\bar{u}}_{\left(\right.\mu }{\bar{\nabla }}_{\nu \left.\right)}\bar{\theta }+{\bar{g}}_{\mu \nu }\left(\frac{1}{2}{\bar{\theta }}^{2}+\bar{u}\cdot \bar{\nabla }\bar{\theta }\right)-{\bar{u}}_{\mu }{\bar{u}}_{\nu }(\bar{u}\cdot \bar{\nabla }\bar{\theta })\right)\end{equation} \tag{ 3.4 }$$

and again indices are raised/lowered with the leading metric ${\bar{g}}_{\mu \nu }$. Since the twist vanishes and the twist term in the action is quadratic it does not backreact. Now contracting with u^μ u^ν yields the condition,

$$\begin{equation}{\bar{\theta }}^{2}=0\end{equation} \tag{ 3.5 }$$

which in fact forces the expansion to vanish for consistency.

Thus we see that we may consistently take the spacetime to be that of vacuum GR at leading order O(⁰) with a 'painted on' aether vector that is twist and expansion free, provided we can find a potential function f which obeys the vanishing expansion condition on this spacetime, given by the non-linear p.d.e.,

$$\begin{equation}\left({\bar{g}}^{\mu \nu }-\frac{1}{{(\partial f)}^{2}}{\partial }^{\mu }f{\partial }^{\nu }f\right){\bar{\nabla }}_{\mu }{\partial }_{\nu }f=0.\end{equation} \tag{ 3.6 }$$

We note the two derivative terms in f are contracted by the induced inverse metric on a constant f surface. This induced inverse metric must be Riemannian, since it is orthogonal to the aether which is constrained to be timelike, and hence this p.d.e. has elliptic character. We may view the aether as locally defining a foliation of the spacetime by the constant f hypersurfaces, with the unit timelike norm and expansion free condition implying that it is a maximal slicing. This leading order Ricci flat spacetime is then perturbed at order O(¹) by ${g}_{\mu \nu }^{(1)}$, and the aether by ${u}_{\mu }^{(1)}$. However for region II we have = O(10⁻⁷) so to one part in 10⁷ the spacetime geometry is simply that of vacuum GR. Thus to a great accuracy the vacuum behaviour of GR can potentially be recovered even though both c_ω and c_θ are O(1).

While for convenience we have excluded matter, including it does not change the picture. With matter, ${\bar{g}}_{\mu \nu }$ and the leading matter behaviour will be a solution to the usual GR Einstein-matter equations. The aether will again be twist and expansion free at leading order, and to high accuracy the usual Einstein-matter dynamics will be recovered for such solutions. Furthermore these conditions on the aether are independent of the values of the large aether couplings c_ω and c_θ . Thus the aether behaviour for these near GR solutions is approximately (i.e. to one part in 10⁷) universal within region II.

An important point to emphasize is that there is no guarantee that solutions to the above twist and expansion free equation (3.6) exist for all GR spacetimes. For example, in a cosmological setting the asymptotic expansion of the aether will not vanish and hence one could not hope to find a global solution to the above conditions. While it was long ago proved that vacuum spacetimes 'close' to Minkowski spacetime admit a maximal foliation, to our knowledge there is no stronger result that could guarantee such a foliation more generally, such as for black hole spacetimes that we will later be interested in ¹¹ . For the case of Kerr, maximal spacelike foliations have been constructed numerically [28]. ¹² We will shortly give numerical evidence that a twist and expansion free aether, associated to such a foliation, indeed describes the behaviour of black holes that we are able to construct associated to region II.

If solutions do exist another question is whether these are unique, and whether there might be a moduli space of such solutions. While the above elliptic p.d.e. for the potential f above locally determines f, such a moduli space may arise if there is global data for solutions once boundary conditions are prescribed. A key point is that if a solution or moduli space of solutions for f exist, they are universal to region II in the sense that they are independent of the O(1) couplings c_ω and c_θ , as we see explicitly from the coupling independence of equation (3.6)

Now let us consider region I which has one large parameter, the twist coupling c_ω . The logic is similar to that above for region II. If the expansion and acceleration couplings vanish exactly, then any solution to the usual GR Einstein equation that admits a twist free aether will be an exact solution to the Einstein-aether equations. Now if the expansion and acceleration couplings are not zero, but are small, we then expect a solution to exist which is a small perturbation of this solution where it does vanish.

We again take the small parameter to perform our expansion to be ≡ c_a = O(10⁻⁵). The aether action is now dominated only by the twist term. The expansion and acceleration terms have small coefficients obeying the phenomenological constraint above, $\frac{{c}_{\theta }}{3}={c}_{a}+O(1{0}^{-7})$. Since we have restricted region I to have 10⁻⁷ ≪ ≲ 10⁻⁵ then we can write,

$$\begin{equation}\frac{{c}_{\theta }}{3}={\epsilon}+{k}_{\theta }{{\epsilon}}^{2}\end{equation} \tag{ 3.7 }$$

so that |k_θ ²| ≪ . Again for simplicity let us consider the spacetime at leading order to be a GR vacuum solution, so ${\bar{g}}_{\mu \nu }$ is Ricci flat and there is no matter. Then 'painting on' a twist free aether, given by the potential function f, we again expand the metric and aether as above in equation (3.1).

Now since at leading order O(⁰) the aether action comprises only the quadratic twist term, and the leading aether is twist free, both the aether and Einstein equations are trivially satisfied at this order. One might then wonder whether the aether potential f can be arbitrary? In fact it cannot, as it is constrained by consistency of the O(¹) equations. At this order the aether equation gives,

$$\begin{equation}{c}_{\omega }\left({\bar{\nabla }}_{\nu }{\omega }^{(1)\nu \mu }+{\bar{a}}_{\nu }{\omega }^{(1)\nu \mu }\right)={\bar{j}}^{\mu }\end{equation} \tag{ 3.8 }$$

where ${\omega }_{\mu \nu }={\epsilon}{\omega }_{\mu \nu }^{(1)}+O({{\epsilon}}^{2})$ and we note that ${\bar{u}}_{\mu }{\omega }^{(1)\mu \nu }=0$, and consequently ${\bar{u}}_{\mu }{\bar{j}}^{\mu }=0$ and,

$$\begin{equation}{\bar{j}}^{\mu }=-\left({\bar{\nabla }}^{\mu }\bar{\theta }+{\bar{u}}^{\mu }\bar{u}\cdot \bar{\nabla }\bar{\theta }\right)+\left(\bar{u}\cdot \bar{\nabla }{\bar{a}}^{\mu }-{\bar{u}}^{\mu }{\bar{a}}^{2}+\frac{2}{3}\bar{\theta }{\bar{a}}^{\mu }-{\bar{\sigma }}^{\mu \nu }{\bar{a}}_{\nu }\right)\end{equation} \tag{ 3.9 }$$

which we emphasize is constructed only from the leading metric and aether, ${\bar{g}}_{\mu \nu }$ and ${\bar{u}}_{\mu }$. We may regard equation (3.8) as determining the aether perturbation ${\omega }_{\mu \nu }^{(1)}$. However, somewhat analogously to source conservation in the Maxwell equation, one also finds the constraint,

$$\begin{equation}\bar{\nabla }\cdot \bar{j}+\bar{a}\cdot \bar{j}=0\end{equation} \tag{ 3.10 }$$

a key point being that this involves only the leading order aether and metric. This appears to be a complicated differential equation for the potential f which is fourth order in derivatives and which we may think of as its equation of motion. Another important point is that solutions for f of this p.d.e. are universal for region I in the sense that they do not depend on the O(1) aether coupling c_ω , as we explicitly see from the above equation.

The origin of this condition may be understood by considering the aether action under an O() variation of the aether, ${u}_{\mu }^{(1)}\to {u}_{\mu }^{(1)}+\delta {u}_{\mu }^{(1)}$; the variation of the action at order O(²) is,

$$\begin{equation}\delta {S}^{(2)}=\frac{1}{8\pi G}\int {\mathrm{d}}^{4}x\sqrt{-\bar{g}}\left(-{c}_{\omega }{\omega }^{(1)\mu \nu }\delta {\omega }_{\mu \nu }^{(1)}-\delta {u}^{(1)\mu }{\bar{j}}_{\mu }\right)\end{equation} \tag{ 3.11 }$$

which gives rise to the equation (3.8) for ω⁽¹⁾ above. Now consider the variation generated by changing the foliation at order O();

$$\begin{equation}f\to f+{\epsilon}\delta {f}^{(1)}\,\Longrightarrow\,\delta {u}_{\mu }^{(1)}=k{\bar{h}}_{\mu }^{\nu }{\bar{\nabla }}_{\nu }\delta {f}^{(1)}.\end{equation} \tag{ 3.12 }$$

In analogy to a gauge transformation for Maxwell leaving the field strength invariant, this variation leaves the twist ω⁽¹⁾ invariant, so $\delta {\omega }_{\mu \nu }^{(1)}=0$. Then taking the remaining variation and integrating by parts,

$$\begin{equation}\delta {S}^{(2)}=\frac{1}{8\pi G}\int {\mathrm{d}}^{4}x\sqrt{-\bar{g}}\delta {f}^{(1)}{\bar{\nabla }}^{\mu }\left(k{\bar{j}}_{\mu }\right)\end{equation} \tag{ 3.13 }$$

and since the action should be stationary for a solution, ${\bar{\nabla }}^{\mu }\left(k{\bar{j}}_{\mu }\right)=0$. Then using ${\bar{a}}_{\mu }={\bar{h}}_{\mu }^{\enspace \nu }{\partial }_{\nu }(\mathrm{ln}\,k)$ implies equation (3.10) above.

As for the region II discussion, we have no argument that guarantees solutions to (3.10) for f must generally exist given an underlying GR spacetime. If they do, given that it is a fourth order differential equation in f, it is also unclear how much global data would characterize them. However, we will shortly give numerical evidence that for stationary black holes a solution does exist and furthermore it appears to be unique given the boundary condition that the aether is asymptotically at rest with respect to the black hole.

A possible solution to (3.10) is $\bar{j}=0$, which is a simpler condition involving only three derivatives of f. For example, if one considers static spherically symmetric black holes where the aether is twist free by symmetry, then indeed $\bar{j}$ must vanish (as ω⁽¹⁾ vanishes in equation (3.8)). Could it be the case that while (3.10) holds true, it only does so because of the more fundamental condition $\bar{j}=0$ which really plays the role of the equation governing the potential f? We later show that this is not the case. For stationary rotating black holes our numerical solutions will show that in region I these take the approximate form of Kerr with the aether 'painted on', and further that equation (3.10) is satisfied by an aether where $\bar{j}\ne 0$. Hence it appears that (3.10) is indeed the true condition governing the leading aether potential. Moreover, we find that $\bar{\theta }$, ${\bar{a}}^{\mu }$, and ${\bar{\sigma }}^{\mu \nu }$ are all nonvanishing, so there is no obvious simplification of the form of ${\bar{j}}^{\mu }$ in (3.9).

Subject to the existence of solutions for f, again we have solutions in region I that are very close to GR solutions, now to one part in $\sim 1{0}^{5}$, and where the aether is close to being twist free and determined by a single function f obeying a somewhat complicated four derivative condition. Interestingly this condition for f, and hence the leading aether solution, is independent of where the theory is in region I, given by the large parameter c_ω . Thus such GR-like solutions will have an approximately universal aether behaviour, to one part in $\sim 1{0}^{5}$, in this phenomenological parameter region ¹³ .

The conditions (2.8) define regions I and II precisely so that the deviation of the metric from a GR response to matter in the weak field regime is small at leading order in a PPN expansion. However our discussion above suggests that not only should the metric behaviour in region I and II be close to GR, but the aether should also be approximately twist free and in the case of region II additionally expansion free. It is instructive to see this emerge in the weak field PPN calculation. In appendix B we review the relevant results from the original Einstein-aether PPN calculation [23] and use them to compute the twist and expansion of the aether.

Interestingly we find that the condition c_σ = 0 is sufficient at leading non-trivial PPN order to give a twist free aether. Thus the aether being twist free in weak field at leading PPN order is generic for c_σ = 0, even if the metric response does not behave as GR (as seen through the preferred frame parameters α_1,2). Indeed this is precisely the reason that c_ω does not enter the expressions for α_1,2 in equation (2.8).

Further restricting to regions I and II the aether takes the expected form ${u}_{\mu }={\bar{u}}_{\mu }+{c}_{a}{u}_{\mu }^{(1)}+O({c}_{a}^{2})$ with universal ${\bar{u}}_{\mu }=k{\partial }_{\mu }f$. As shown in appendix B, in region II one finds the potential function simply goes as time, f = t, at leading non-trivial PPN order, which indeed results in ${\bar{u}}_{\mu }$ being expansion free. For region I the potential involves a contribution from the matter, and consequently the expansion does not vanish, again as expected. Since c_σ = 0 is sufficient to give a twist free aether at leading PPN order, in this weak field limit the leading correction to ${\bar{u}}_{\mu }$, given by ${u}_{\mu }^{(1)}$, and all subsequent terms in the c_a expansion of u_μ must also be twist free. More generally for strong fields the aether will only be twist free in the c_a → 0 limit, so these corrections will have non-vanishing twist as we explicitly see later for stationary black holes.

An example of this GR limit for black holes in region II can be given analytically in the static case via the exact solutions in [16] which were found for c_a = 0. Taking c_σ = 0 as we do here, then for any c_ω and c_θ , an exact solution for the metric and aether is simply given by Schwarzschild,

$$\begin{equation}\mathrm{d}{s}^{2}=-\left(1-\frac{{r}_{0}}{r}\right)\mathrm{d}{v}^{2}+2\,\mathrm{d}v\,\mathrm{d}r+{r}^{2}\,\mathrm{d}{{\Omega}}^{2}\end{equation} \tag{ 3.14 }$$

with the aether taking the form,

$$\begin{equation}u=\alpha (r)\frac{\partial }{\partial v}+\beta (r)\frac{\partial }{\partial r},\qquad \beta (r)=-\frac{{r}_{\text{\AE}}^{2}}{{r}^{2}}\end{equation} \tag{ 3.15 }$$

with α(r) determined from β(r) by the norm condition, and r_æ being a parameter. In spherical symmetry any aether is twist free. Writing the aether as u_μ = k∂_μ f we see the potential takes the form, f = v + ϕ(r). One can verify that this aether has zero expansion. Hence this represents the leading aether behaviour as c_a → 0 for the parameter region II. From the discussion above, a static black hole in region II would approximate Schwarzschild with this aether to an accuracy of order O(10⁻⁷).

The parameter r_æ above reflects the fact that the expansion free condition is a differential equation with global data. However, as shown in [16], if one requires the existence of a universal horizon, then ${r}_{\text{\AE}}=\frac{{3}^{3/4}}{4}{r}_{0}$. For region II in the limit c_a → 0 both the spin-0 and spin-1 wavespeeds diverge, and hence existence of a smooth universal horizon is the same statement as existence of smooth horizons for these degrees of freedom. Thus in the static spherically symmetric case for region II it seems that, for a given mass, determined here by r₀, there is a unique aether independent of the aether coupling parameters in the limit c_a → 0.

The universal horizon is located at $r=\frac{3}{4}{r}_{0}$. While the aether vector field is smooth there, the function ϕ, and hence f, diverges as $\sim -\mathrm{log}\vert r-\frac{3}{4}{r}_{0}\vert$, being explicitly given as,

$$\begin{equation}\phi (r)=-r-{r}_{0}\,\mathrm{log}\left[\left(\sqrt{3}(20r+7{r}_{0})+9{\Delta}\right){\left(\frac{\left\vert r-\frac{3}{4}{r}_{0}\right\vert }{4r+\frac{3}{2}{r}_{0}+\frac{3\sqrt{2}}{4}{\Delta}}\right)}^{\sqrt{\frac{27}{32}}}\right]\end{equation} \tag{ 3.16 }$$

where we have defined ${\Delta}=\sqrt{16{r}^{2}+8r{r}_{0}+3{r}_{0}^{2}}$. Thus the twist potential f is not globally smooth, but is smooth separately in the interior and exterior of the universal horizon where it defines maximal foliations. The exterior foliation is one discussed in [31, 32]. In the case of Kerr again one expects that maximal foliations with the correct asymptotics only penetrate a certain distance inside the Kerr horizon [28].

For region I we expect an analogous solution, but with a different function β corresponding to an aether potential function that satisfies the condition $\bar{j}=0$. (As discussed above, since any static spherically symmetric aether is twist free, the only solution to (3.10) is that with vanishing $\bar{j}$.) We note this solution for the aether is implicit in numerical solutions in [17] where static spherical black holes in region I were found, and indeed seen to be approximately Schwarzschild.

The setting under discussion is a spacetime with various tensor fields, including a metric, which are all invariant under the flow of a Killing vector χ that possesses a Killing horizon, i.e. a null surface generated by χ. We divide the discussion into two parts. In the first part, we suppose that the Killing horizon has a bifurcation surface where χ vanishes and all the fields are regular, and explain, following an argument in [33], why this implies that the Killing horizon is a Killing horizon with respect to the effective metrics for all wave modes. We also highlight why the Einstein-aether theory cannot have such solutions. In the second part we review the results of [34] that establish conditions under which the existence of such a bifurcation surface is guaranteed.

We label the different linearized modes in the theory by an index i, and we suppose that the mode equations are all hyperbolic, with local characteristic surfaces (a.k.a. 'causal cones') determined by metric tensors ${g}_{\mu \nu }^{i}$ that are constructed from the fields of the theory. Since by assumption all the fields of the configuration are invariant under the flow of χ, we have ${\mathcal{L}}_{\chi }{g}_{\mu \nu }^{i}=0$. That is, χ is a Killing vector for all of the metrics. Moreover, this implies that the scalar quantities ${g}_{\mu \nu }^{i}{\chi }^{\mu }{\chi }^{\nu }$ are constant on the flow lines of χ, and in particular on the null curves that generate the Killing horizon. Since these generators all pass through the bifurcation surface where χ = 0, it follows that these scalars all vanish everywhere on the Killing horizon, which means that it is a Killing horizon for all of the metrics ${g}_{\mu \nu }^{i}$.

This is a very powerful argument, but it relies on a strong assumption: the existence of a bifurcation surface at which all of the fields are regular. The Einstein-aether theory cannot have such solutions, since a regular unit timelike vector field cannot be invariant under the Killing flow at the bifurcation surface, because the flow acts there as a boost in the tangent space at each point and a nonzero timelike vector is not invariant under any boost. Instead the Einstein-aether theory has black hole solutions which have different future horizons for each wave mode, as one can see for the earlier example in equations (3.14) and (3.15), but these cannot be smoothly extended back to a bifurcation surface where the aether is regular. Furthermore for rotating black holes we will later see these horizons are not even Killing horizons.

Let us briefly return to theories that admit such black holes with a common Killing horizon for all wave modes. The above argument for a common Killing horizon requires a bifurcation surface and fields to be smooth there. In a maximally extended, analytic solution that requirement might be manifestly met, but this is not an adequate criterion in the more phenomenological setting in which a black hole forms via collapse of matter, and asymptotically approaches a stationary black hole solution in the future. What, if anything, can be said in that case? Some strong results due to Rácz and Wald [34] are useful here. What they showed (among other things) is, roughly speaking, that if (1) a spacetime possesses a neighbourhood of a future portion of a Killing horizon; (2) any fields on the spacetime other than the metric are also invariant under the Killing flow; (3) the metric and other fields are either static (and therefore time-reflection symmetric) or stationary-axisymmetric and invariant under a t–ϕ reflection isometry, then the spacetime metric and other fields are all extendable to a neighbourhood of a regular bifurcation surface. Field equations play no role in the argument. What this means is that while a physical spacetime in which a black hole forms by collapse has no real bifurcation surface, it has what might be called a 'virtual bifurcation surface'. The existence of such a virtual bifurcation surface suffices for running the argument given above, so we may conclude that, under the hypotheses of the Rácz–Wald theorem, a Killing horizon is a causal horizon for all of the fields in a gravity theory.

The first two assumptions are routine, but that of the reflection symmetry is more subtle. For a nonrotating black hole, it seems plausible that the metric would generically admit a static Killing vector, whose time reflection isometry would be shared by any scalar fields. But a unit timelike vector, as occurs in Einstein-aether theory, would reverse its time orientation under the reflection, and hence would not satisfy the requirement. Furthermore we will later show that rotating Einstein-aether black holes do not possess the t–ϕ orthogonality property, which presumably means they do not admit a t–ϕ reflection isometry.

Above we have assumed the local characteristic surfaces are governed by metric tensors. However as recently discussed in [35] one may have more complicated situations where the principle symbol may not be decomposed in such a simple manner. Interestingly in the explicit scalar-tensor theories considered there, which all had second order equations of motion, it was shown by direct analysis of the principle symbol, that a metric Killing horizon is a characteristic surface for all degrees of freedom in the theory, without direct reference to the bifurcation surface of the metric horizon.

In order to numerically construct the stationary black hole spacetimes we will employ the harmonic method outlined in [10, 36, 37]. The harmonic Einstein equation takes the form (in four dimensions),

$$\begin{equation}{R}_{\mu \nu }^{H}[g]\equiv {R}_{\mu \nu }+{\nabla }_{\left(\right.\mu }{\xi }_{\nu \left.\right)}={T}_{\mu \nu }-\frac{1}{2}T{g}_{\mu \nu }\end{equation} \tag{ 4.1 }$$

where T_μν is the stress tensor, and the vector ${\xi }^{\mu }={g}^{\alpha \beta }\left({{\Gamma}}_{\alpha \beta }^{\mu }-{\bar{{\Gamma}}}_{\alpha \beta }^{\mu }\right)$, with $\bar{{\Gamma}}$ a smooth reference connection on the manifold which we take to be the Levi-Civita connection of a smooth reference metric $\bar{g}$. In the context of Riemannian geometry this modification of the Ricci tensor was due to DeTurck who used it in the analysis of Ricci flow [38]. The harmonic Einstein equation is then solved in conjunction with the matter equations of motion. This formulation removes the coordinate invariance of the equations, and in the vacuum gravity case (i.e. T_μν = 0) gives a principle symbol controlled by the metric itself. The vector ξ should vanish if we wish to find a solution to the Einstein equation, rather than a solution with non-vanishing ξ, which we term a 'soliton' as in vacuum equation (4.1) is known as the Ricci soliton equation. We may view the vanishing of ξ as a gauge condition, with the resulting coordinates determined by the prescribed reference metric $\bar{g}$. This takes the form of a generalized harmonic gauge condition [39, 40] with the coordinate functions locally obeying ${\nabla }_{S}^{2}{x}^{\mu }=-{g}^{\alpha \beta }{\bar{{\Gamma}}}_{\alpha \beta }^{\mu }$, where ${\nabla }_{S}^{2}$ is the scalar Laplacian. Clearly in order to find solutions of the Einstein equation we must ensure boundary/asymptotic conditions such that ξ may vanish there.

The existence of solutions with non-vanishing ξ is constrained. We may understand this as the vector satisfies the linear equation

$$\begin{equation}O{\xi }_{\mu }={\nabla }^{2}{\xi }_{\mu }+{R}_{\mu }^{\enspace \nu }{\xi }_{\nu }=0.\end{equation} \tag{ 4.2 }$$

Clearly ξ = 0 is a solution. However for ξ ≠ 0 to be a solution, the linear operator O must have a non-trivial kernel. If boundary/asymptotic conditions ensure ξ → 0, then the kernel may be forced to be trivial. This has been proven to occur in the vacuum GR setting for stationary black holes in [41, 42]. In practice one may simply check that solutions obtained by solving the harmonic Einstein equation are not solitons.

We now briefly review the use of this harmonic formulation in the conventional stationary black hole setting where there is a single Killing horizon. In the rotating case we assume rigidity holds so that the black hole rotates in the direction of an asymptotically spatial Killing vector (e.g. ∂_ϕ for Kerr). Then the problem of finding the exterior to the horizon may be phrased as an elliptic boundary value problem on a spatial slice that extends from infinity and terminates at the bifurcation surface [36]. The horizon and asymptotic regions form the boundaries, and the lack of dependence of the metric components on time or the rotation direction implies that the principle symbol is elliptic. The boundary conditions at the horizon ensure its regularity, and also fix the physical data of the black hole; its surface gravity and the velocities of the horizon.

In this stationary setting, with a reasonable initial guess, one can hope to solve the resulting elliptic p.d.e.'s to find the desired ξ = 0 solution. This is typically achieved either as a flow (such as DeTurck flow in the vacuum case) or more directly as we will do here, using a Newton method. While one usually thinks of using the Newton method after finite differencing in order to solve the resulting coupled non-linear algebraic equations, we may formally think of it prior to this in the continuum. In the case of the vacuum equations it is simply,

$$\begin{equation}{g}_{i}={g}_{i-1}-{\Delta}{[{g}_{i-1}]}^{-1}{R}^{H}[{g}_{i-1}]\end{equation} \tag{ 4.3 }$$

where g_i are metrics with g₀ an initial guess and Δ is the linearization, R^H [g + h] = R^H [g] + Δ[g]h + O(²), and the aim is that given a good initial guess g₀ then g_i tends to a solution as i → ∞. Considering the Newton method in the continuum we see that provided g_i−1 is a smooth metric, and the problem is well posed so that the inverse Δ⁻¹ exists, then we expect the new metric g_i to remain smooth.

The approach described above is based implicitly on the black hole of interest having a smooth Killing horizon, and hence a bifurcation surface that forms a boundary of the domain. This will not work if the theory does not admit black holes with a common horizon for all wave modes, as for the Einstein-aether theory. If there are other fields with effective metrics which have horizons inside the metric one, as is the situation for Einstein-aether due to the Cherenkov constraints ${s}_{0,1}^{2}\geqslant 1$, such an approach cannot correctly impose regularity of these horizons and the problem will be ill-posed.

Instead here we will introduce a new approach using in-going coordinates, taking the coordinate domain to pierce inside the horizon of the problem, or horizons if there are multiple effective metrics. This in-going approach has previously been used in [37, 43] where certain exotic black holes with 'non-Killing horizons' were found that naturally occur in the holographic context with AdS asymptotics. Here we will use a similar approach, although the reason is somewhat different. In that holographic setting there was a unique horizon, but it failed to be Killing. Here the issue we must address is the presence of multiple effective horizons for different degrees of freedom inside the metric horizon. Our aim is to extend the domain far enough inside the metric horizon to capture all these effective horizons. An early version of this method applied to Einstein-aether black holes can be found in [44].

Let us firstly consider the case where the metric controls all wave speeds using such in-going horizon piercing coordinates. The problem no longer has an elliptic character and becomes hyperbolic inside the Killing horizon. Smoothness of the horizon is imposed simply by the metric functions being smooth in the interior of the domain where the horizon lies. As discussed above, considering the Newton method in the continuum, then if the metric at some step is smooth we expect the updated metric will also be smooth, and hence provided the method converges and one starts with a smooth initial guess, the resulting solution should have a smooth horizon. Inside the horizon the problem is hyperbolic with the solution determined by the initial data set by the metric functions at the horizon. Thus the innermost boundary of the domain should be regarded as the last slice of this hyperbolic evolution, and no boundary condition should be imposed there. Now consider the case where there are multiple wave modes with different effective metrics controlling their propagation. We may apply the method in the same way, solving the Einstein and matter equations on an in-going slice that pierces all the horizons associated to these wave modes. Provided the metric functions are smooth, regularity of all the wave mode horizons will be ensured. Starting with a smooth initial guess close enough to a solution we expect the Newton method to maintain smoothness as it iterates to the solution. Outside all the horizons we expect the problem to be elliptic in character, and inside all the horizons we expect it to be hyperbolic. It will have a mixed character in the region between the various wave mode horizons. Again since the problem should be hyperbolic inside all the horizons, no boundary condition should be imposed on the innermost points of the coordinate domain, and instead the equations of motion should be solved there. This is schematically illustrated in figure 2.

Zoom In Zoom Out Reset image size

In the previously described elliptic setting we expect a solution given the boundary data which fixes the particular black hole of interest by specifying the moduli of the solution, the surface gravity and angular velocities. In our ingoing setting we are imposing only smoothness at the horizon. Thus there is no unique solution to the problem, but rather any stationary black hole will be a solution. In order to have a tractable numerical problem we must have a unique solution (or at least a discrete set of solutions) and thus must fix the surface gravity, angular velocities and any other moduli. Without fixing these physical data a method such as the Newton method will fail as the linearized operator Δ above will not be invertible. An iterative method will also typically fail to produce a solution, and will likely drift around in the space of solutions not settling on any one, or if it does settle it will be arbitrary which one it has chosen. Here we will employ a simple and numerically stable way to fix the physical data. Suppose we have n moduli that should be fixed, we fix the values of n metric functions at some specific point in the domain. Thus the values of certain metric functions are mapped to physical data, assuming solutions exist with those values.

While this in-going approach is certainly more complicated than the usual method outlined above for a single Killing horizon, its strength is that the generalization to multiple wave mode horizons is straightforward. One must simply ensure the coordinate domain pierces sufficiently into the interior of the black hole to capture all horizons of the effective propagation metrics of the various wave modes. Before applying this approach to black holes in Einstein-aether we firstly illustrate this method explicitly with the toy example of recovering the Kerr solution in the simplest setting of vacuum Einstein gravity.

Suppose we wish to numerically 'find' the Kerr solution in vacuum gravity using this ingoing approach. Then we wish to solve the harmonic Einstein equation (4.1) with no stress tensor. We begin by taking an in-going coordinate chart that will cover the exterior of the black hole and penetrate the future horizon so that smoothness is imposed there. We do this by taking the most general metric ansatz compatible with stationarity and axisymmetry generated by ∂_v and ∂_ϕ respectively;

$$\begin{align}\hfill g& =-\left({T}_{0}+T\right)\mathrm{d}{v}^{2}-\frac{2\left({V}_{0}+V\right)}{{z}^{2}}\,\mathrm{d}v\,\mathrm{d}z+\frac{2U}{z}\,\mathrm{sin}\,\theta \,\mathrm{cos}\,\theta \,\mathrm{d}v\,\mathrm{d}\theta \hfill \\ \hfill & \quad -2\left({W}_{0}+W\right){\mathrm{sin}}^{2}\,\theta \,\mathrm{d}v\,\mathrm{d}\phi +\frac{A}{{z}^{4}}\,\mathrm{d}{z}^{2}+\frac{2F}{{z}^{2}}\,\mathrm{sin}\,\theta \,\mathrm{cos}\,\theta \,\mathrm{d}z\,\mathrm{d}\theta \hfill \\ \hfill & \quad +\frac{2\left({P}_{0}+P\right)}{{z}^{2}}{\mathrm{sin}}^{2}\,\theta \,\mathrm{d}z\,\mathrm{d}\phi +\frac{{B}_{0}}{{z}^{2}}\,{\mathrm{e}}^{S+{\mathrm{sin}}^{2}\,\theta B}\,\mathrm{d}{\theta }^{2}+\frac{2Q}{{z}^{2}}\,\mathrm{sin}\,\theta \,\mathrm{cos}\,\theta \,\mathrm{d}\theta \,\mathrm{d}\phi \hfill \\ \hfill & \quad +\frac{{S}_{0}}{{z}^{2}}\,{\mathrm{sin}}^{2}\,\theta \,{\mathrm{e}}^{S}\,\mathrm{d}{\phi }^{2}\hfill \end{align} \tag{ 4.4 }$$

where ${\mathcal{F}}_{0}=\left\{{T}_{0},{V}_{0},{W}_{0},{P}_{0},{S}_{0},{B}_{0}\right\}$ are fixed functions of z and θ and the unknown functions to be solved for are then the 10 functions $\mathcal{F}=\left\{T,V,W,U,A,P,F,S,Q,B\right\}$, again depending on z and θ. We choose the reference metric $\bar{g}$ to be that above, g, with the unknown functions $\mathcal{F}$ vanishing, so that it is given by the fixed functions ${\mathcal{F}}_{0}$. Here we will restrict our solutions to have a reflection symmetry in the equatorial plane of the black holes. We emphasize that the form of the metric above does not impose the 't–ϕ' orthogonality property, i.e. that the two planes orthogonal to the Killing vectors ∂_v and ∂_ϕ are integrable. As we discuss later, this implies that the above ansatz does not assume horizons will take the form of a Killing horizon.

Here we think of v, θ and ϕ as analogous to those in ingoing Eddington–Finklestein coordinates for Schwarzschild, with M/z being analogous to the usual radial coordinate r in that coordinate system, where M is the black hole mass. Hence we take z = 0 to be the asymptotic infinity ${\mathcal{I}}^{-}$, and the computation domain (given reflection symmetry in the equatorial plane) is then rectangular with z ∈ [0, z_max] and θ ∈ [0, π/2] with θ = 0 the rotation axis and θ = π/2 the equatorial plane. The value z_max should be large enough that this boundary of the domain is entirely contained within the horizon so that smoothness is imposed there provided the functions $\mathcal{F}$ are smooth. If the value z_max is too small so that the boundary is outside (any) horizon then the problem will not be well-posed—it will be analogous to an elliptic problem lacking data on one boundary.

The various powers of z and factors of sin θ and cos θ are included in (4.4) to simplify the various boundary and asymptotic conditions. Regularity of the axis of symmetry and equatorial plane requires the functions ${\mathcal{F}}_{0}$ and $\mathcal{F}$ to be even in θ and in (π/2 − θ). Taking the ${\mathcal{F}}_{0}$ to be smooth functions at the boundary z = 0 with Dirichlet boundary conditions,

$$\begin{equation}{T}_{0}\to 1,\qquad {W}_{0}\to 0,\qquad {V}_{0}\to \mu ,\qquad {P}_{0}\to \mu a,\qquad {B}_{0},{S}_{0}\to {\mu }^{2}\end{equation} \tag{ 4.5 }$$

where μ and a are constants then imposes the requirement that the reference metric be asymptotically flat. We then impose analogous conditions on the metric g by taking $\mathcal{F}$ to be smooth functions at z = 0 that all vanish there.

We now must make a choice for ${\mathcal{F}}_{0}$, which determines the metric ansatz and the reference metric, subject to these boundary conditions. Later we will take the reference metric to be the Kerr solution in appropriate ingoing coordinates by choosing,

$$\begin{align}\hfill \begin{aligned}\hfill {T}_{0}& =1-\frac{2{\mu }^{2} z}{{\Sigma}},\qquad {W}_{0}=\frac{2a{\mu }^{2}z}{{\Sigma}},\qquad {V}_{0}=\mu ,\qquad {P}_{0}=\mu a\hfill \\ \hfill {B}_{0}& ={\Sigma},\qquad {S}_{0}=\frac{{\left({\mu }^{2}+{a}^{2}{z}^{2}\right)}^{2}-{z}^{2}{\Delta}{a}^{2}\,{\mathrm{sin}}^{2}\,\theta }{{\Sigma}}\hfill \end{aligned}\end{align} \tag{ 4.6 }$$

with Δ = μ² − 2μ² z + a² z² and Σ = μ² + a² z² cos²  θ. Here μ is the mass of the reference metric Kerr spacetime and its angular momentum is aμ. However, in the example of finding Kerr itself in vacuum gravity it would be 'cheating' to take Kerr as the reference already! Thus we choose the simpler ${\mathcal{F}}_{0}$;

$$\begin{equation}\begin{gathered}{T}_{0}=1-2z,\qquad {W}_{0}=2az,\qquad {V}_{0}=\mu ,\qquad {P}_{0}=\mu a\\ {B}_{0}={\mu }^{2},\qquad {S}_{0}={\mu }^{2}\end{gathered}\end{equation} \tag{ 4.7 }$$

so that the reference metric is not Kerr, and we have not built in Kerr to our metric ansatz. In either case of ${\mathcal{F}}_{0}$ above we compute the mass and angular momentum of the metric g to be;

$$\begin{align}\hfill \begin{aligned}\hfill M& ={\left.\frac{\mu }{2}{\partial }_{z}{g}_{vv}\right\vert }_{z=0}\hfill \\ \hfill J& ={\left.-\frac{\mu }{2\,{\mathrm{sin}}^{2}\,\theta }{\partial }_{z}{g}_{v\phi }\right\vert }_{z=0}\hfill \end{aligned}\end{align} \tag{ 4.8 }$$

and hence these are determined by the normal gradients of the functions T and W at the z = 0 boundary. For a solution these quantities should be independent of θ, and we indeed see this in all numerical solutions found. We also define the dimensionless 'spin' of the black hole,

$$\begin{equation}j=\frac{J}{{M}^{2}}\end{equation} \tag{ 4.9 }$$

which for Kerr equals a/μ and obeys |j| ⩽ 1. An important point is that since the metric and reference metric agree asymptotically as z → 0, this implies that the vector ξ vanishes there. As discussed above, it is crucial to ensure that ξ may vanish on boundaries of the domain where boundary conditions are imposed.

In our ingoing method we construct a portion of the spacetime interior to the horizons. Smoothness at these and elsewhere is imposed when solving by the Newton method provided our initial guess is smooth since, as discussed above, we expect each update of the metric will preserve smoothness. We generally take this initial guess simply to be the reference metric. While the reference metric above has a particular mass and angular momentum, determined by μ and a, the solutions (which in this simple example should be Kerr) still admit a two parameter moduli space given by their mass M and angular momentum J. The metric and reference metric need not have the same values for these. As discussed above, in order for the problem to be well posed we must constrain the physical data M and J. There are many such ways to do this, but we have found a convenient and reasonably stable one is to simply fix the value of the metric functions T and W at the innermost point of the equatorial plane, z = z_max and θ = π/2. We call the values there T_data and W_data respectively. Using the Newton method to solve the system we have found this much more stable than trying to fix the point in moduli space by constraining the actual mass and angular momentum. Having found solutions for given T_data and W_data one can then tune these to obtain the desired mass and angular momentum by using a Newton method wrapped around the Newton method that solves the harmonic Einstein equation.

Finally one must choose a sufficiently large value of z_max that the coordinates penetrate inside the horizon, but not too near the singularity that the numerical system becomes destabilized. Of course one does not know where this coordinate position of the horizon will be a priori, but in practice it is straightforward by trial and error to find values of z_max.

This method is pragmatic and works well, and accurately reproduces the Kerr family of metrics using the Newton method starting with an initial guess where the metric g is equal to the reference metric $\bar{g}$. We use sixth order finite differencing to compute derivatives of the metric functions in their rectangular domain. Interestingly we have found higher order pseudospectral methods render our scheme numerically unstable, presumably due to their non-local nature. At the asymptotic boundary the metric functions $\mathcal{F}$ are fixed to zero. We choose to include the boundary points θ = 0, π/2 in the domain (one could choose not to) and there we impose Neumann boundary conditions (one could also choose to impose the equations at these points instead) with the exception of the functions T, W at the innermost equatorial point whose values we fix to T_data and W_data. We then impose the harmonic Einstein equation at all interior points and also at the remaining boundary of the domain, z = z_max. Having found a solution we check that it is indeed consistent with having ξ = 0. In this simple vacuum gravity setting 'solitons' with ξ ≠ 0 cannot exist [42]. While there is no proof for general matter that solitons cannot exist, later when we construct black holes in Einstein-aether we will only find solutions with ξ = 0 in the continuum.

We wish to prescribe the dimensionless spin j = J/M² of the black hole solution. However, due to scale invariance of vacuum GR we do not need to fix the overall mass as we will be interested in only dimensionless measures of the solution. We arbitrarily fix the scale choosing μ = 1 in the reference metric and ansatz. We take a to be equal to the desired ratio a = J/M. As discussed above, the metric solution found will generally not have the same values of mass and angular momentum as the reference metric. Instead the values of T_data and W_data will determine the physical data of the solution. Due to the scale invariance we choose to fix T_data to be zero. Then we find solutions varying W_data using a Newton method until the one with the desired dimensionless spin j is found. See appendix D for more details.

We provide details of convergence to the continuum for these Kerr solutions in appendix C, but in summary here using 40 × 40 points for the given reference metric allows a wide family of Kerr solutions to be found with a pointwise metric accuracy of $\sim 1{0}^{-6}$.

Before turning to the phenomenologically allowed couplings we first verify our numerical code by reproducing results for static spherically symmetric black holes found previously in [14, 15] where c₁ is varied taking c₃ = c₄ = 0 and ${s}_{0}^{2}=1$, with the last condition fixing c₂. While in the rest of this paper we take vanishing c_σ motivated by the strong LIGO constraints, in order to compare with these older results we take non-vanishing c_σ . In terms of the irreducible couplings c_ω,σ,θ,a this corresponds to varying c_ω while fixing c_ω = c_σ = c_a , again with ${s}_{0}^{2}=1$ which determines c_θ .

For that parameter range the static black holes tend to Schwarzschild as c_ω → 0, as then all the aether couplings become small. On the other hand taking c_ω → 1⁻ gives strong deviations from GR. It is interesting to then consider intermediate c_ω , so that deviations from GR behaviour are quite small in the static case. Do rotating solutions deviate more from GR behaviour and hence from Kerr?

We show data here for the intermediate value c_ω = c₁ = 0.5. We find qualitatively similar results for other intermediate values where we are able to find solutions. Varying the spin j up to $\sim 0.8$ it is straightforward to find rotating black hole solutions. In figure 3 we show various quantities computed for these Einstein-aether black holes and compared to Kerr for the same mass and spin. We show the ratio of the equatorial radius of the horizon (defined as the circumference divided by 2π) of Kerr to that of the metric horizon in the aether theory, r_GR/r_AE. For zero spin this value is given in [15] as r_GR/r_AE = 0.933 04. We reproduce this value for zero spin, and see that at least out to spins of ≃0.8 this ratio is remarkably independent of the spin, indicating a larger equatorial radius for the aether black holes than for Kerr of the same mass and spin. We also plot the quantity $\sqrt{{A}_{H}/4\pi {r}_{H}^{2}}$, with A_H the horizon area and r_H the equatorial radius, which measures the deviation from sphericity, both for Kerr and the aether black holes. Again we see little variation between these as spin is added, the aether black holes being marginally more spherical than Kerr. Interestingly, taken together these results suggest the aether carries a net negative energy for these solutions, allowing a larger black hole for the same mass and spin than for Kerr. Finally we show the ratio of the ISCO frequency of Kerr to that of the aether black hole with the same mass and spin. Here we see that a greater deviation from GR develops as we increase the spin of the solutions, although it is not a dramatic change in behaviour, at least out to spin ≃0.8. Since the ISCO is closer to the horizon for higher spin, it is natural that the deviation in the frequency is greater for higher spin, since it is more sensitive to the change of horizon radius. Note also that the larger horizon for the aether black holes relative to Kerr naturally leads to a lower ISCO frequency.

Zoom In Zoom Out Reset image size

An interesting possibility is that closed timelike curves may be present in the effective metrics for the various degrees of freedom behind their horizons. A simple class of such curves are those generated by the ϕ circle, such as occur inside the inner horizon of Kerr. We have checked whether such curves occur by examining the component ${({g}_{\text{eff}})}_{\phi \phi }$ which if negative would yield such a closed timelike curve, and found it to be positive within the coordinate domain of all obtained solutions. Of course this does not preclude closed timelike curves of this form further in the interior of the spacetimes than we have constructed, or the existence of more general closed timelike curves than those purely generated by ∂_ϕ .

A black hole horizon is usually defined as the boundary of the causal past of future null infinity, where the causal structure is determined by the spacetime metric. In Einstein-aether theory, each of the wave modes is associated with an effective metric ${g}_{\mu \nu }^{\text{eff}}$ (2.3), which defines a corresponding notion of event horizon ¹⁴ . These horizons are all hypersurfaces that are null with respect to the appropriate metric. For stationary axisymmetric solutions, as considered here, the Killing vectors ∂_v and ∂_ϕ are tangent to any horizon, so the horizon is the hypersurface which, at least locally, can be defined by z = z_H (θ). Then defining the normal one-form,

$$\begin{equation}n\equiv d(z-{z}_{H}(\theta )),\end{equation} \tag{ 5.4 }$$

the condition that the horizon is a null hypersurface is the condition that n is null with respect to g_eff, i.e. ${({g}_{\text{eff}}^{-1})}^{\mu \nu }{n}_{\mu }{n}_{\nu }=0$. ¹⁵

Stationary axisymmetric black hole horizons in GR are typically Killing horizons, that is, there is a Killing vector tangent to their null generators. This is automatically the case for static, spherically symmetric black holes, but for rotating black holes it is a nontrivial property. For the Kerr horizon, for example, the null generators are tangent to the Killing field ∂_v + Ω_H ∂_ϕ , where Ω_H , the angular velocity of the horizon, is a constant. A theorem of Hawking [45] proved that stationary black hole horizons in vacuum or electrovac analytic solutions to Einstein's equation are Killing horizons. That theorem is not relevant to us here, however, since the Einstein-aether field equations are different, and since we have no reason to assume analyticity. However, a different theorem does apply: without invoking field equations, Carter [46] showed that an event horizon must be a Killing horizon if the spacetime is stationary and axisymmetric, with two commuting Killing fields whose integral two-surfaces are themselves orthogonal to two-surfaces. This integrability condition is called the t–ϕ orthogonality property. We find that in our Einstein-aether theory solutions the t–ϕ orthogonality property fails to hold, and the horizons are in fact not Killing horizons! On each wave mode horizon the surface spanned by the two Killing vectors is spacelike with respect to the wavemode metric except at the poles and the equator, and hence does not include the null direction. The null generators are therefore not parallel to a linear combination of the Killing vectors; rather, they are parallel to

$$\begin{equation}\chi ={\partial }_{v}+{{\Omega}}_{\phi }{\partial }_{\phi }+{{\Omega}}_{\theta }{\partial }_{\theta }+{\chi }^{z}{\partial }_{z},\end{equation} \tag{ 5.5 }$$

with nonzero Ω_θ and χ^z .

We expect that in any solution found by our numerical method all the wave mode horizons are captured within the computational domain, since without imposing smoothness at the horizons the system is not expected to be well posed with a unique solution (or possibly a finite set of solutions). Indeed this is borne out in our computations. We generally find solutions only when the coordinate domain extends deep enough to capture all the horizons. If the coordinate domain ends too far out, so that it does not contain all the horizons, the Newton method does not converge ¹⁶ .

The t–ϕ—here actually v–ϕ—orthogonality property discussed above corresponds (according to the Frobenius theorem) to the vanishing of the two four-forms, ${\partial }_{v}^{{\flat}}\wedge {\partial }_{\phi }^{{\flat}}\wedge \mathrm{d}{\partial }_{v}^{{\flat}}$ and ${\partial }_{v}^{{\flat}}\wedge {\partial }_{\phi }^{{\flat}}\wedge \mathrm{d}{\partial }_{\phi }^{{\flat}}$ (where ${({\partial }_{v}^{{\flat}})}_{a}{:=}{g}_{ab}^{\text{eff}}{\partial }_{v}^{b}$). In figure 4 we plot the norm of the Hodge dual of the first of these (which is a scalar), for the metric horizon in the lefthand of the two solutions shown in figure 5. And for the righthand solution in that figure we plot the same quantity for the spin-1 horizon, where we note the quantity is very small for the metric horizon since that solution is close to Kerr. As seen in the plots, this integrability condition fails to hold, so there is no reason the horizon must be a Killing horizon. We see similar non-vanishing behaviour for the second integrability condition. In appendix C both these conditions are checked for the Kerr solution, where it is shown they hold to very good numerical accuracy, despite the use of a non-adapted coordinate system.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The function z_H (θ) is determined by the condition that the normal one-form (5.4) be null, ${({g}_{\text{eff}}^{-1})}^{\mu \nu }{n}_{\mu }{n}_{\nu }=0$, which defines a non-linear first order ordinary differential equation for z_H (θ). It is a quadratic in ${z}_{H}^{\prime }(\theta )$ with two solutions,

$$\begin{equation}{z}_{H}^{\prime }(\theta )=\frac{{({g}_{\text{eff}}^{-1})}^{z\theta }\pm \sqrt{{\left({({g}_{\text{eff}}^{-1})}^{z\theta }\right)}^{2}-{({g}_{\text{eff}}^{-1})}^{zz}{({g}_{\text{eff}}^{-1})}^{\theta \theta }}}{{({g}_{\text{eff}}^{-1})}^{\theta \theta }}.\end{equation} \tag{ 5.6 }$$

At θ = 0 and θ = π/2 the normal one-form reduces to n = dz, so this equation reduces to the condition ${({g}_{\text{eff}}^{-1})}^{zz}=0$, whose solution serves as the initial data for integration of z_H (θ) from pole to equator or vice versa. In principle multiple solutions might exist, in which case the outermost, i.e. the one with the smallest z, would correspond to the event horizon. In practice we find only one root within our coordinate domains. Away from the poles and equator the plane spanned by ∂_v and ∂_ϕ is spacelike with respect to ${({g}_{\text{eff}})}_{\mu \nu }$ so has two orthogonal null directions which correspond to the two roots in the solution for ${z}_{H}^{\prime }(\theta )$. Only one of these null directions limits to the unique null normal at the poles and equator ¹⁷ .

In figure 5 we plot the various horizons within the coordinate domain for an example solution, the spin j = 0.77 solution considered in the previous section. For that solution the chart was taken to extend to z_max = 0.78. The spin-0, 1 and 2 horizons are distinct, and the curves defining them are θ dependent. Note that for those parameters the spin-0 horizon coincides with the metric horizon. In the figure we also present a similar plot for a typical solution taken from a family we construct in the next section 5.3, where the ordering of the wave mode horizons is different, and now the metric horizon coincides with the spin-2 horizon as c_σ is taken to vanish. As discussed in detail later, the spacetime geometry for this solution is close to the Kerr metric with the same spin, and the spin-0 and spin-2 horizons are close to each other and lie close to the corresponding Kerr value 0.625. Interestingly for our choice of reference metric the coordinate positions, z_H (θ), of the various wave mode horizons appear to have a surprisingly weak dependence on the angle θ for all the Einstein-aether black holes we have found in this work.

In figure 6 we plot the angular velocities Ω_ϕ and Ω_θ , corresponding to the horizons shown in figure 5. The presence of a nonzero Ω_θ , and the explicit dependence of Ω_ϕ on θ, indicate that the horizons are not Killing horizons. It goes linearly to zero at the pole and the equator, and in each case has one sign everywhere else. Thus on any given horizon the generators spiral either from the pole to the equator, or vice versa, over an infinite range of Killing time, with exponential behaviour of θ with v at the two ends.

Zoom In Zoom Out Reset image size

Since the rotating black hole horizons are not Killing, the velocity Ω_ϕ is not constant on them. For small spins we expect the variation in Ω_ϕ across the horizon to be quadratic in the spin, rather than linear. This is because the degree of variation should have an expansion in powers of the spin and yet be independent of the sense of the black hole rotation, hence independent of the sign of the spin. We confirm this by computing the spin dependence of the ratio of the value of Ω_ϕ at the pole to that at the equator, for the various wave mode horizons. This quantity is shown in the previous figure 3 for the same solutions that are discussed there. It deviates from 1 with precisely j² dependence for small spins. We see a similar quadratic dependence on spin for the maximum absolute value of the angular velocity Ω_θ , as the above argument also predicts.

While black holes with non-Killing horizons have been numerically constructed in asymptotically AdS spacetimes [37, 43, 47], we believe our solutions are the first instance of black holes which are asymptotically flat and lack t–ϕ orthogonality, and further have non-Killing horizons of the effective metrics for the physical degrees of freedom, including in this case the usual spacetime metric horizon. This was not previously evident for Einstein-aether black holes since a static, spherically symmetric horizon is always a Killing horizon. Presumably it is also not evident in the slow rotation approximation, working to first order in the spin as in [19], since the non-constancy of the angular velocity Ω_ϕ and non-vanishing Ω_θ appear only at second order.

Following our earlier discussion showing that in the phenomenological parameter regions I and II cf figure 1 a solution may take the approximate form of a GR solution for the metric, with a twist free aether 'painted' on top we might expect to find that the metric is close to Kerr. Indeed we will show numerically that this does occur, by examining three families of aether parameters that tend towards region I when taking |c_a | small, and three that tend towards region II. We will call these the region I and region II families respectively. We write

$$\begin{equation}{c}_{a}={\epsilon}\end{equation} \tag{ 5.7 }$$

and then explore the black hole behaviours in the c_a = → 0 limit, to understand whether the behaviour matches the expected behaviour in equation (3.1). In particular we will demonstrate that Kerr with a twist free aether 'painted' on emerges in the → 0 limit, and also we will study the leading O() approach to this limit. For both families I and II we make the choice

$$\begin{equation}{c}_{\omega }=\lambda +{\epsilon},\end{equation} \tag{ 5.8 }$$

so that as → 0 we have c_ω → λ. For the region I families we take

• Family IA: λ = 0.25 and $\frac{{c}_{\theta }}{3}={\epsilon}+5{{\epsilon}}^{2}$
• Family IB: λ = 1.25 and $\frac{{c}_{\theta }}{3}={\epsilon}+5{{\epsilon}}^{2}$
• Family IC: λ = 0.5 and $\frac{{c}_{\theta }}{3}={\epsilon}+3{{\epsilon}}^{2}$

so that for || = O(10⁻⁵) these parameters are within the viable region I. For region II families we take

• Family IIA: λ = 1 and $\frac{{c}_{\theta }}{3}=0.7-{\epsilon}$
• Family IIB: λ = 0.2 and $\frac{{c}_{\theta }}{3}=0.475-{\epsilon}$
• Family IIC: λ = 0.5 and $\frac{{c}_{\theta }}{3}=1+0.5{\epsilon}$

so that for || = O(10⁻⁷) these lie within the viable region II. The families above are representative of a more general behaviour. We have computed black holes for other parameter choices that limit to region I and II for small , and these give the same behaviours in the limit → 0.

For these parameter families we have constructed numerical black hole solution with various spins j starting with = O(1) and reducing it to see the limiting behaviour as → 0. The largest spin we were able to comfortably find for a good range of was j = 0.8. With more effort it would presumably be possible to push to greater spins. The results we present in what follows are for this j = 0.8 case, but we emphasize that we see exactly the same behaviour for the smaller spins as well.

As we construct the solutions we must make sure z_max is sufficiently large to capture all the horizons. In the appendix D we discuss details of the construction of these solutions, and in particular the values of z_max chosen for them. An important point is that, for both the region I and II families, as → 0 the spin-1 speed diverges. Thus the spin-1 horizon becomes a universal horizon in this limit. We find that constructing such black holes becomes numerically harder for small , precisely as to reach a given one must ensure the extent of the coordinate chart, z_max, is large enough to still capture this spin-1 horizon. Typically we find it difficult to extend the families to smaller than $\sim O(1{0}^{-2})$. However, as we shall see, this is certainly sufficient to see the limiting → 0 behaviour emerging. For O(1) values of , and hence away from the phenomenological regime, it also becomes hard to find solutions. In both limits, i.e. small and large , the difficulties seem to arise from the fact that the speeds of the various modes differ significantly, and consequently the coordinate locations of the corresponding horizons become widely separated in our grid. Then the computational domain necessarily covers a large part of the interior of some of the horizons, and we observe it becomes harder to obtain convergence to a solution.

We wish to show the metric (4.4) tends to Kerr in the limit → 0. Due to our choice of the reference metric being Kerr with a = J/M, and since we further fix T_data = 0, then conveniently if the metric indeed tends to Kerr it does so by becoming equal to the reference metric so we may see this simply by showing that the various metric functions $\mathcal{F}$ all tend to zero ¹⁸ . The L²-norm of the metric function T over the coordinate domain is displayed in figure 7 for black holes with spin j = J/M² = 0.8, where

$$\begin{equation}{\Vert}T{{\Vert}}_{2}={\left[{\int }_{0}^{{z}_{\mathrm{c}\mathrm{u}\mathrm{t}}}\mathrm{d}z{\int }_{0}^{\pi /2}\mathrm{d}\phi \vert T{\vert }^{2}\right]}^{\frac{1}{2}}.\end{equation} \tag{ 5.9 }$$

To compute this norm, we first interpolate the function T to sixth order (in accordance with the order of our differencing scheme) and then calculate the integral using Mathematica's NIntegrate function. Figure 7 clearly shows that ||T||₂ → 0 as → 0 for all the families. Furthermore the gradient of the linear behaviour in this log-log plot is precisely consistent with ||T||₂ ∝ for small , agreeing with the 'painting' picture encapsulated in equation (3.1). While we must take different values of z_max for the various solutions in a family as is varied, we choose z_cut ⩽ z_max to be the largest value common to all solutions found so that the norms of these are taken over the same coordinate range and hence can be compared. For the data in the figure we have taken z_cut = 0.6 which is the smallest value of z_max we used in constructing that family of solutions.

Zoom In Zoom Out Reset image size

We see the same behaviour for all the other metric functions for each of the families of solutions with various spin values of j for which we have constructed solutions. Thus we indeed confirm that Kerr emerges in the limit → 0 of the families of black holes we study. We then expect the aether to be 'painted' on to the Kerr spacetime in a twist free configuration as discussed in sections 3.1 and 3.2. We confirm this in figure 8 by plotting the L²-norm of the twist field,

$$\begin{equation}{\Vert}\omega {{\Vert}}_{2}={\left[{\int }_{0}^{{z}_{\mathrm{c}\mathrm{u}\mathrm{t}}}\mathrm{d}z{\int }_{0}^{\pi /2}\mathrm{d}\phi \vert {\omega }_{\mu \nu }{\omega }^{\mu \nu }\vert \right]}^{\frac{1}{2}},\end{equation} \tag{ 5.10 }$$

for the different families for a particular spin j = 0.8. This figure shows the twist vanishing linearly with as we take → 0, again compatible with the 'painting' picture in equation (3.1), and shows that ${u}_{\mu }^{(1)}$, which gives the O() correction to the leading twist free aether ${\bar{u}}_{\mu }$, has non-vanishing twist. The same behaviour is seen for all the spins we have constructed. Note that in this figure we see the expected deviation from the small linear behaviour, corresponding to the leading term in the expansion (3.1), as is increased towards O(1) values. This is seen also for the other quantities we shortly discuss.

Zoom In Zoom Out Reset image size

From our earlier discussion in sections 3.1 and 3.2, as → 0 we expect the aether to obey equation (3.10) in region I and (3.6) in region II. These are different relations, the latter being the expansion free condition, the former allowing non-trivial expansion. Indeed in figure 9 we see the L²-norm of the expansion θ vanishing for the region II families linearly in as → 0, whereas for the region I families it does not vanish.

Zoom In Zoom Out Reset image size

We now proceed to confirm another important aspect of our 'painted' on aether picture. We expect the aether (5.1) in the limit → 0, when ${u}_{\mu }\to {\bar{u}}_{\mu }$, to be universal in the sense that it depends only on the GR spacetime that it is 'painted' onto, ${\bar{g}}_{\mu \nu }$, and whether one is in region I or II, but does not depend on the O(1) couplings that define the theory in that region in the → 0 limit. Recall this is because the aether potential is determined by the equation (3.10) for region I in the limit → 0, or equation (3.6) for region II in that limit, and these are explicitly independent of all the couplings. Furthermore for these stationary black holes we can also determine the answer to our earlier question of whether there is a moduli space of solutions for the limiting aether ${\bar{u}}_{\mu }$ or not. Our numerical solutions demonstrate that as → 0 there exists a unique and universal aether for both region I and II families given our asymptotic boundary conditions. For a given mass and spin, the region I families give the same aether vector field as → 0. Likewise one obtains the same aether for the region II families in this limit, although this is a different aether from that for the region I families. This is suggested in figure 9 in the case of region I families where the non-zero value of the norm of the expansion is the same in the limit → 0 for the different families shown. In figure 10 we confirm this in detail by giving the L²-norm (with the same z_cut as previously) of the difference of the aether vector function H between different region I families. We see these norms vanish as → 0. In that figure we plot the same for different region II families, and again see that the function H tends to the same form as → 0. The same holds true for the other aether functions K, X and Y. The expansion free condition associated to region II given in equation (3.6), and likewise equation (3.10) for region I, locally govern the leading twist free aether behaviour in a manner independent of the large aether couplings. However there could have been global data which would be reflected in the different families giving different limiting forms of the aether. Since we see the same limiting aether for all the families within region I, and correspondingly for region II, in fact it seems that regularity of the various aether horizons is sufficient to uniquely fix the solution. Thus, the aether behaviour in the limit → 0 is universal, i.e. different for region I or II but not depending on the large aether couplings.

Zoom In Zoom Out Reset image size

Finally, we compare with some analytical results available in the literature. Recall that for the region II families in the static spherically symmetric case we expect that as → 0 the solution tends to the exact analytic solution with Schwarzschild metric as in equation (3.14) and aether as in equation (3.15). If we compute these static family II solutions, and rescale them to have mass M = 1, so that r₀ = 2, then in the limit → 0 the aether function H should take the form,

$$\begin{equation}{H}_{\text{exact}}=-\sqrt{1-2z+{r}_{\text{\ae}}^{4}{z}^{4}}\end{equation} \tag{ 5.11 }$$

in our coordinate system for the resulting Schwarzschild metric, which is related to that in equation (3.14) by z = 1/r. As discussed earlier, the parameter r_æ is determined as ${r}_{\text{\ae}}=\frac{{3}^{3/4}}{2}$ (for r₀ = 2) to ensure existence of a regular spin-0 and spin-1 horizon (or equivalently ensure existence of a universal horizon). We indeed see precisely this behaviour emerge as → 0. For family IIA we plot in figure 11 the L²-norm of the difference H − H_exact for solutions as we vary , where we scale the z coordinate for each solution so that it has mass M = 1, and we use a cut-off z_cut = 0.6 for the norm in this rescaled coordinate (note that as → 0 the metric horizon is at z = 0.5 so the norm is computed including some interior to this horizon).

Zoom In Zoom Out Reset image size

The comparison with analytic solutions can be extended a bit further, to first order in the spin of rotating solutions. In [19] the linear in spin correction to a static solution was considered analytically in the case c_a = 0. For c_σ = 0 this is precisely the slowly spinning correction to our limiting family II static solution. In this case the metric found in [19] is the Kerr metric, to linear order in the spin. As discussed above, we indeed find Kerr in the → 0 limit for the family II solutions for all spins (as expected due to the vanishing of the aether stress tensor in this limit for an aether orthogonal to a maximal foliation), so our numerical solutions for the metric match the approximate analytic ones for slow rotation. As for the aether covector u_a , reference [19] found that it is not perturbed at linear order in the spin. This is also consistent with our numerical limiting family II solution. As discussed in the next section (see figure 12), the u_ϕ aether component vanishes as → 0 for any spin. This is a consequence of the aether being orthogonal to an axisymmetric maximal foliation, hence in particular orthogonal to ∂_ϕ . The remaining components, u_v , u_z and u_θ , should be even functions of spin (as they should be invariant under flipping the spin j → −j), and hence have no linear deviation from their static form, which is indeed what we see numerically.

Zoom In Zoom Out Reset image size

As discussed above we have seen the limiting behaviour we expect for small , namely, an aether 'painted on' to a near-Kerr spacetime, independent of the coupling constants for regions I and II. In this subsection we shall characterize the configuration of these universal aether fields.

For reasonable spins, such as j = 0.8 as presented above, it becomes difficult numerically to find the black hole solutions when ≲ 0.01, as one must solve the system far inside the metric horizon as the spin-1 wavespeed diverges (as does the spin-0 wavespeed for region II). We believe our failure to find solutions is a technical problem due to stability of the Newton solver, rather than a reflection that such solutions do not exist. To find solutions truly in region I we require ∼ 10⁻⁵ or less, and for region II we need ∼ 10⁻⁷, and this is not possible with our numerical code. However already for ∼ 10⁻² we clearly see the solutions take the form of the expansion in in equation (3.1), and we can accurately extrapolate the leading behaviours ${\bar{g}}_{\mu \nu }$ and ${\bar{u}}_{\mu }$. We have seen the leading metric ${\bar{g}}_{\mu \nu }$ is Kerr. By linear extrapolation of our solutions we can deduce the leading aether ${\bar{u}}_{\mu }$ for both the region I and region II families as a function of the Kerr parameters. We obtain the same result, up to expected numerical accuracy, for all the region I families, confirming the universal behaviour. Likewise for the region II families such extrapolation yields the same aether ${\bar{u}}_{\mu }$.

This twist free leading aether takes the form $\bar{u}=k\partial f$. Recalling our aether ansatz in equation (5.1) and the asymptotic boundary conditions the potential must take the form

$$\begin{equation}f=v+{\Phi}(z,\theta )\end{equation} \tag{ 5.12 }$$

and so we should have,

$$\begin{equation}k=H,\qquad {\partial }_{z}{\Phi}=-\frac{X}{{z}^{2}H},\qquad {\partial }_{\theta }{\Phi}=\frac{Y\,\mathrm{sin}\,\theta \,\mathrm{cos}\,\theta }{zH}\end{equation} \tag{ 5.13 }$$

in the → 0 limit. We note that the requirement that the aether asymptotically tends to ∂_v implies there is no linear term in ϕ, so that the aether congruence consists of zero angular momentum worldlines. By linearly extrapolating the functions H, X and Y for our small solutions we can accurately deduce this limiting form for → 0. Comparison of the aether ansatz with the form $\bar{u}=k\partial f$ also implies that in the limit → 0 the aether function K should vanish, and we indeed see this in the numerical solutions. In figure 12 we display the L²-norm of the function K in the aether co-vector, confirming that indeed K vanishes in the → 0 limit.

We have also checked from the numerical solutions that the integrability condition, ${\partial }_{\theta }\left(\frac{X}{{z}^{2}H}\right)+{\partial }_{z}\left(\frac{Y\,\mathrm{sin}\,\theta \,\mathrm{cos}\,\theta }{zH}\right)=0$, holds in the limit → 0, see figure 13. Using our small solutions to extrapolate the potential, in figure 14 we show the z and θ derivatives of the potential, ∂_z Φ and ∂_θ Φ, over the coordinate domain for solutions with j = 0.8 for region I and for region II. We see that for both region I and region II the radial derivatives, ∂_z Φ, are far greater than the angular ones, ∂_θ Φ, in magnitude. One can see by eye that, as expected, the radial dependence is different in detail for region I and region II. We note that ${\bar{g}}_{\mu \nu }$ is the same Kerr metric in the same coordinates for both region I and II solutions, and hence this difference is physical and not related to a difference of coordinates. There is a small θ dependence in the aether as we see from the plots of ∂_θ Φ, and it is again quite different for the two regions.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

It is not easy to see from the aether potential how the aether vector behaves in the → 0 limit. Aether worldlines are determined by the 'upstairs' components of the aether, u^μ , and while the metric approaches Kerr for → 0, the off-diagonal components of the metric entail a complicated relationship of these 'upstairs' aether components with the derivatives of the aether potential. Thus in figure 15 we plot these upstairs aether vector components for the region I family black holes with spin j = 0.8 extrapolated to → 0 using our small solutions. We note that in this → 0 Kerr limit our coordinate z = 1/r where r is the corresponding Boyer–Lindquist radial coordinate, and θ is the same as the Boyer–Lindquist coordinate. We see that u^v and u^z are everywhere positive, so dz/dv > 0 along the aether worldlines, and hence they 'fall into' the black hole spin-0 and spin-2 horizons. It is very likely they also fall into the spin-1 horizon too, although our smallest solution, with = 0.01, is constructed with a coordinate domain that does not extend as far as is required to cover the extrapolated position of the spin-1 horizon in the → 0 limit, and hence we cannot say for sure. We have computed the aether worldlines and do not display them on these diagrams as their motion in the (z, θ) plane appears indistinguishable by eye from a horizontal line, due to the relative smallness of u^θ . We find very similar behaviour for the region II family in the → 0 limit, although the aether components are quantitatively different.

Zoom In Zoom Out Reset image size

In figure 16 we show the limiting shear tensor (squared), σ_μν σ^μν , and the limiting acceleration (squared), a_μ a^μ , for families IA and IIB for black holes with j = 0.8. We have obtained the limiting shear and acceleration by pointwise quadratic extrapolation of the same objects calculated for finite values of . As the figure shows, both σ_μν σ^μν and a_μ a^μ are smooth and look quite similar in regions I and II. In particular, these objects are O(1) and develop non-trivial θ-dependence near the horizons, while they decay quite fast near infinity.

Zoom In Zoom Out Reset image size

Recall that for region I the limiting aether obeys the equation (3.10). As discussed earlier this may be simply satisfied if ${\bar{j}}^{\mu }=0$. Here we show that the limiting region I aether has a non-vanishing ${\bar{j}}^{\mu }$, confirming that the fourth order equation (3.10) for the twist potential is the true equation of motion, rather than the simpler condition ${\bar{j}}^{\mu }=0$. From the earlier discussion we may compute ${\bar{j}}^{\mu }$ two different ways. Firstly we may compute it using equation (3.8) as;

$$\begin{equation}{\bar{j}}^{\mu }=\underset{{\epsilon}\to 0}{\mathrm{lim}}\,\frac{{c}_{\omega }\left({\nabla }_{\nu }{\omega }^{\nu \mu }+{a}_{\nu }{\omega }^{\nu \mu }\right)}{{\epsilon}}\end{equation} \tag{ 5.14 }$$

and secondly we may compute it using equation (3.9) as,

$$\begin{equation}{\bar{j}}^{\mu }=\underset{{\epsilon}\to 0}{\mathrm{lim}}\left(-\left({\nabla }^{\mu }\theta +{u}^{\mu }u\cdot \nabla \theta \right)+\left(u\cdot \nabla {a}^{\mu }-{u}^{\mu }{a}^{2}+\frac{2}{3}\theta {a}^{\mu }-{\sigma }^{\mu \nu }{a}_{\nu }\right)\right).\end{equation} \tag{ 5.15 }$$

We compute $\bar{j}$ both ways, using multiple small solutions to extrapolate to compute the → 0 limits. Both methods agree, providing a further check on our numerical code. Since the form of the aether in this limit, $\bar{u}$, appears to be universal and hence the same for the different families IA, B and C we see the same function $\bar{j}$ for each. In figure 17 we show ${\bar{j}}^{2}$ for these families and see clearly that it is non-vanishing and O(1).

Zoom In Zoom Out Reset image size

A hypersurface everywhere orthogonal to the aether, and lying inside the fastest wave mode horizon for a black hole, is called a universal horizon, because it is a causal barrier for influences of any speed with respect to the aether rest frame. Universal horizons would be relevant if higher spatial derivatives are present in the field equations, producing unbounded propagation speeds. Such higher derivative terms are present for example in Hořava–Lifshitz gravity, and in that theory the aether is always hypersurface orthogonal, so black holes in that theory presumably always possess universal horizons. The presence of universal horizons in static spherically symmetric black holes in Einstein-aether theory is also generic, since in spherical symmetry the aether is always hypersurface orthogonal, and indeed they were found in numerical solutions of this sort [15]. For rotating black holes, however, their existence is not guaranteed. In fact, already in [19] it was shown perturbatively that slowly rotating Einstein-aether black holes have no universal horizons. Given our numerical solutions for rotating black holes, what can be said at the nonperturbative level?

If a universal horizon does exist then, since the aether is orthogonal to that hypersurface, its twist must vanish there. To rule out the presence of a universal horizon, it thus suffices to show that the twist is nowhere vanishing. This is the method that was used in the perturbative analysis of [19]. For the various rotating black holes we have constructed here, this criterion allows us to conclude that there is no universal horizon within our coordinate domains, which do cover at least the spin-0, 1, and 2 horizons. If a universal horizon were to be present it would of course lie inside these horizons, however, and we cannot presently rule out the possibility that one might be lurking deeper in the interior. Nevertheless, we do have some evidence against that possibility.

Recall that for the 'painted' on aether the twist vanishes as → 0, going as, ${\omega }_{\mu \nu }={\epsilon}{\omega }_{\mu \nu }^{(1)}+\cdots \,$. In figure 18 is plotted the squared norm of the leading twist ${\omega }_{\mu \nu }^{(1)}$ for the region I and region II families, extrapolated from our small solutions. This quantity takes a slightly different form for each family of solutions, but common to all of them is firstly that it does not vanish in the interior of our coordinate domain, and secondly this leading twist grows larger at points deeper in the black hole interior. If it continues to grow, there can be no universal horizon. In the limit c_a = → 0, on the other hand, the aether becomes twist free. We expect that in that limit a universal horizon arises, because the spin-1 (and spin-0 for region II) wave speed diverges, so the horizon for that wave mode becomes a universal horizon.

Zoom In Zoom Out Reset image size

In section 4.5 we demonstrated our in-going method using the toy problem of 'finding' the Kerr solution in vacuum GR. Here we present several convergence tests that provide evidence that our numerical solutions indeed converge to Kerr.

In the first test, we show that each of the functions that we solve for in our code converges, pointwise, to the continuum limit with a convergence order which is close to 6, as it should be given our differencing stencils. To see this, for each of the metric functions, for example A, we consider the differences between solutions computed at some resolution (N_z , N_θ ) = (N, N), at (2N, 2N) and at (3N, 3N). We denote these solutions A_N , A_2N and A_3N respectively. Let us define,

$$\begin{equation}{L}_{N}^{(A)}={\Vert} {A}_{N}-{A}_{2N}{\Vert} ,\qquad {R}_{N}^{(A)}={\Vert} {A}_{2N}-{A}_{3N}{\Vert} \times \frac{{h}_{N}^{Q}-{h}_{2N}^{Q}}{{h}_{2N}^{Q}-{h}_{3N}^{Q}},\end{equation} \tag{ C.1 }$$

where the norm above is taken to be the L₁-norm of the function values at the grid points at the resolution (N, N) (which are common to the higher resolutions since the grids are uniform), and where h_N , h_2N , h_3N are the grid spacing for those resolutions. Then if the numerical scheme converges with an order of convergence Q, asymptotically for increasing resolution we would expect,

$$\begin{equation}\underset{N\to \infty }{\mathrm{lim}}\left({L}_{N}^{(A)}-{R}_{N}^{(A)}\right)\to 0\end{equation} \tag{ C.2 }$$

and similarly for the other nine metric functions $\mathcal{F}=\left\{T,V,W,\dots \right\}$. To estimate the order of convergence from our solutions in figure 19 we show L₂₀ and R₂₀ for all the metric functions for convergence orders Q = 2, 4, 6, 8 for the solution with spin j = 0.70 and z_max = 0.65. This compares the resolutions (N_z , N_θ ) = (20, 20), (40, 40) and (60, 60). As this figure shows, the order of convergence for all variables is very close to 6. We see similar behaviours for other spins.

Zoom In Zoom Out Reset image size

Another relevant test that we can carry out to check that the geometry constructed with our method approaches the Kerr spacetime is to verify that it enjoys the v–ϕ orthogonality property, i.e. that the planes spanned by the Killing directions v–ϕ are integrable. Recall that our metric ansatz does not have this property built in so recovering it in the continuum limit is non-trivial. Following our discussion in section 5.2, in figure 20 we show that $\star ({\partial }_{v}^{{\flat}}\wedge {\partial }_{\phi }^{{\flat}}\wedge \mathrm{d}{\partial }_{v}^{{\flat}})$ approaches zero in the continuum limit with a slope that is consistent with our differencing scheme. We observe a similar behaviour for $\star ({\partial }_{v}^{{\flat}}\wedge {\partial }_{\phi }^{{\flat}}\wedge \mathrm{d}{\partial }_{\phi }^{{\flat}})$. The vanishing of these two quantities implies that indeed the continuum spacetime has the v–ϕ orthogonality property.

Zoom In Zoom Out Reset image size

In figure 21 we display the L₁-norm of $\sqrt{\vert {\xi }^{\mu }{\xi }_{\mu }\vert }$ evaluated in our computational domain for resolutions (N_z , N_θ ) = (N, N) as N is varied. We compute the norm as,

$$\begin{equation}{\Vert}\sqrt{\vert {\xi }^{\mu }{\xi }_{\mu }\vert }{{\Vert}}_{1}={\int }_{0}^{{z}_{\mathrm{max}}}\mathrm{d}z{\int }_{0}^{\pi /2}\,\mathrm{d}\phi \sqrt{\vert {\xi }^{\mu }{\xi }_{\mu }\vert }\end{equation} \tag{ C.3 }$$

first constructing the function $\sqrt{\vert {\xi }^{\mu }{\xi }_{\mu }\vert }$ from our numerical solution using sixth order interpolation, and then computing the integral using Mathematica's NIntegrate function. Recall that we are solving the Einstein–DeTurck equation (4.1) and we need ξ^a → 0 in the continuum limit to recover an actual solution of the Einstein equations. This figure indicates that indeed $\sqrt{\vert {\xi }^{\mu }{\xi }_{\mu }\vert }\to 0$ in the continuum limit. Moreover, the slope of the fit turns out to be $\sim -5.93$, and hence close to the order of our finite differencing scheme.

Zoom In Zoom Out Reset image size

Finally in figure 22 since we know the numerical solution should tend to Kerr in the continuum, we study the convergence of the radius of the horizon S² at its equator to its true value for the Kerr black hole. This is a geometric quantity which is gauge invariant in the symmetry class of spacetimes that we consider in this paper. For the actual Kerr black hole, this quantity turns out to be 2M, where M is the mass parameter. The figure shows data for spin j = 0.70, but we observe similar behaviour for other spins. In order to compute this quantity the equatorial coordinate position of the horizon must be found, and we do this by interpolating the metric functions to sixth order (in accordance with our differencing order) in Mathematica and then using the FindRoot function to locate the horizon. As this plot indicates, our numerical solutions converge to the analytic solution. The fit to our numerical data has slope $\sim -3.38$, which would suggest a convergence order less than four (but still greater than two, which is the highest derivative order in the Einstein equations). Note however that this convergence reflects not only that of the underlying numerical solution of the harmonic Einstein equation, but also the subsequent post-processing to locate the horizon. It is therefore not surprising that the convergence order is less than that of the sixth order differencing of the p.d.e.s.

Zoom In Zoom Out Reset image size

In this subsection we present results of the convergence tests for some representative rotating Einstein-aether black holes. We show results for a black hole with j = 0.8 in family IA with = 0.10 and z_max = 0.75 and also one with = 0.61 and z_max = 0.63 as two representative examples. The latter is the largest value of for which we could construct black holes with this spin in this family. We have considered other families, other values of , and also different spins and find similar results.

We first consider the pointwise convergence of the 15 metric and aether functions, $\mathcal{J}$, by computing the L₁ norms of their differences for various resolutions as for the Kerr toy example using the resolutions (N_z , N_θ ) = (20, 20), (40, 40) and (60, 60). This is shown in figure 23. For the solution with = 0.10 (top plot) we observe a sixth order or better convergence. For the = 0.61 solution we see most metric functions give sixth order convergence, although some seem to display slightly worse fourth order convergence, although we note that since this is the largest where solutions were found for this family and spin this is perhaps not surprising. Also as mentioned in the main text, we observe that the order of convergence is sensitive to the relative location of the various horizons and also the extent of our computational domain in the radial direction, z_max.

Zoom In Zoom Out Reset image size

In figure 24 we show the L₁ norm of $\sqrt{\vert {\xi }_{\mu }{\xi }^{\mu }\vert }$ for different resolutions for the same Einstein-aether black holes in family IA. For the solutions with = 0.1, the convergence order of this quantity is better than six, while for the solutions with = 0.61, the convergence order is slightly less than four. We note for the black holes with = 0.1, the quantity $\sqrt{\vert {\xi }_{\mu }{\xi }^{\mu }\vert }$ obtained from solutions with 50 × 50 grid points or more is already affected by finite machine precision (for our C implementation we represent real numbers using double precision).

Zoom In Zoom Out Reset image size