Black hole spacetimes in Gkeyll

In principle, the Gkeyll Moment app provides the requisite algorithmic infrastructure for solving the equations of general relativistic hydrodynamics within an arbitrary curved spacetime background. However, Gkeyll provides specific in-built support for stationary black hole spacetimes in the Cartesian Kerr-Schild coordinate system \(\left\lbrace t, x, y, z \right\rbrace\), since such spacetimes and coordinate systems prove particularly useful when performing hydrodynamic simulations of black hole accretion and astrophysical jet formation. In the absence of electromagnetic fields, the no-hair theorem implies that any stationary black hole in general relativity may be characterized by a mass \(M\) and a non-dimensional spin parameter:

\[a = \frac{J}{M} \in \left[ 0, 1 \right],\]

where \(J\) is a dimensional angular momentum parameter. Gkeyll uses the formulation of the Kerr-Schild coordinate system presented in [Moreno2002], based upon the earlier work of [Debney1969]. In this short technical note, we will follow the notational and terminological conventions introduced within the description of the GRHD equations in Gkeyll. In particular, we assume a \({3 + 1}\) “ADM” decomposition of our 4-dimensional spacetime into 3-dimensional spacelike hypersurfaces with induced metric tensor \(\gamma_{i j}\), lapse function \(\alpha\), shift vector \(\beta^i\), and extrinsic curvature tensor \(K_{i j}\). Einstein summation convention (in which repeated tensor indices are implicitly summed over) is assumed throughout. The Greek indices \(\mu, \nu\) range over the full spacetime coordinate directions \(\left\lbrace 0, \dots 3 \right\rbrace\), while the Latin indices \(i, j, k\) range over the spatial coordinate directions \(\left\lbrace 1, \dots 3 \right\rbrace\) only.

See also this note on the complete eigensystem for the GRHD equations (as implemented as part of Gkeyll’s stable time-step calculation), and this note on Gkeyll’s “robustified” conservative to primitive variable reconstruction algorithm for both special and general relativity.

Cartesian Kerr-Schild coordinates

By default, Gkeyll represents all black hole spacetimes as generalized Kerr-Schild perturbations of flat (Minkowski) spacetime, of the form:

\[g_{\mu \nu} = \eta_{\mu \nu} - V l_{\mu} l_{\nu},\]

where \(\eta_{\mu \nu}\) denotes the standard Minkowski metric:

\[\begin{split}\eta_{\mu \nu} = \begin{bmatrix} -1 & 0 & 0 & 0\\ 0 & 1 & 0 & 0\\ 0 & 0 & 1 & 0\\ 0 & 0 & 0 & 1 \end{bmatrix},\end{split}\]

\(V\) is a scalar (referred to as the Kerr-Schild scalar within the Gkeyll code), and \(l_{\mu} = l^{\mu}\) is a null (co)vector with respect to the background spacetime (referred to as the Kerr-Schild vector within the Gkeyll code). For a black hole with mass \(M\) and dimensionless spin \(a \in \left[ 0, 1 \right]\) in Cartesian coordinates \(\left\lbrace t, x, y, z \right\rbrace\), it is helpful to introduce a generalized radial quantity \(R\), defined implicitly by the following algebraic equation:

\[\frac{x^2 + y^2}{R^2 + M^2 a^2} + \frac{z^2}{R^2} = 1,\]

with the explicit (positive) solution:

\[R = \left( \frac{1}{\sqrt{2}} \right) \sqrt{x^2 + y^2 + z^2 - M^2 a^2 + \sqrt{\left( x^2 + y^2 + z^2 - M^2 a^2 \right)^2 + 4 M^2 a^2 z^2}},\]

with respect to which the Kerr-Schild scalar \(V\) and the components of the null Kerr-Schild (co)vector \(l_{\mu} = l^{\mu}\) may be straightforwardly expressed as:

\[V = - \frac{M R^3}{R^4 + M^2 a^2 z^2},\]
\[l_0 = l^0 = -1,\]
\[l_1 = l^1 = - \frac{R x + a M y}{R^2 + M^2 a^2},\]
\[l_2 = l^2 = - \frac{R y - a M x}{R^2 + M^2 a^2},\]

and:

\[l_3 = l^3 = - \frac{z}{R},\]

respectively. From here, we are able to perform a \({3 + 1}\)/ADM decomposition of the spacetime into spacelike hypersurfaces, which are themselves generalized Kerr-Schild perturbations of flat (Euclidean) space, of the form:

\[\gamma_{i j} = \delta_{i j} - 2 V l_i l_j,\]

where \(\delta_{i j}\) denotes the standard Euclidean metric:

\[\begin{split}\delta_{i j} = \begin{bmatrix} 1 & 0 & 0\\ 0 & 1 & 0\\ 0 & 0 & 1 \end{bmatrix}.\end{split}\]

Within such a foliation of a Kerr-Schild spacetime, the ADM lapse function \(\alpha\) and the components of the ADM shift vector \(\beta^i\) may be written in terms of the Kerr-Schild scalar \(V\) and (co)vector components \(l_i = l^i\) as:

\[\alpha = \frac{1}{\sqrt{1 - 2V}},\]

and:

\[\beta^i = \frac{2V}{1 - 2V} l^i,\]

respectively. The expression for the extrinsic curvature tensor \(K_{i j}\) is somewhat more involved:

\[\begin{split}K_{i j} = \alpha \left( - l_i \frac{\partial}{\partial x^j} \left( V \right) - l_j \frac{\partial}{\partial x^i} \left( V \right) - V \frac{\partial}{\partial x^j} \left( l_i \right) - V \frac{\partial}{\partial x^i} \left( l_j \right) \right.\\ \left. + 2 V^2 \left( l_i l^k \frac{\partial}{\partial x^k} \left( l_j \right) + l_j l^k \frac{\partial}{\partial x^k} \left( l_i \right) \right) + 2 V l_i l_j l^k \frac{\partial}{\partial x^k} \left( V \right) \right).\end{split}\]

Transformation to standard Boyer-Lindquist coordinates

A more common representation of the Kerr metric for a spinning black hole is in the Boyer-Lindquist spherical (or, more precisely, oblate spheroidal) coordinate system \(\left\lbrace t, r, \theta, \phi \right\rbrace\) proposed by [Boyer1967], which can be straightforwardly related to the Cartesian coordinate system \(\left\lbrace t, x, y, z \right\rbrace\) in the usual way:

\[x = \sqrt{r^2 + a^2} \sin \left( \theta \right) \cos \left( \phi \right),\]
\[y = \sqrt{r^2 + a^2} \sin \left( \theta \right) \sin \left( \phi \right),\]
\[z = r \cos \left( \theta \right).\]

The \({3 + 1}\)/ADM decomposition of the Kerr metric in Boyer-Lindquist coordinates then yields a family of spacelike hypersurfaces with a purely diagonal spatial metric tensor \(\gamma_{i j}\), of the form:

\[\begin{split}\gamma_{i j} = \begin{bmatrix} \frac{R^2}{r^2 + a^2 - 2 M r} & 0 & 0\\ 0 & R^2 & 0\\ 0 & 0 & \frac{R^2 \left( r^2 + a^2 \right) + 2 M^3 a^2 r \sin^2 \left( \theta \right)}{R^2} \sin^2 \left( \theta \right), \end{bmatrix}\end{split}\]

with the ADM lapse function \(\alpha\) and shift vector components \(\beta^i\) now given by:

\[\alpha = \sqrt{\frac{R^2 \left( r^2 + a^2 - 2 M r \right)}{R^2 \left( r^2 + a^2 \right) + 2 M^3 a^2 r \sin^2 \left( \theta \right)}},\]

and:

\[\begin{split}\boldsymbol\beta = \begin{bmatrix} 0\\ 0\\ - \frac{2 a M r}{R^2 \left( r^2 + a^2 \right) + 2 M^3 a^2 r \sin^2 \left( \theta \right)} \end{bmatrix},\end{split}\]

respectively. Although expressing the Kerr metric in Boyer-Lindquist coordinates has certain aesthetic advantages (e.g. the diagonality of the spatial metric \(\gamma_{i j}\), the non-zero angular component \(\beta^{\phi} \neq 0\) of the shift vector explicitly indicating frame-dragging effects, etc.), it is less appropriate for numerical simulations than the Kerr-Schild coordinate system due to the appearance of coordinate singularities wherever:

\[r^2 + a^2 - 2 M r = 0,\]

whose solutions \(r = r_{-}\) and \(r = r_{+}\) correspond to the interior and exterior black hole horizons:

\[r_{\pm} = M \left( 1 \pm \sqrt{1 - a^2} \right),\]

respectively. By contrast, the Kerr-Schild coordinate system is horizon-adapted, meaning that the solution may be smoothly extended across both the interior and exterior horizons, guaranteeing greater numerical stability in the near-horizon region; it is for this reason that Gkeyll represents black hole geometries in the Kerr-Schild coordinate system by default. In order to see how to transform from the Cartesian Kerr-Schild coordinate system into the spherical (oblate spheroidal) Boyer-Lindquist coordinate system, we begin by transforming to a spherical (oblate spheroidal) form of the Kerr-Schild coordinate system \(\left\lbrace t, r, \theta, \phi \right\rbrace\), related to the old Cartesian form \(\left\lbrace t, x, y, z \right\rbrace\) by the transformation:

\[x = \left( r \cos \left( \phi \right) - a \sin \left( \phi \right) \right) \sin \left( \theta \right) = \sqrt{r^2 + a^2} \sin \left( \theta \right) \cos \left( \phi - \arctan \left( \frac{a}{r} \right) \right),\]
\[y = \left( r \sin \left( \phi \right) + a \cos \left( \phi \right) \right) \sin \left( \theta \right) = \sqrt{r^2 + a^2} \sin \left( \theta \right) \sin \left( \phi - \arctan \left( \frac{a}{r} \right) \right),\]
\[z = r \cos \left( \theta \right).\]

Note that this is almost the same as the Boyer-Lindquist coordinate system, but with a small modification in the \(\phi\) coordinate that causes the new system to remain regular across the black hole horizon(s). The \({3 + 1}\)/ADM decomposition of the Kerr metric in spherical (oblate spheroidal) Kerr-Schild coordinates now yields a family of spacelike hypersurfaces with a single off-diagonal term in the spatial metric tensor \(\gamma_{i j}\):

\[\begin{split}\gamma_{i j} = \begin{bmatrix} 1 + \frac{2 M r}{R^2} & 0 & - M a \sin^2 \left( \theta \right) \left( 1 + \frac{2 M r}{R^2} \right)\\ 0 & R^2 & 0\\ - M a \sin^2 \left( \theta \right) \left( 1 + \frac{2 M r}{R^2} \right) & 0 & \left( r^2 + M^2 a^2 + \frac{2 M^3 a^2 r}{R^2} \sin^2 \left( \theta \right) \right) \sin^2 \left( \theta \right) \end{bmatrix},\end{split}\]

coupling the radial and angular directions, with the ADM lapse function \(\alpha\) and shift vector components \(\beta^i\) now given by:

\[\alpha = \frac{1}{\sqrt{1 + \frac{2 M r}{R^2}}},\]

and:

\[\begin{split}\boldsymbol\beta = \begin{bmatrix} \frac{2 M r}{R^2 \left( 1 + \frac{2 M r}{R^2} \right)}\\ 0\\ 0 \end{bmatrix},\end{split}\]

