Excision and avoiding the use of boundary conditions in numerical relativity

TL;DR

This paper introduces a boundary-condition-free method for evolving hyperbolic systems in numerical relativity by excising the domain along a spacelike direction, demonstrated through scalar field simulations.

Contribution

It proposes a novel excision technique that eliminates the need for boundary conditions in numerical relativity simulations, improving stability and simplicity.

Findings

The excision method successfully evolves scalar fields without boundary conditions.

Comparison shows the excision approach is effective and comparable to existing methods.

Numerical experiments demonstrate the method's stability and accuracy.

Abstract

A procedure for evolving hyperbolic systems of equations on compact computational domains with no boundary conditions was recently described in [arXiv:1905.08657]. In that proposal, the computational grid is expanded in spacelike directions with respect to the outermost characteristic and initial data is imposed on the expanded grid boundary. We discuss a related method that removes the need for imposing boundary conditions: the computational domain is excised along a direction spacelike with respect to the innermost going characteristic. We compare the two methods, and provide example evolutions from a code that implements the excision method: evolution of a massless self-gravitating scalar field in spherical symmetry.

Equations

$$g_{\mu \nu }d{x}^{\mu }d{x}^{\nu }={\alpha }_{ab}d{x}^{a}d{x}^b+{\gamma }_{AB}\left(d{x}^A+{\beta }_a^{A}d{x}^a\right)\left(d{x}^B+{\beta }_b^{B}d{x}^b\right),$$

$$g^{\mu \nu }{\xi }_{\mu }{\xi }_{\nu }={\alpha }^{ab}\left({\xi }_a-{\beta }_a^A{\xi }_A\right)\left({\xi }_b-{\beta }_b^B{\xi }_B\right)+{\gamma }^{AB}{\xi }_A{\xi }_B.$$

$$g_{\mu \nu }d{x}^{\mu }d{x}^{\nu }=2{\alpha }_{pw}dpdw+{\alpha }_{ww}d{w}^2+{\gamma }_{AB}\left(d{x}^A+{\beta }_a^{A}d{x}^a\right)\left(d{x}^B+{\beta }_b^{B}d{x}^b\right).$$

$$g_{\mu \nu }d{x}^{\mu }d{x}^{\nu }=2{\alpha }_{uv}dudv+{\gamma }_{AB}\left(d{x}^A+{\beta }_a^{A}d{x}^a\right)\left(d{x}^B+{\beta }_b^{B}d{x}^b\right).$$

$$R_{\mu \nu }-\frac{1}{2}{g}_{\mu \nu }R=T_{\mu \nu },$$

$$T_{\mu \nu }={\nabla }_{\mu }\varphi {\nabla }_{\nu }\varphi -\frac{1}{2}{g}_{\mu \nu }(\nabla \varphi {)}^2,$$

$$\square \varphi =0.$$

$$d{s}^2=-\alpha (t,r{)}^{2}d{t}^2+{\left(dr+\alpha (t,r)\zeta (t,r)dt\right)}^2+r^2\left(d{\vartheta }^2+{\sin }^2\vartheta d{\phi }^2\right).$$

$$Q\equiv {\partial }_r\varphi ,$$

$$P\equiv \frac{1}{\alpha }{\partial }_t\varphi -\zeta Q,$$

$${\partial }_{t}Q-{\partial }_r\left(\alpha \left[P+\zeta Q\right]\right)=0,$$

$${\partial }_{t}P-\frac{1}{r^2}{\partial }_r\left(r^2\alpha \left[Q+\zeta P\right]\right)=0.$$

$${\partial }_r\left(r{\zeta }^2\right)+2r\frac{{\partial }_r\alpha }{\alpha }{\zeta }^2-r^2\rho =0,$$

$$\frac{{\partial }_r\alpha }{\alpha }-\frac{1}{2}\frac{r}{\zeta }{j}_r=0,$$

$$\rho \equiv n^{\mu }{n}^{\nu }{T}_{\mu \nu }=\frac{1}{2}\left(P^2+Q^2\right),$$

