Toward General Relativistic Magnetohydrodynamics Simulations in Stationary Nonvacuum Spacetimes

General relativistic magnetohydrodynamic (GRMHD) simulations have emerged as an indispensable tool in modern astrophysical research, providing a robust framework for investigating the complex dynamics of magnetized gas within extreme gravitational environments (e.g., Banyuls et al. 1997; Koide et al. 1999; De Villiers et al. 2003; Gammie et al. 2003). By virtue of their ability to simulate diverse accretion regimes spanning from sub- to super-Eddington rates, GRMHD simulations offer invaluable insights into the structure and evolution of a variety of black hole (BH) and neutron star systems (e.g., Fragile et al. 2007; Narayan et al. 2012; Sadowski & Narayan 2015; Parfrey & Tchekhovskoy 2017; Liska et al. 2018; Porth et al. 2019; Begelman et al. 2022; Chatterjee & Narayan 2022), as well as astrophysical phenomena like gamma-ray bursts (e.g., Gottlieb et al. 2022), tidal disruption events (e.g., Curd & Narayan 2019; Andalman et al. 2022), and high-energy flares from the centers of galaxies (e.g., Chatterjee et al. 2021; Porth et al. 2021; Ripperda et al. 2022). Furthermore, GRMHD simulations self-consistently launch powerful jets that can accelerate to relativistic speeds and extend to galactic scales and thus are crucial for understanding galaxy evolution (e.g., McKinney 2006; Tchekhovskoy et al. 2011; Chatterjee et al. 2019; Narayan et al. 2022; Ricarte et al. 2023). Also, the predictive capabilities of GRMHD simulations facilitate direct comparisons with observational data, enhancing their applicability in interpreting real astrophysical systems (e.g., Mościbrodzka et al. 2009; Dexter et al. 2012; Davelaar et al. 2018; Chael et al. 2019; Chatterjee et al. 2020; Ricarte et al. 2020; Cruz-Osorio et al. 2022; Ressler et al. 2023), and have helped motivate experimental tests of gravitational effects such as frame dragging (Ricarte et al. 2022) and electromagnetic extraction of energy from spinning BHs (Chael et al. 2023).

The Event Horizon Telescope (EHT) collaboration recently demonstrated that it is now possible to "measure" the spacetime metric of astrophysical BHs. More concretely, we are able to measure the level of agreement of the spacetime geometry of the supermassive BHs M87* (The EHT Collaboration et al. 2019a, 2019b; Psaltis et al. 2020) and Sgr A* (The EHT Collaboration et al. 2022a, 2022b) with that of a Kerr BH (Kerr 1963), which is a vacuum, stationary solution of general relativity (GR; see, e.g., Wiltshire et al. 2009). Such investigations allow us, in principle, to test fundamental aspects of GR such as the no-hair hypothesis (see, e.g., Carter 1971) and whether astrophysical BHs truly have no surfaces and are devoid of matter (The EHT Collaboration et al. 2022b).

It has also been possible to test the level of agreement of EHT data with non-BH spacetime geometries (such as boson stars, wormholes, classical naked singularities, etc.) within GR, as well as both BH and non-BH spacetimes in non-GR theories (Völkel et al. 2021; Kocherlakota et al. 2021; The EHT Collaboration et al. 2022b; Vagnozzi et al. 2023). This program has seen intense activity, and current approaches use a combination of GRMHD simulations in Kerr BH spacetimes (see, e.g., The EHT Collaboration et al. 2019b, 2022b) and non-GRMHD (xGRMHD) simulations in nonrotating spacetimes (Mizuno et al. 2018; Olivares et al. 2020; Fromm et al. 2021; Röder et al. 2023). In addition, a diverse range of instructive semianalytic accretion flow models have been considered (see, e.g., Broderick et al. 2014; Shaikh et al. 2019; Shaikh & Joshi 2019; Narayan et al. 2019; Paul et al. 2020; Bauer et al. 2022; Özel et al. 2022; Younsi et al. 2023; Kocherlakota & Rezzolla 2022; Ayzenberg 2022; The EHT Collaboration et al. 2022b) to model the physics of accretion flows and their emission.

Nearly all of the xGRMHD work to date has been restricted to spherically symmetric "static" spacetime metrics. This is because, apart from the spinning Kerr metric (Kerr 1963) and Kerr–Newman metric (Newman et al. 1965), the only metrics available in the horizon-penetrating Kerr–Schild coordinates required by current GRMHD codes (see, e.g., Font et al. 1998; McKinney & Gammie 2004; Sadowski et al. 2014; Porth et al. 2017; for a discussion of representative GRMHD codes) are static metrics. The one other exception is Nampalliwar et al. (2022), who used the Kerr–Schild form of the Johannsen (2013a) stationary metric, which is derived from the well-known Johannsen & Psaltis (2011) parameterized BH metric.

Astrophysical objects are expected to possess angular momentum. Hence, their spacetimes will generally be axisymmetric (not spherically symmetric), and their metrics should be "stationary" (not static). For a more realistic confrontation of the underlying effective theory of gravity and fields in the vicinity of ultracompact objects to current observations, and in anticipation of future higher-quality data, it is necessary to be able to model accretion and emission processes in generic stationary spacetimes. However, other than the aforementioned Kerr, Kerr–Newman, and Johannsen metrics, no other stationary spacetime metrics have been written down in the horizon-penetrating coordinates needed for computer simulations. The goal of the present paper is to eliminate this roadblock.

We introduce a general family of stationary and axisymmetric metrics (Azreg-Aïnou 2014a, 2014b) that are derived from a corresponding family of "seed" spherically symmetric metrics. We then cast these stationary metrics in horizon-penetrating Kerr–Schild form. This is the main contribution of the present paper. Using the results described in this paper, we will report elsewhere the first high-resolution, 3D GRMHD simulations in spinning non-Kerr spacetimes (K. Chatterjee et al. 2023a, 2023b, in preparation).

We begin in Section 2 by introducing the general form of a (spherically symmetric) static metric. This form includes the Schwarzschild metric as a special case but covers a wide range of other models as well. We pick out and discuss a few well-known examples of the latter.

In Section 3, we introduce the Azreg–Aïnou (AA; Azreg-Aïnou 2014b) metric, which is a stationary, axisymmetric generalization of the static metric described in Section 2. The AA metric was proposed as an ansatz, inspired by the Newman & Janis (1965) algorithm, to describe generic spinning spacetimes. This metric has several attractive features, which we touch on in Section 1.1 and consider in greater detail in Section 3.

In Section 4, we report the form of the AA metric in horizon-penetrating coordinates, as needed for most GRMHD codes, and in Section 5, we write down the 3+1 decomposition of the metric required for certain other codes. The spherical polar Kerr–Schild coordinates we use to write the horizon-penetrating AA metric are adapted to the ingoing principal null congruence (PNC) of spacetime.

Finally, in Section 6, we summarize our findings and present a representative list of stationary and axisymmetric metrics that are obtained from the AA metric for specific choices of its metric functions. These metrics may be directly used in future GRMHD simulations. We also describe the steps needed to generate other stationary metrics in the future. Throughout, we use geometrized units, G = c = 1, and the metric signature is (−, +, +, +).

Of necessity, this paper includes a fair amount of technical discussion about the AA metric and its various forms. A reader who is interested in utilizing various BH and non-BH metrics to explore with GRMHD simulations can obtain the main results directly by reading only Sections 2 and 6, along with Tables 1 and 2.

Table 1. Metric Functions of Selected Spherically Symmetric and Static Spacetimes

Object	Theory	Spacetime	f(r)	R2(r)
BH	GR	Schwarzschild	$1-\tfrac{2M}{r}$	r2
BH	GR	RN	$1-\tfrac{2M}{r}+\tfrac{{Q}^{2}}{{r}^{2}}$	r2
BH	GR	Modified Hayward	$1-\tfrac{2{{Mr}}^{3}}{{r}^{4}+2{L}^{4}}$	r2
BH	String	GMGHS	$1-\tfrac{2M}{r}$	$r\left(r-\tfrac{{Q}^{2}}{M}\right)$
Naked singularity	GR	JNW	${\left(1-\tfrac{{r}_{\star }}{r}\right)}^{1-\hat{\nu }}$	${r}^{2}{\left(1-\tfrac{{r}_{\star }}{r}\right)}^{\hat{\nu }}$
Naked singularity	GR	JMN-1	$\left(1-\tfrac{2M}{{R}_{{\rm{b}}}}\right){\left(\tfrac{r}{{r}_{{\rm{b}}}}\right)}^{2\sigma }$	${R}_{{\rm{b}}}^{2}{\left(\tfrac{r}{{r}_{{\rm{b}}}}\right)}^{2-2\sigma }$

