Inflation followed by Belinskii-Khalatnikov-Lifshitz collapse inside accreting, rotating black holes

TL;DR

This paper provides numerical evidence that mass inflation inside accreting, rotating black holes leads to a Belinskii-Khalatnikov-Lifshitz oscillatory collapse, resulting in a spacelike singularity, by simulating all gravitational degrees of freedom.

Contribution

It is the first to numerically demonstrate the transition from mass inflation to BKL collapse inside rotating black holes using full gravitational dynamics.

Findings

Mass inflation is followed by BKL oscillatory collapse.

Collapse results in a spacelike singularity.

Full gravitational degrees of freedom are involved in the simulation.

Abstract

Numerical evidence is presented that the Poisson-Israel mass inflation instability at the inner horizon of an accreting, rotating black hole is generically followed by Belinskii-Khalatnikov-Lifshitz oscillatory collapse to a spacelike singularity. The computation involves following all 6 degrees of freedom of the gravitational field. To simplify the problem, the computation takes as initial conditions the conformally separable solutions of arXiv:1010.1269, arXiv:1010.1271 just above the inner horizon of a slowly accreting, rotating black hole, and integrates the equations inward along single latitudes.

Figures

Equations

$$d{s}^2=-d{T}^2+\sum _{a=1}^3(a_{a}d{X}^a{)}^2,$$

$$a_a\propto |T{|}^{q_a},$$

$$\sum _{a=1}^3{q}_a=1,\sum _{a=1}^3(q_a{)}^2=1.$$

$$q_1=\frac{-u}{1+u+u^2},q_2=\frac{1+u}{1+u+u^2},q_3=\frac{u(u+1)}{1+u+u^2}.$$

$$-\frac{1}{3}\le q_1\le 0\le q_2\le \frac{2}{3}\le q_3\le 1.$$

$$\mathbf{e}$$

$$\mathbf{\Gamma }$$

$$\mathbf{\pi }\equiv -\mathbf{e}\wedge \mathbf{\Gamma }.$$

$$\mathbf{S}\equiv \mathrm{d}\mathbf{e}+\frac{1}{2}[\mathbf{\Gamma },\mathbf{e}]$$

$$\mathbf{\Pi }\equiv \mathrm{d}\mathbf{\pi }+\frac{1}{2}[\mathbf{\Gamma },\mathbf{\pi }]-\frac{1}{4}\mathbf{e}\wedge [\mathbf{\Gamma },\mathbf{\Gamma }]$$

$$\mathbf{e}={\mathbf{e}}_{\bar{t}}+{\mathbf{e}}_{\bar{\alpha }},\mathbf{\pi }={\mathbf{\pi }}_{\bar{t}}+{\mathbf{\pi }}_{\bar{\alpha }}.$$

$${\mathbf{S}}_{\bar{t}}$$

$${\mathbf{\Pi }}_{\bar{t}}$$

$${\mathbf{S}}_{\bar{\alpha }}$$

$${\mathbf{\Pi }}_{\bar{\alpha }}$$

$$0<v<u\ll 1.$$

$$\displaystyle ds^{2}=\rho_{\rm s}^{2}e^{2vt}\Biggl{\{}\left[{e^{-5\xi}dr^{2}\over R^{4}\Delta_{r}}-{e^{\xi}\Delta_{r}\over(1-\omega_{r}\omega_{\theta})^{2}}\left(dt-\omega_{\theta}d\phi\right)^{2}\right]$$

$$\displaystyle\mbox{}+e^{-2\xi}\left[d\theta^{2}+{\Delta_{\theta}\over(1-\omega_{r}\omega_{\theta})^{2}}\left(d\phi-\omega_{r}dt\right)^{2}\right]\Biggr{\}}\ ,$$

$${\rho }_{\mathrm{s}}=\sqrt{r^2+a^2{\cos }^2\theta },$$

$${\omega }_r=\frac{a}{R^2},R\equiv \sqrt{r^2+a^2},{\omega }_{\theta }=a{\sin }^2\theta ,$$

$${\Delta }_r\stackrel{r\gg r_{-}}{⟶}{\mathring{\Delta }}_r=\frac{1}{R^2}\left(1-\frac{2Mr}{R^2}\right),{\Delta }_{\theta }={\sin }^2\theta .$$

$$d{s}^2=e^k{}_{\mu }{e}^l{}_{\nu }{\mathbf{\gamma }}_k\cdot {\mathbf{\gamma }}_{l}d{x}^{\mu }d{x}^{\nu },$$

$${\mathbf{\gamma }}_{\stackrel{v}{u}}\equiv \frac{1}{\sqrt{2}}({\mathbf{\gamma }}_r\pm {\mathbf{\gamma }}_t),{\mathbf{\gamma }}_{\pm }\equiv \frac{1}{\sqrt{2}}({\mathbf{\gamma }}_{\theta }\pm i{\mathbf{\gamma }}_{\varphi }).$$

$$d{r}^{\ast }=\frac{e^{-3\xi }dr}{R^2{\Delta }_r},$$

$$\frac{d\xi }{d{r}^{\ast }}$$

$$\frac{dU}{d{r}^{\ast }}$$

$$\frac{d\ln |R^4{\Delta }_r/r|}{d{r}^{\ast }}$$

$$\xi =\frac{1}{4}\ln \left(\frac{U^2-v^2}{u^2-v^2}\right),$$

$$\frac{(U-v)(u+v)}{(U+v)(u-v)}=e^{4v{r}^{\ast }}.$$

$$T_{kl}=n_k^{+}{p}_l^{+}+n_k^{-}{p}_l^{-},$$

$$n_k^{\pm }=N^{\pm }{p}_k^{\pm },$$

$$p_k^{\pm }{D}^k{p}_l^{\pm }=0,$$

$$D^k{n}_k^{\pm }=0.$$

$$p_1^{\pm }=\frac{P_1^{\pm }}{{\rho }_{\mathrm{s}}\sqrt{{\Delta }_r}{e}^{\xi /2}},p_k^{\pm }=\frac{P_k^{\pm }}{{\rho }_{\mathrm{s}}\sqrt{{\Delta }_{\theta }}{e}^{-\xi }}(k=2,3).$$

$$P_1^{\pm }$$

Full text

I Introduction

During the 1970s, Belinskii, Khalatnikov, and Lifshitz Belinskii et al. (1970); Belinskii and Khalatnikov (1971); Belinskii et al. (1972, 1982) (BKL) developed arguments that the generic outcome of general relativistic collapse to a spacelike singularity would be complicated and oscillatory. BKL’s arguments were asymptotic in nature. They assumed that the spacetime was already deep into collapse towards a spacelike singularity, and they analyzed the behavior in the asymptotic limit as the spacetime approached the singularity.

