Introduction

In recent years, there has been a great amount of interest in precisely controlling the electromagnetic response of artificial materials. By introducing subwavelength structural features, the permittivity and permeability tensors of the medium can be tuned to exhibit a wide range of interesting and useful phenomena, such as cloaking1,2,3,4,5,6,7, negative refraction1,8,9, and subwavelength microscopy with superlenses10,11,12,13.

Analogue spacetimes1,2,14,15,16,17,18,19 use optical materials to implement coordinate transformations between a physical space and a virtual “electromagnetic space”, via the formal equivalence between Maxwell’s equations in curved spacetime and those in flat spacetime within a corresponding bianisotropic medium19,20,21,22,23. This allows one to build optical analogues to gravitational systems24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44. In particular, there has been a fair amount of interest in reproducing the metrics of black holes45,46,47,48,49,50. The null geodesics and polarizations of light moving in the spacetime metric can be reproduced exactly within a fully bianisotropic material; if one simply wishes to reproduce the null geodesics of the metric, however, it is much simpler to use an appropriately designed gradient-index material that is easier to construct experimentally.

In this paper, we discuss the bianisotropic and gradient-index materials that imitate the exterior equatorial Kerr–Newman black hole solution. We first carry out the analysis for optical systems reproducing the null geodesics of the Schwarzschild black hole. We recover the familiar results for the permittivity and permeability tensors and scalar refractive index reproducing the metric in isotropic coordinates, as well as the permittivity and permeability tensors reproducing the metric in the Schwarzschild coordinates20,46,47,51. We then present the scalar index that reproduces the null geodesics for Schwarzschild coordinates, which, by comparison with the isotropic result, has the significant experimental benefit of remaining finite all the way to the horizon. We then carry out these same analyses for the equatorial Kerr–Newman metric in Boyer–Lindquist coordinates, reproducing the metric within a fully bianisotropic material23, and finding the scalar index required to reproduce the null geodesics. We use finite-difference frequency-domain simulations of systems that approximate the gradient-index solutions of the Schwarzschild and Kerr–Newman black holes with concentric circular shells of constant index, and use ray tracing to perform an analysis of the error sensitivity of such systems. These analyses demonstrate that these approximate gradient-index systems, which are far simpler to construct than true gradient-index systems or full bianisotropic media, can adequately reproduce null geodesics and are forgiving to fabrication and experimental error for reasonable geodesics. As such, they are practical tabletop analogues for charged and/or rotating black holes.

Results

Throughout this paper we use Gaussian Planck units, with c = ℏ = G = 4πϵ0 = 1. Greek indices range over temporal and spatial coordinates, e.g., μ = 0, …, 3, while Roman indices range over only spatial coordinates, e.g., i = 1, …, 3. We use uppercase Greek and Roman letters to indicate variables related to the optical system, while we use lowercase letters to indicate variables related to the spacetime metric that it is replicating. We refer to these respectively as “real space” and “spacetime” variables. We typically use hats to indicate the dimensionless versions of variables. When we map spacetime coordinates onto real space coordinates, we always do so by equating the dimensionless coordinates. Spacetime variables are dedimensionalized via multiplication by the appropriate power of the black hole mass M. Real space dimensionless variables are then dimensionalized by a convenient length scale for construction. Using this matching of coordinates allows one to more easily keep track of the relationship between real space coordinates and the spacetime coordinates they represent.

The Schwarzschild black hole

We will begin by studying the Schwarzschild black hole and various optical analogues thereof. The Schwarzschild metric describes the spacetime geometry of a static, uncharged black hole of mass M, and is given in dimensionless Schwarzschild coordinates \(\hat{s},\hat{t},\rho ,\theta ,\phi\) (related to the usual dimensionful quantities via \(s=M\hat{s}\), \(t=M\hat{t}\), r = Mρ) by ref. 52

$${\rm{d}}{\hat{s}}^{2}=-\left(1-\frac{2}{\rho }\right){\rm{d}}{\hat{t}}^{2}+{\left(1-\frac{2}{\rho }\right)}^{-1}{\rm{d}}{\rho }^{2}+{\rho }^{2}\left({\rm{d}}{\theta }^{2}+{\sin }^{2} \, \theta \,\,{\rm{d}}{\phi }^{2}\right)\ .$$
(1)

Making the coordinate transformation \(\rho =\tilde{\rho }\,{(1+\frac{1}{2\tilde{\rho }})}^{2}\), the Schwarzschild metric (1) can be written in the form52

$${\rm{d}}{\hat{s}}^{2}=-\frac{{\left(1-\frac{1}{2\tilde{\rho }}\right)}^{2}}{{\left(1+\frac{1}{2\tilde{\rho }}\right)}^{2}}{\rm{d}}{\hat{t}}^{2}+{\left(1+\frac{1}{2\tilde{\rho }}\right)}^{4}\left({\rm{d}}{\hat{x}}^{2}+{\rm{d}}{\hat{y}}^{2}+{\rm{d}}{\hat{z}}^{2}\right)\ ,$$
(2)

where the spacetime isotropic coordinates \(\left(\hat{x},\hat{y},\hat{z}\right)\) are related to the transformed Schwarzschild coordinates \(\left(\tilde{\rho },\theta ,\phi \right)\) via the transformation from Cartesian to spherical coordinates.

We first replicate the metric in isotropic coordinates, given in Eq. (2), in order to make contact with existing literature. As discussed in refs. 15,19, there is a formal equivalence between the equations of electrodynamics in a curved spacetime and those in flat space in a macroscopic medium. Specifically, the behavior of light in a curved spacetime background described by metric gμν is reproduced in flat space within an impedance-matched bianisotropic medium with permittivity ϵij, permeability μij, and magnetoelectric coupling αi given by

$${\epsilon }^{ij}={\mu }^{ij}=-\frac{\sqrt{-\det g}}{{g}_{00}\sqrt{\det \gamma }}{g}^{ij}\ ,\quad {\alpha }_{i}=\frac{{g}_{0i}}{{g}_{00}\sqrt{\det \gamma }}\ ,$$
(3)

where γij is the three-dimensional metric tensor of the real space coordinate system in which we construct the medium, onto which we map the spatial components gij. Here, g and γ denote the determinants of gμν and γij, respectively. The macroscopic fields D, H are related to the microscopic fields E, B via

$${\bf{D}}={\boldsymbol{\epsilon }}{\bf{E}}+{\boldsymbol{\alpha }}\;\times\, {\bf{H}}\ ,\quad {\bf{B}}={\boldsymbol{\mu }}{\bf{H}}-{\boldsymbol{\alpha }}\;\times\, {\bf{E}}\ .$$
(4)