Note. We have used the coordinate freedom to set $-{\hat{g}}_{{tt}}{\hat{g}}_{{rr}}=g(r)=1$ in the general metric (Equation (1)); hence, we have only two metric functions, f(r) and R(r). The Schwarzschild and RN spacetimes describe the vacuum and electrovacuum BH solutions of GR, respectively. The modified Hayward BH spacetime contains a regular BH and is generated by an anisotropic fluid. The GMGHS BH metric is a solution to the low-energy theory of the heterotic string and contains a scalar field (a dilaton, which vanishes asymptotically), as well as an electromagnetic field. The JNW naked singularity spacetime is generated by a scalar field. The JMN-1 naked singularity spacetime is an interior solution to the Schwarzschild spacetime and is generated by an anisotropic fluid. For easier comparison with the other five models, the JMN-1 metric is presented here with g(r) = 1 (see Appendix A for details), thus differing from the original form given in Equation (29) of Joshi et al. (2011). The parameter Q denotes an electromagnetic charge, L denotes a (de Sitter) length scale, R_b ≥ 2.5M denotes a matching radius (σ = M/(R_b − M) and r_b = (1 − σ)R_b ), and $\hat{\nu }$ denotes a scalar charge (${r}_{\star }=2M/(1-\hat{\nu })$ locates the curvature singularity).

Download table as:  ASCIITypeset image

1.1. Background on the AA Metric

The line element ${\rm{d}}{s}^{2}={\hat{g}}_{\mu \nu }{\rm{d}}{x}^{\mu }{\rm{d}}{x}^{\nu }$ of an arbitrary static and spherically symmetric metric ${\hat{g}}_{\mu \nu }$ can be written in spherical polar coordinates, x^μ = (t, r, ϑ, φ), in the form

$$\begin{eqnarray}&&{\rm{d}}{s}^{2}=-f(r){\rm{d}}{t}^{2}+\displaystyle \frac{g(r)}{f(r)}\,{\rm{d}}{r}^{2}+{R}^{2}(r)\,{\rm{d}}{{\rm{\Omega }}}_{2}^{2},\end{eqnarray} \tag{ 1 }$$

where ${\rm{d}}{{\rm{\Omega }}}_{2}^{2}={\rm{d}}{\vartheta }^{2}+{\sin }^{2}\vartheta \,{\rm{d}}{\varphi }^{2}$ is the standard line element on a unit 2-sphere. This form of the metric involves three functions of the coordinate r: f(r), g(r), and R(r). However, an arbitrary static metric can be defined using only two metric functions; i.e., any one of the metric functions, f, g, and R, can be eliminated by an appropriate change of coordinates (see, e.g., Chapter 14 of Plebanski & Krasinski 2012). For example, it is possible to find a different radial coordinate r in which $-{\hat{g}}_{{tt}}{\hat{g}}_{{rr}}=g(r)=1$ (as we do in Table 1 and discuss further in Appendix A). However, the required coordinate transformation does not always lead to analytically convenient metric functions. We therefore keep an arbitrary g(r) > 0 in Equation (1)).

The function R(r) gives the curvature radius or areal radius as a function of the radial coordinate r. In simple models, R(r) = r (see Table 1). When this is the case, the coordinate radius r of a 2-sphere is its curvature radius, ${\kappa }_{{\rm{G}}}(r)={{\mathscr{R}}}^{(2)}(r)/2\,=\,1/{r}^{2}$, as well as its areal radius, ${ \mathcal A }(r)=4\pi {r}^{2}$. Here we have used κ_G, ${{\mathscr{R}}}^{(2)}$, and ${ \mathcal A }$ to denote the Gaussian curvature, the Ricci scalar, and the area of a 2-sphere, respectively. While it is always possible to find a coordinate transformation to curvature coordinates in which R(r) = r, this again may only be tractable numerically for some models. Finding simple analytic forms for all of the metric functions in such coordinates is not generally possible (see, e.g., Section IV D of Kocherlakota & Rezzolla 2020). Thus, for general convenience, we will leave R(r) free.

When the metric (Equation (1)) describes a BH spacetime, its event horizon is a null, stationary surface. The normal to such a constant-r surface satisfies ${\hat{g}}^{\mu \nu }{\partial }_{\mu }r\cdot {\partial }_{\nu }r=0$, from which we find that the horizon is located at the outermost root of ${\hat{g}}^{{rr}}=0$, i.e., at the outermost root of f(r).

In Table 1, we list the metric functions corresponding to six well-known and representative spherically symmetric static solutions in GR and string theory. For uniformity, we have written all of the metrics in coordinates such that g(r) = 1. The first two solutions listed are the canonical nonrotating BH solutions of GR, namely, the Schwarzschild BH and the charged Reissner–Nordström (RN) BH, which contain, respectively, a space-like and a time-like curvature singularity (see, e.g., Poisson 2004). Several regular BH spacetime models have been proposed to mimic the desired effect of singularity-resolving physics (see, e.g., Bardeen 1968; Ayón-Beato & García 1998; Bronnikov 2001; Dymnikova 2004; Hayward 2006). The modified Hayward BH (Zhou & Modesto 2023a; n = 4) spacetime included in Table 1 is generated by an anisotropic fluid with a radial equation of state, ω_r = p_r /ρ = − 1, where ρ is the total energy density, and p_r is the radial pressure in the fluid rest frame. The fluid is dark energy–like, producing the pressure necessary to avoid a singularity. The charged Gibbons–Maeda–Garfinkle–Horowitz–Strominger (GMGHS) BH (Gibbons & Maeda 1988; Garfinkle et al. 1991) is different from the RN BH of GR because of the presence of a dilaton (a scalar field) that mediates the interaction between electromagnetism and gravity in the low-energy effective action of string theory. In particular, the central singularity of the GMGHS BH remains space-like.

In addition to the above four BH solutions, Table 1 lists two other non-BH models. Over the years, non-BH solutions to various theories have been discussed in the literature, and significant progress has been made to determine ways in which such objects can be distinguished from BHs using observations (see, e.g., Shaikh et al. 2019; Kocherlakota et al. 2021; The EHT Collaboration et al. 2022b). The Janis–Newman–Winicour naked singularity (JNW; Janis et al. 1968) spacetime is a static solution generated by a minimally coupled, massless scalar field in GR. Indeed, this is also a solution to string theory when an electromagnetic field is absent (Virbhadra 1997), and, in addition, it is a solution to the Brans–Dicke theory with the parameter ω = −1 (Kar 1997). Finally, the Joshi–Malafarina–Narayan-1 naked singularity (JMN-1; Joshi et al. 2011) is constructed in GR using an anisotropic fluid with vanishing radial pressure, p_r = 0.

Inspired by the Newman & Janis (1965) algorithm, Azreg-Aïnou (2014b) proposed a stationary and axisymmetric generalization of the spherically symmetric metric (Equation (1)) described in Section 2. We will refer to this stationary metric as the AA metric and denote it by g_{μ ν} (to distinguish it from the ${\hat{g}}_{\mu \nu }$ of the spherically symmetric metric). In Boyer–Lindquist (BL; Boyer & Lindquist 1967) coordinates, x^μ = (t, r, ϑ, φ), ⁶ the AA metric is given by

$$\begin{eqnarray}\begin{array}{rcl}{\rm{d}}{s}^{2} & = & \displaystyle \frac{X}{{\rm{\Sigma }}}\left[-\left(1-\displaystyle \frac{2F}{{\rm{\Sigma }}}\right){\rm{d}}{t}^{2}-2\displaystyle \frac{2F}{{\rm{\Sigma }}}a{\sin }^{2}\vartheta \,{\rm{d}}t{\rm{d}}\varphi \right.\\ & & \left.+\displaystyle \frac{{\rm{\Sigma }}}{{\rm{\Delta }}}{\rm{d}}{r}^{2}+{\rm{\Sigma }}\,{\rm{d}}{\vartheta }^{2}+\displaystyle \frac{{\rm{\Pi }}}{{\rm{\Sigma }}}{\sin }^{2}\vartheta \,{\rm{d}}{\varphi }^{2}\right].\end{array}\end{eqnarray} \tag{ 2 }$$

Here the parameter a is the spin of the central object, and the stationary metric functions {F, Δ, Σ, Π} can be related to the static metric functions {f, g, R} in Equation (1), as discussed below. This then leaves only the metric function X = X(r, ϑ) to be fixed by the field equations (Azreg-Aïnou 2014b).

While we focus in this paper on the AA metric, we note that different generalizations of the static metric (Equation (1)) are possible, in principle, and can be achieved by modifying the complex coordinate transformation involved in the Newman & Janis–type (1965) solution-generating technique (see, e.g., Azreg-Aïnou 2014a). It is important to note, however, that the particular transformations used in Azreg-Aïnou (2014a, 2014b) send the Schwarzschild, RN, and GMGHS metrics to their appropriate spinning generalizations, namely, the Kerr, Kerr–Newman, and Kerr–Sen metrics, respectively. We direct the reader to Erbin (2017) for a review of such solution-generating techniques and Section 7.1 of Wald (1984) for a general discussion of stationary, axisymmetric spacetimes and the construction of BL-like coordinates.

Expressions for the metric functions {F, Δ, Σ, Π} in the AA metric (Equation (2)) are most conveniently written by first defining two auxiliary functions that are specific combinations of the original static metric functions, f(r), g(r), and R(r),

$$\begin{eqnarray*}&&A(r)={R}^{2}/\sqrt{g},\ B(r)=(f/g){R}^{2},\end{eqnarray*}$$