The BKL scenario was brought into question Ori (1991, 1992, 1999) by Poisson & Israel’s Poisson and Israel (1989, 1990); Barrabès et al. (1990) discovery of the mass inflation instability at the inner horizons of black holes. Mass inflation is the nonlinear consequence of the infinite blueshift at the inner horizon first pointed out by Penrose Penrose (1968). Poisson & Israel argued that cross flow between outgoing and ingoing streams just above the inner horizon would drive exponential growth of the interior mass.

It has been commonly but incorrectly asserted that the generic outcome of inflation is a weak null singularity on the Cauchy horizon Dafermos:2012ny; Luk (2018), not a spacelike singularity. The conclusion that the singularity is weak and null is premised on a black hole that collapses and thereafter remains isolated, in which case the outgoing and ingoing streams that drive inflation are provided by Price tails of gravitational radiation generated during collapse, as originally suggested by Poisson and Israel (1990). However, this premise never holds for an astronomical black hole. The energy density in the most slowly decaying mode (the quadrupole) of Price radiation decays as $t^{-12}$. It is straightforward to estimate that the energy density of a Price tail in the collapse of a stellar-mass black hole falls below the energy density of, for example, the cosmic microwave background in of the order of $10^{4}$ black hole crossing times, or of order 1 second.

As argued by Hamilton and Avelino (2010), during inflation the tidal force grows exponentially with a radial scale length that is inversely proportional to the accretion rate. The reason the singularity is weak and null in the “standard” picture of an isolated black hole is that the accretion rate goes to zero, so the growth rate of the tidal force becomes infinite. The growth takes place over such a small proper time that the tidal force becomes infinite before it has caused the spacetime to deform significantly Ori (1999). This is the weak singularity. But in any real black hole, the accretion rate, however tiny, is never zero. For any finite accretion rate, the outcome of the growing tidal force is deformation leading to collapse, not a weak null singularity.

Even if the accretion rate were vanishingly small, the diverging tidal force associated with the weak null singularity would surely result in diverging pair creation, and such pairs would surely act as an effective source of accretion, again precipitating collapse. Frolov et al. (2006) have shown that in the simplified case of a 1+1-dimensional charged black hole, if the effects of pair creation of charged particles are taken into account, then the result is collapse to a spacelike singularity rather than a null singularity on the Cauchy horizon. The effects of pair creation on mass inflation in black holes in 4 spacetime dimensions has yet to be explored in the literature.

There is a complication to the argument that the outgoing and ingoing streams that drive inflation are dominated by accretion, not by primeval Price tails. Accreted particles must necessarily be ingoing at the outer horizon. Thanks to centrifugal force, freely falling particles that are sufficiently prograde and sufficiently close to the equator become outgoing at the inner horizon, while other particles remain ingoing at the inner horizon. Thus at sufficiently low latitudes, both outgoing and ingoing streams near the inner horizon can be fueled by direct accretion of freely falling particles from outside the outer horizon. However, for rotating black holes, there is a latitude above which all particles that free fall from outside the outer horizon remain ingoing at the inner horizon. The critical latitude is highest for massless particles. For massless particles that free fall from infinity, the latitude above which all particles are ingoing at the inner horizon ranges from $90^{\circ}$ (i.e. all latitudes are accessible) for a non-rotating black hole, to $\sin^{-1}\!\sqrt{-3+2\sqrt{3}}\approx 43^{\circ}$ for a maximally-rotating (extremal) black hole.

Although there is a latitude above which particles freely falling from outside the outer horizon cannot become outgoing, in real astronomical black holes collisional process such as electron-photon scattering will inevitably result in outgoing particles accreting on to the inner horizon at any latitude. Even absent collisional processes, accretion at different rates over the inner horizon will perturb the geometry of the black hole away from a no-hair configuration. The perturbation will result in gravitational waves that relax the black hole back towards a no-hair configuration. These gravitational waves will provide a source of both outgoing and ingoing gravitational radiation at all latitudes on the inner horizon. One way or another, directly or indirectly, in any real astronomical black hole, ongoing accretion will provide the dominant source of outgoing and ingoing energy at all latitudes on the inner horizon.

The purpose of the present paper is to explore numerically inflation and collapse in an accreting, rotating black hole, in fully nonlinear general relativity. Generically this is a challenging problem. One of the challenges is that the numerical method must be robust enough to follow the physical inflationary instability without allowing unphysical gauge instabilities to overwhelm the calculation. Part of the intent of this paper is to test a new covariant Hamiltonian tetrad formalism for general relativity coupled to matter fields, developed by the author Hamilton (2017) to meet the challenges.

To give confidence that the numerical results can be trusted, it makes sense to start with the simplest computation consistent with the goal of following all 6 physical degrees of freedom of the gravitational field through inflation and collapse. The simplest computation is for a slowly accreting black hole, where the geometry is almost Kerr down to near the inner horizon.

In the asymptotic limit of small accretion rates, there exist conformally separable solutions in which inflation occurs Hamilton and Polhemus (2011); Hamilton (2011a, b). The Kerr-Newman geometry has the property, discovered by Carter Carter (1968), that the equations of motion of particles are Hamilton-Jacobi separable. Conformal separability imposes the weaker condition that the equations of motion are Hamilton-Jacobi separable only for massless, not massive, particles. The conformally separable solutions found by Hamilton and Polhemus (2011); Hamilton (2011a, b) are axisymmetric and conformally time-translation invariant (that is, self-similar). Whereas strictly separable solutions (Kerr-Newman) are time-translation invariant and therefore admit no accretion flow, conformally separable solutions allow a conformally growing geometry, and thus admit accretion, hence inflation.

The conformally separable solutions of Hamilton and Polhemus (2011); Hamilton (2011a, b) require a special “monopole” accretion rate, in which the outgoing and ingoing accretion flows on to the inner horizon are uniform in latitude. This might seem to limit the applicability of the solutions. However, because the mass inflation takes place over a proper time that is short for small accretion rates, and becomes shorter as inflation develops, what happens at one point on the inner horizon becomes progressively causally disconnected from what happens at other points. During inflation, outgoing and ingoing accretion flows, regardless of their initial orbital parameters, focus ever more narrowly along the principal outgoing and ingoing null directions. For small accretion rates, the counterstreaming outgoing and ingoing beams become hyper-relativistic while the geometry is still scarcely perturbed from the Kerr geometry. By the time the streams start to back-react on the geometry, the transverse size of the causal patch of an inflating region is tiny. Thus transverse gradients in the accretion flow on to the inner horizon have a subdominant influence. The argument that transverse gradients are subdominant echoes a similar argument by BKL, that spatial gradients become dominated by time gradients during BKL collapse.

The conformally separable solutions predict that inflation is followed by collapse, but the solutions fail at some point during collapse because of growing rotational motions. A motivation for the present paper was to explore numerically what happens then. As will be seen, the result is BKL collapse. The inflationary and collapse phases of the conformally separable solution prove to be the first and second Kasner epochs of BKL collapse.

Units in this paper are geometric, $c=\kappa\,(\equiv 8\pi G)=1$.

II BKL bounces and Kasner epochs