As discussed in ref. 53, this choice of identification between the spacetime geometry and the electromagnetic analogue, elaborated first in ref. 19, is not unique, and cannot reproduce all measurable properties of light moving in the spacetime metric. However, it is sufficient to reproduce both the null geodesic trajectory and the polarizations of light moving along these geodesics, which makes analogues produced with this identification worthy subjects of study.

Using Eq. (3) to map the dimensionless spacetime isotropic coordinates \(\left(\hat{x},\hat{y},\hat{z}\right)\) onto the corresponding dimensionless real space Cartesian coordinates \(\left(\hat{X},\hat{Y},\hat{Z}\right)\) (and thus mapping the dimensionless spacetime isotropic radial coordinate \(\tilde{\rho }\) onto the dimensionless real space radial coordinate P), we find that the behavior of light in the Schwarzschild metric (2) is reproduced in flat space within a medium described by

$${\epsilon }^{ij}={\mu }^{ij}=\frac{{(2P+1)}^{3}}{4{P}^{2}(2P-1)}{{\mathbb{1}}}^{ij}\ ,\quad i,j\in \left\{\hat{X},\hat{Y},\hat{Z}\right\}.$$
(5)

In this case, the medium is isotropic, and the scalar index can be read off immediately from Eq. (5) as

$$n(P)=\frac{{(2P+1)}^{3}}{4{P}^{2}(2P-1)}\ .$$
(6)

Note that the results of Eqs. (5) and (6) are well-established in the literature20,46,47,51. Equation (6) has the benefit that there is a single scalar index that reproduces all null geodesics and the polarization of light moving along these geodesics, and as the material is isotropic (though still inhomogeneous) it is thus easier to construct experimentally. However, this refractive index diverges approaching the horizon, i.e., as \(P\to \frac{1}{2}\), so it is not useful for investigating geodesics in the vicinity of the horizon.

Another approach is to instead use the Schwarzschild coordinates (1), which produces an anisotropic medium distinct from Eq. (5). As before, we use Eq. (3) to map the dimensionless spacetime Schwarzschild coordinates \(\left(\rho ,\theta ,\phi \right)\) (now with the Schwarzschild radial coordinate ρ rather than the isotropic radial coordinate \(\tilde{\rho }\)) onto the corresponding dimensionless real space spherical coordinates \(\left(P,{\Theta} ,{\Phi} \right)\), yielding

$${\epsilon }^{ij}={\mu }^{ij}=\left(\begin{array}{ccc}1&0&0\\ 0&\frac{1}{P(P-2)}&0\\ 0&0&\frac{{\mathrm{csc}}^{2}\Theta}{P(P-2)}\end{array}\right)\ ,\quad i,j\in \left\{P,\Theta ,\Phi \right\}.$$
(7)

This same system is described in dimensionless real space Cartesian coordinates \(\left({\hat{X}}, {\hat{Y}}, {\hat{Z}}\right)\) by

$${\epsilon }^{ij}={\mu }^{ij} =\frac{1}{{P}^{2}(2-P)}\left(\begin{array}{ccc}2{\hat{X}}^{2}-{P}^{3}&2\hat{X}\hat{Y}&2\hat{X}\hat{Z}\\ 2\hat{X}\hat{Y}&2{\hat{Y}}^{2}-{P}^{3}&2\hat{Y}\hat{Z}\\ 2\hat{X}\hat{Z}&2\hat{Y}\hat{Z}&2{\hat{Z}}^{2}-{P}^{3}\end{array}\right)\\ =\frac{2{P}^{i}{P}^{j}-{P}^{3}{{\mathbb{1}}}^{ij}}{{P}^{2}(2-P)}\ ,\quad i,j\in \left\{\hat{X},\hat{Y},\hat{Z} \right\},$$
(8)

where \({P}^{i}=\left(\hat{X},\hat{Y},\hat{Z}\right)\) and the dimensionless real space Cartesian coordinates \(\left(\hat{X},\hat{Y},\hat{Z}\right)\) are related to the real space spherical coordinates \(\left(P,\Theta ,\Phi \right)\) in the usual way. This result matches those presented in refs. 46,47.

If we only wish to reproduce the trajectories of light in the Schwarzschild metric (and not the proper polarizations), then a radially varying scalar index n(P) is sufficient. In Schwarzschild coordinates, we will find that the radial profile depends on the initial conditions defining the geodesic. We consider null geodesics of the metric (1); all such geodesics are planar, and so the spherical symmetry allows us to take θ = π/2 without loss of generality. Such null geodesics of the Schwarzschild metric are parametrized by a conserved energy at infinity, \(\varepsilon =\left(1-\frac{2M}{r}\right)\frac{{\rm{d}}t}{{\rm{d}}\sigma }\), and the conserved angular momentum, \(\ell ={r}^{2}\frac{{\rm{d}}\phi }{{\rm{d}}\sigma }\), with σ the affine parameter of the geodesic. Dedimensionalizing these parameters via \(\ell =M\hat{\ell }\) and \(\sigma =M\hat{\sigma }\) (note that the energy is already dimensionless, \(\hat{\varepsilon }=\varepsilon\)), null geodesics satisfy the geodesic equation52

$$-{\hat{\varepsilon }}^{2}+{\left(\frac{{\rm{d}}\rho }{{\rm{d}}\hat{\sigma }}\right)}^{2}+\left(1-\frac{2}{\rho }\right)\frac{{\hat{\ell }}^{2}}{{\rho }^{2}}=0\ .$$
(9)

Combining this equation with \({\left(\frac{{\rm{d}}\phi }{{\rm{d}}\hat{\sigma }}\right)}^{2}={\hat{\ell }}^{2}/{\rho }^{4}\) yields

$$\frac{{\rm{d}}\phi }{{\rm{d}}\rho }=\pm {\left(\frac{{\hat{\varepsilon }}^{2}}{{\hat{\ell }}^{2}}{\rho }^{4}-{\rho }^{2}+2\rho \right)}^{-1/2}\ .$$
(10)

We then make use of the spacetime impact parameter \(\hat{b}(\rho )=\rho \sin \beta\), where β is defined by the relation

$$\rho \frac{{\rm{d}}\phi }{{\rm{d}}\rho }=-\tan \beta \ .$$
(11)

Plugging this relation into Eq. (10) and making the sign choice consistent with our definition of β, we find that

$$\hat{b}(\rho )={\left(\hat{b}_{\infty }^{-2}+2{\rho }^{-3}\right)}^{-1/2}\ ,$$
(12)

where we have defined \(\hat{b}_\infty =\hat{\ell }/\hat{\varepsilon }\). Fermat’s principle relates the real space impact parameter and index of refraction by