which can be written even more transparently as

$$\begin{eqnarray}&&A=\displaystyle \frac{{\hat{g}}_{\vartheta \vartheta }}{\sqrt{-{\hat{g}}_{{tt}}{\hat{g}}_{{rr}}}};\ B=\displaystyle \frac{{\hat{g}}_{\vartheta \vartheta }}{{\hat{g}}_{{rr}}}={\hat{g}}^{{rr}}{\hat{g}}_{\vartheta \vartheta }.\end{eqnarray} \tag{ 3 }$$

Here A and B are each functions of r alone, and in terms of them, the functions {F, Δ, Σ, Π} in the AA metric (Equation (2)) take the remarkably simple form

$$\begin{eqnarray}\begin{array}{rcl}F(r) & = & (A-B)/2,\\ {\rm{\Delta }}(r) & = & B+{a}^{2},\\ {\rm{\Sigma }}(r,\vartheta ) & = & A+{a}^{2}{\cos }^{2}\vartheta ,\\ {\rm{\Pi }}(r,\vartheta ) & = & {\left(A+{a}^{2}\right)}^{2}-{\rm{\Delta }}{a}^{2}{\sin }^{2}\vartheta .\end{array}\end{eqnarray} \tag{ 4 }$$

Note that F and Δ depend only on r, while Σ and Π are functions of r and ϑ. The spinning AA metric g_{μ ν} (Equation (2)) reduces to the nonspinning metric ${\hat{g}}_{\mu \nu }$ (Equation (1)) in the limit of vanishing spin (a = 0) for appropriate choices of X (see Azreg-Aïnou 2014a for a discussion of conformal versus normal fluids).

The determinant $\det [{g}_{\mu \nu }]$ of the AA metric and its t φ − sector $\det [{g}_{t\varphi }]:= {g}_{{tt}}{g}_{\varphi \varphi }-{g}_{t\varphi }^{2}$ are given in BL coordinates by

$$\begin{eqnarray}&&\det [{g}_{\mu \nu }]=-({X}^{4}/{{\rm{\Sigma }}}^{2}){\sin }^{2}\vartheta ,\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&\det [{g}_{t\varphi }]=-({X}^{2}/{{\rm{\Sigma }}}^{2}){\rm{\Delta }}{\sin }^{2}\vartheta .\end{eqnarray} \tag{ 6 }$$

The AA metric can also be expressed in the closely related quasi-isotropic (QI) coordinates (see, e.g., Section 2.3.2 of Gourgoulhon 2010), ${x}^{\tilde{\mu }}=(t,\tilde{r},\vartheta ,\varphi )$, as

$$\begin{eqnarray}\begin{array}{rcl}{\rm{d}}{s}^{2} & = & \displaystyle \frac{X}{{\rm{\Sigma }}}\left[-\displaystyle \frac{X{\rm{\Delta }}}{{\rm{\Pi }}}{\rm{d}}{t}^{2}+\displaystyle \frac{{\rm{\Sigma }}}{{\tilde{r}}^{2}}\left({\rm{d}}{\tilde{r}}^{2}+{\tilde{r}}^{2}{\rm{d}}{\vartheta }^{2}\right)\right.\\ & & \left.+\,\displaystyle \frac{{\rm{\Pi }}}{{\rm{\Sigma }}}{\sin }^{2}\vartheta {\left({\rm{d}}\varphi -{{\rm{\Omega }}}_{{\rm{Z}}}{\rm{d}}t\right)}^{2}\right],\end{array}\end{eqnarray} \tag{ 7 }$$

where we have introduced the angular velocity of the zero angular momentum observer (ZAMO), Ω_Z = − g_{t φ} /g_{φ φ} = 2aF/Π (see Equations (2.4) and (2.5) of Bardeen et al. 1972). The form above (Equation (7)) can be obtained from Equation (2) by transforming the radial coordinate, $r\mapsto \tilde{r}$, via $\tilde{r}(r)=\exp [\int ({\rm{d}}r/\sqrt{{\rm{\Delta }}})]$. The r ϑ − sector of the AA metric in these QI coordinates is clearly conformally flat, as desired. The related cylindrical coordinates, ${x}^{\mu ^{\prime} }=(t,\rho ,z,\varphi )$, with $\rho := \tilde{r}\sin \vartheta$ and $z:= \tilde{r}\cos \vartheta$, are called the Lewis–Papapetrou coordinates (Lewis 1932; Papapetrou 1966).

We now briefly illustrate how the remaining metric function X can be fixed by considering, as an example, the case when the background matter is an anisotropic fluid in GR that flows around the spin axis. Since, in its own rest frame, $\{{e}_{(a)}^{\mu }\ (a=0-3)\}$, the fluid stress–energy–momentum tensor is diagonal, the Einstein equations imply that the Einstein tensor in this frame, ${{\mathscr{G}}}_{(a)(b)}={{\mathscr{G}}}_{\mu \nu }{e}_{(a)}^{\mu }{e}_{(b)}^{\nu }$, must also be diagonal. However, it can be shown that two of the off-diagonal elements of the Einstein tensor, ${{\mathscr{G}}}_{(r)(\vartheta )}$ and ${{\mathscr{G}}}_{(t)(\varphi )}$, will generally be nonzero. Demanding that these two terms vanish yields one nonlinear partial differential equation (PDE) and one linear PDE, each of which involves the free metric function X, the auxiliary function A(r), and the fluid angular velocity Ω. These PDEs are given in Equations (15) and (18) of Azreg-Aïnou (2014a), respectively, where the fluid angular velocity was chosen to be Ω = a/(A + a²) (see ${e}_{t}^{\mu }$ in their Equation (16)). It is interesting to note that this angular velocity does not correspond to the ZAMO angular velocity introduced above, Ω_Z = 2aF/Π. It matches, however, the angular velocity of the PNCs of the spacetime introduced below (see Equation (8)). For nonfluid matter models, additional equations of motion for the matter fields must be solved. For example, if the matter is a scalar field, the associated Klein–Gordon equations must additionally be solved. Similarly, if the matter is an electromagnetic field, the Maxwell equations have to be accounted for as well. We will not enter into a discussion of this topic but refer the reader to Erbin (2017).

For the purposes of the present paper, it suffices to note that the subclass of AA metrics with X = Σ is particularly interesting, as it corresponds to spacetimes that are asymptotically flat (see Appendix B for further details). Furthermore, it can be checked that this choice for the conformal factor X always solves the ${{\mathscr{G}}}_{(r)(\vartheta )}=0$ equation for axially spinning matter. Finally, all of the BH solutions discussed in this paper have X = Σ, whether they are vacuum or contain scalar, electromagnetic, and/or axion fields. Even though these spacetimes arise as "solution metrics" to a variety of field equations, $X={\rm{\Sigma }}=A(r)+{a}^{2}{\cos }^{2}\vartheta$ consistently remains a valid choice. For this subclass of AA metrics, the possibly divergent behavior of the Ricci ${\mathscr{R}}$ and Kretschmann ${\mathscr{K}}$ scalars can be seen from the behavior of their denominators, which go roughly as ∼Σ⁻³ and ∼Σ⁻⁶, respectively. Thus, when this metric describes a singular spacetime, the curvature singularity is located at r = r_⋆, such that Σ(r_⋆, π/2) = 0, which corresponds to a ring singularity located in the equatorial plane. This last equation is equivalent to R(r_⋆) = 0.

When the general AA metric (Equation (2)) describes a BH spacetime, as above, a horizon is present at every location, r = r_H > r_⋆, where g^rr = Δ/Σ = 0. Equivalently, horizons are located at the real, positive (R(r_H) > 0) roots of Δ(r). Conversely, if no such roots exist, then the AA metric does not correspond to a BH spacetime (e.g., the a > M Kerr metric). The event horizon in particular is located at the largest such root. Since Δ is a function of r alone, it is reassuring to find that the event horizon is indeed a round sphere in BL coordinates. Furthermore, the horizon angular velocity, Ω_H $:=$ Ω_Z(r_H), is simply given as ${{\rm{\Omega }}}_{{\rm{H}}}\,=2aF({r}_{{\rm{H}}})/(A({r}_{{\rm{H}}})+{a}^{2}{)}^{2}=a/(A({r}_{{\rm{H}}})+{a}^{2})$. Clearly, Ω_H is independent of ϑ, i.e., an AA BH rotates like a rigid body, as a consequence of the weak rigidity theorem (see, e.g., Section 8.4.4 of Straumann 2013). These desirable properties, which are well known for the Kerr and Kerr–Newman metrics, are now seen to be true for the very wide class of AA stationary BH metrics.

It is interesting to note that the general AA metric (Equation (2)) can also describe spacetimes that contain regions that admit closed time-like curves. In such regions, g_{φ φ} < 0 (see, e.g., Section V.B of Johannsen 2013b), i.e., Π < 0 (e.g., the a > M Kerr metric).

