Optical analogues to the equatorial Kerr–Newman black hole

Optical analogues to black holes allow the investigation of general relativity in a laboratory setting. Previous works have considered analogues to Schwarzschild black holes in an isotropic coordinate system; the major drawback is that required material properties diverge at the horizon. We present the dielectric permittivity and permeability tensors that exactly reproduce the equatorial Kerr–Newman metric, as well as the gradient-index material that reproduces equatorial Kerr–Newman null geodesics. Importantly, the radial profile of the scalar refractive index is finite along all trajectories except at the point of rotation reversal for counter-rotating geodesics. Construction of these analogues is feasible with available ordinary materials. A finite-difference frequency-domain solver of Maxwell’s equations is used to simulate light trajectories around a variety of Kerr–Newman black holes. For reasonably sized experimental systems, ray tracing confirms that null geodesics can be well-approximated in the lab, even when allowing for imperfect construction and experimental error. 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 cloaking [1–7], negative refraction [1, 8, 9], and subwavelength microscopy with superlenses [10–13]. Analogue spacetimes [1, 2, 14–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 medium [19–23]. This allows one to build optical analogues to gravitational systems [24–44]. In particular, there has been a fair amount of interest in reproducing the metrics of black holes [45–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 †rating at mit.edu, ‡apturner at mit.edu

Analogue spacetimes [1,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 medium [19][20][21][22][23]. This allows one to build optical analogues to gravitational systems . In particular, there has been a fair amount of interest in reproducing the metrics of black holes [45][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 coordinates [20,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 material [23], 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 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ŝ,t, ρ, θ, φ (related to the usual dimensionful quantities via s = Mŝ, t = Mt, r = M ρ) by [52] Making the coordinate transformation ρ =ρ 1 + 1 2ρ 2 , the Schwarzschild metric (1) can be written in the form [52] where the spacetime isotropic coordinates (x,ŷ,ẑ) are related to the transformed Schwarzschild coordinates (ρ, θ, φ) 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 [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 impedancematched bianisotropic medium with permittivity ij , permeability µ ij , and magnetoelectric coupling α i given by 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 g ij . 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 As discussed in [53], this choice of identification between the spacetime geometry and the electromagnetic analogue, elaborated first in [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 (x,ŷ,ẑ) onto the corresponding dimensionless real space Cartesian coordinates X ,Ŷ ,Ẑ (and thus mapping the dimensionless spacetime isotropic radial coordinateρ 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 In this case, the medium is isotropic, and the scalar index can be read off immediately from Eq. (5) as .
Note that the results of Eqs. (5) and (6) are well-established in the literature [20,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 → 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 (ρ, θ, φ) (now with the Schwarzschild radial coordinate ρ rather than the isotropic radial coordinateρ) onto the corresponding dimensionless real space spherical coordinates (P, Θ, Φ), yielding This same system is described in dimensionless real space Cartesian coordinates X ,Ŷ ,Ẑ by where P i = X ,Ŷ ,Ẑ and the dimensionless real space Cartesian coordinates X ,Ŷ ,Ẑ are related to the real space spherical coordinates (P, Θ, Φ) in the usual way. This result matches those presented in [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, ε = 1 − 2M r dt dσ , and the conserved angular momentum, = r 2 dφ dσ , with σ the affine parameter of the geodesic. Dedimensionalizing these parameters via = Mˆ and σ = Mσ (note that the energy is already dimensionless,ε = ε), null geodesics satisfy the geodesic equation [52] −ε 2 + dρ dσ Combining this equation with dφ We then make use of the spacetime impact parameterb(ρ) = ρ sin β, where β is defined by the relation Plugging this relation into Eq. (10) and making the sign choice consistent with our definition of β, we find thatb where we have definedb ∞ =ˆ /ε. Fermat's principle relates the real space impact parameter and index of refraction by n(P ) ∝B(P ) −1 .
Equation (12) is then taken as input to Eq. (13) by equating the spacetime coordinates (ρ, φ) with the real space coordinates (P, Φ), which also equates the dimensionless spacetime impact parameterb with the dimensionless real space impact parameterB. This yields 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 conditionb ∞ , which is related to initial angle β and initial radius P 0 byb 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, wherê b ∞ 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 asb ∞ = 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 P 0 to the most feasible value for construction. Finally, note that we can relate conserved quantities ε, in spacetime to E, L in real space in the following way: in flat space, a photon with frequency f and wavelength λ has energy E = 2πf and angular momentum L = 2πB/λ, with B the dimensionful real space impact parameter. Thus, L/ER S = nB/2 = constant. Setting n = 1 at P → ∞ in Eq. (14) and equating the real space and spacetime dimensionless impact parameters,B =b, yields L/ER S = /εr S . Here, r S = 2M , and R S is the real space radius onto which r S 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 Q m . Dedimensionalizing the quantities via a = Mâ, Q = MQ, Q m = MQ m , the metric is given in dimensionless Boyer-Lindquist coordinates by [52] whereΣ = ρ 2 +â 2 cos 2 θ , Here, M is the total mass-equivalent, which contains contributions from the irreducible mass, the rotational energy, and the Coulomb energy of the black hole [54]. 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 where∆ 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 the''se properties [23]. There is a subtlety here-although the radial and azimuthal components P, Θ appear to diverge at the ergosphere∆ =â 2 , this is a spurious divergence. As discussed in [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∆ = 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 bŷ The geodesic equations describing the equatorial motion are whereV As before, we find the impact parameterb(ρ) = ρ sin β by plugging Eq. (11) into dφ dρ = dφ/dσ dρ/dσ , which yieldsb We proceed by equating spacetime coordinates (ρ, θ, φ) and real space coordinates (P, Θ, Φ), which setŝ B(P ) =b(ρ = P ). The scalar index for an optical Kerr-Newman black hole is again given by 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ˆ of opposite sign toâ) must turn around and become co-rotating before crossing into the ergosphere; such a reversal of the sign of dφ dρ 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 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.,B ≤ P 0 , it is found that n(P ) ≤ c 0 n 0 P 0 , where c 0 = 31/108 ≈ 0.54 and n 0 = n(P 0 ). Thus, the construction of an optical Schwarzschild black hole with n 0 = 1 and moderate P 0 ≤ 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., λ n/ ∇n , 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 [28,29,31,32,38,41,46,47,57].
In this study, all optical black holes are modeled with dimensionless outer radius P 0 = 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 Figures 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.
It is important to note here that this geometry was chosen for simplicity in the finite-difference frequencydomain 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 dn dR 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) solver [58,59] of Maxwell's equations. Simulation details are provided in the Methods section. Figure 1 shows the profiles n(R) used when modeling light incident on an optical Schwarzschild black hole with four different impact parameters,b ∞ = 2, 3, 4, and 5; the resulting FDFD simulations are shown in Figure 3, with the wavelength of light λ = 0.5 µm and optical Schwarzschild radius Consider the simulation shown in Figure 3a, for which the dedimensionalized impact parameter at infinity isb ∞ = 2. The peak of the electric field normalized to its maximum, |E|/max|E|, follows the path of the geodesic quite closely. Here, |E| = √ EE * , with E * the complex conjugate of E. Time-averaged Poynting vectors are calculated as Re 1 2 E × H * 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 Figure 3b, for whichb ∞ = 3. Note in Figure 1 how the profile of scalar index n(R) increases in amplitude as the initial impact parameter increases, in order to further bend light toward the horizon.
Forb ∞ = 4 and 5, seen in Figures 3c and 3d, respectively, the brightest regions of |E|/max|E| (and longest Poynting vectors) predominantly follow the geodesics. This is actually seen more clearly in the energy contained in the electric field (∝ |E| 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 Figure 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, b ∞ = 3 √ 3 ≈ 5.2. Only traces of the photon "ring" are resolved in Figure 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(R) in Figure 2 corresponding to the FDFD solutions in Figure 4. For each case, the impact parameter isb ∞ = 3, and the outer edge of the optical black hole is again at radius P 0 = 6. These can be compared to the optical Schwarzschild black hole of Figure 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.
An extremal Kerr black hole (â = 1, ρ Q = 0), with beam trajectory co-rotating with the black hole spin, is shown in Figure 4a. Here, n(P ) diverges at P * = 4/3; however, the true geodesic escapes the black hole with dr dσ = 0 at P = 2. Therefore, though the horizon is at P h = 1, the system is modeled with innermost radius P = 1.4. Comparing to the Schwarzschild case in Figure 3b, we see that light is "dragged" further around the co-rotating black hole, as expected. Additionally, 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 Figure 4a, an extremal Reissner-Nordström black hole (â = 0, ρ Q = 1) is simulated, with FDFD results depicted in Figure 4b. Here, the optical black hole could be modeled completely to the horizon at R h = 2.5 µm. In general, the peak of |E|/max|E| follows the geodesic to the horizon. Little difference is seen when comparing to the Schwarzschild case of Figure 3b, except that light now propagates within R S = 5 µm. Two non-extremal Kerr-Newman black holes, with the same charge (ρ Q = 4/5) but opposite spins (â = ±2/5), are also simulated and shown in Figures 4c and 4d. The co-rotating black hole is modeled to the horizon at P h = 1 + 1/5 ≈ 1.45. Compared to the extremal Kerr black hole in Figure 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 dφ dσ = 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

Ray tracing calculations
In principle, the optical black holes of the previous section could be scaled in size from µm to cm 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 Figures 1 and 2 is n/ ∇n ∼ 0.6 µm, so λ < n/ ∇n 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 10 3 , i.e., from µm to mm, 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 ∆B 0 from the desired initial impact parameter B 0 . We explore the impacts of these below for light incident on an optical Schwarzschild black hole with outer radius P 0 = 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 [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 Figure 5 forb ∞ ∈ [0, 5] and the number of annuli ranging 1 to 50. We see that only 25 annuli are needed to reproduce trajectories witĥ b ∞ ≤ 3 to within ∆Φ = 3°. As expected, ∆Φ increases rapidly for largeb ∞ and fixed annulus number. However, even the trajectory withb ∞ = 5 can achieve ∆Φ ≤ 3°with 1000 annuli.  Figure 5: Impact of annulus number on ray trajectories. The angular deviation of the ray trajectory (Φ ray ) from the geodesic (Φ geo ) at the horizon for an optical Schwarzschild black hole with outer radius P 0 = 6, as a function of the initial impact parameterb ∞ and number of annuli used in the construction.
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 parameterb ∞ = 3 to connect with the FDFD simulations of the previous section. The range spans ∆n ∈ [0, 0.5] in Figure 6; this is a significant percent change compared to profile b in Figure 1. In Figure 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 Figure 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 grey 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.
In addition, a scan in initial impact parameter is performed to assess how experimental error would affect the ray trajectory. The ratio ∆B 0 /B 0 is varied within ±10%, with results shown in Figure 6. The ray trajectories (Figure 6c) vary as expected: as |∆B 0 | 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 Figure 6d. For large |∆B 0 |/B 0 , ∆Φ increases rapidly as the trajectory approaches the horizon. The deviation can be as large as ∆Φ = 30°at P = 2 when ∆B 0 /B 0 ≈ 10%. Interestingly, the contours of ∆Φ versus P and ∆B 0 /B 0 are not symmetric about ∆B 0 /B 0 = 0 in Figure 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 (∆B 0 /B 0 < 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 finite-difference frequency-domain (FDFD) solver of Maxwell's equations was used to simulate the path of light incident on a Schwarzschild black hole with varying impact parameterb ∞ = b ∞ /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 parameterŝ b ∞ = 2-3, but the discrepancy grew forb ∞ = 4-5. Interestingly, forb ∞ = 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.

Numerical simulations
The trajectories of light around an optical black holes are modeled using a finite-difference frequency-domain (FDFD) solver of Maxwell's equations [58,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δ 2 ). Here, B 0 is the dimensionful real space impact parameter at P 0 , and δ = λ/2 so that the beam satisfies the paraxial approximation [60]. 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 R S = 10λ = 5 µm (M = 2.5 µm) and edge at R 0 = 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 R h ) is modeled as a disc with dielectric permittivity = in − iπ. Here, in is the scalar permittivity ( = n 2 ) 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 R i < R i−1 , so that the annuli are numbered 1, 2, . . . , N from the outside in, and the outer edge of the system is at R 0 . The scalar index of each annulus is n(R i < R ≤ R i−1 ) = n i , which monotonically increases from annulus 1 → N , so n i < n i+1 . Let the scalar index for R > R 0 be n 0 . For a light ray incident on annulus (i + 1) (propagating in the region R i ≤ R ≤ R i−1 ), let the impact parameter be B i = R i sin Φ i , where Φ i is the azimuthal angle at which the ray intersects the annulus at R i . Then, the azimuthal angle at which the light ray intersects the next annulus (i + 2) at R i+1 is given by provided that B i+1 ≤ R i+1 . Note that the impact parameter always satisfies n i B i = constant. Thus, given an optical system with a well-defined profile n(R) and an initial impact parameter B 0 , 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 Expanding this, we find where µ, ν run overt, ρ, θ, φ. This metric has inverse and determinant det g = −ρ 4 . We map the curved spacetime coordinates (t, ρ, θ, φ) onto the flat spacetime spherical coordinates (T , P, Θ, Φ), so the flat space coordinate metric is in this case with determinant det γ = P 4 sin 2 Θ. Because we have restricted our attention to θ = π/2, we similarly have Θ = π/2, and so this simply becomes det γ = P 4 . After this coordinate matching, we have where∆ is now interpreted as a function of P , as opposed to ρ. Plugging these values into Eq.

Data availability
Data is available from the corresponding author upon request.

Code availability
Code for the simulations shown here is available from the corresponding author upon request. The finitedifference frequency-domain solver used in this work is available at https://github.com/wsshin/maxwellfdfd.