$$j_r\equiv -{\gamma }_r{}^{\mu }{n}^{\nu }{T}_{\mu \nu }=-PQ.$$

$$\alpha =1,\zeta =\sqrt{\frac{2m}{r}},$$

$$m_{MS}(t,r)\equiv \frac{r}{2}\left(1-(\nabla r{)}^2\right)=\frac{r}{2}\zeta (t,r{)}^2.$$

$$c_{\pm }=\alpha \left(\pm 1-\zeta \right).$$

$${\partial }_t\zeta -\frac{\alpha }{2r}{\partial }_r\left(r{\zeta }^2\right)-\frac{r}{2\zeta }{T}_{tr}=0.$$

$$\displaystyle\phi\big{|}_{t=0}=a_{0}\left(\frac{r}{w_{0}}\right)^{4}\mathrm{exp}\left(-\frac{(r-r_{0})^{2}}{w_{0}^{2}}\right),$$

Full text

I Introduction

Many physically relevant solutions to the Einstein equations are asymptotically flat (or asymptote to flat/open Friedman-Lemaitre-Robertson-Walker cosmologies), and so are infinite in spatial extent. To numerically generate these solutions on finite computational domains, several methods are presently in use. One approach is to use compactification on the spatial slices so that the boundary of the computational domain is spatial infinity Garfinkle and Duncan (2001); Pretorius (2005) or (future) null infinity Husa (2002); Frauendiener (2004); Winicour (2012) (for asymptotically flat spacetimes, there is additionally a proposal to extend the computational domain past future null infinity by employing an “artificial” cosmological constant; see e.g. van Meter et al. (2006)). A common approach is to evolve a finite subregion of each spatial slice. One then solves an initial boundary value problem for the Einstein equations; well posed formulations include Friedrich and Nagy (1999); Babiuc et al. (2006). One can also use approximate boundary conditions and hope that they do not spoil the constraints or the well posedness of the system of equations. Further discussion of these procedures may be found in Sarbach and Tiglio (2012) and references therein.

Recently, the authors in Bieri et al. (2019) introduced a simple way to avoid the mathematical complications involved in finding a well posed, constraint preserving initial boundary value problem for the Einstein equations. Instead of evolving the computational boundary along a timelike hypersurface, they propose expanding the computational domain along a spacelike hypersurface at each time step. As all characteristics are ingoing on the expansion surface, this allows for the imposition of initial data instead of boundary conditions along the boundary of the computational domain.

Here we discuss another simple method to avoid the use of imposing boundary conditions or spatial compactification on a finite computational domain. For each time step in the simulation, excise inwards along the innermost characteristic (or on a surface spacelike with respect to the innermost characteristic) so that the computational domain remains within the domain of influence of the initial data. As no characteristics are ingoing on this surface, there is no need to impose boundary conditions or initial data on the boundary (see e.g. Kreiss and Lorenz (1989) and references therein). We refer to the proposal discussed in Bieri et al. (2019) as an “expansion method”, and the idea discussed here as an “excision method”. Figs. (1a)-(1b) provide an illustration of the two ideas.

Excision methods have long been used in numerical relativity Seidel and Suen (1992) for excising the interiors of trapped regions (although excision outside of trapped regions have been applied in, e.g. Boyle et al. (2007), and Pretorius and Choptuik (2000) discusses a “singularity excision” method that could in principle work outside a trapped surface). In this work we discuss numerical and coordinate conditions to implement excision in computational boundaries exterior to trapped regions, and to discuss how the excision and expansion methods relate to one another.

We follow the sign conventions of Misner, Thorne, and Wheeler Misner et al. (1973), and set $c=1$, $8\pi G=1$.

II The excision method in more detail

II.1 CFL condition and

Implementation of excision method with finite difference methods with $3+1$ evolution