$$n(P)\propto \hat{B}{(P)}^{-1}\ .$$
(13)

Equation (12) is then taken as input to Eq. (13) by equating the spacetime coordinates \(\left(\rho ,\phi \right)\) with the real space coordinates \(\left(P,\Phi \right)\), which also equates the dimensionless spacetime impact parameter \(\hat{b}\) with the dimensionless real space impact parameter \(\hat{B}\). This yields

$$n(P)\propto \sqrt{\hat{b}_{\infty }^{-2}+2{P}^{-3}}\ .$$
(14)

This solution has a number of noteworthy features. First, we reiterate that Eq. (14) only reproduces the geodesic trajectories of light moving in the Schwarzschild metric (1), but does not faithfully reproduce its polarizations. The radial profile depends on the initial condition \(\hat{b}_\infty\), which is related to initial angle β and initial radius P0 by

$$\hat{b}_{\infty }^{2}=\frac{{P}_{0}^{2}}{{\csc }^{2}{\beta }_{0}-2{P}_{0}^{-1}}\ .$$
(15)

This is somewhat inconvenient for experimental application, as it means that a different apparatus must be constructed for each family of geodesics; to address this, one could in principle construct a cylinder, where \(\hat{b}_\infty\) varies along the cylinder axis and each 2D slice recreates the corresponding family of null geodesics. A significant benefit of this coordinate system, however, is that n(P) approaches a finite value as P → 2 so long as \(\hat{b}_\infty \, {\ne} \, 0\), which means that geodesics can be studied in the vicinity of the event horizon, in contrast to the solution (6). The constant of proportionality in Eq. (14) allows us to tune the scalar index at the initial P0 to the most feasible value for construction.

Finally, note that we can relate conserved quantities ε, \(\ell\) in spacetime to \({\mathcal{E}},L\) in real space in the following way: in flat space, a photon with frequency f and wavelength λ has energy \({\mathcal{E}}=2\pi f\) and angular momentum L = 2πB/λ, with B the dimensionful real space impact parameter. Thus, \(L/{\mathcal{E}}R_{\mathrm{S}}=n\hat{B}/2={\mathrm{constant}}\). Setting n = 1 at P → ∞ in Eq. (14) and equating the real space and spacetime dimensionless impact parameters, \(\hat{B}=\hat{b}\), yields \(L/{\mathcal{E}}R_{\mathrm{S}}=\ell /\varepsilon r_{\mathrm{S}}\). Here, rS = 2M, and RS is the real space radius onto which rS is mapped.

The Kerr–Newman black hole

We now apply the same approaches to investigate optical analogues of the Kerr–Newman black hole, of which the Kerr, Reissner–Nordström, and Schwarzschild results are special cases. We will restrict our attention to equatorial null geodesics.

The Kerr–Newman metric describes the spacetime geometry surrounding a black hole of mass M, angular momentum per unit mass a = J/M, electric charge Q, and magnetic charge Qm. Dedimensionalizing the quantities via \(a=M\hat{a}\), \(Q=M\hat{Q}\), \(Q_{\rm{m}}=M\hat{Q}_{\rm{m}}\), the metric is given in dimensionless Boyer–Lindquist coordinates by52

$${\rm{d}}{\hat{s}}^{2}= \,\hat{\Sigma }\left(\frac{{\rm{d}}{\rho }^{2}}{\hat{\Delta }}+{\rm{d}}{\theta }^{2}\right)-\frac{\hat{\Delta }}{\hat{\Sigma }}{\left({\rm{d}}\hat{t}-\hat{a}\,{\sin }^{2}\,\theta \,{\rm{d}}\phi \right)}^{2}\\ +\;\frac{{\sin }^{2}\theta }{\hat{\Sigma }}{\left[\left({\rho }^{2}+{\hat{a}}^{2}\right){\rm{d}}\phi -\hat{a}\,{\rm{d}}\hat{t}\right]}^{2}\ ,$$
(16)

where

$$\begin{array}{lll}\hat{\Sigma }\,=\,{\rho }^{2}+{\hat{a}}^{2}{\cos }^{2}\theta \ ,\\ \hat{\Delta }\,=\,{\rho }^{2}-2\rho +{\hat{a}}^{2}+\rho _{Q}^{2}\ ,\\ \rho _{Q}^{2}\,=\,{\hat{Q}}^{2}+\hat{Q}_{{\rm{m}}}^{2}\ .\end{array}$$
(17)

Here, M is the total mass-equivalent, which contains contributions from the irreducible mass, the rotational energy, and the Coulomb energy of the black hole54.

After setting θ = π/2 and dθ = 0 to restrict to the equatorial case, we use Eq. (3) to map the dimensionless spacetime coordinates (ρ, θ, ϕ) onto the dimensionless real space coordinates (P, Θ, Φ), as before, yielding

$${\epsilon }^{ij}={\mu }^{ij}=\left(\begin{array}{ccc}\frac{\hat{\Delta }}{\hat{\Delta }-{\hat{a}}^{2}}&0&0\\ 0&\frac{1}{\hat{\Delta }-{\hat{a}}^{2}}&0\\ 0&0&\frac{1}{\hat{\Delta }}\end{array}\right)\ ,\quad {\alpha }_{i}=\left(\begin{array}{c}0\\ 0\\ \hat{a}\left(\frac{1}{\hat{\Delta }-{\hat{a}}^{2}}-\frac{1}{{P}^{2}}\right)\end{array}\right)\ ,$$
(18)

where \(\hat{\Delta }\) should now be interpreted as a function of P. The derivation is given in full in the “Methods” section. The equatorial geodesics and polarizations of the Kerr–Newman metric are exactly reproduced in flat space within a medium with these properties23. There is a subtlety here—although the radial and azimuthal components P, Θ appear to diverge at the ergosphere \(\hat{\Delta }={\hat{a}}^{2}\), this is a spurious divergence. As discussed in refs. 23,55,56, the physically relevant covariant quantity is the tensor χ defined therein, which relates the macroscopic and microscopic fields. This quantity diverges only at the horizon \(\hat{\Delta }=0\).

As before, we can also replicate equatorial null geodesics of the Kerr–Newman metric using only a scalar index. As in the Schwarzschild case, these geodesics are parametrized by the dimensionless conserved energy at infinity and conserved angular momentum, given in this case by