respectively. The spatial velocity components \(v_{KS}^{r}\), \(v_{KS}^{\theta}\), and \(v_{KS}^{\phi}\) within this spherical/oblate spheroidal Kerr-Schild coordinate system may then be transformed to the spherical/oblate spheroidal Boyer-Lindquist coordinate system, yielding \(v_{BL}^{r}\), \(v_{BL}^{\theta}\), and \(v_{BL}^{\phi}\), respectively, by means of the following prescription due to [Font1999]:

\[v_{BL}^{r} = \frac{1}{\Psi} \left( v_{KS}^{r} - \frac{\beta_{KS}^{r}}{\alpha_{KS}} \right),\]
\[v_{BL}^{\theta} = \frac{1}{\Psi} v_{KS}^{\theta},\]

and:

\[\begin{split}v_{BL}^{\phi} = \frac{1}{\Psi} \left( v_{KS}^{\phi} - \left( \frac{M a}{r^2 + a^2 - 2 M r} \right) v_{KS}^{r} \right)\\ - \frac{1}{\Psi} \left( \frac{\beta_{KS}^{\phi}}{\alpha_{KS}} - \left( \frac{M a}{r^2 + a^2 - 2 M r} \right) \left( \frac{\beta_{KS}^{r}}{\alpha_{KS}} \right) \right) + \frac{\beta_{BL}^{\phi}}{\alpha_{BL}},\end{split}\]

where we have introduced the scalar quantity \(\Psi\), defined by:

\[\Psi = \frac{\alpha_{BL}}{\alpha_{KS}} - \left( \frac{2 M r}{r^2 + a^2 - 2 M r} \right) \alpha_{BL} \left( v_{KS}^{r} - \frac{\beta_{KS}}{\alpha_{KS}} \right),\]

and where \(\alpha_{KS}\)/\(\beta_{KS}^{i}\) and \(\alpha_{BL}\)/ \(\beta_{BL}^{i}\) denote the lapse function and shift vector components in the spherical/oblate spheroidal Kerr-Schild coordinate system and the spherical/oblate spheroidal Boyer-Lindquist coordinate system, respectively.

Excision boundary conditions

Although no physical information may propagate from the interior region of a black hole’s event horizon out to the exterior region, the type of finite-volume numerical methods employed within Gkeyll’s moment app may make use of cells from the black hole’s interior region when updating near-horizon cells in the exterior region. Thus, unless one is particularly careful about how this finite-volume update is performed, this can lead to unphysical wave propagation and numerical instabilities in the near-horizon region. In order to avoid this, Gkeyll imposes a variant of the excision boundary conditions of [Alcubierre2001], in which the black hole interior region, defined by:

\[r < r_{-} = M \left( 1 - \sqrt{1 - a^2} \right),\]

i.e. the region whose boundary corresponds to the inner event horizon of the Kerr geometry, corresponds to a region of zero flux. The advantage of this approach is that it can be implemented using only very minor modifications to the overall finite-volume update algorithm.

References

[Moreno2002]

C. Moreno, D. Núñez and O. Sarbach, “Kerr-Schild type initial data for black holes with angular momenta”, Classical and Quantum Gravity 19 (23): 6059-6073. 2002.

[Debney1969]

G. C. Debney, R. P. Kerr and A. Schild, “Solutions of the Einstein and Einstein-Maxwell Equations”, Journal of Mathematical Physics 10 (10): 1842-1854. 1969.

[Boyer1967]

R. H. Boyer and R. W. Lindquist, “Maximal Analytic Extension of the Kerr Metric”, Journal of Mathematical Physics 8 (2): 265-281. 1967.

[Font1999]

J. A. Font, J. M. Ibáñez and P. Papadopoulos, “Non-axisymmetric relativistic Bondi-Hoyle accretion on to a Kerr black hole”, Monthly Notices of the Royal Astronomical Society 305 (4): 920-936. 1999.

[Alcubierre2001]

M. Alcubierre and B. Brügmann, “Simple excision of a black hole in 3 + 1 numerical relativity”, Physical Review D 63 (10): 104006. 2001.