Consider an outer grid boundary, and denote the innermost characteristic speed orthogonal to the grid boundary by $c_{-}$. For each time step $\Delta t$, we must excise by a value $|\delta x|\geq|c_{-}\Delta t|$ so the domain of dependence of the grid at $t+\Delta t$ is a subset of the grid at time $t$. For explicit finite difference methods, the Courant-Friedrichs-Lewy (CFL) condition states that the numerical domain of dependence of the solution method must contain by the mathematical domain of dependence of the underlying partial differential equation. This condition sets $c\Delta t\leq\lambda_{max}\Delta x$, where $\lambda_{max}$ is the maximum CFL number and $c$ is the characteristic speed Courant et al. (1967). We see that when using excision on the outer boundary with a CFL number $\lambda_{max}\leq 1$ (which is typically the case for explicit time solving finite difference methods), we can only integrate for a time at most equal to $t\leq N_{x}\Delta t\leq\lambda_{max}T$, where $N_{x}$ is the number of initial spatial grid points and $T=(N_{x}\Delta x)/c$ is the light crossing time of the initial data, as illustrated in Fig. (2a). This condition may be relaxed if one is willing to incur a small amount of violation of the CFL condition on the outer boundary. In particular, to excise directly along the ingoing characteristic one should excise one spatial grid point every $1/(c_{-}\lambda)$ time steps, as is illustrated in Fig. (2b). Provided the characteristic speeds of any potential errors incurred at the boundary by violating the CFL condition are bounded by $c_{-}$, the error incurred by this approach will be contained in a region near the excision boundary, and the size of that region will converge to zero as the resolution increases. See Pretorius and Choptuik (2000) for an example of a stable and convergent code that excises along an ingoing null ray, and also the discussion in Sec. (III.4). Along the excision surface one may use, e.g. upwind difference stencils as is done in the example code described in this paper; see Sec. (III.2).

By contrast, with the expansion method one may evolve in principle for an indefinite amount of time by specifying a larger and larger spacelike expansion region Bieri et al. (2019). Conversely, the computational resources to evolve to another time step increases with the growth in the computational domain in the expansion method, while with the excision method they decrease as the domain shrinks in size.

II.2 Null coordinates and excision condition

As discussed earlier, if the computational domain boundary was tangent to the innermost characteristic, there would be no ingoing characteristics into the domain. There would neither be a need to apply boundary conditions on that boundary nor to excise along that boundary. Here we describe coordinate conditions that automatically enforce that setup. We assume that all the characteristics are contained within the null cone, which is the case for many physical fields (e.g. Geroch (1996) and references therein). Gravitational and electromagnetic waves travel along null characteristics, so we look for coordinates $\{x^{\mu}\}$ such that, e.g. a $x^{1}=const.$ hypersurface is tangent to the ingoing null ray at the computational boundary. As is well-known (e.g. Hawking and Ellis (1975)), at each point of a Lorentzian manifold there are two linearly independent real null vectors, thus we may choose at most two null coordinate surfaces.

To organize this in a more concrete form, we write the metric as (e.g. Ripley and Yagi (2018))

[TABLE]

where lower-case Latin indices range over $\{0,1\}$, upper case Latin indices range over $\{2,3\}$, $\alpha_{ab}$ has Lorentzian signature, and $\gamma_{AB}$ is positive definite. The matrices $\alpha_{ab}/\alpha^{ab}$ lower/raise in indices $\{a\}$, while $\gamma_{AB}/\gamma^{AB}$ lower/raise indices $\{A\}$. We then have

[TABLE]

We consider the case of one null coordinate, $x^{1}\equiv w$. We then have $g^{\mu\nu}\partial_{\mu}w\partial_{\nu}w=0\implies\alpha^{ww}=0$. Labeling $p\equiv x^{0}$, the metric then reads

[TABLE]

There remain three gauge degrees of freedom. We can view $p$ as a time coordinate if $g^{\mu\nu}\partial_{\mu}p\partial_{\nu}p<0$ (i.e. $a_{ww}>0$). We could then evolve with $p$ with a $w=const.$ surface acting as a boundary where no boundary conditions needs to be imposed, as is illustrated in Fig. (3a). If $a_{ww}<0$, then $p$ is a spatial coordinate. We would then need to treat $w$ as the timelike variable. As an aside, in that case one can exhaust the remaining gauge freedom to set $\beta^{A}_{p}=0$ and $\mathrm{det}\gamma_{AB}=p^{4}\mathfrak{q}$, where $\mathfrak{q}$ is the determinant of the unit 2-sphere, to obtain Bondi-Sachs coordinates Bondi (1960); Bondi et al. (1962); Sachs (1962). For a discussion of a numerical implementation of gravitational collapse in spherical symmetry with Bondi-Sachs coordinates see e.g. Pretorius and Lehner (2004).