$$\begin{array}{lll}\hat{\varepsilon }\,=\,\left(1-\frac{2}{\rho }+\frac{\rho _{Q}^{2}}{{\rho }^{2}}\right)\frac{{\rm{d}}\hat{t}}{{\rm{d}}\hat{\sigma }}+\left(\frac{2\hat{a}}{\rho }-\frac{\rho_{Q}^{2}\hat{a}}{{\rho }^{2}}\right)\frac{{\rm{d}}\phi }{{\rm{d}}\hat{\sigma }}\ ,\\ \hat{\ell }\,=\,-\left(\frac{2\hat{a}}{\rho }-\frac{\rho_{Q}^{2}\hat{a}}{{\rho }^{2}}\right)\frac{{\rm{d}}\hat{t}}{{\rm{d}}\hat{\sigma }}+\left({\rho }^{2}+{\hat{a}}^{2}+\frac{2{\hat{a}}^{2}}{\rho }-\frac{\rho_{Q}^{2}{\hat{a}}^{2}}{{\rho }^{2}}\right)\frac{{\rm{d}}\phi }{{\rm{d}}\hat{\sigma }}\ .\end{array}$$
(19)

The geodesic equations describing the equatorial motion are

$$\begin{array}{rcl}\frac{{\rm{d}}\phi }{{\rm{d}}\hat{\sigma }}\,=\,\frac{1}{\hat{\Delta }}\left[\left(1-\frac{2}{\rho }+\frac{\rho_{Q}^{2}}{{\rho }^{2}}\right)\hat{\ell }+\left(\frac{2\hat{a}}{\rho }-\frac{\rho_ {Q}^{2}\hat{a}}{{\rho }^{2}}\right)\varepsilon \right]\ ,\\ {\left(\frac{{\rm{d}}\rho }{{\rm{d}}\hat{\sigma }}\right)}^{2}\,=\,\frac{\left[{\left({\rho }^{2}+{\hat{a}}^{2}\right)}^{2}-{\hat{a}}^{2}\hat{\Delta }\right]{\hat{\ell }}^{2}}{{\rho }^{4}}\left(\hat{b}_{\infty }^{-1}-{\hat{V}}_{+}\right)\left(\hat{b}_{\infty }^{-1}-{\hat{V}}_{-}\right)\ ,\end{array}$$
(20)

where

$${\hat{V}}_{\pm }=\frac{\hat{a}\left(2\rho -\rho_{Q}^{2}\right)\pm {\rm{sgn}}(\hat{\ell }){\rho }^{2}\sqrt{\hat{\Delta }}}{{\left({\rho }^{2}+{\hat{a}}^{2}\right)}^{2}-{\hat{a}}^{2}\hat{\Delta }}\ .$$
(21)

As before, we find the impact parameter \(\hat{b}(\rho )=\rho \sin \beta\) by plugging Eq. (11) into \(\frac{{\rm{d}}\phi }{{\rm{d}}\rho }=\frac{{\rm{d}}\phi /{\rm{d}}\hat{\sigma }}{{\rm{d}}\rho /{\rm{d}}\hat{\sigma }}\), which yields

$$\hat{b}(\rho )=\frac{{\rho }^{2}\left[\left(\hat{\Delta }-{\hat{a}}^{2}\right)+\left(2\rho -\rho_{Q}^{2}\right)\hat{a}\hat{b}_{\infty }^{-1}\right]}{\sqrt{{\rho }^{2}{\left[\left(\hat{\Delta }-{\hat{a}}^{2}\right)+\left(2\rho -\rho_{Q}^{2}\right)\hat{a}\hat{b}_{\infty }^{-1}\right]}^{2}+{\hat{\Delta }}^{2}\left[{\left({\rho }^{2}+{\hat{a}}^{2}\right)}^{2}-{\hat{a}}^{2}\hat{\Delta }\right]\left(\hat{b}_{\infty }^{-1}-{\hat{V}}_{+}\right)\left(\hat{b}_{\infty }^{-1}-{\hat{V}}_{-}\right)}}\ .$$
(22)

We proceed by equating spacetime coordinates (ρ, θ, ϕ) and real space coordinates (P, Θ, Φ), which sets \(\hat{B}(P)=\hat{b}(\rho =P)\). The scalar index for an optical Kerr–Newman black hole is again given by

$$n(P)\propto \hat{B}{(P)}^{-1}\ .$$
(23)

An optical system with this scalar index reproduces the equatorial null geodesic trajectories of the Kerr–Newman metric.

Unlike the Schwarzschild case, this scalar index is not always sufficient to fully reproduce the given family of Kerr–Newman geodesics. This can be seen immediately by noting that initially counter-rotating geodesics (those with \(\hat{\ell }\) of opposite sign to \(\hat{a}\)) must turn around and become co-rotating before crossing into the ergosphere; such a reversal of the sign of \(\frac{{\rm{d}}\phi }{{\rm{d}}\rho }\) is not possible with a finite (and positive) scalar index. This shortcoming manifests itself as a divergence of the scalar index; the outermost divergence occurs at radius

$${P}_{* }=1-\hat{a}\hat{b}_{\infty }^{-1}+\sqrt{\left(1-\hat{a}\hat{b}_{\infty }^{-1}\right)\left(1-\hat{a}\hat{b}_{\infty }^{-1}-\rho_ {Q}^{2}\right)}\ .$$
(24)

This is a removable pole in the Schwarzschild and Reissner–Nordström cases. For rotating black holes, the divergence occurs at the point in the trajectory where the direction of rotation reverses, consistent with the above observation that a finite radially varying scalar index is insufficient to implement such a reversal. Thus, the pole only affects initially counter-rotating geodesics that enter the ergosphere.

Simulations of constructible optical black holes

Optical analogues to black holes are particularly useful if their constructions are realizable. In the following sections, we model optical black holes with radially varying scalar refractive indices n(P), as given by Eqs. (14) and (23). For Schwarzschild (and many Kerr–Newman) black holes, n(P) is maximal at the horizon. Because the impact parameter of light on the optical black hole must be less than or equal to the radius of the “edge” of the system, i.e., \(\hat{B}\le {P}_{0}\), it is found that n(P) ≤ c0n0P0, where \({c}_{0}=\sqrt{31/108}\approx 0.54\) and n0 = n(P0). Thus, the construction of an optical Schwarzschild black hole with n0 = 1 and moderate P0 ≤ 6 is plausible and achievable with indices of refraction in the range of ordinary materials, such as water, glass, and plastic. (As will be seen, many optical Kerr–Newman black holes are also constructible.) True gradient-index profiles of the form (14) could perhaps be achieved with metamaterials; however, it is not clear how easily realizable such systems are, so in this work we approximate the profiles with concentric annuli of constant scalar index.