The 3+1 (ADM Arnowitt et al. (1959, 1963)) formalism shows that the 6 physical degrees of freedom of gravity are encoded in the 6 components of the spatial metric, a $3\times 3$ symmetric matrix. The spatial metric can be pictured as an ellipsoid, characterized by the lengths of its three axes, and by three rotation angles. In BKL collapse, the volume of the ellipsoid (determinant of the spatial metric) decreases monotonically to zero in a finite time, but one axis always expands while the other two collapse. When one of the collapsing axes has collapsed to sufficiently small size, it “bounces,” turning from collapse into expansion, while the previously expanding axis turns around and starts collapsing. The directions of the three axes also change during a bounce. The reason for BKL bounces is that the Einstein equations (eq. (8b)) involve a potential energy term proportional to squares of connection coefficients, which provides a repulsive potential that diverges as any axis becomes small. The sensitivity to initial conditions makes the behavior chaotic.

BKL refer to the epochs between BKL bounces as “Kasner epochs” since between bounces the spacetime evolves in accordance with a vacuum solution of Einstein’s equations discovered by Kasner in 1921 Kasner (1921),

[TABLE]

in which the scale factors $a_{a}$ evolve as power laws with time $T$

[TABLE]

with exponents $q_{a}$ satisfying

[TABLE]

A parametric solution for the exponents $q_{a}$ is

[TABLE]

If the exponenents $q_{a}$ are ordered such that $q_{1}\leq q_{2}\leq q_{3}$, then

[TABLE]

III Method

III.1 Numerical method

The numerical method is described by Hamilton (2017), to which the reader is referred for details. The degrees of freedom of the gravitational field are contained in the line interval $\bm{e}$, a vector 1-form with $4\times 4=16$ degrees of freedom, and the Lorentz connection $\bm{\Gamma}$, a bivector 1-form with $6\times 4=24$ degrees of freedom,

[TABLE]

The coefficients $e_{k\kappa}$ of the line interval are commonly called the vierbein coefficients, while the coefficients $\Gamma_{kl{\kappa}}$ of the Lorentz connection are commonly called Ricci rotation coefficients, or, especially when referred to a Newman-Penrose double-null tetrad, spin coefficients. The momentum canonically conjugate to the line interval $\bm{e}$ is the 24-component pseudovector 2-form $\bm{\pi}$ defined by

[TABLE]

The momentum $\bm{\pi}$ is invertibly related to the Lorentz connection $\bm{\Gamma}$. The $16+24=40$ gravitational coordinates and momenta $\bm{e}$ and $\bm{\pi}$ are governed by 40 Hamilton’s equations (eqs. (19) of Hamilton (2017)),

[TABLE]

where $\tilde{\bm{\Sigma}}$ and $\tilde{\bm{T}}$ are the (modified) spin angular-momentum and energy-momentum of matter. Equations (8b) are the Einstein equations. The $\tfrac{1}{4}\bm{e}\operatorname{\wedge}[\bm{\Gamma},\bm{\Gamma}]$ term in equation (8b) is the potential energy term that leads to BKL bounces in BKL collapse.

The gravitational coordinates $\bm{e}$ and momenta $\bm{\pi}$ contain excess degrees of freedom, for two reasons. First, there are gauge degrees of freedom arising from the symmetries of general relativity, namely Lorentz transformations and coordinate transformations; and second, there are redundant degrees of freedom arising from the 4-dimensional covariant description of the coordinates and momenta.

The standard way to remove the excess degrees of freedom from the equations of motion is to perform a 3+1 space+time split of spacetime. The split decomposes the 16-component line interval into 4 time components $\bm{e}_{\bar{t}}$ (subscripted $\bar{t}$), and 12 spatial components $\bm{e}_{\bar{\alpha}}$ (subscripted $\bar{\alpha}$), and the conjugate momentum into 12 time components $\bm{\pi}_{\bar{t}}$ and 12 spatial components $\bm{\pi}_{\bar{\alpha}}$,

[TABLE]

The $\bar{t}$ and $\bar{\alpha}$ subscripts should be interpreted as labels, not indices. After the 3+1 space+time split, the physical degrees of freedom of the gravitational field comprise the 12 spatial components $\bm{e}_{\bar{\alpha}}$ of the line interval and the 12 spatial components $\bm{\pi}_{\bar{\alpha}}$ of their conjugate momenta. The 4 time components $\bm{e}_{\bar{t}}$ of the line interval (the lapse and shift) are gauge degrees of freedom arising from symmetry under coordinate transformations, while the 12 time components $\bm{\pi}_{\bar{t}}$ of the momentum comprise 6 gauge degrees of freedom associated with symmetry under Lorentz transformations, and another 6 degrees of freedom that are simply redundant (see §II D.E of Hamilton (2017)).

Thus Hamilton’s equations (8) split into $12+12=24$ equations of motion involving time derivatives, together with 6 Gaussian constraints, 6 identities, and 4 Hamiltonian constraints,

[TABLE]

Equations (10a) and (10b) are equations of motion in the sense that they determine the time derivatives ${\rm d}_{\bar{t}}\bm{e}_{\bar{\alpha}}$ and ${\rm d}_{\bar{t}}\bm{\pi}_{\bar{\alpha}}$ of the gravitational coordinates $\bm{e}_{\bar{\alpha}}$ and momenta $\bm{\pi}_{\bar{\alpha}}$. The time derivative here is the 1-form ${\rm d}_{\bar{t}}\equiv(\partial/\partial T)dT$ where $T$ is a suitable time coordinate. Hamilton (2017) shows how to isolate the 6 Gaussian constraints and 6 identities in equations (10c). The 6 identities define a 6-component traceless spatial tensor that Hamilton (2017) calls the “gravitational magnetic field” since its role in the gravitational equations is similar to that of the magnetic field in the equations of electromagnetism.

III.2 Conformally separable initial conditions

The initial conditions adopted in this paper are those of the conformally separable solution Hamilton and Polhemus (2011); Hamilton (2011a) for a rotating, slowly accreting black hole. The conformally separable solution is not exact, but rather holds in the vicinity of the inner horizon in the asymptotic limit of small accretion rates.

The conformally separable solution is sourced by two collisionless null streams, one outgoing and one ingoing. The solution is parameterized by three dimensionless parameters, which are the spin parameter $a$ of the rotating black hole, and two constants $u$ and $v$. The accretion rates of the outgoing and ingoing streams are proportional to $u\pm v$ (see eq. 27). Physically realistic solutions have positive accretion rates; small positive accretion rates require that the constants $v$ and $u$ satisfy

[TABLE]

The advantage of choosing the initial conditions to be those of the conformally separable solution is that it is possible to start the integration inward into the black hole near the inner horizon, at a stage of evolution where inflation is already well advanced. The behavior of the accretion flow simplifies greatly near the inner horizon, because outgoing and ingoing accretion streams focus along just two directions, the outgoing and ingoing principal null directions, regardless of their initial orbital parameters.