We show in Appendix C that the null geodesic equations are Liouville-separable. Thus, the AA metric always admits a Carter constant for null geodesics. Furthermore, we also identify a subclass of AA metrics that additionally admit separable time-like geodesic equations. These require the conformal factor to take the form X(r, ϑ) = X_r (r) + X_ϑ (ϑ). For this subclass, a Carter constant exists for all geodesics. Concomitantly, arbitrary geodesics of all AA metrics with X = Σ possess Carter constants.

The rr-component of the AA metric (Equation (2)) in BL coordinates diverges at the horizon (Δ(r) = 0). The Ricci and Kretschmann scalars, however, reveal this to be merely a coordinate singularity, an artifact of the choice of coordinates. The singularity can be eliminated by a change of coordinates, and this is the topic we turn to next.

In this section, we convert the AA metric, which is written in BL coordinates in Equation (2), to a horizon-penetrating form in which no coordinate singularity is present at the horizon. This is the form that is most useful for GRMHD and xGRMHD simulations of astrophysical accretion flows.

The ingoing (−) and outgoing (+) PNCs of the AA spacetime consist of null geodesics x^μ (λ) whose tangents ${{\ell }}_{\pm }^{\mu }={\dot{x}}^{\mu }={\rm{d}}{x}^{\mu }/{\rm{d}}\lambda$ satisfy ${{\ell }}_{\pm }^{\vartheta }=\dot{\vartheta }=0$ and $\ddot{\vartheta }=0$ (see, e.g., Misner et al. 1973; Hioki & Miyamoto 2008 for further details). The 1 − form fields ${({l}_{\pm })}_{\mu }$ associated with the PNCs in an AA spacetime can be expressed elegantly as ⁷

$$\begin{eqnarray}&&{\left({{\ell }}_{\pm }\right)}_{\mu }=E\left[-1,\pm \displaystyle \frac{{\rm{\Sigma }}}{{\rm{\Delta }}},0,a{\sin }^{2}\vartheta \right],\end{eqnarray} \tag{ 9 }$$

where E is some constant.

We now define the spherical polar ingoing Kerr–Schild (siKS) coordinates, ${x}^{\bar{\mu }}=(\tau ,r,\vartheta ,\phi )$, as those in which the ingoing PNC takes the form

$$\begin{eqnarray}&&{\left({{\ell }}_{-}\right)}_{\bar{\mu }}=E\left[-1,-1,0,a{\sin }^{2}\vartheta \right].\end{eqnarray} \tag{ 10 }$$

To distinguish between the BL (x^μ ) and siKS (${x}^{\bar{\mu }}$) coordinate systems, we will use a bar for the indices of the latter. We are interested in the ingoing PNC because our coordinate system should penetrate the future horizon ${{\mathscr{H}}}^{+}$ when one exists and cover patches I and II of the Kruskal (Kruskal 1960) or Penrose–Carter (Penrose 1963; Carter 1966) diagrams. The Jacobian ${{\rm{\Lambda }}}_{\ \mu }^{\bar{\mu }}={\partial }_{\mu }{x}^{\bar{\mu }}$ for the coordinate transformation from BL to siKS coordinates (i.e., ${\rm{d}}{x}^{\bar{\mu }}={{\rm{\Lambda }}}_{\ \mu }^{\bar{\mu }}{\rm{d}}{x}^{\mu }$) can be inferred from Equations (9) and (10) to be

$$\begin{eqnarray}{{\rm{\Lambda }}}_{\ \mu }^{\bar{\mu }}=\left[\begin{array}{cccc}1 & 2F/{\rm{\Delta }} & 0 & 0\\ 0 & 1 & 0 & 0\\ 0 & 0 & 1 & 0\\ 0 & a/{\rm{\Delta }} & 0 & 1\end{array}\right].\end{eqnarray} \tag{ 11 }$$

The time coordinate τ is then clearly a "tortoise" time coordinate (see, e.g., Blau 2023),

$$\begin{eqnarray}&&\tau =t+\int (2F/{\rm{\Delta }})\,{\rm{d}}r,\end{eqnarray} \tag{ 12 }$$

and we recognize immediately that, in the limit of vanishing spin (a → 0), the siKS coordinates are analogous to the coordinates presented in Eddington (1924) and Finkelstein (1958).

Using the above coordinate transformation, we can write the AA metric in siKS coordinates as

$$\begin{eqnarray}{g}_{\bar{\mu }\bar{\nu }}=\displaystyle \frac{X}{{\rm{\Sigma }}}\left[\begin{array}{cccc}-\left(1-\displaystyle \frac{2F}{{\rm{\Sigma }}}\right) & \displaystyle \frac{2F}{{\rm{\Sigma }}} & 0 & -\left(\displaystyle \frac{2F}{{\rm{\Sigma }}}\right)a{\sin }^{2}\vartheta \\ * & \left(1+\displaystyle \frac{2F}{{\rm{\Sigma }}}\right) & 0 & -\left(1+\displaystyle \frac{2F}{{\rm{\Sigma }}}\right)a{\sin }^{2}\vartheta \\ * & * & {\rm{\Sigma }} & 0\\ * & * & * & \left(\displaystyle \frac{{\rm{\Pi }}}{{\rm{\Sigma }}}\right){\sin }^{2}\vartheta \end{array}\right],\end{eqnarray} \tag{ 13 }$$

where the asterisks denote components that are fixed by symmetry, ${g}_{\bar{\mu }\bar{\nu }}={g}_{\bar{\nu }\bar{\mu }}$. Notice that the metric components now diverge only at the spacetime/curvature singularity, present at Σ = 0. The inverse AA metric in siKS coordinates is

$$\begin{eqnarray}{g}^{\bar{\mu }\bar{\nu }}=\displaystyle \frac{{\rm{\Sigma }}}{X}\left[\begin{array}{cccc}-\left(1+\displaystyle \frac{2F}{{\rm{\Sigma }}}\right) & \displaystyle \frac{2F}{{\rm{\Sigma }}} & 0 & 0\\ * & \displaystyle \frac{{\rm{\Delta }}}{{\rm{\Sigma }}} & 0 & \displaystyle \frac{a}{{\rm{\Sigma }}}\\ * & * & \displaystyle \frac{1}{{\rm{\Sigma }}} & 0\\ * & * & * & \displaystyle \frac{1}{{\rm{\Sigma }}{\sin }^{2}\vartheta }\end{array}\right].\end{eqnarray} \tag{ 14 }$$

For completeness, we cast the general AA metric into its "classic" Kerr–Schild form in Equation (B2) of Appendix B.

In siKS coordinates, the event horizon is still located at the outermost root of ${g}^{\bar{1}\bar{1}}={g}^{{rr}}=0$, i.e., Δ(r) = 0 (Equation (14)), but as we can see from Equation (13), the metric components ${g}_{\bar{\mu }\bar{\nu }}$ no longer diverge at Δ = 0. Thus, the coordinate singularity at the horizon has been removed by changing from BL to siKS coordinates, showing that the latter are horizon-penetrating coordinates. Furthermore, a comparison of the Kerr metric in particular, in these coordinates, with that reported in Font et al. (1998), with k = 1 there, and in McKinney & Gammie (2004), shows that these are the horizon-penetrating coordinates canonically used in the context of GRMHD simulations.

Note that the original set of Kerr–Schild coordinates used to describe the Kerr spacetime (Kerr 1963; Kerr & Schild 2009) are adapted to the outgoing PNC in which ${l}_{+}^{\breve{\mu }}\propto {\delta }_{r}^{\breve{\mu }}$ (see, e.g., Wiltshire et al. 2009; Azreg-Aïnou 2014b). The corresponding time coordinate u is the retarded time, as well as a null affine parameter that parameterizes the outgoing PNC on future null infinity ${{\mathscr{I}}}^{+}$. Equivalently, the ingoing version of the original set of Kerr–Schild coordinates would be adapted to the ingoing PNC in which ${l}_{-}^{\breve{\mu }}\propto {\delta }_{r}^{\breve{\mu }}$, with the corresponding time coordinate v being the advanced time, as well as a null affine parameter that parameterizes the ingoing PNC on past null infinity ${{\mathscr{I}}}^{-}$. Clearly, the "original" ingoing Kerr–Schild system is quite different from the spherical ingoing Kerr–Schild coordinates that we use in this section and that are matched to current GRMHD codes. From Equation (10), the tangent to the ingoing PNC in our coordinates is given as

$$\begin{eqnarray}&&{{\ell }}_{-}^{\bar{\mu }}=E\displaystyle \frac{{\rm{\Sigma }}}{X}\left[1,-1,0,0\right].\end{eqnarray} \tag{ 15 }$$

We can cast the AA metric when written in horizon-penetrating siKS coordinates into its associated 3+1 form (Arnowitt et al. 2008) by introducing a foliation into space-like hypersurfaces Σ_τ that are isosurfaces of the scalar time function τ (see, e.g., Poisson 2004; Gourgoulhon 2007; Alcubierre 2008; Rezzolla & Zanotti 2013). The unit time-like normal to these hypersurfaces is

$$\begin{eqnarray}&&{n}_{\bar{\mu }}:= -\alpha {{\rm{\nabla }}}_{\bar{\mu }}\tau =-\alpha {\delta }_{\bar{\mu }}^{\tau },\end{eqnarray} \tag{ 16 }$$