For a system size in which the wavelength of the source light is much smaller than the gradient length scale of the scalar-index profile, i.e., \(\lambda \ll n/\left\Vert \nabla n\right\Vert\), a highly localized and highly directional light source, like a laser, would nearly approximate the geodesics of Eqs. (10) and (20). Simulating these trajectories amounts to ray tracing, which we pursue in the following section. Specifically, we investigate the number of annuli needed to sufficiently mimic the true scalar-index profile and explore the impact of imperfect construction and experimental error on the deviation of the ray trajectory from the geodesic. However, in the following section, we will first consider the case in which the source wavelength is similar to the size of the optical black hole, i.e., M/λ ~ O(10). This is done to demonstrate the strengths and limitations of this study’s approach, as well as to be consistent with previous studies such as refs. 28,29,31,32,38,41,46,47,57.

In this study, all optical black holes are modeled with dimensionless outer radius P0 = 6. The system comprises either 16 or 21 concentric annuli, with the number depending on acceptable annulus thickness (i.e., greater than the wavelength) and the minimum modeled radius \({P}_{\min }\). The innermost and outermost annuli each have half the width of each interior annulus. The scalar index of each annulus is uniform, so that the simulated n(P) profiles are piecewise functions, as shown in Figs. 1 and 2. The values of n for the innermost and outermost annuli are taken as n(P) at the minimum and maximum radii, respectively; the refractive index of an inner annulus is taken as the value of n(P) at its center.

Fig. 1: Scalar refractive index of simulated optical Schwarzschild black holes.
figure 1

Radial profiles of the scalar refractive index used for simulations of optical Schwarzschild black holes with impact parameters \(\hat{b}_\infty =\) a 2, b 3, c 4, and d 5. The outer radius is R0/M = 6 with M the black hole mass. Note the logarithmic scale of the vertical axis.

Full size image
Fig. 2: Scalar refractive index of simulated optical Kerr–Newman black holes.
figure 2

Radial profiles of the scalar refractive index n used for simulations of optical Kerr–Newman black holes: a maximally co-rotating (\(\hat{a}=1,\rho_Q=0\)); b maximally charged (\(\hat{a}=0,\rho_Q=1\)); c charged and co-rotating \((\hat{a}=2/5,\rho_Q=4/5)\); and d charged and counter-rotating \((\hat{a}=-2/5,\rho_Q=4/5)\). The impact parameter is \(\hat{b}_\infty =3\), and outer radius is R0/M = 6, with M the black hole mass. Note the logarithmic scales and limits of the vertical axes.

Full size image

It is important to note here that this geometry was chosen for simplicity in the finite-difference frequency-domain simulations of the next section, in which dimensions are constrained by the wavelength. For consistency, the same geometry is scaled linearly for the ray tracing analyses that follow. In practice, non-uniform annulus thicknesses could be used to minimize steps Δn in regions of high \(\frac{{\rm{d}}n}{{\rm{d}}R}\) and to reduce light scattering at each boundary, but this is left for future work.

Finite-difference frequency-domain simulations

In this section, the trajectory of light around an optical black hole is modeled using a finite-difference frequency-domain (FDFD) solver58,59 of Maxwell’s equations. Simulation details are provided in the “Methods” section. Figure 1 shows the profiles n(P) used when modeling light incident on an optical Schwarzschild black hole with four different impact parameters, \(\hat{b}_\infty =2\), 3, 4, and 5; the resulting FDFD simulations are shown in Fig. 3, with the wavelength of light λ = 0.5 μm and optical Schwarzschild radius RS = 2M = 5 μm.

Fig. 3: Numerical simulations of optical Schwarzschild black holes.
figure 3

Finite-difference frequency-domain simulations of light incident on an optical Schwarzschild black hole with impact parameters \(\hat{b}_\infty =\) a 2, b 3, c 4, and d 5. True geodesics are plotted as thick lines. Poynting vectors (white arrows) are scaled by ∝1/R. Edge radii and Schwarzschild radii (RS) are solid circles. Each interior annulus's edge is marked. Color scales for the normalized electric field amplitude \(\left|E\right|/\!\max \left|E\right|\) are the same for each subplot.

Full size image

Consider the simulation shown in Fig. 3a, for which the dedimensionalized impact parameter at infinity is \(\hat{b}_\infty =2\). The peak of the electric field normalized to its maximum, \(\left|E\right|/\max \left|E\right|\), follows the path of the geodesic quite closely. Here, \(\left|E\right|=\sqrt{E{E}^{* }}\), with E* the complex conjugate of E. Time-averaged Poynting vectors are calculated as \(Re\left(\frac{1}{2}{\bf{E}}\times {{\bf{H}}}^{* }\right)\) and scaled by ∝ 1/R in the figures. Those with largest magnitude point mostly along the geodesic, and much of the energy flux is directed into the optical black hole. The same spatial trend is seen in Fig. 3b, for which \(\hat{b}_\infty =3\). Note in Fig. 1 how the profile of scalar index n(P) increases in amplitude as the initial impact parameter increases, in order to further bend light toward the horizon.

For \(\hat{b}_\infty =4\) and 5, seen in Fig. 3c, d, respectively, the brightest regions of \(\left|E\right|/\!\max \left|E\right|\) (and longest Poynting vectors) predominantly follow the geodesics. This is actually seen more clearly in the energy contained in the electric field (\(\propto\!\!{\left|E\right|}^{2}\)); however, only the electric field amplitude is shown here for better visualization of both small and large amplitude features. Agreement between the simulated light path and actual geodesic is expected to improve as the wavelength and beam width decrease relative to the size of the optical black hole, as described in the following section.

Another interesting effect is observed in Fig. 3d: the FDFD simulation does not show light following the geodesic all the way to the horizon. Instead, light begins to orbit at the photon sphere, R = 3M = 7.5 μm. This results because the impact parameter is nearly equal to that at which light becomes trapped, \(\hat{b}_\infty =3\sqrt{3}\approx 5.2\). Only traces of the photon “ring” are resolved in Fig. 3d. Higher fidelity simulations, with the optical black hole comprising many more annuli, would likely be required to properly simulate and study this phenomenon. This is left to future work.

Several optical Kerr–Newman black holes are also simulated, with profiles n(P) in Fig. 2 corresponding to the FDFD solutions in Fig. 4. For each case, the impact parameter is \(\hat{b}_\infty =3\), and the outer edge of the optical black hole is again at radius P0 = 6. These can be compared to the optical Schwarzschild black hole of Fig. 3b. The innermost modeled radius varies for each simulation, depending on whether n(P) diverges outside of the horizon. Each simulation is described in detail below.

Fig. 4: Numerical simulations of optical Kerr–Newman black holes.
figure 4