Choosing the initial conditions to be those of the conformally separable solution has the additional merit that agreement between the solution and numerical calculations gives confidence that both are reliable.

With respect to Boyer-Lindquist coordinates $x^{\mu}\equiv\{r,t,\theta,\phi\}$, the conformally separable Kerr line element may be written

[TABLE]

where $\rho_{\rm s}$ is the separable conformal factor

[TABLE]

$\omega_{r}$ and $\omega_{\theta}$ are functions respectively only of $r$ and $\theta$,

[TABLE]

and the horizon function $\Delta_{r}$ and polar function $\Delta_{\theta}$ are similarly functions respectively only of $r$ and $\theta$. In the conformally separable solution, the polar function $\Delta_{\theta}$ always takes its Kerr value, but the horizon function $\Delta_{r}$, which satisfies the equation of motion (19c), equals its Kerr value $\mathring{\Delta}_{r}$ only at radii sufficiently above the inner horizon, $r\gg r_{-}$,

[TABLE]

Inside the outer horizon, the horizon function $\Delta_{r}$ is negative, the radial coordinate $r$ is timelike, and the time coordinate $t$ is spacelike. The horizon function in Hamilton and Polhemus (2011); Hamilton (2011a), written there as $\Delta_{x}$, is $\Delta_{x}=e^{3\xi}\Delta_{r}$ where $\xi$ is the inflationary exponent given by equations (20). The convention for the horizon function adopted by Hamilton and Polhemus (2011); Hamilton (2011a) is more natural from the perspective of separating the Einstein equations, but the present definition has the advantage that the horizon function $\Delta_{r}$ varies slowly (“freezes”) in the collapse regime where the inflationary exponent $\xi$ grows large.

The geometry described by the line element (III.2) is conformally time-translation symmetric (i.e. self-similar), expanding as $e^{vt}$. The radial coordinate $r$ is a conformal (self-similar) coordinate, and likewise $M$ is the conformal mass of the black hole. The proper mass and radius of the black hole increase as $e^{vt}$ as seen by a distant observer.

The conformally separable Boyer-Lindquist line element (III.2) defines not only a metric, but also, through

[TABLE]

a vierbein $e^{k}{}_{\mu}$, and an associated locally inertial tetrad $\bm{\gamma}_{k}\equiv\{\bm{\gamma}_{r},\bm{\gamma}_{t},\bm{\gamma}_{\theta},\bm{\gamma}_{\phi}\}$, whose scalar products form the Minkowski metric, $\bm{\gamma}_{k}\cdot\bm{\gamma}_{l}=\eta_{kl}$. Corresponding to the orthonormal tetrad is a double-null Newman-Penrose tetrad $\{\bm{\gamma}_{v},\bm{\gamma}_{u},\bm{\gamma}_{+},\bm{\gamma}_{-}\}$ defined by

[TABLE]

The nonvanishing scalar products of the Newman-Penrose tetrad are $\bm{\gamma}_{v}\cdot\bm{\gamma}_{u}=-1$ and $\bm{\gamma}_{+}\cdot\bm{\gamma}_{-}=1$. The Boyer-Lindquist tetrad is constructed precisely so that the null directions $\bm{\gamma}_{v}$ and $\bm{\gamma}_{u}$ point along respectively the principle outgoing and ingoing null directions of the black hole.

Outgoing and ingoing principal null geodesics follow $\theta=\mbox{constant}$, $d\phi/dt=\omega_{r}$, and $r^{\ast}{\pm}\,t=\mbox{constant}$, where $r^{\ast}$ is the tortoise coordinate. The tortoise coordinate $r^{\ast}$ satisfies

[TABLE]

which is a function only of conformal radius $r$. The tortoise coordinate $r^{\ast}$ increases inward towards the inner horizon.

Separation of Einstein’s equations Hamilton and Polhemus (2011); Hamilton (2011a) sourced by outgoing and ingoing collisionless null streams leads to the following equations governing the evolution of the inflationary exponent $\xi$ and the horizon function $\Delta_{r}$ in the line element (III.2),

[TABLE]

with initial conditions $\xi=0$, $U=u$, and $\Delta_{r}=\mathring{\Delta}_{r}$. Equations (19a) and (19b) solve to

[TABLE]

Hamilton and Polhemus (2011); Hamilton (2011a) give a solution of equation (19c) for the horizon function $\Delta_{r}$ valid in the approximation that the radius $r$ is frozen at its inner horizon value $r_{-}$, which is valid in the asymptotic limit of small accretion rates. In the present paper, where accretion rates are finite, the horizon function $\Delta_{r}$ is solved instead from its differential equation (19c) with the Kerr initial condition (15), so that the integration can be started from a (small but) finite distance above the inner horizon, and consistency between integrations from different starting radii can be tested.

III.3 Energy-momenta

The accretion flow consists of two null, pressureless, collisionless streams, one outgoing ($+$) and one ingoing ($-$). The tetrad-frame energy-momentum tensor $T_{kl}$ is a sum over the two collisionless streams (see §VII of Hamilton (2011a)),

[TABLE]

where $p^{\pm}_{k}$ and $n^{\pm}_{k}$ are the outgoing and ingoing momenta and number currents

[TABLE]

with $N^{\pm}$ being outgoing and ingoing scalar densities. The equations of motion for each of the two streams are the geodesic equation

[TABLE]

and number conservation

[TABLE]

In the conformally separable initial conditions, the tetrad-frame outgoing and ingoing momenta $p^{\pm}_{k}$ are related to Hamilton-Jacobi parameters $P^{\pm}_{k}$ by

[TABLE]

The Hamilton-Jacobi parameters $P^{\pm}_{k}$ of the outgoing and ingoing collisionless streams are, equation (22) of Hamilton and Polhemus (2011),

[TABLE]

where $\Delta_{r}^{\prime}\equiv-R^{2}d\Delta_{r}/dr$. The outgoing and ingoing densities $N^{\pm}$ in the conformally separable initial conditions are

[TABLE]

In practice, because the conformally separable solution is not exact, the densities and momenta of the collisionless streams in the initial conditions are adjusted slightly so that the energy-momentum tensor agrees as well as possible with the Einstein tensor deduced from the line element (III.2).

Initially the outgoing and ingoing momenta are closely (though not exactly) aligned with the outgoing and ingoing principal null directions; the momenta become misaligned with the principal directions at later times.

The conformally separable solution assumes zero spin-angular momentum $\tilde{\bm{\Sigma}}$, hence vanishing torsion $\bm{S}$, as is the common assumption in general relativity.

III.4 Integrate numerically along a single radial direction

For simplicity, this paper adopts the approach of integrating inward into the black hole along a single radial direction, at a definite latitude. Gradients in the spatial directions $t,\theta,\phi$ are taken to be given by the conformally separable solution. The assumption of axisymmetry means that gradients in the azimuthal direction $\phi$ vanish identically, while the assumption of conformal time symmetry means that derivatives $\partial\ln e_{k\mu}/\partial t=v$ of logarithmic vierbein elements with respect to $t$ are all equal to the constant $v$. Gradients in the latitude direction $\theta$ are nontrivial.