where $\alpha := \sqrt{-1/{g}^{\tau \tau }}$ is called the lapse function and, from Equation (14), is given by

$$\begin{eqnarray}&&\alpha ={[X/({\rm{\Sigma }}+2F)]}^{1/2}.\end{eqnarray} \tag{ 17 }$$

The four-vector that is dual to the normal, i.e., ${n}^{\bar{\mu }}={g}^{\bar{\mu }\bar{\nu }}{n}_{\bar{\nu }}\,=-\alpha {g}^{\bar{\mu }\tau }=[1/\alpha ,-\alpha (2F/X),0,0]$, corresponds to the four-velocity of the local time-like Eulerian observer. We note that, due to Frobenius' theorem, the rotation or vorticity tensor for this congruence of Eulerian observers vanishes identically since it is hypersurface-orthogonal or hypersurface-forming (see, e.g., Section 2.3.3. of Poisson 2004). The four-velocity of this observer in BL coordinates can be obtained as ${n}^{\mu }={{\rm{\Lambda }}}_{\ \bar{\mu }}^{\mu }{n}^{\bar{\mu }}$, where ${{\rm{\Lambda }}}_{\bar{\mu }}^{\mu }=({{\rm{\Lambda }}}_{\mu }^{\bar{\mu }}{)}^{-1}$. It is straightforward to then check that the angular velocity of the Eulerian observer in BL coordinates Ω_E $:=$ n^φ /n^t matches the angular velocity of the ZAMO introduced above, Ω_Z = − g_{t φ} /g_{φ φ} = 2aF/Π, exactly. However, the Eulerian observer has nonzero BL radial velocity, i.e., ${n}^{r}=-2\alpha F/X(={n}^{\bar{r}})$.

With Equation (16), the induced (Riemannian) metric γ on the space-like hypersurfaces can now be introduced as

$$\begin{eqnarray}&&{\gamma }_{\bar{\mu }\bar{\nu }}:= {g}_{\bar{\mu }\bar{\nu }}+{n}_{\bar{\mu }}{n}_{\bar{\nu }}.\end{eqnarray} \tag{ 18 }$$

The projection tensor, which is used to obtain the spatial components of any four-vector or tensor, is simply ${\gamma }_{\ \bar{\nu }}^{\bar{\mu }}={g}^{\bar{\mu }\bar{\alpha }}{\gamma }_{\bar{\alpha }\bar{\nu }}={\delta }_{\ \bar{\nu }}^{\bar{\mu }}+{n}^{\bar{\mu }}{n}_{\bar{\nu }}$, since ${\gamma }_{\ \bar{\nu }}^{\bar{\mu }}{n}^{\bar{\nu }}=0$.

Although the 1-form ${n}_{\bar{\mu }}{\rm{d}}{x}^{\bar{\mu }}$ is collinear with dτ (see Equation (16)), the dual four-vector ${n}^{\bar{\mu }}{\partial }_{\bar{\mu }}$ is not, in general, collinear with ∂_τ (since $-\alpha {g}^{\bar{\mu }\tau }{/}\!\!\!\!{\propto }\ {\delta }_{\tau }^{\bar{\mu }}$). Indeed, this noncollinearity is typically captured via a four-vector ${\beta }^{\bar{\mu }}$, which is defined such that ${n}_{\bar{\mu }}{\beta }^{\bar{\mu }}=0$:

$$\begin{eqnarray}&&{\beta }^{\bar{\mu }}:= {\left({\partial }_{\tau }\right)}^{\bar{\mu }}-\alpha {n}^{\bar{\mu }}={\delta }_{\tau }^{\bar{\mu }}-\alpha {n}^{\bar{\mu }}.\end{eqnarray} \tag{ 19 }$$

Since, by construction, β^τ = 0, it is convenient to work instead with the space-like three-vector (see, e.g., McKinney & Gammie 2004; Porth et al. 2017),

$$\begin{eqnarray}&&{\beta }^{i}:= {\alpha }^{2}{g}^{i\tau }=-{g}^{i\tau }/{g}^{\tau \tau }=[2F/(2F+{\rm{\Sigma }})]{\delta }_{r}^{i},\end{eqnarray} \tag{ 20 }$$

where the index i takes values i = 1, 2, and 3. This three-vector measures the shift of the isolines of the three spatial coordinates with respect to the normal to the hypersurfaces and is therefore called the shift vector (see, e.g., Section 1.3 of Gourgoulhon 2010).

The desired 3+1 form of the AA metric can now be written using the lapse function α (Equation (17)), the shift vector βⁱ (Equation (20)), and the components γ_ij (Equation (18)) of the induced metric on the hypersurfaces Σ_τ as

$$\begin{eqnarray}&&{\rm{d}}{s}^{2}=-{\alpha }^{2}{\rm{d}}{\tau }^{2}+{\gamma }_{{ij}}({\rm{d}}{x}^{i}+{\beta }^{i}{\rm{d}}\tau )({\rm{d}}{x}^{j}+{\beta }^{j}{\rm{d}}\tau ).\end{eqnarray} \tag{ 21 }$$

For completeness, we demonstrate how the temporal u^τ and spatial uⁱ components of a vector ${u}^{\bar{\mu }}=({u}^{\tau },{u}^{i})$ can be obtained using the unit normal and the projection tensor, respectively. First, with ${\rm{\Gamma }}=-{n}_{\bar{\mu }}{u}^{\bar{\mu }}$, we have u^τ = Γ/α. Second, introducing the purely spatial vectors (v⁰ = 0) ${v}^{i}:= ({\gamma }_{\ \bar{\mu }}^{i}{u}^{\bar{\mu }})/{\rm{\Gamma }}$, we obtain uⁱ = Γ(vⁱ − βⁱ /α). One also finds that ${\rm{\Gamma }}=1/\sqrt{1-{v}^{2}}$, where v² = v_i vⁱ . When u corresponds to the four-velocity of a fluid element, Γ has the significance of being the fluid Lorentz factor with respect to the Eulerian observer. An excellent summary of the 3+1 treatment relevant for GRMHD can be found in Gammie et al. (2003).

We have presented in Equation (13) in Section 4 a horizon-penetrating Kerr–Schild form of the rather general stationary and axisymmetric Azreg–Aïnou (AA) metric. This form of the metric can be used as an input when performing general relativistic (GR), as well as non-GR, magnetohydrodynamics simulations of accretion flows onto compact objects in arbitrary metric theories of gravity. The AA metric can be used to generate a number of popular metrics that describe the stationary spacetimes in GR corresponding to electrovacuum black holes (BHs; Kerr 1963; Newman et al. 1965), regular BHs such as the Kerr–Hayward models (Hayward 2006; Bambi & Modesto 2013; Zhou & Modesto 2023a, 2023b), naked singularities such as the spinning Janis–Newman–Winicour spacetime (Solanki et al. 2022), etc. It can also be used to describe BH solutions arising in alternative theories of gravity (e.g., Sen 1992).

To illustrate the ease of generating the Kerr–Schild forms of new stationary metrics, we now summarize the steps involved. We start with an arbitrary static and spherically symmetric seed metric, e.g., the ones listed in Table 1, in standard spherical polar coordinates, x^μ = (t, r, ϑ, φ),

$$\begin{eqnarray}&&{\rm{d}}{s}^{2}=-f(r){\rm{d}}{t}^{2}+\displaystyle \frac{g(r)}{f(r)}{\rm{d}}{r}^{2}+{R}^{2}(r){\rm{d}}{{\rm{\Omega }}}_{2}^{2}.\end{eqnarray} \tag{ 22 }$$

As we noted in Section 2, it is always possible to find a radial coordinate r in which g(r) = 1. Since this reduces the number of free functions, we used such coordinates for the static models listed in Table 1. However, the required transformation may not always lead to simple analytical forms for f(r) and R(r). When it does not, it is preferable to use the more general form (Equation (22)).

The AA stationary generalization of the metric (Equation (22)) in spherical polar ingoing Kerr–Schild coordinates, ${x}^{\bar{\mu }}=(\tau ,r,\vartheta ,\phi )$, takes the form

$$\begin{eqnarray}\begin{array}{rcl}{\rm{d}}{s}^{2} & = & \displaystyle \frac{X}{{\rm{\Sigma }}}\left[-\left(1-\displaystyle \frac{2F}{{\rm{\Sigma }}}\right){\rm{d}}{\tau }^{2}+\left(1+\displaystyle \frac{2F}{{\rm{\Sigma }}}\right){\rm{d}}{r}^{2}\right.\\ & & \left.+\,{\rm{\Sigma }}\,{\rm{d}}{\vartheta }^{2}+\displaystyle \frac{{\rm{\Pi }}}{{\rm{\Sigma }}}{\sin }^{2}\vartheta \,{\rm{d}}{\phi }^{2}-2\displaystyle \frac{2F}{{\rm{\Sigma }}}a{\sin }^{2}\vartheta {\rm{d}}\tau {\rm{d}}\phi\right. \\ & & \left.+\,2\displaystyle \frac{2F}{{\rm{\Sigma }}}{\rm{d}}\tau {\rm{d}}r-2\left(1+\displaystyle \frac{2F}{{\rm{\Sigma }}}\right){\rm{d}}r{\rm{d}}\phi \right],\end{array}\end{eqnarray} \tag{ 23 }$$