We next consider two null coordinates. We label $u\equiv x^{0}$, $v\equiv x^{1}$. The conditions $g^{\mu\nu}\partial_{\mu}u\partial_{\nu}u=0$ and $g^{\mu\nu}\partial_{\mu}v\partial_{\nu}v=0$ give $\alpha^{uu}=0$ and $\alpha^{vv}=0$, respectively. The metric is

[TABLE]

There remain three gauge degrees of freedom when we include the simultaneous null rescaling degree of freedom $u\to e^{\kappa}u$, $v\to e^{-\kappa}v$. As is well known (e.g. Winicour (2012) and references therein), with characteristic initial data we do not need to impose boundary conditions on $u=const.$ and $v=const.$ boundaries, as illustrated in Fig. (3b). For a discussion of a numerical implementation of gravitational collapse in double null coordinates, see e.g. Garfinkle (1995) and references therein.

III Example evolution implementing excision method with a finite

difference code: self gravitating scalar in spherical symmetry

As an example implementation of an excision method, we consider a massless scalar field coupled to Einstein gravity.

III.1 Equations of motion

The massless self gravitating scalar field equations are

[TABLE]

We evolve this system using the following coordinate system

[TABLE]

These are called Painlevé-Gullstrand (PG) coordinates as they reduce to those named coordinates for Schwarzschild black hole solutions; earlier analytic and numerical studies with this type of coordinates in four dimensional spacetime include Adler et al. (2005); Ziprick and Kunstatter (2009); Kanai et al. (2011).

Writing

[TABLE]

the evolution equation for $\phi$, Eqn. (5b) can be written as the following system of equations

[TABLE]

The Hamiltonian and momentum constraints give ordinary differential equations (ODEs) for the metric fields

[TABLE]

where (here $n_{\mu}\equiv(-\alpha,0,0,0)$)

[TABLE]

From the PG coordinate solution for the Schwarzschild black hole,

[TABLE]

we see that these coordinates are horizon penetrating, but not singularity avoiding.

As PG coordinates are spatially flat the ADM mass is always zero, we instead use the Misner-Sharp mass Misner and Sharp (1964) to characterize the mass of our solutions 111From the Hamiltonian and momentum constraints Eqs. (9a) and (9a), we see that in vacuum ($\rho=j_{r}=0$) generically the lapse $\alpha=const.$ and $\zeta$ falls off as $\zeta\propto r^{-1/2}$, which violates the asymptotically flat condition used in deriving the ADM mass Arnowitt et al. (2008); as is well known there is no contradiction with a zero ADM mass and nonzero Misner-Sharp mass.

[TABLE]

The radial speeds of propagation of the scalar field $\phi$ are given by $c_{\pm}=\mp\xi_{t}/\xi_{r}$, where $\xi_{\mu}=(\xi_{t},\xi_{r},0,0)$ solve the characteristic equation of the scalar field equation of motion (Eq. (5b)): $g^{\mu\nu}\xi_{\mu}\xi_{\nu}=0$. We see the null characteristics define the domain of influence for $\phi$. Solving for $c_{\pm}$ gives us

[TABLE]

The condition $\zeta(t,r)=1$ signals the formation of a marginally outer trapped surface.

III.2 Description of code