The approximation that angular gradients are those of the conformally separable solution is at least partially tested by the degree to which the Hamiltonian and Gaussian constraints are satisfied.

III.5 Factor vierbein into dynamical and fixed parts

The vierbein $e_{k\mu}$ defined by the conformally separable line element (16) can be written as the product of a dynamical vierbein $\tilde{e}_{k\kappa}$ that is a function only of the timelike coordinate $r$, and a fixed vierbein $\mathring{e}^{\kappa}{}_{\mu}$ whose elements are fixed to those of the parent Kerr black hole,

[TABLE]

The fixed vierbein $\mathring{e}^{\kappa}{}_{\mu}$ with respect to Boyer-Lindquist coordinates $x^{\mu}\equiv\{r,t,\theta,\phi\}$ is

[TABLE]

which serves to align the tetrad with the principle null directions of the black hole. In the conformally separable initial conditions, the dynamical vierbein is the diagonal matrix

[TABLE]

As the numerical integration proceeds inward, the dynamical vierbein $\tilde{e}_{k\kappa}$ ceases to be diagonal.

The numerical integration works with a time variable $T$ related to the timelike radial coordinate $r$ by

[TABLE]

where $\alpha$ is the dynamical lapse, a function only of $r$. With respect to coordinates $x^{\mu}\equiv\{T,t,\theta,\phi\}$, the conformally separable dynamical line interval (30) is

[TABLE]

where the spatial scale factors $a_{a}$ are

[TABLE]

By assumption, the only part of the line element that depends on spatial variables $t,\theta,\phi$ is the fixed part $\mathring{e}^{\kappa}{}_{\mu}$, equation (29). The spatial variation, which is assumed to continue to be given by equation (29) even after the conformally separable solution breaks down, in turn depends on the radius $r$. The radius $r$ is taken to be given by a plausible extrapolation of equation (31),

[TABLE]

with the scale factor $a_{a}$ in each tetrad direction $a=1,2,3$ being taken to be determined by the expansion $\Gamma^{a}_{0a}$ (no sum over $a$) in that direction $a$,

[TABLE]

Equations (34) and (35) are not part of the system of equations (8) governing the evolution of the gravitational field. Rather, equations (34) and (35) are guesses necessitated by the simplifying approximation adopted in this paper that gradients in angular directions may be approximated adequately by (a plausible extrapolation of) the conformally separable solution. The justification for this approximation is the BKL argument that transverse gradients are subdominant to gradients in the time direction during BKL collapse. The validity of the approximation can be tested, at least in part, by the extent to which the (Hamiltonian and Gaussian) constraints are satisfied (see §IV.3).

III.6 First two Kasner epochs

The conformally separable dynamical vierbein (32) has the form of two successive Kasner epochs. As discussed by Hamilton and Polhemus (2011), the traditional inflationary regime introduced by Poisson & Israel Poisson and Israel (1990) corresponds to the situation where the first and second radial gradients $\partial\xi/\partial r$ and $\partial^{2}\xi/\partial r^{2}$ of the inflationary exponent grow large, while the inflationary exponent $\xi$ itself scarcely budges from zero. In this inflationary regime where the inflationary exponent remains close to zero, $\xi\approx 0$, the dynamical vierbein (32) near the inner horizon looks like a Kasner spacetime (1) with exponents

[TABLE]

meaning that the spacetime is collapsing along the radial direction only. For a black hole that continues to accrete, as is always true in an astronomically realistic black hole, the focusing of accretion streams along the principal outgoing and ingoing null directions produces a tidal force that eventually causes inflation to stall, and the geometry to start collapsing along the transverse directions and expanding along the radial direction. This corresponds to the collapse phase of the conformally separable solution, where the inflationary exponent $\xi$ grows large while the horizon function $\Delta_{r}$ freezes to a constant. In this collapse regime, the dynamical vierbein (32) looks like a second Kasner spacetime (1) with exponents

[TABLE]

Thus the conformally separable solution appears to describe the first two Kasner epochs of BKL-like collapse.

The Schwarzschild solution for a non-rotating black hole has no inner horizon. Near its singular surface (the singularity is a spacelike surface, not a point), the Schwarzschild solution resembles a Kasner spacetime with exponents (37), coinciding with those of a rotating black hole in the second Kasner epoch.

III.7 Aligning the Lorentz frame with the principal null frame

The tetrad defined by the conformally separable line element (III.2) is aligned with the principle null directions. An observer at rest in the tetrad frame sees outgoing and ingoing principal null streams to be $180^{\circ}$ apart on the sky (outgoing appears below, ingoing above), equally blueshifted, and equally rotated. The condition that the principle null directions be geodesic (their transverse momenta vanish, $p_{2}=p_{3}=0$) imposes the following 4 conditions on the tetrad-frame Lorentz connections $\Gamma_{klm}$ (related to the components of the bivector 1-form Lorentz connection by $\Gamma_{klm}=e_{m}{}^{\mu}\Gamma_{kl\mu}$),

[TABLE]

The condition that outgoing and ingoing streams along the principle null directions appear equally blueshifted (their energies $p^{0}$ are the same) is

[TABLE]

Equation (39) holds for any spacetime that is conformally time-translation invariant. The condition that the outgoing and ingoing streams along the principle null directions appear equally rotated about the null directions is

[TABLE]

In the conformally separable solution, the principal null directions are geodesic. It is natural to define a “principal frame” by the geodesic continuation of the principal null directions.

III.8 Gauge choices

There are 10 gauge choices to be made, 4 associated with symmetry under coordinate transformations, and 6 associated with symmetry under Lorentz transformations.

The 4 components of the time $T$ component $\tilde{e}_{kT}$ of the dynamical vierbein, commonly called the (dynamical) lapse $\alpha$ and shift $\beta_{a}$, can be treated as gauge variables arbitrarily adjustable under a coordinate transformation. This paper chooses the dynamical lapse to be 1 and the dynamical shift to be zero,

[TABLE]

For the remaining 6 gauge choices, this paper tries 2 different sets of gauge choices, referred to here as principal gauge, and BSSN (Baumgarte-Shapiro-Shibata-Nakamura) Baumgarte and Shapiro (2010); Brown et al. (2012) gauge. BSSN gauge should perhaps be called BSSN-like, because the numerical method followed in this paper is still the covariant Hamiltonian tetrad approach of Hamilton (2017), as opposed to the coordinate-based, Lorentz-gauge-invariant approach commonly referred to as BSSN. The key distinction between principal and BSSN gauges is that whereas principal gauge imposes 6 conditions on the Lorentz connections, BSSN gauge replaces 3 of the 6 gauge conditions by the 3 ADM gauge conditions (42) on the vierbein.

III.8.1 Principal gauge