Finite-difference frequency-domain simulations of light incident \((\hat{b}_\infty =3)\) on four optical Kerr–Newman black holes, which are a maximally co-rotating (\(\hat{a}=1,\rho_Q=0\)); b maximally charged (\(\hat{a}=0,\rho_Q=1\)); c charged and co-rotating (\(\hat{a}=2/5,\rho_Q=4/5\)); and d charged and counter-rotating (\(\hat{a}=-2/5,\rho_Q=4/5\)). True geodesics are plotted as thick lines. Poynting vectors (white arrows) are scaled by ∝1/R. Maximum and minimum radii are solid circles; radii of interest, such as the horizon radius (Rh) or Schwarzschild radius (RS), are also plotted and labeled. Each annulus's edge is marked. Color scales for the normalized electric field amplitude \(\left|E\right|/\max \left|E\right|\) are the same for each subplot.

Full size image

An extremal Kerr black hole (\(\hat{a}=1,\rho_Q=0\)), with beam trajectory co-rotating with the black hole spin, is shown in Fig. 4a. Here, n(P) diverges at P* = 4/3; however, the true geodesic escapes the black hole with \(\frac{{\rm{d}}r}{{\rm{d}}\sigma }=0\) at P = 2. Therefore, though the horizon is at Ph = 1, the system is modeled with innermost radius P = 1.4. Comparing to the Schwarzschild case in Fig. 3b, we see that light is “dragged” further around the co-rotating black hole, as expected. In addition, more light “escapes”, although not all is directed along the geodesic. Inevitably, some energy flux is directed into the optical black hole, as indicated by the Poynting vectors; this is partially due to the finite width of the beam, and partially to the discrete annular approximation of the true gradient-index profile.

In contrast to Fig. 4a, an extremal Reissner–Nordström black hole (\(\hat{a}=0,\rho_Q=1\)) is simulated, with FDFD results depicted in Fig. 4b. Here, the optical black hole could be modeled completely to the horizon at Rh = 2.5 μm. In general, the peak of \(\left|E\right|/{\max} \left|E\right|\) follows the geodesic to the horizon. Little difference is seen when comparing to the Schwarzschild case of Fig. 3b, except that light now propagates within RS = 5 μm.

Two non-extremal Kerr–Newman black holes, with the same charge (ρQ = 4/5) but opposite spins (\(\hat{a}=\pm \!2/5\)), are also simulated and shown in Fig. 4c, d. The co-rotating black hole is modeled to the horizon at \(P_{\rm{h}}=1+\sqrt{1/5}\approx 1.45\). Compared to the extremal Kerr black hole in Fig. 4a, light is not dragged as far around the black hole.

For the counter-rotating Kerr–Newman black hole, the profile n(P) diverges at \(P_*\) ≈ 1.88; this is the radius at which the geodesic begins co-rotating with the black hole spin, i.e., where \(\frac{{\rm{d}}\phi }{{\rm{d}}\sigma }=0\). Thus, the system is modeled only to P = 1.96, where n(P = 1.96) ≈ 6. Comparing the co- and counter-rotating black holes, we see that light travels further in the Φ-direction for the former system, as expected.

As described in this section, a variety of optical Schwarzschild and Kerr–Newman black holes can be constructed feasibly with low indices of refraction. If such systems are built at a small scale, FDFD simulations show that the trajectories of light behave as expected, mostly following the true geodesics despite the discrete approximation to the proper gradient-index profile. The benefits of building larger systems are discussed in the next section.

Ray tracing calculations

In principle, the optical black holes of the previous section could be scaled in size from micrometer to centimeter or larger. This would simplify not only the construction of the optical black hole but also the calculation of light propagation, since the wavelength and width of the light source would be much smaller than the system size and related gradient length scales. The minimum gradient scale length of the scalar-index profiles in Figs. 1 and 2 is \(n/\left\Vert \nabla n\right\Vert \sim 0.6 \, \upmu \mathrm{m}\), so \(\lambda \, < \, n/\left\Vert \nabla n\right\Vert\) is valid for the above FDFD simulations. If visible light, λ ≈ 0.3 − 0.7 μm, is used, scaling the system size by even a factor of 103, i.e., from micrometer to millimeter, or greater would be appropriate for the validity of the ray tracing approximation made in this section.

It is of interest to calculate the deviation of a ray trajectory around the optical black hole from the true geodesic. These deviations could occur for a number of reasons: for instance, the discretization of n(R) due to the finite number of annuli; manufacturing error, leading to an offset Δn of the desired scalar index; or experimental error, resulting in a deviation ΔB0 from the desired initial impact parameter B0. We explore the impacts of these below for light incident on an optical Schwarzschild black hole with outer radius P0 = 6.

First, we investigate the number of annuli (with uniform thicknesses) needed to sufficiently approximate the scalar-index profile for a range of initial impact parameters. We define our performance metric as the deviation of the ray trajectory from the geodesic, quantified by the difference in azimuthal angle ΔΦ = Φray − Φgeo. Note that this performance metric is design-specific and does not account for scattering, whereas the semi-classical calculations of refs. 29,32 do. However, as there is no analytic solution to the wave equation for the system under consideration, the pursuit of a more appropriate metric is left to future work. Here, we are concerned with the deviation at the horizon. This value is shown in Fig. 5 for \(\hat{b}_\infty \in [0,5]\) and the number of annuli ranging 1–50. We see that only 25 annuli are needed to reproduce trajectories with \(\hat{b}_\infty \le 3\) to within ΔΦ = 3∘. As expected, ΔΦ increases rapidly for large \(\hat{b}_\infty\) and fixed annulus number. However, even the trajectory with \(\hat{b}_\infty =5\) can achieve ΔΦ ≤ 3∘ with 1000 annuli.

Fig. 5: Impact of annulus number on ray trajectories.
figure 5

The angular deviation of the ray trajectory (Φray) from the geodesic (Φgeo) at the horizon for an optical Schwarzschild black hole with outer radius P0 = 6, as a function of the initial impact parameter \(\hat{b}_\infty\) and number of annuli used in the construction.

Full size image

Next, we consider the scenario in which the scalar-index profile is imperfect, offset by a constant Δn due to some manufacturing error. We choose a specific trajectory with impact parameter \(\hat{b}_\infty =3\) to connect with the FDFD simulations of the previous section. The range spans Δn ∈ [0, 0.5] in Fig. 6; this is a significant percent change compared to profile b in Fig. 1. In Fig. 6a, we see that the ray trajectory skews radially outward as Δn increases. We are again interested in the deviation of the ray trajectory from the geodesic, ΔΦ  = Φray − Φgeo, shown in Fig. 6b as a function of radius P = R/M. Most trajectories follow the geodesic closely, within ΔΦ ≤ 2∘, for P > 3; however, within P < 3, ΔΦ grows rapidly. The small gray region, near P ≈ 2 and Δn ≈ 0.5, indicates that the ray trajectory escapes the black hole, so that ΔΦ diverges. In this case, if errors of ΔΦ ≤ 5∘ were allowable, then n(R) must be constrained with Δn ≤ 0.1.