We work with a unigrid code; the initial computational domain covers $r\in[0,R_{max}]$. We set initial data at $t=0$ by specifying the values of $P$ and $Q$, and then solve for $\alpha$ and $\zeta$ using momentum and Hamiltonian constraints, respectively. The ODEs for these constraints are discretized using the trapezoid rule and solved with a Newton relaxation method. At every new time step we solve for $\alpha$, $\zeta$, $P$, and $Q$ by alternating between an iterative Crank-Nicolson solver for $P$ and $Q$ and the ODE solvers for $\alpha$ and $\zeta$ until the discrete infinity norm of all the residuals are below a pre-defined tolerance. Iterative Crank-Nicolson being a two level scheme, initial data only need be specified on the $t=0$ grid. Regularity at the origin sets $Q|_{r=0}=\zeta|_{r=0}=\partial_{r}\alpha|_{r=0}=\partial_{r}P|_{r=0}=0$. We integrate $\zeta$ outwards from $r=0$ using the regularity condition $\zeta=0$. Solving $\alpha$ using Eq. (9b) from $r=0$ automatically enforces its regularity condition. We use $Q=0$ and $\partial_{r}P=0$ in lieu of Eqs. (8a) and (8b) respectively at the grid point $r=0$. There is a residual gauge symmetry $\alpha(t,r)\to c(t)\times\alpha(t,r)$, which we use to rescale $\alpha$ so that $\alpha=1$ at the outermost grid point.

Formation of a marginally outer trapped surface is signaled by $\zeta=1$; see Eq. (13). If $\zeta>1$ for some connected buffer region interior to the trapped surface, we excise all grid points interior to that buffer region. At this excision surface, which we call the inner excision surface to distinguish it from the excision we apply at the outer computational boundary, we evolve $\zeta$ using the $tr$ component of the Einstein equations with upwind stencils

[TABLE]

This provides the boundary condition for $\zeta$ at the excision surface. We then integrate outwards in $r$ using the Hamiltonian constraint as described above to solve for $\zeta$. The lapse $\alpha$ is held fixed at the excision surface and integrated outwards using the momentum constraint; the value of $\alpha$ at the excision surface is arbitrary due to the $\alpha(t,r)\to c(t)\alpha(t,r)$ residual gauge symmetry. We evolve the $Q$ and $P$ fields at the inner excision surface using Eqs. (8a) and (8b) respectively with upwind stencils.

Following the discussion in Sec. (II.1), as PG coordinates are not adapted to the characteristics of the scalar field, we must excise one grid point at the exterior boundary for each time step we take. We refer to this boundary as the outer excision boundary. As the computational domain decreases by one grid point every time step, we can evolve at most for a time $\lambda T$, where $\lambda\leq 1$ is the CFL number and $T$ is the light crossing time of the initial time slice.

III.3 Results and convergence

We consider initial data for $\phi$ of the following form

[TABLE]

where $\{a_{0},w_{0},r_{0}\}$ are constant. We set $Q|_{t=0}=\partial_{r}\phi|_{t=0}$, and $P|_{t=0}=Q|_{t=0}$, which gives approximately ingoing scalar field pulses.

Evolution of an initial scalar pulse that does not form a black hole formation are shown in Figs. (4), (5), and (6). For these runs, we use an initial grid size of $R_{max}=100$, $N_{r}=N_{t}=2^{15}+1$, and CFL number $\lambda=0.5$. The initial conditions are $a_{0}=1\times 10^{-3}$, $w_{0}=5$, $r_{0}=10$; the Misner-Sharp mass at the outer grid point on the initial slice is $m\approx 3.1\times 10^{-2}$. With this setup we can evolve the simulation for $t\approx 1.6\times 10^{3}m$ before the grid shrinks to zero size.

We show results from evolution of an initial scalar pulse that does result black hole in Figs. (7) and (8). For these runs, we use an initial grid size of $R_{max}=800$, $N_{r}=N_{t}=2^{13}+1$, and CFL number $\lambda=0.5$. The initial conditions are $a_{0}=1\times 10^{-2}$, $w_{0}=4$, $r_{0}=8$; the Misner-Sharp mass at the outer grid point on the initial slice is $m\approx 3.6$. With this setup we can evolve the simulation for $t\approx 110m$ before the grid shrinks to zero size.

III.4 Numerical investigation of relaxation of excision CFL condition