where the stationary metric functions {F, Δ, Σ, Π} are fixed by the static metric functions {f, g, R} in the seed metric (Equation (22)) as

$$\begin{eqnarray}\begin{array}{rcl}2F(r) & = & {R}^{2}/\sqrt{g}-(f/g){R}^{2},\\ {\rm{\Delta }}(r) & = & (f/g){R}^{2}+{a}^{2},\\ {\rm{\Sigma }}(r,\vartheta ) & = & {R}^{2}/\sqrt{g}+{a}^{2}{\cos }^{2}\vartheta ,\\ {\rm{\Pi }}(r,\vartheta ) & = & {\left({R}^{2}/\sqrt{g}+{a}^{2}\right)}^{2}-{\rm{\Delta }}{a}^{2}{\sin }^{2}\vartheta .\end{array}\end{eqnarray} \tag{ 24 }$$

The parameter a describes the spin of the spacetime. The only remaining freedom is the stationary metric function X(r, ϑ).

If the background matter in the stationary spacetime is a fluid flowing around the spin axis (ϑ = 0) in GR, the unknown function X = X(r, ϑ) is fixed by solving two Einstein equations, ${{\mathscr{G}}}_{(r)(\vartheta )}=0$ and ${{\mathscr{G}}}_{(t)(\varphi )}=0$, in the fluid rest frame (the energy–momentum–stress components T_(r)(ϑ) and T_(t)(φ) of the fluid vanish in this frame). We will not enter into the detailed implications of these equations here but note that X = Σ is always a solution of the first equation, ${{\mathscr{G}}}_{(r)(\vartheta )}=0$. The choice X = Σ is additionally appealing because spacetimes with this property are asymptotically flat (provided, of course, that the original seed static metric is also asymptotically flat). Verifying that X = Σ also solves the second Einstein equation, ${{\mathscr{G}}}_{(t)(\varphi )}=0$, in the fluid frame, and that the resulting energy–momentum–stress tensor is physically valid, i.e., satisfies various energy conditions, is an involved calculation and beyond the scope of this paper. However, if one merely wants a novel stationary metric in which to carry out GRMHD simulations, and if one is not necessarily too concerned about satisfying energy conditions, then the metric (Equation (23)) with X = Σ and a suitable choice of the seed functions f(r), g(r), and R(r) would be a reasonable choice.

We should note that, for nonfluid models, additional equations of motion for the fields (e.g., the Klein–Gordon or Maxwell equations) must be solved. The status and successes of such metric solution–generating techniques in being able to also yield solutions to the nongravitational field equations is discussed in Erbin (2017).

Table 2 presents a compilation of the metric functions of the stationary AA spacetimes corresponding to the six static seed models listed in Table 1. These are all models for which the choice X = Σ is valid. The stationary Kerr (1963) metric arises naturally as the AA generalization of the static Schwarzschild metric. It describes vacuum, spinning BHs in GR and has been used extensively in GRMHD simulations. Similarly, the AA transformation automatically generates the Kerr–Newman metric (Newman et al. 1965), which describes electromagnetically charged, spinning BHs in GR, as the spinning generalization of the nonspinning Reissner–Nordström BH. The Kerr–Hayward model proposed here (the third model in Table 2) is expected to describe spinning regular (i.e., no singularities) BHs in GR (see, e.g., Bambi & Modesto 2013; Zhou & Modesto 2023b). The spinning equivalent of the GMGHS BH solution in string theory is given by the Kerr–Sen metric (Sen 1992). The form of the Kerr–Sen metric shown in Table 2 can be put into the form presented in Equation (20) of Xavier et al. (2020) by replacing r with r + Q²/M. It is worth noting that the spacetime of Kerr–Sen BHs contains not just electromagnetic fields but also a scalar field (dilaton), as well as an axion field. All of the BH solutions in Table 2 reduce to the Kerr BH in appropriate limits (Q → 0, L → 0).

Table 2. Metric Functions for the Stationary and Axisymmetric Generalizations of the Static and Spherically Symmetric Spacetimes Given in Table 1

Spacetime	2F = (1 − f)R2	Δ = fR2 + a2	${\rm{\Sigma }}=X={R}^{2}+{a}^{2}{\cos }^{2}\vartheta$
Kerr	2Mr	r2 − 2Mr + a2	${r}^{2}+{a}^{2}{\cos }^{2}\vartheta$
Kerr–Newman	2Mr − Q2	r2 − 2Mr + Q2 + a2	${r}^{2}+{a}^{2}{\cos }^{2}\vartheta$
Kerr–Hayward	$\tfrac{2{{Mr}}^{5}}{{r}^{4}+2{L}^{4}}$	${r}^{2}-\tfrac{2{{Mr}}^{5}}{{r}^{4}+2{L}^{4}}+{a}^{2}$	${r}^{2}+{a}^{2}{\cos }^{2}\vartheta$
Kerr–Sen	2Mr − 2Q2	${r}^{2}-\left(2M+\tfrac{{Q}^{2}}{M}\right)r+2{Q}^{2}+{a}^{2}$	$r\left(r-\tfrac{{Q}^{2}}{M}\right)+{a}^{2}{\cos }^{2}\vartheta$
Spinning JNW	${r}^{2}\left[{\left(1-\tfrac{{r}_{\star }}{r}\right)}^{\hat{\nu }}-\left(1-\tfrac{{r}_{\star }}{r}\right)\right]$	${r}^{2}-\tfrac{2{Mr}}{1-\hat{\nu }}+{a}^{2}$	${r}^{2}{\left(1-\tfrac{{r}_{\star }}{r}\right)}^{\hat{\nu }}+{a}^{2}{\cos }^{2}\vartheta$
Spinning JMN-1	$\tfrac{{r}^{2}}{{(1-\sigma )}^{2}}\left[{\left(\tfrac{r}{{r}_{{\rm{b}}}}\right)}^{-2\sigma }-\left(1-\tfrac{2M}{{r}_{{\rm{b}}}}\right)\right]$	$\left(1-\tfrac{2M}{{R}_{{\rm{b}}}}\right)\tfrac{{r}^{2}}{{(1-\sigma )}^{2}}+{a}^{2}$	${R}_{{\rm{b}}}^{2}{\left(\tfrac{r}{{r}_{{\rm{b}}}}\right)}^{2-2\sigma }+{a}^{2}{\cos }^{2}\vartheta$

Note. The spinning generalizations of a nonspinning seed metric (Equation (1)) are given in BL coordinates in Equation (2) and horizon-penetrating spherical ingoing Kerr–Schild coordinates in Equation (23). These involve five stationary metric functions {F, Δ, Σ, Π, X}, the first four of which are related to f and R via Equation (24). We show above the spinning generalizations of the specific static solutions listed in Table 1. Since we used coordinates for the static metrics in which g(r) = 1, g does not appear in the expressions below. For brevity, we only show the four metric functions 2F, Δ, Σ, and X above. The fifth metric function is given by ${\rm{\Pi }}={({R}^{2}+{a}^{2})}^{2}-{\rm{\Delta }}{a}^{2}{\sin }^{2}\vartheta$. All of the spacetime models listed here are asymptotically flat, and a denotes the usual spin parameter. A variety of physical scenarios are captured by the exemplary set listed here: the vacuum Kerr BH, the electromagnetically charged Kerr–Newman BH, the (modified) Kerr–Hayward spinning BH, the spinning Janis–Newman–Winicour (JNW) naked singularity, and the spinning JMN-1 naked singularity are solutions of GR, whereas the Kerr–Sen spacetime describes the electromagnetically charged BHs of string theory, which are also charged under a dilaton field, as well as an axion field.

Download table as:  ASCIITypeset image

The spinning JNW spacetime (the fifth model in Table 2) was proposed in Solanki et al. (2022) and is expected to describe a spinning naked singularity spacetime containing scalar field matter. Finally, the spinning JMN-1 metric proposed here is expected to describe a spinning naked singularity spacetime containing an anisotropic fluid in GR. Careful explorations of these spacetimes, especially with a focus on the physicality of the underlying matter (e.g., how well energy conditions are satisfied), is left for future work.

Vagnozzi et al. (2023) gave an immense compilation of spacetime metrics, including those corresponding to static and spherically symmetric spacetimes. These metric functions describe alternative nonspinning BHs, which are of considerable interest and can be studied using simulations. More importantly, each static model could potentially be used as a seed metric to obtain the horizon-penetrating form of the stationary metric that it may belong to via Equations (23) and (24). The only open issue, which needs to be checked in each case, is whether the choice of X = Σ is physically valid.