Principal gauge imposes the 6 conditions (38)–(40) on the Lorentz connections. The conditions ensure that, in rest frame of the tetrad, the (geodesic continuation of the) outgoing and ingoing principal null streams remain $180^{\circ}$ apart on the sky, and equally blueshifted, and equally rotated.

III.8.2 BSSN gauge

BSSN Baumgarte and Shapiro (2010); Brown et al. (2012) gauge imposes the 3 ADM gauge conditions

[TABLE]

on the vierbein. The equations of motion for the 3 components $e_{0\alpha}$ are then reinterpreted as identities, equation (45) of Hamilton (2017). The ADM conditions (42) use up 3 of the 6 Lorentz gauge freedoms on the vierbein. For the remaining 3 Lorentz gauge freedoms, BSSN gauge in this paper imposes the conditions (38b) and (40).

III.9 Integration scheme

The integrator is a 4th-order predictor-corrector (Adams-Bashforth). The predictor-corrector routine calls two subroutines in succession, a subroutine that sets time derivatives of the 12 spatial coordinates $\bm{e}$ and 12 spatial momenta $\bm{\pi}$ in accordance with the equations of motion (10a) and (10b), and a subroutine that applies gauge choices and identities to determine the remaining coordinates and momenta, namely the time components $\bm{e}_{\bar{t}}$ and $\bm{\pi}_{\bar{t}}$.

The routine that sets time derivatives has an option to compute the evolution of the conformal parts of the line interval $\bm{e}$ and its conjugate momentum $\bm{\pi}$ separately from the nonconformal parts. The nonconformal (dimensionless) parts of the spatial line interval and its conjugate momentum are $\bm{e}/|e|^{1/3}$ and $\bm{\pi}/|e|^{1/3}$, where $|e|$ is the determinant of the spatial vierbein. The equations of motion for the conformal parts of $\bm{e}$ and $\bm{\pi}$ are equations (33) and (34) of Hamilton (2017) for the spatial volume element $\bm{e}^{3}$ and the spatial expansion $\bm{\vartheta}\equiv\tfrac{1}{2}\bm{e}\operatorname{\wedge}\bm{\pi}$. The equation for the expansion $\bm{\vartheta}$ has an option to adjust it by adding an arbitrary factor of $\bm{e}_{\bar{t}}$ wedged with the Hamiltonian constraints, equation (35) of Hamilton (2017). In practice, adjusting the equation of motion for the expansion proves to have no appreciable effect on the results. By default, therefore, the results reported in this paper do not treat the conformal parts in any special way: the line interval and conjugate momentum are left unscaled, and the evolution of the conformal parts of $\bm{e}$ and $\bm{\pi}$ is computed from equations (10a) and (10b) just like the nonconformal parts.

After calling the subroutine that sets time derivatives of the spatial coordinates and momenta $\bm{e}$ and $\bm{\pi}$, the predictor-corrector routine calls a subroutine that sets the 4 time components $\bm{e}_{\bar{t}}$ and 12 time components $\bm{\pi}_{\bar{t}}$.

The 4 time components $\bm{e}_{\bar{t}}$, the lapse and shift, are gauge choices, equations (41).

The 12 time components $\bm{\pi}_{\bar{t}}$ and 12 spatial components $\bm{\pi}$ of the conjugate momentum together satisfy 30 equations, including the 12 equations of motion that determine $\bm{\pi}$. The 30 equations may expressed as the matrix equation

[TABLE]

where $C$ is a vector with 18 entries, consisting of

[TABLE]

Given $\bm{\pi}$ from their equations of motion, equation (43) reduces to

[TABLE]

Equation (45) is a set of 18 equations for the 12 time components $\bm{\pi}_{\bar{t}}$ of the conjugate momentum. The 18 equations include 6 redundant equations, the 6 Gaussian constraints. The 6 redundant equations in equations (45) could potentially be eliminated by dropping the constraint equations. That this may not be the best strategy is evidenced by the fact that the constraint equation associated with the ADM gauge condition (42) is formally the BSSN momentum equation, equation (43) of Hamilton (2017), and keeping the BSSN momentum equation is precisely what, apparently, makes the BSSN approach superior to the traditional ADM approach.

Rather than attempt to determine which constraints to drop a priori, this paper adopts the numerical strategy of solving equation (45) for $\bm{\pi}_{\bar{t}}$ by singular value decomposition (SVD).

III.10 Numerical precision

The results shown in this paper are shown up to the point that the numerical integration grinds to a halt. The numerical integration stalls when the matrix equation (45) solved by singular value decomposition becomes ill-conditioned, that is, the ratio of the largest to smallest singular values exceeds the numerical precision.

To mitigate loss of precision, the code, written in c, is implemented in long double. On x86 architecture, a long double uses 80 bits, comprising a sign, 64 bits in the mantissa, and 15 bits in the exponent, whereas an ordinary double uses 64 bits, comprising a sign, 52 bits in the mantissa, and 11 bits in the exponent. The numerical integration stalls when the condition number of the singular-value matrix $A_{\bar{t}}$ exceeds the precision $\sim 2^{64}\approx 10^{19}$.

III.11 Decomposition of the line interval

Only 6 of the 16 components $e_{k\kappa}$ of the vierbein represent physical, coordinate and Lorentz gauge-invariant degrees of freedom of the gravitational field. The 4 time components $e_{kT}$ of the vierbein, the lapse and shift, are gauge variables, adjustable by a coordinate transformation. The remaining 12 components of the vierbein are the spatial components $e_{k\alpha}$. Of these, the physical, Lorentz gauge-invariant degrees of freedom are the 6 components of the symmetric spatial metric $g_{\alpha\beta}=e_{k\alpha}\eta^{kl}e_{l\beta}$.

Because of the block diagonal form (in time and space) of the fixed vierbein $\mathring{e}^{\kappa}{}_{\mu}$ in the decomposition (28) of the vierbein into dynamical and fixed factors, the degrees of freedom of the spatial vierbein $e_{k\alpha}$ are the same as those of the spatial dynamical vierbein $\tilde{e}_{k\alpha}$,

[TABLE]

The physical, Lorentz gauge-invariant degrees of freedom of the dynamical spatial vierbein are those of the dynamical spatial metric $\tilde{g}_{\alpha\beta}=\tilde{e}_{k\alpha}\eta^{kl}\tilde{e}_{l\beta}$. In abbreviated notation

[TABLE]

The $3\times 3$ dynamical spatial metric $\tilde{g}$ can be decomposed as

[TABLE]

where $a$ is a $3\times 3$ diagonal matrix with all positive entries, and $U$ is a $3\times 3$ matrix that can be chosen in a variety of ways: $U$ could be upper triangular with unit diagonals, or lower triangular with unit diagonals, or orthogonal. The $4\times 3$ spatial dynamical vierbein $\tilde{e}$ can then be written

[TABLE]