In Sec. (II.1) we showed that in order to obey the CFL condition with the excision method with a finite difference code with a CFL number $\lambda\leq 1$, one could only evolve for a time $\lambda T$, where $T$ is the light crossing time of the initial data surface. Here we compare our earlier excision results with an excision method where we excise every $1/(c_{-}\lambda)$ time steps (i.e. directly along the ingoing null characteristic) and examine the form of the solution on the excision surface. We refer to the excision method where we excise to maintain the CFL condition on the boundary as excision method $I$ (illustrated in Fig. (2a)), while the method where we excise along the null ray we call excision method $II$ (illustrated in Fig. (2b)). Note the methods only differ if the CFL number is less than one. For this investigation we rescaled $\alpha$ every time step so that $c_{-}=-1$ (instead of rescaling $\alpha=1$) at the outer excision boundary. As in Sec. (III.2) we use a CFL $\lambda=0.5$, so to implement method $I$ we excised one grid point at every time step, while to implement method $II$ we excised a grid point every other time step. For initial data, we chose an outgoing Gaussian pulse of scalar field (Eq. (15)) with $a=0.003$, $w=4$, and $r_{0}=8$, so the Misner-Sharp mass on the initial outer boundary is $m\sim 0.19$.

We find that methods $I$ and $II$ both produce convergent and stable evolution in the region interior to the excision surface. The methods though produce different results near the outer excision boundary. We find that errors begin to accumulate near the excision surface when using method $II$, which very slowly spreads (e.g. only by a few grid points for resolution $\Delta r\approx 0.1$ over evolution of time $t\sim 50m$) to the interior solution with respect to the moving excision boundary. The width of this error-filled region converges to zero with increasing resolution for any fixed time. This region is not present for runs using method $I$. Figs. (9) provide examples of the behavior of the $P$ variable near the excision boundary for each method. To capture approximately the same point in time to compare the two methods, for these plots we began the run using method $I$ with a domain that extended to radius $R=200(\approx 1000m)$, while with method $II$ the initial domain extended to $R=100(\approx 500m)$. We have also ran simulations where we began with the same $R$ for the initial data and have found a similar accumulation of errors near the excision boundary for method $II$, and the lack of accumulation of error near the excision boundary for method $I$. The difference in $P$ in the solution regions interior to the excision boundary in each plot is due to the fact that we normalize $\alpha$ such that $c_{-}=-1$ on the excision boundaries, and the excision boundaries are at different places at any give time as the CFL obeying excision moves at half the speed as the null excision boundary. This difference in $P$ disappears (to within truncation error) if we do not rescale $\alpha$ (note then $c_{-}\neq-1$ generically on the boundary, so we do not excise precisely on the null ray with excision method $II$ in that case).

We caution that these results may not capture the range of differences that may occur when using methods $I$ and $II$ for different coordinate systems or in axisymmetric/full $3+1$ evolution codes. In particular, there may exist gauges/numerical setups where the errors propagate with characteristic speeds smaller or larger than light speed. As discussed earlier, if the CFL number $\lambda\geq 1$, then there is no difference between methods $I$ and $II$. It is outside the scope of this note (and the outside the capabilities of the $1+1$ evolution code used here) to investigate whether or not the errors incurred by a CFL violating excision method remain near the outer boundary and converge to zero for the gauge conditions more commonly used in symmetry unconstrained $1+3$ dimensional numerical relativity.

IV Discussion

Considerable mathematical and computational challenges accompany developing and implementing well posed, constraint preserving initial boundary value problems for the Einstein equations (e.g. Sarbach and Tiglio (2012) and references therein). We discussed an excision method that removes the need for boundary conditions when evolving the Einstein equations on compact spatial domains, and compared the method to the recently proposed expansion method of Bieri et al. (2019). If one uses coordinates not adapted to the characteristics of the hyperbolic degrees of freedom with an explicit numerical time integrator, using excision that obeys the CFL condition restricts ones evolution to to be comparable to $\lambda T$, where $T$ is the light crossing time of the initial data surface and $\lambda\leq 1$ is the CFL number. If one uses an scheme that is stable with $\lambda=1$, then it is possible to excise directly along the ingoing null characteristic while still obeying the CFL condition on the boundary.