We conclude on an optimistic note, envisioning the widespread adoption of GRMHD simulations facilitated by our easy-to-use metric formulation. A Kerr simulation library has recently been used successfully to infer the properties of plasma present in the close vicinity of Sgr A*, as well as the spacetime geometry of this supermassive BH (The EHT Collaboration et al. 2022a, 2022b). While it might be a challenging endeavor in the near future to construct similar libraries of fully 3D simulations in the other stationary spacetimes described in this paper (because of the many additional spacetime parameters in these models, such as Q, L, $\hat{\nu }$, and σ), it should at least be possible to build extensive 2D GRMHD simulation libraries for these spacetimes. This will enable valuable exploration of magnetized relativistic gas dynamics in diverse spacetimes.

We thank the referee for helpful suggestions. P.K., R.N., and K.C. acknowledge support in part from grants from the Gordon and Betty Moore Foundation and the John Templeton Foundation to the Black Hole Initiative at Harvard University and from NSF award OISE-1743747. Y.M. is supported by the National Natural Science Foundation of China (grant No. 12273022) and the Shanghai target program of basic research for international scientists (grant No. 22JC1410600).

Given an arbitrary spherically symmetric and static metric ${\hat{g}}_{\mu \nu }$ in arbitrary spherical polar coordinates of the form (Equation (1))

$$\begin{eqnarray}&&{\rm{d}}{s}^{2}=\ -f(r){\rm{d}}{t}^{2}+\displaystyle \frac{g(r)}{f(r)}{\rm{d}}{r}^{2}+{R}^{2}(r){\rm{d}}{{\rm{\Omega }}}_{2}^{2},\end{eqnarray} \tag{ A1 }$$

we can put it into the form where g(r) = 1 as

$$\begin{eqnarray}\begin{array}{rcl}{\rm{d}}{s}^{2} & = & -f(r){\rm{d}}{t}^{2}+\displaystyle \frac{g(r)}{f(r)}{\left({\partial }_{\rho }r\right)}^{2}{\rm{d}}{\rho }^{2}+{R}^{2}(r){\rm{d}}{{\rm{\Omega }}}_{2}^{2},\\ & = & -f(\rho ){\rm{d}}{t}^{2}+\displaystyle \frac{{\rm{d}}{\rho }^{2}}{f(\rho )}+{R}^{2}(\rho ){\rm{d}}{{\rm{\Omega }}}_{2}^{2},\end{array}\end{eqnarray} \tag{ A2 }$$

where, in writing the last equality, we have imposed the condition that r(ρ) satisfies the ordinary differential equation,

$$\begin{eqnarray}&&\sqrt{g(r)}\,{\partial }_{\rho }r=1,\end{eqnarray} \tag{ A3 }$$

and f(ρ) = f(r(ρ)). Note that the above equation can be rewritten as

$$\begin{eqnarray}&&\sqrt{-{\hat{g}}_{{tt}}(r){\hat{g}}_{{rr}}(r)}\,{\partial }_{\rho }r=1.\end{eqnarray} \tag{ A4 }$$

Relabeling ρ by r in Equation (A2), we have

$$\begin{eqnarray}&&{\rm{d}}{s}^{2}=\ -f(r){\rm{d}}{t}^{2}+\displaystyle \frac{{\rm{d}}{r}^{2}}{f(r)}+{R}^{2}(r){\rm{d}}{{\rm{\Omega }}}_{2}^{2}.\end{eqnarray} \tag{ A5 }$$

We now demonstrate how we can straightforwardly obtain the JMN-1 naked singularity metric in the form given above (Equation (A5)), where $-{\hat{g}}_{{tt}}{\hat{g}}_{{rr}}=1$. Its original form, as given in Equation (29) of Joshi et al. (2011), uses areal radial coordinates, x^μ = (t, R, ϑ, φ), in which $-{\hat{g}}_{{tt}}{\hat{g}}_{{RR}}\ne 1$. From the original form,

$$\begin{eqnarray*}&&{\rm{d}}{s}^{2}=-\left(1-{M}_{0}\right){\left(\displaystyle \frac{R}{{R}_{{\rm{b}}}}\right)}^{\tfrac{{M}_{0}}{1-{M}_{0}}}{\rm{d}}{t}^{2}+\displaystyle \frac{{\rm{d}}{R}^{2}}{\left(1-{M}_{0}\right)}+{R}^{2}\,{\rm{d}}{{\rm{\Omega }}}_{2}^{2},\end{eqnarray*}$$

where the compactness parameter M₀ is given in terms of the Arnowitt–Deser–Misner (ADM) mass M and the physical boundary or matching areal radius R_b is M₀ = 2M/R_b, we can write the desired equation (Equation (A4)) to put it in the g(r) = 1 coordinates (Equation (A5)) used in Table 1 as

$$\begin{eqnarray}&&{\left(\displaystyle \frac{R}{{R}_{{\rm{b}}}}\right)}^{\alpha }{\rm{d}}R={\rm{d}}r,\end{eqnarray} \tag{ A6 }$$

where α = M₀/(2 − 2M₀), not to be confused with the lapse function, which yields the solution

$$\begin{eqnarray}&&\displaystyle \frac{{R}_{{\rm{b}}}}{\alpha +1}{\left(\displaystyle \frac{R}{{R}_{{\rm{b}}}}\right)}^{\alpha +1}=r+k.\end{eqnarray} \tag{ A7 }$$

The integration constant k above can be set to zero so that at R = 0, we also have r = 0. The above equation is easily invertible as

$$\begin{eqnarray}&&R(r)={R}_{{\rm{b}}}{\left[\left(\displaystyle \frac{\alpha +1}{{R}_{{\rm{b}}}}\right)r\right]}^{\tfrac{1}{\alpha +1}}.\end{eqnarray} \tag{ A8 }$$

Now, with M₀ = 2α/(1 + 2α) and

$$\begin{eqnarray}&&{\partial }_{r}R(r)={\left[\left(\displaystyle \frac{\alpha +1}{{R}_{{\rm{b}}}}\right)r\right]}^{\tfrac{-\alpha }{\alpha +1}},\end{eqnarray} \tag{ A9 }$$

we can rewrite the metric given in Equation (A6) as

$$\begin{eqnarray}\begin{array}{rcl}{\rm{d}}{s}^{2} & = & -\left(\displaystyle \frac{1}{1+2\alpha }\right){\left[\left(\displaystyle \frac{\alpha +1}{{R}_{{\rm{b}}}}\right)r\right]}^{\tfrac{2\alpha }{\alpha +1}}{\rm{d}}{t}^{2}\\ & & +\,\left(1+2\alpha \right){\left[\left(\displaystyle \frac{\alpha +1}{{R}_{{\rm{b}}}}\right)r\right]}^{\tfrac{-2\alpha }{\alpha +1}}{\rm{d}}{r}^{2}\\ & & +\,{R}_{{\rm{b}}}^{2}{\left[\left(\displaystyle \frac{\alpha +1}{{R}_{{\rm{b}}}}\right)r\right]}^{\tfrac{2}{\alpha +1}}\,{\rm{d}}{{\rm{\Omega }}}_{2}^{2}.\end{array}\end{eqnarray} \tag{ A10 }$$

Introducing σ and r_b,

$$\begin{eqnarray}\begin{array}{rcl}\sigma & = & \displaystyle \frac{\alpha }{\alpha +1}=\displaystyle \frac{M}{{R}_{{\rm{b}}}-M},\\ {r}_{{\rm{b}}} & = & (1-\sigma ){R}_{{\rm{b}}}=\left(\displaystyle \frac{{R}_{{\rm{b}}}-2M}{{R}_{{\rm{b}}}-M}\right){R}_{{\rm{b}}},\end{array}\end{eqnarray} \tag{ A11 }$$

so that α = σ/(1 − σ), and we can simplify Equation (A10) and instead write

$$\begin{eqnarray}\begin{array}{rcl}{\rm{d}}{s}^{2} & = & -\left(1-\displaystyle \frac{2M}{{R}_{{\rm{b}}}}\right){\left(\displaystyle \frac{r}{{r}_{{\rm{b}}}}\right)}^{2\sigma }{\rm{d}}{t}^{2}\\ & & +\,{\left(1-\displaystyle \frac{2M}{{R}_{{\rm{b}}}}\right)}^{-1}{\left(\displaystyle \frac{r}{{r}_{{\rm{b}}}}\right)}^{-2\sigma }{\rm{d}}{r}^{2}\\ & & +\,{R}_{{\rm{b}}}^{2}{\left(\displaystyle \frac{r}{{r}_{{\rm{b}}}}\right)}^{2-2\sigma }{\rm{d}}{{\rm{\Omega }}}_{2}^{2}.\end{array}\end{eqnarray} \tag{ A12 }$$

This is the form of the JMN-1 metric reported in Table 1. We note that while inverting r(R), given in Equation (A7), to R(r), given in Equation (A8), was trivial in this case, this step might not be analytically feasible in general (assuming Equation (A4) admits a closed-form solution in the first place). Therefore, it is convenient to use three metric functions, f, g, and R, in general, as in Equation (1).

Since we are focused here on asymptotically flat spacetimes (see, e.g., Adamo et al. 2009), let us note first that for the seed metric (Equation (1)) to be asymptotically flat, its metric functions must become R(r) = r, f(r) = 1, and g(r) = 1 in the limit r → ∞ . Then, if we define ${\tilde{\eta }}_{\bar{\mu }\bar{\nu }}$ as