where $L$ is a Lorentz transformation, satisfying $L^{\top}\!\eta L=\eta$. The right hand side of equation (49) is a product of the $4\times 4$ matrix $L$, the $4\times 3$ matrix \left(\begin{array}[]{c}0\\ a\end{array}\right), and the $3\times 3$ matrix $U$. The Lorentz transformation $L$ contains the 6 redundant Lorentz degrees of freedom of the spatial dynamical vierbein. In BSSN gauge, where $\tilde{e}_{0\alpha}=0$, the Lorentz transformation $L$ reduces to a purely spatial rotation.

The choice of whether $U$ is upper or lower triangular, or orthogonal, has no effect on the numerical integration; the decomposition is needed only to project out the physical degrees of freedom. The choice adopted in this paper is $U$ upper triangular. The 6 physical degrees of freedom of the gravitational field are then the 3 components $a_{a}$ of the spatial diagonal matrix $a$, and the 3 above-diagonal components $U_{ab}$ of the upper-triangular matrix $U$,

[TABLE]

IV Results

IV.1 Parameters

The conformally separable solutions hold in the asymptotic limit of small (but nonzero) outgoing and ingoing accretion rates. But the smaller the accretion rates, the more rapidly inflation exponentiates, and the sooner the integration stalls because the singular-value matrix $A_{\bar{t}}$ in equation (45) becomes ill-conditioned. The accretion rates adopted in this paper are a compromise,

[TABLE]

small enough that the conformally separable solution is a satisfactory approximation (a proposition that can be tested by the degree to which the constraint equations are satisfied), but large enough that the numerical integration continues through several BKL bounces. The outgoing and ingoing accretion rates in the initial conditions are, equation (27), proportional to $u\pm v$,

[TABLE]

The black hole spin adopted in this paper is

[TABLE]

which is chosen to be large, but short of extremal. As noted by Hamilton (2011a), the conformally separable solutions do not admit an extremal black hole, since they require that the inner horizon be separate from the outer horizon.

Numerical experiment indicates that the most “difficult” computation (in the sense that the Hamiltonian constraints are least well satisfied) is at mid latitudes, where gradients in the latitude direction are greatest. To illustrate the results, this paper therefore adopts

[TABLE]

Results at a variety of latitudes are shown later, in §IV.6.

IV.2 BKL behavior

Figure 1 shows the lengths $a_{a}$ and spatial rotation $U_{ab}$, equations (50), of the 3 axes of the dynamical spatial vierbein as a function of the mean scale factor $\bar{a}\equiv(a_{1}a_{2}a_{3})^{1/3}$. The parameters are those of equations (51) and (53). The Figure shows results computed using respectively the BSSN and principal gauges. The two gauges yield similar but slightly different results.

The behavior illustrated in Figure 1 is characteristic of BKL collapse (see §II). The evolution encompasses 4 Kasner epochs during which the scale factors $a_{a}$ evolve as powers of the mean scale factor $\bar{a}$, punctuated by 3 BKL bounces during which the power-law behavior changes abruptly, and at the same time the rotation of the axes changes abruptly.

Figure 2 shows the Kasner exponents $q_{a}$ inferred from the evolution of the scale factors shown in Figure 1,

[TABLE]

By definition of the mean scale factor $\bar{a}\equiv(a_{1}a_{2}a_{3})^{1/3}$, the sum of the exponents is unity, $\sum q_{a}=1$. Figure 2 shows in addition the sum $\sum q_{a}^{2}$ of the squares of the Kasner exponents, which the BKL model predicts should also equal unity, equations (3). The numerically calculated Kasner exponents $q_{a}$ generally conform to the BKL model during Kasner epochs, but deviate during BKL bounces.

The dashed lines shown in Figures 1 and 2 at small mean scale factor $\bar{a}$ are the predictions of the conformally separable solution, continued up to the point where the conformally separable solution is predicted to fail because of growing rotational motions Hamilton and Polhemus (2011); Hamilton (2011a). The numerically computed results are indeed consistent with the predictions of the conformally separable solution approximately up to the expected point of failure.

IV.3 Hamiltonian constraint

A figure of merit for how well the numerical computation is performing is provided by the Hamiltonian constraint, which is the all-spatial component ($123\alpha\beta\gamma$) of the difference between the curvature $\bm{\Pi}$ and ($\kappa\equiv 8\pi G$ times) the (modified) energy-momentum tensor $\tilde{\bm{T}}$, and which should equal zero if the computation is accurate,

[TABLE]

There are of course other constraints; but the Hamiltonian constraint (56) is symptomatic.

Figure 3 shows the Hamiltonian constraint (56) for the same model as in Figures 1 and 2, for the BSSN and principal gauges. Figure 3 also shows for comparison the (modified) energy-momentum tensor $\tilde{T}_{123\alpha\beta\gamma}$, and the potential energy term $\left(\tfrac{1}{4}\bm{e}\operatorname{\wedge}[\bm{\Gamma},\bm{\Gamma}]\right)_{123\alpha\beta\gamma}$ in equation (8b). Figure 3 shows that BSSN gauge performs better than principal gauge, in the sense that the Hamiltonian constraint is relatively smaller than the the energy-momentum.

The Hamiltonian constraint shown in Figure 3 is not zero either in the initial conditions or subsequently. The reason the Hamiltonian constraint is not zero in the initial conditions is that the initial conditions are taken to be those of the conformally separable solution, which is not exact for finite accretion rates. The Hamiltonian constraint gets worse as the integration proceeds because, as described in §§III.4 and III.5, the integration proceeds along a single radial direction at fixed latitude, with gradients in angular directions being taken from the conformally separable solution as opposed to being computed numerically. Thus the growing departure of the Hamiltonian constraint from zero can be attributed at least in part to a failure of the conformally separable solution to yield a good approximation to angular gradients.

Numerical experiment shows that, up to a certain point, the closer to the inner horizon the numerical integration is started, the better the Hamiltonian and Gaussian constraints are satisfied in the conformally separable initial conditions. The constraints are best satisfied when the integration is started around the end of inflation, at approximately the boundary between the first two Kasner epochs. This accounts for the location of the starting point of the numerical integration indicated by the solid lines in Figures 1 and 2. In practice, the radius at which the integration is started is

[TABLE]

which is $r_{\rm init}-r_{-}=0.00002\,M$ above the radius $r_{-}=0.72\,M$ of the inner horizon for spin parameter $a$ given by equation (53). If the integration is started at larger radius, then the Hamiltonian constraints are less well satisfied, and the integration breaks down earlier than shown in Figures 1 and 2.

I experimented, unsuccessfully, with adjusting the initial conditions in an attempt to improve the Hamiltonian constraints one way or another. The conformally separable solutions are approximate, not exact, and within the scope of these solutions it is impossible to satisfy all constraints (Hamiltonian and Gaussian) exactly.