As the computational domain becomes smaller as one evolves in time, the excision method is not useful for simulations where one needs to, for example extract gravitational wave information near future null infinity. The method may be useful though in simulations where one is more interested in the nearby/local physics during gravitational collapse. For example, in critical gravitational collapse one is interested in understanding the nature of the (discrete/continuous) self-similarity of the collapse, the time and spatial scales of which are decreasing exponentially in time Choptuik (1993). Thus while with the excision method one can only evolve for a time comparable to the light-crossing time of the initial data surface, this may not be a serious impediment if one uses an initial data surface that is sufficiently large and with initial data sufficiently close to the critical collapse regime.

Another potential application of the method would be in studies of the interior structure of black hole spacetimes. There has been renewed interest in understanding the stability of the Cauchy horizon of rotating and charged black holes to small perturbations caused either by the infall of matter or gravitational waves Franzen (2016); Luk and Oh (2017); Dafermos and Luk (2017). The stability of the Cauchy horizon of these spacetimes is related to the cosmic censorship conjecture Penrose (1979); Christodoulou (2008), as one consequence of that the conjecture (if true) is that extendability across the Cauchy horizon of the Kerr and Reissner-Nordstrom black holes is not generically possible. An extensive amount of analysis suggests (but does not prove) that small perturbations inside the black hole will seed curvature blowup on or near the Cauchy horizon of those black holes (e.g. McNamara and Penrose (1978); Poisson and Israel (1990)), and that this curvature blowup would make the spacetime inextendible past the Cauchy horizon. Interestingly, recent work suggests that regardless of a possible blowup in curvature on the Cauchy horizon, the metric of Kerr black holes may be generically extendible by a continuous ($C^{0}$) metric that solves the vacuum Einstein equations Dafermos and Luk (2017) (provided the black hole solution exterior to the event horizon is stable). Numerical studies of black hole interiors may provide further insight on the dynamics of the Cauchy horizon of Kerr and Reissner-Nordstrom black hole interiors (for a recent such study see Chesler et al. (2019)). The excision method could be useful for the following setup: begin with an initial data surfaces that is interior to a slightly perturbed Reissner-Nordstrom/Kerr black hole, and have one end of the initial data surface terminate “close” to the putative Cauchy horizon. Then, evolve and excise along the outer surface-if curvature blowup begins to occur during the course of evolution on the outer boundary, this would suggest a spacelike curvature singularity is forming starting from the Cauchy horizon. This approach closely mirrors the reasoning used in previous numerical studies of the interior of black holes using double null coordinates, see e.g. Brady (1999); Nakonieczna et al. (2019) and references therein. Excising along the outer boundary allows one to approach this problem using $1+3$ time evolution, although the CFL condition (Sec. (II.1)) would restrict one to excising on a spacelike, instead of an exactly null, surface if the CFL number is less than one. Using for example implicit time stepping routines, it may be possible to devise stable, CFL preserving, and convergent $3+1$ codes that implement excision directly along the null ray. For CFL numbers less than one, if one is willing to violate the CFL condition along the excision surface one can also excise directly along the null ray by excising only one point every $1/(c_{-}\lambda)$ time steps (see Sec. (III.4) or Pretorius and Choptuik (2000)). Our numerical investigations suggest the error this method incurs on the excision surface converges to zero with higher resolution. It would be interesting to investigate this issue further in a non-spherically symmetric evolution code, with a gauge condition (such as generalized harmonic or BSSN) more commonly used in the numerical relativity literature.

Finally, we note that the expansion and excision methods could be profitably combined: first the expansion method would be used, then the excision method would be employed once the grid reached some maximum size. This would allow for the computational domain to not grow to an impractically large size, and for longer time evolution than the pure excision scheme would allow for a given fixed initial grid size.

Acknowledgements.

I am grateful to David Garfinkle and Frans Pretorius for helpful discussions regarding earlier drafts of this work, Alex Pandya for a helpful discussion on an earlier application of an excision method in Pretorius and Choptuik (2000), Ted Jacobson for alerting me to van Meter et al. (2006), and to the anonymous referees for their useful comments and references.