$$\begin{eqnarray}{\tilde{\eta }}_{\bar{\mu }\bar{\nu }}=\left[\begin{array}{cccc}-1 & 0 & 0 & 0\\ 0 & 1 & 0 & -a{\sin }^{2}\vartheta \\ 0 & 0 & A+{a}^{2}{\cos }^{2}\vartheta & 0\\ 0 & -a{\sin }^{2}\vartheta & 0 & (A+{a}^{2}){\sin }^{2}\vartheta \end{array}\right],\end{eqnarray} \tag{ B1 }$$

we will find that as r → ∞ , the Riemann tensor associated with $\tilde{\eta }$ vanishes identically for such asymptotically flat seeds. Thus, asymptotically, $\tilde{\eta }$ is the flat Minkowski metric η in disguise, i.e., ${\mathrm{lim}}_{r\to \infty }\tilde{\eta }=\eta$. We stress that the Riemann tensor does not generally vanish at finite coordinate radii r for the metric tensor $\tilde{\eta }$.

Now, the AA metric in the siKS coordinates (Equation (13)) can be expressed in terms of $\tilde{\eta }$ and the tangent to the ingoing PNC ℓ_- (Equation (10)) everywhere as

$$\begin{eqnarray}&&{g}_{\bar{\mu }\bar{\nu }}=\displaystyle \frac{X}{{\rm{\Sigma }}}\left[{\tilde{\eta }}_{\bar{\mu }\bar{\nu }}+\displaystyle \frac{2F}{{\rm{\Sigma }}}{({{\ell }}_{-})}_{\bar{\mu }}{({{\ell }}_{-})}_{\bar{\nu }}\right].\end{eqnarray} \tag{ B2 }$$

We note further that ℓ_- is null with respect to both g and $\tilde{\eta }$. We refer to the form of the metric above as the generalized (spherical ingoing) Kerr–Schild form of the AA metric (see Section 32.5 of Stephani et al. 2009 for related discussion). This can be seen to be conformally related to the classic Kerr–Schild form, as given in Equation (1.1) of Kerr & Schild (2009).

Therefore, for the AA metric itself to be asymptotically flat, we require (a) ${\mathrm{lim}}_{r\to \infty }(2F/{\rm{\Sigma }})=0$, and (b) ${\mathrm{lim}}_{r\to \infty }(X/{\rm{\Sigma }})=1$. The first of these conditions is met due to the properties of the asymptotically flat seed metric functions,

$$\begin{eqnarray}&&\mathop{\mathrm{lim}}\limits_{r\to \infty }\displaystyle \frac{2F}{{\rm{\Sigma }}}=\mathop{\mathrm{lim}}\limits_{r\to \infty }\displaystyle \frac{(1/\sqrt{g}-f/g){R}^{2}}{(f/g){R}^{2}+{a}^{2}{\cos }^{2}\vartheta }=0,\end{eqnarray} \tag{ B3 }$$

whereas the second is trivially met when X = Σ. Thus, asymptotically flat seed metrics admit asymptotically flat stationary generalizations when X = Σ. Note that this latter condition (X = Σ) is a sufficient but not necessary condition. It is instructive to compare the asymptotically Minkowski metric in our siKS coordinates (Equation (B1)) with the analogous form of the Minkowski metric in the original (outgoing) Kerr–Schild coordinates for the Kerr metric, as given, e.g., in Equations (1.7) and (1.13) of Wiltshire et al. (2009).

The Lagrangian ${\mathscr{L}}$ describing a geodesic orbit x^μ (λ) is given as $2{\mathscr{L}}={u}_{\mu }{u}^{\mu }$, where ${u}^{\mu }={\rm{d}}{x}^{\mu }/{\rm{d}}\lambda ={\dot{x}}^{\mu }$ is the four-velocity along the geodesic, and λ is an affine parameter along it. Working now in BL coordinates (Equation (2)) due to the two Killing symmetries of the spacetime that are generated by T = ∂_t and Φ = ∂_φ , we can find momenta, ${p}_{\mu }={\partial }_{{\dot{x}}^{\mu }}{\mathscr{L}}={\partial }_{{u}^{\mu }}{\mathscr{L}}={u}_{\mu }$, that are conserved, corresponding to the two cyclic variables. These are identified as being the energy E = −u_μ T^μ and azimuthal angular momentum L = u_μ Φ^μ of the orbit, respectively, ⁸ which can be used to obtain

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{\dot{t}}{E} & = & -\displaystyle \frac{X}{{{\rm{\Sigma }}}^{2}\det [{g}_{t\varphi }]}\left[{\rm{\Pi }}-2{Fa}\xi \right]{\sin }^{2}\vartheta ,\\ \displaystyle \frac{\dot{\varphi }}{E} & = & -\displaystyle \frac{X}{{{\rm{\Sigma }}}^{2}\det [{g}_{t\varphi }]}\left[2{Fa}{\sin }^{2}\vartheta +({\rm{\Sigma }}-2F)\xi \right],\end{array}\end{eqnarray} \tag{ C1 }$$

where we have introduced the first impact parameter, ξ $:=$ L/E, and the determinant of the t φ − sector of the AA metric tensor in BL coordinates (Equation (2)), $\det [{g}_{t\varphi }]:= {g}_{{tt}}{g}_{\varphi \varphi }-{g}_{t\varphi }^{2}=-({X}^{2}/{{\rm{\Sigma }}}^{2}){\rm{\Delta }}{\sin }^{2}\vartheta$.

Now, with Equation (C1), we can write

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{2{\mathscr{L}}X}{{E}^{2}}{\rm{\Delta }} & = & \left[{X}^{2}\displaystyle \frac{{\dot{r}}^{2}}{{E}^{2}}-{\left({\rm{\Delta }}+2F-a\xi \right)}^{2}\right]\\ & & +\,{\rm{\Delta }}\left[{X}^{2}\displaystyle \frac{{\dot{\vartheta }}^{2}}{{E}^{2}}+{\left(a\sin \vartheta -\xi \csc \vartheta \right)}^{2}\right].\end{array}\end{eqnarray} \tag{ C2 }$$

It is easy to see then that the geodesic equation for null geodesics ($2{\mathscr{L}}=0$) is fully separable,

$$\begin{eqnarray}&&{X}^{2}\displaystyle \frac{{\dot{\vartheta }}^{2}}{{E}^{2}}={\eta }^{2}-{\left(a\sin \vartheta -\xi \csc \vartheta \right)}^{2}=: {{\rm{\Theta }}}_{0}(\vartheta ),\end{eqnarray} \tag{ C3 }$$

$$\begin{eqnarray}&&{X}^{2}\displaystyle \frac{{\dot{r}}^{2}}{{E}^{2}}={\left({\rm{\Delta }}+2F-a\xi \right)}^{2}-{\rm{\Delta }}{\eta }^{2}=: {{\mathscr{R}}}_{0}(r),\end{eqnarray} \tag{ C4 }$$

where we have introduced a separation constant η², which is related to the Carter constant C through η² = C/E². The Carter constant, in turn, can be used to demonstrate the existence of a Killing–Yano tensor and an associated hidden symmetry of the motion.

The fundamental PNCs of the AA spacetime consist of null geodesics that satisfy $\dot{\vartheta }=0$ and $\ddot{\vartheta }=0$ (see, e.g., Misner et al. 1973; Hioki & Miyamoto 2008; see also Section 2.3 of Adamo et al. 2009). This is equivalent to requiring that Θ₀ = 0 and ∂_ϑ Θ₀ = 0, which yields a solution η = 0 and $\xi =a{\sin }^{2}\vartheta$. On a related note, the Newman & Penrose (1962) complex null tetrad adapted to the outgoing ($\dot{r}\gt 0$) PNC for the AA metric can be found in Azreg-Aïnou (2014a).

If X(r, ϑ) is of the form X(r, ϑ) = X_r (r) + X_ϑ (ϑ), then the geodesic equation is separable even for non-null orbits,

$$\begin{eqnarray}\begin{array}{rcl}{X}^{2}\displaystyle \frac{{\dot{\vartheta }}^{2}}{{E}^{2}} & = & {\eta }^{2}-{\left(a\sin \vartheta -\xi \csc \vartheta \right)}^{2}+\displaystyle \frac{2{\mathscr{L}}}{{E}^{2}}{X}_{\vartheta }=: {\rm{\Theta }}(\vartheta ),\\ {X}^{2}\displaystyle \frac{{\dot{r}}^{2}}{{E}^{2}} & = & {\left({\rm{\Delta }}+2{mr}-a\xi \right)}^{2}-{\rm{\Delta }}{\eta }^{2}+\displaystyle \frac{2{\mathscr{L}}}{{E}^{2}}{\rm{\Delta }}{X}_{r}=: {\mathscr{R}}(r).\end{array}\end{eqnarray} \tag{ C5 }$$

For the special case where X = Σ, these are given simply as X_r (r) = A(r) and ${X}_{\vartheta }(\vartheta )={a}^{2}{\cos }^{2}\vartheta$.