Fig. 6: Effects of construction and experimenter errors on ray trajectories.
figure 6

a, b A uniform offset Δn from the true scalar refractive index profile n(R). c, d A deviation ΔB0 from the desired impact parameter B0. a, c Ray trajectories in real space, computed from Eq. (25), compared to the true geodesic (gray dashed). b, d Angular deviation of the ray trajectory (Φray) from the geodesic (Φgeo) versus radius R normalized to the black hole mass M. Both scans use the optical Schwarzschild black hole with \(\hat{b}_\infty =3\). The scale of the color bar of subplot a is the same as the scale of the vertical axis of subplot b; the same is true for subplots c and d.

Full size image

In addition, a scan in initial impact parameter is performed to assess how experimental error would affect the ray trajectory. The ratio ΔB0/B0 is varied within  ±10%, with results shown in Fig. 6. The ray trajectories (Fig. 6c) vary as expected: as \(\left|\Delta {B}_{0}\right|\) increases, the ray path moves farther from the true geodesic, but keeps the same general shape. Again, the deviation in azimuthal angle is shown in Fig. 6d. For large \(\left|\Delta {B}_{0}\right|/{B}_{0}\), ΔΦ increases rapidly as the trajectory approaches the horizon. The deviation can be as large as ΔΦ = 30∘ at P = 2 when ΔB0/B0 ≈ 10%. Interestingly, the contours of ΔΦ versus P and ΔB0/B0 are not symmetric about ΔB0/B0 = 0 in Fig. 6d. This results from the discretization of n(R). Therefore, if the number of annuli cannot be increased, it could actually be beneficial to purposefully shift the impact parameter (ΔB0/B0 < 0, in this case) to better match the light trajectory with the true geodesic.

Discussion

The application of analogue spacetimes to the study of general relativity has seen a resurgence in theory, simulation, and experiment in the past two decades. Many recent works have focused on optical analogues to static, uncharged (Schwarzschild) black holes in an isotropic coordinate system. In this paper, we have calculated the dielectric permittivity and permeability tensors ϵ, μ that reproduce the equatorial null geodesics and polarizations of light moving in the metric of spinning, charged (Kerr–Newman) black holes. Furthermore, we have conceived, for the first time, a gradient-index material that exactly reproduces families of equatorial Kerr–Newman null geodesics in almost all cases. Importantly, the radial profile of the scalar refractive index n(R) is finite along the entire trajectory (even to the horizon, if applicable), except at the point of rotation reversal for initially counter-rotating null geodesics. Values of n ≲ 6 can be achieved for many trajectories of interest, meaning that such gradient-index optical analogues could be constructed with conventional materials and metamaterials.

Simulations of a variety of optical black holes were performed, each with n(R) approximated by concentric circular annuli of constant scalar index. First, a FDFD solver of Maxwell’s equations was used to simulate the path of light incident on a Schwarzschild black hole with varying impact parameter \(\hat{b}_\infty =b_\infty /\!M\). Good agreement was observed between the light trajectory (indicated by maximum values of the electric field and Poynting vectors) and geodesic for low impact parameters \(\hat{b}_\infty\) = 2–3, but the discrepancy grew for \(\hat{b}_\infty\) = 4–5. Interestingly, for \(\hat{b}_\infty =5\), some features of light orbiting at the photon sphere were observed. Utilizing the same FDFD framework, several optical Kerr–Newman black holes were simulated: extremal Kerr, extremal Reissner–Nordström, and non-extremal Kerr–Newman with initially co- and counter-rotating trajectories. Each of these optical systems was simulated within the Schwarzschild radius, some even to the horizon. While there exist some discrepancies between the simulated light trajectories and true geodesics, the qualitative feature of light “dragged” in the direction of the black hole’s spin was observed. The three co-rotating cases require n ≲ 3, meaning that constructions of these optical Kerr–Newman black holes are feasible; the counter-rotating case requires n ≲ 6, which might be realized with more exotic materials like metamaterials.

Finally, we have investigated the number of annuli used in construction as well as the effects of fabrication and experimental errors on these optical black holes. The results demonstrate that with a modest number of annuli, the approximate gradient-index systems adequately reproduce null geodesics and are robust to small variations in refractive index and impact parameter. As these systems are far easier to manufacture than true gradient-index or bianisotropic media, they are thus practical tabletop analogues for equatorial Kerr–Newman black holes.

Methods

Numerical simulations

The trajectories of light around optical black holes are modeled using a FDFD solver of Maxwell’s equations58,59. The wavelength of light is chosen to be λ = 0.5 μm. The 2D simulation domain is modeled as a vacuum, with scalar properties ϵ = μ = n = 1 and size 60λ × 60λ; a perfectly matching layer of width λ/5 is applied at its boundary. A Gaussian beam of light is approximated as an array of line sources, each of width λ/25 = 20 nm and electric field amplitude calculated from a Gaussian envelope of the form \(\exp (-{(X-{B}_{0})}^{2}/2{\delta }^{2})\). Here, B0 is the dimensionful real space impact parameter at P0, and δ = λ/2 so that the beam satisfies the paraxial approximation60. The total width of the beam is truncated at 2λ by imposing two absorbing (ϵ = 1 − iπ) boundaries as vertically aligned “waveguides” of the light from the edge of the domain to the edge of the optical black hole. These restrict the beam to travel along a straight path in free space, as a directional light source would in the laboratory. Note that the factor of  −π is arbitrarily chosen for the imaginary (damping) component.

Each simulated optical black hole is centered in the domain, with the Schwarzschild radius always RS = 10λ  = 5 μm (M = 2.5 μm) and edge at R0 = 30λ = 15 μm. The Gaussian light source propagates in the vertical direction toward the black hole. For all simulations, the region within the minimum radius (oftentimes the horizon radius Rh) is modeled as a disc with dielectric permittivity ϵ = ϵin − iπ. Here, ϵin is the scalar permittivity (ϵ = n2) of the innermost annulus, and a factor of  −π is used for the imaginary (damping) component, as with the aforementioned “waveguides.”

Ray tracing algorithm