In Principal gauge, the Hamiltonian constraint is least well satisfied around $\bar{a}\sim 10^{-3}$. This coincides with the proper acceleration $\Gamma_{a00}$ of the tetrad frame becoming large, especially in the $a=2$ latitude direction (large compared for example to the acceleration in BSSN gauge). A related circumstance is that in BSSN gauge also around $\bar{a}\sim 10^{-3}$, the principal null directions cease to point mainly along the tetrad 1-direction, as seen in Figure 5, the direction along which the null directions are forced to point in Principal gauge.

IV.4 Radius and horizon function

Figure 4 shows the radius $r$ and horizon function $\Delta_{r}$ in the same model as in Figures 1 and 2. The evolution of the radius and horizon function are governed by equations (34) and (19c).

The radius $r$ plays a role in setting up the conformally separable initial conditions, but, as alluded to after equation (35), thereafter does not enter directly into the either the equations of motion (8) of the gravitational field or (23) and (24) of the collisionless streams. Rather, the radius $r$ is needed only to determine angular derivatives of the fixed vierbein $\mathring{e}^{\kappa}{}_{\mu}$, equation (29) (specifically, the latitude gradient $\partial/\partial\theta$ of the separable conformal factor $\rho_{\rm s}$, equation (13), depends on $r$). In other words, the behavior of the radius $r$ is not a central element of the evolution. Figure 4 serves merely to demonstrate that $r$ does not go crazy.

Similarly, the horizon function $\Delta_{r}$ plays a role in setting up the conformally separable initial conditions, but after the initial conditions are set, nothing in the equations of motion of either the gravitational field or the collisionless streams depends on the horizon function. Figure 4 serves merely to demonstrate that the horizon function $\Delta_{r}$ “freezes” following inflation, as commented after equation (15).

IV.5 Collisionless streams

Figure 5 shows the spatial components $p^{\pm}_{a}$ of the tetrad-frame momenta, and densities $N^{\pm}$, of the outgoing and ingoing collisionless null streams that provide the source of energy-momentum in the model. Also shown are the momenta and densities of principal outgoing and ingoing null streams, designated (PN). There is no actual energy-momentum in the principal streams; the behavior of the densities $N^{\pm}\mbox{\scriptsize{(PN)}}$ is shown only for comparison. The collisionless streams that constitute the source of energy-momentum are initially almost aligned with the principal null directions, but then deviate.

Principal gauge is defined by the requirement that the tetrad frame remains aligned with the principal null directions, §III.8.1, the principal null momenta satisfying $p^{+}_{1}=-p^{-}_{1}$ and $p^{\pm}_{2}=p^{\pm}_{3}=0$. The top-right panel of Figure 5 confirms that the principal null momenta indeed satisfy these conditions.

IV.6 Results at various latitudes

Figure 6 shows the evolution of the scale factors $a_{a}$ at a variety of latitudes, for the same accretion rates (51) and black hole spin (53) as before. The gauge is BSSN; results in principal gauge are similar but not identical. Results are shown for latitudes from the equator almost to the pole. The highest latitude, $89^{\circ}$, is not quite at the pole in order to avoid the coordinate singularity at the pole that occurs when polar coordinates $\theta$, $\phi$ are used, as here. In polar coordinates, the determinant of the vierbein (29) is zero at the pole, and the inverse vierbein is singular. The singularity could presumably be removed by using different coordinates, but that is not done in this paper.

Figure 6 shows BKL-like behavior at all latitudes. At all latitudes except the equator, $0^{\circ}$, the numerical integration continues over 4 Kasner epochs and 3 BKL bounces before grinding to a halt. At the equator, the integration covers just 2 Kasner epochs and 1 BKL bounce before stalling.

Figure 7 shows the Kasner exponents $q_{a}$ corresponding to the models shown in Figure 6. Again, the behavior is consistent with the BKL model, where the exponents $q_{a}$ satisfy the Kasner relations (3) during power-law Kasner epochs.

V Summary and conclusions

The inner horizons of accreting, rotating black holes are subject to the mass inflation instability discovered by Poisson & Israel Poisson and Israel (1989, 1990); Barrabès et al. (1990). During inflation, collisionless accretion streams focus along the outgoing and ingoing principal null directions, their momenta and densities growing exponentially huge.

This paper has explored numerically the behavior of collisionless outgoing and ingoing accretion streams, and of all 6 physical degrees of freedom of the gravitational field, during and following inflation in the vicinity of the inner horizon of an accreting, rotating black hole. The generic outcome appears to be BKL collapse, as originally proposed in the 1970s by Belinskii, Khalatnikov, and Lifshitz Belinskii et al. (1970); Belinskii and Khalatnikov (1971); Belinskii et al. (1972, 1982). The result differs from the common claim that the generic outcome of inflation is a weak null singularity on the Cauchy horizon Ori (1992, 1999); Dafermos:2012ny; Luk (2018). As argued by Hamilton and Avelino (2010), the outcome of a weak null singularity is an artifact of the assumption that a black hole remains isolated for ever after it first collapses, whereas real astronomical black holes are never isolated, but rather accrete.

Part of the aim of this paper was to test numerically the conformally separable solutions for accreting, rotating black holes discovered by Hamilton and Polhemus (2011); Hamilton (2011a). The conformally separable solutions are not exact, but are claimed by Hamilton and Polhemus (2011); Hamilton (2011a) to hold asymptotically in the limit of small but nonzero accretion rates. In the conformally separable solutions, inflation is followed by collapse, but the solutions are expected to break down at small scales when rotational motions become important. The present paper takes the conformally separable solutions as initial conditions for numerical integration. The numerical computations confirm that the conformally separable solutions are valid over the expected range of validity. In fact, the conformally separable solutions appear to describe the first two Kasner epochs of BKL collapse. The numerical integrations in the present paper continue through 4 Kasner epochs punctuated by 3 BKL bounces, before the integrations grind to a halt.

The numerical computations use a covariant Hamiltonian tetrad approach devised by the author Hamilton (2017). Computations have been carried out in two gauges, principal gauge, in which the tetrad is aligned with the principal null directions, and BSSN (strictly, BSSN-like) gauge, which imposes the traditional ADM gauge choice (42). The two gauges yield similar results, although BSSN appears to be superior in the sense that the Hamiltonian constraints are better satisfied by BSSN.

A major simplifying assumption made in this paper, §III.4, is to assume that gradients in spatial directions $t$, $\theta$, $\phi$ continue to be given by the conformally separable solution even after the latter breaks down. The assumption is motivated in part by the BKL argument that spatial gradients are subdominant to gradients in the time direction during BKL collapse. The simplification allows the numerical integration to be carried inward into the black hole along a single radial direction at a single latitude. The approximation may hold approximately in the conformally separable regime (first two Kasner epochs), but certainly fails thereafter, as evidenced both by the growing failure of the Hamiltonian constraint, Figure 3, and by the differences in evolution at different latitudes, Figure 6.

Full exploration of the outcome of the mass inflation instability in accreting, rotating black holes will require calculating spatial gradients in full numerical general relativity. This will be a challenging task.

Acknowledgements.

This research was supported in part by FQXI mini-grant FQXI-MGB-1626.