Consider an optical system consisting of N concentric annuli. Let the radii bounding each annulus i be Ri < Ri−1, so that the annuli are numbered 1, 2, …, N from the outside in, and the outer edge of the system is at R0. The scalar index of each annulus is n(Ri < R ≤ Ri−1) = ni, which monotonically increases from annulus 1 → N, so ni < ni+1. Let the scalar index for R > R0 be n0. For a light ray incident on annulus (i + 1) (propagating in the region Ri ≤ R ≤ Ri−1), let the impact parameter be \({B}_{i}={R}_{i}\sin {\Phi }_{i}\), where Φi is the azimuthal angle at which the ray intersects the annulus at Ri. Then, the azimuthal angle at which the light ray intersects the next annulus (i + 2) at Ri+1 is given by

$${\Phi }_{i+1}-{\Phi }_{i}=\arcsin \left(\frac{{B}_{i+1}}{{R}_{i+1}}\right)-\arcsin \left(\frac{{B}_{i+1}}{{R}_{i}}\right)\ ,$$
(25)

provided that Bi+1 ≤ Ri+1. Note that the impact parameter always satisfies \({n}_{i}{B}_{i}={\rm{constant}}\). Thus, given an optical system with a well-defined profile n(R) and an initial impact parameter B0, the trajectory of a light ray can be iteratively computed via Eq. (25) until the ray reaches its minimum radius. Note that only in-going trajectories are considered here, so light escaping the optical black hole is not modeled. Furthermore, it is assumed that all light is transmitted at each boundary; absorption and reflection are left for future work.

Derivation of Kerr–Newman analogue material properties

Here, we derive Eq. (18), beginning with Eqs. (3) and (16). Restricting our attention to equatorial geodesics, we have θ = π/2 and dθ = 0 (as the motion will always remain equatorial). With this, the metric simplifies to

$${\rm{d}}{\hat{s}}^{2}={\rho }^{2}\left(\frac{{\rm{d}}{\rho }^{2}}{\hat{\Delta }}+{\rm{d}}{\theta }^{2}\right)-\frac{\hat{\Delta }}{{\rho }^{2}}{\left({\rm{d}}\hat{t}-\hat{a}\,{\rm{d}}\phi \right)}^{2}+\frac{1}{{\rho }^{2}}{\left[\left({\rho }^{2}+{\hat{a}}^{2}\right){\rm{d}}\phi -\hat{a}\,{\rm{d}}\hat{t}\right]}^{2}\ .$$
(26)

Expanding this, we find

$${g}_{\mu \nu }=\left(\begin{array}{cccc}-\frac{\hat{\Delta }-{\hat{a}}^{2}}{{\rho }^{2}}&0&0&\frac{\hat{a}\left(\hat{\Delta }-{\rho }^{2}-{\hat{a}}^{2}\right)}{{\rho }^{2}}\\ 0&\frac{{\rho }^{2}}{\hat{\Delta }}&0&0\\ 0&0&{\rho }^{2}&0\\ \frac{\hat{a}\left(\hat{\Delta }-{\rho }^{2}-{\hat{a}}^{2}\right)}{{\rho }^{2}}&0&0&\frac{{\left({\rho }^{2}+{\hat{a}}^{2}\right)}^{2}-{\hat{a}}^{2}\hat{\Delta }}{{\rho }^{2}}\end{array}\right)\ ,$$
(27)

where μ, ν run over \(\hat{t},\rho ,\theta ,\phi\). This metric has inverse

$${g}^{\mu \nu }=\left(\begin{array}{cccc}\frac{{\hat{a}}^{2}\hat{\Delta }-{\left({\rho }^{2}+{\hat{a}}^{2}\right)}^{2}}{{\rho }^{2}\hat{\Delta }}&0&0&\frac{\hat{a}\left(\hat{\Delta }-{\rho }^{2}-{\hat{a}}^{2}\right)}{{\rho }^{2}\hat{\Delta }}\\ 0&\frac{\hat{\Delta }}{{\rho }^{2}}&0&0\\ 0&0&\frac{1}{{\rho }^{2}}&0\\ \frac{\hat{a}\left(\hat{\Delta }-{\rho }^{2}-{\hat{a}}^{2}\right)}{{\rho }^{2}\hat{\Delta }}&0&0&\frac{\hat{\Delta }-{\hat{a}}^{2}}{{\rho }^{2}\hat{\Delta }}\end{array}\right)$$
(28)

and determinant \(\det g=-{\rho }^{4}\). We map the curved spacetime coordinates \((\hat{t},\rho ,\theta ,\phi )\) onto the flat spacetime spherical coordinates \((\hat{T},P,\Theta ,\Phi )\), so the flat space coordinate metric is in this case

$${\gamma }_{ij}=\left(\begin{array}{ccc}1&0&0\\ 0&{P}^{2}&0\\ 0&0&{P}^{2}{\sin }^{2}\Theta \end{array}\right)$$
(29)

with determinant \(\det \gamma ={P}^{4}{\sin }^{2}\Theta\). Because we have restricted our attention to θ = π/2, we similarly have Θ = π/2, and so this simply becomes \(\det \gamma ={P}^{4}\). After this coordinate matching, we have

$${g}^{ij} = \left(\begin{array}{ccc}\frac{\hat{\Delta }}{{P}^{2}}&0&0\\ 0&\frac{1}{{P}^{2}}&0\\ 0&0&\frac{\hat{\Delta }-{\hat{a}}^{2}}{{P}^{2}\hat{\Delta }}\end{array}\right)\ ,\\ {g}_{00} = -\frac{\hat{\Delta }-{\hat{a}}^{2}}{{P}^{2}}\ ,\\ {g}_{0i} = \left(\begin{array}{ccc}0&0&\frac{\hat{a}\left(\hat{\Delta }-{P}^{2}-{\hat{a}}^{2}\right)}{{P}^{2}}\end{array}\right), \\ {{\mathrm{det}}\, {g} = {-{P}^{4}}},$$
(30)

where \(\hat{\Delta }\) is now interpreted as a function of P, as opposed to ρ. Plugging these values into Eq. (3), we arrive at

$${\epsilon }^{ij}={\mu }^{ij}=\left(\begin{array}{ccc}\frac{\hat{\Delta }}{\hat{\Delta }-{\hat{a}}^{2}}&0&0\\ 0&\frac{1}{\hat{\Delta }-{\hat{a}}^{2}}&0\\ 0&0&\frac{1}{\hat{\Delta }}\end{array}\right)\ ,\quad {\alpha }_{i}=\left(\begin{array}{c}0\\ 0\\ \hat{a}\left(\frac{1}{\hat{\Delta }-{\hat{a}}^{2}}-\frac{1}{{P}^{2}}\right)\end{array}\right)\ .$$
(31)