Metamaterials mimic the black holes: the effects of charge and rotation on the optical properties

Motivated by investigation of black hole properties in the lab, some interesting subjects such as analogue gravity and transformation optics are generated. In this paper, we look for the analogies between the geometry of a gravitating system and the optical medium. In addition, we recognize that appropriate metamaterials can be used to mimic the propagation of light in the curved spacetimes and behave like black holes. The resemblance of metamaterials with Kerr and Reissner-Nordstr\"{o}m spacetimes is studied. At last, we compare the full-wave numerical calculation of light with its optical limit of geometry.


I. INTRODUCTION
One of the interesting subjects in gravitational physics is analogue gravity in the optical frameworks which is attributed to the pioneer work of Gordon [1]. Some of the Gordon's activities are devoted to describing dielectric media by an effective metric or strictly speaking, mimicking a gravitational field with a dielectric medium. With the advent of transformation optics and analogue gravity, one can create systems which closely resemble to general relativity objects such as black holes. Regarding the analogue gravity, one can investigate the resemblances of general relativistic gravitational fields with the equation of motion of other physical systems [2], such as mimicking some aspects of black holes with super fluids in the Bose-Einstein condensation [3]. Although such systems may not have exact analogues of any real spacetimes with corresponding metric tensors, transformation optics, where an optical material appears to perform a transformation of space, may use a formalism of general relativity to build gradient index static media in order to control the light paths in these media [4,5]. Theoretical concept of transformation optics is similar to equations of general relativity that describe how gravity warps space and time. However, instead of space and time, these equations show how light can be directed in a chosen manner, analogous to warp space [6]. In addition, transformation optics applies metamaterials to produce spatial variations, derived from the coordinate transformation. Therefore, it is possible to construct new artificial composite devices with a desired permittivity and permeability. Moreover, interactions of light and matter with spacetime, as predicted by the general relativity, can be studied by using a new type of artificial optical materials that feature extraordinary abilities to bend light. This research creates an interesting link between the newly emerging field of artificial optical metamaterials and that of general relativity, and therefore, it opens a new possibility to investigate various general relativity phenomena in a laboratory setting: chaotic motions observed in celestial objects that have been subjected to gravitational fields [7,8], reproduce an accurate laboratory model of the physical multiverse [8,9], optical analogues of black holes [10][11][12][13][14], Schwarzschild spacetime [15][16][17] and Hawking radiation [18].
The main purpose of this paper is to generate an optical analogue of the known gravitational black hole and its characteristic phenomena, such as the existence of an event horizon (a one-way membrane) and the bending of light rays. Optical materials can mimic the geometry of light near black holes, both theoretically and experimentally [6,10,14]. It has been shown that the permittivity and permeability tensors can be used instead of refractive index in order to build the equivalent optical medium mimicking black holes [6,17]. In Ref. [17] the propagation of light waves outside the Schwarzschild black hole were simulated and the results were compared with ray paths obtained from the Hamiltonian method. Interestingly, the phenomenon of "photon sphere" which is an important feature of the black hole systems, was observed.
Here as a first step, we generalize the results of Ref. [17] to the case of charged black hole. We consider a charged black hole and its optical simulation, and discuss the effect of electric charge on the photon sphere. Motivated by the lacking of a definitive answer on the existence of astrophysical non-rotating black holes, we are encouraged to look at Kerr spacetime and look for its closely optical simulation as the second step. In addition, we were interested in the effect of spin parameter on the photon sphere. The basic structure of this paper is twofold. First, we determine the permittivity and permeability mimicking the desired spacetime. Then, the Maxwell equations are numerically solved by using constitutive relations. Second, the Lagrangian for each spacetime is obtained. Then, by using the equations of motion, one can find the ray paths in each spacetime. Finally, in order to check the accuracy of simulations, we compare the propagation of waves in optical materials with the ray paths driven from the equations of motion.

II. SPACETIME GEOMETRY AND MEDIA
Using the mathematical machinery of differential geometry, one can write Maxwell equations in arbitrary coordinates and geometry. It has been shown that the source-free Maxwell equations in arbitrary right-handed spacetime coordinates can be written as the macroscopic Maxwell equations in the right-handed Cartesian coordinates with Plebanski's constitutive equations [6]. It means that Maxwell equations in curved coordinates are equivalent to their standard form for the flat space but in the presence of an effective medium. The Plebanski's constitutive relations have been found in the form [19]: where permittivity and permeability tensors and a vector Γ are given by: Here, g ij is the inverse of g ij , the spatial part of full spacetime metric g µν , and g is the determinant of g µν . Designing artificial materials with the permittivity and permeability introduced in Eq. (2) enables us to mimic different spacetimes. It is notable that all the information about the gravitational field is embedded in the material properties of the effective medium. Anisotropic or isotropic materials in which the constitutive relations described by Eq. (1) are called bianisotropic or bi-isotropic materials [20][21][22][23][24]. The vector Γ is the magnetoelectric coupling parameter coupling the electric and magnetic fields. Metamaterials are artificial materials made of sub-wavelength metallic constituents which are randomly or periodically distributed in a dielectric background. In some metamaterials the cross polarization effect (an electric polarization results from an applied magnetic field and vice versa) occurs which are called bianisotropic or isotropic metamaterials [20]. By choosing suitable distribution of the metallic inclusions in these materials, it is possible to obtain desired electric permittivity, magnetic permeability and magnetoelectric coupling parameter following Eq. (2). Fabrication of bianisotropic and biisotropic metamaterials operating in the visible wavelengths is possible using different techniques such as lithography and laser writing [20]. To understand how light wave propagates near a black hole in different space-times, we can experimentally and theoretically study the propagation of electromagnetic waves in the corresponding metamaterials. Then, to prove the correctness of this equivalence, we compare our results with trajectories obtained using Lagrangian formalism. In the other words, since the trajectories near black holes in different spacetimes are well studied, we compare them with the trajectories obtained from the equivalence optical medium to prove the correctness of our method. It should be noted that the electromagnetic (optical) black hole [25] is fabricated for the first time in 2010 based on metamaterials composed of nonresonant I-shaped inclusions and electric field coupled resonators [11]. This structure can absorb the electromagnetic wave incident from every direction like a black hole or a black body.
In other words, the mixing of electric and magnetic fields is addressed by Γ with the physical dimension of velocity. It is also shown that for a slow-moving medium, u/c << 1, Γ is proportional to the velocity of medium [26].

III. MAXWELL EQUATIONS IN MEDIUM
Regarding to the first approach, Maxwell equations written in a source-free form [27] ∇. D = 0, should be supplemented by the constitutive relations Eq. (1) and medium tensors expressed spacetime, in order to be solved. Since the constitutive equations couple the magnetic and electric fields, Eq. (3) become too complicated to be solved, we have to seek for a method to decouple equations. In this regard, we apply the time-harmonic Maxwell equations, assuming a time dependence e −iωt with t and ω as time variable and angular wave frequency, respectively, where E, H, D and B are, respectively, electric field, magnetic field, electric displacement and magnetic flux density. We have considered that the constitutive equations obey Eq. (1), and therefore, substituting this equation to Eq. (4) we can find the following Maxwell-Tellegen equations [28]: We have assumed that the material parameters, , µ and Γ, and the electromagnetic fields are invariant in one direction, here z direction is considered. Since , µ and Γ are z−anisotropic tensors, in general, we can write: Applying Eq. (6) to Eq. (5) leads to: It is patently obvious from Eq. (7-12) that electric and magnetic fields are coupled because of the constitutive equations. In order to decouple these equations, we follow the approach introduced in Ref. [28]. At the first step, we rewrite Eqs. (7) and (8) as the following unified relation and also Eqs. (10) and (11) as where Next, one can drive H and E from Eqs. (13) and (14), respectively and We can also rewrite Eqs. (9) and (12) in the following forms and Substituting Eqs. (17) and (16) into Eqs. (18) and (19), one finds and It is notable that all Eqs. (7)(8)(9)(10)(11)(12) are resolved to two Eqs. (20) and (21), which are decoupled equations with the longitudinal fields E 3 and H 3 as unknowns.
We have considered TE polarization for an electromagnetic wave for which E is perpendicular to the xy plane. If E is determined in the z direction, the only non-zero component is E 3 . Hence, as the electric and magnetic fields are perpendicular, H 3 is zero. By applying these assumptions to Eqs. (20) and (21), one concludes that only Eq. (21) is required to solve, since Eq. (20) is always satisfied.

IV. LAGRANGIAN FORMALISM
In the previous section, we used wave optics in order to determine the propagation of light in optical materials mimicking spacetime. In this section, we apply geometrical optics in order to determine the propagation of ray in that materials to be assure that our model works well. In this regard, the Lagrangian is introduced, then based on the equation of motion the ray path is obtained.
It was shown [30] that the equations governing the geodesics in a spacetime with the line element ds 2 = g αβ dx α dx β can be derived from the Lagrangian where τ is the affine parameter along the geodesics. From the calculus of variations, one can write the following Euler-Lagrange equation where q α = (t, r, θ, φ) is the generalized coordinate, and the generalized momentum p α is defined as It is easy to determineṗ i from the Euler-Lagrange equation, Eq. (23), as beloẇ where dot denotes the derivative over the affine parameter, τ . Based on Eqs. (22) and (24), we can rewrite the Lagrangian as a function of generalized momentum, p i . Since we are seeking the ray path, the Lagrangian, Eq. (22), is zero.
V. REISSNER-NORDSTRÖM SPACETIME First, we apply the analogy to the Reissner-Nordström spacetime. The Reissner-Nordström metric is a static solution of the Einstein-Maxwell field equations, which corresponds to the gravitational field of a charged, nonrotating, spherically symmetric body of mass M [29]. This metric in the Cartesian coordinates can be written as Here Q is the charge of black hole, m is the geometrical mass of the source of gravitation and r = x 2 + y 2 + z 2 . The anisotropic equivalent medium tensors can be obtained from Eq. (2). It is notable that Γ i = 0 because the Reissner-Nordström metric is static. Therefore, permittivity and permeability tensors in the xy plane can be obtained as Since we are working in the xy plane, here r = x 2 + y 2 . Having permittivity and permeability tensors, Eq. (21) can be solved. We have solved the wave equation Eq. (21) numerically for a beam with the appropriate boundary conditions.
It is known that having a non-trivial spherically symmetric electromagnetic wave is not possible. Its reason comes from the fact that the polarization vectors of such a field configuration would introduce a continuous nowhere-vanishing vector field tangent to the 2-sphere, thereby contradicting the well-known fact that the 2-sphere is not parallelizable, and therefore, we motivate to use other coordinate systems. The optical parameters of the material equivalent to the black hole depend on the coordinate system. Because working with Cartesian coordinate in COMSOL is easier than other coordinate systems, we use it in this paper. It is also notable that although the (non-physical) coordinate singularity can be removed in the Cartesian coordinate, its one-way membrane property is kept for the equivalent metamaterials.
The Reissner-Nordström metric in spherical coordinates reads [29] Substituting this metric into Eq. (22) and considering the xy plane, we find Considering the obtained Lagrangian and using Eqs. (24) and (25), p i andṗ i can be easily calculated aṡ We define p t = 2E and p φ = 2L that are constant sinceṗ t = 0 andṗ φ = 0, which just state the conservation of energy and angular momentum for the static spherically symmetric system. Hence, the Lagrangian Eq. (29) can be rewritten as Considering L = 0, one findsṙ andφ can be determined from Eq. (33), therefore Since photons are massless, the ray trajectories can be determined by only one character which is called the impact parameter D = L E [30]. The ray trajectories can be obtained by solving Eq. (36) [30,31]. The results for Reissner-Nordström spacetime, both for wave optics and ray optics are shown in Figs. 1-4. Interestingly, depending on the roots of f (r), we have three types of trajectories: "terminating orbit", "bound orbit" and "flyby orbit". Those roots depend only on the impact parameter, D, since both m and Q are constants. Therefore, for D = D c we have "bound orbit", if D > D c there is "flyby orbit" and if D < D c there is "terminating orbit" [31]. The critical value of impact parameter, D c , can be driven by considering f (r c ) = 0 and f (r c ) = 0 It is obvious that for different values of charge, Q, the critical value of impact parameter differs.
Regarding the Lagrangian formalism, first, we write the Kerr metric in the spherical coordinates [29] where ρ 2 = r 2 + a 2 cos 2 θ and ∆ = r 2 − 2mr + a 2 . Since we are working in the xy plane, we set θ = π 2 . Substituting this metric into Eq. (22), one can show that Using Eqs. (24) and (25), the following relations are obtaineḋ where, we call p t = 2E and p φ = −2L, as before. By solving Eqs. (45) and (46) simultaneously,ṫ andφ can be calculated asṫ Moreover, substituting these parameters into the Lagrangian Eq. (42), one finds Since for light beam we have L = 0, It is noteworthy to mention that only the impact parameter affects the ray trajectories which its effects lead to three types of orbits mentioned in the previous section. At this point, since the equations related to Kerr spacetime are too complicated, we just focus on critical value of impact parameter. In this regard, by calculating f (r c ) = 0 and f (r c ) = 0, one finds [30] It is clear that the spin parameter, a, controls the critical value of the impact parameter. By considering u = 1 r and considering relations in Eqs. (50) and (51), we can rewrite Eq. (50) aṡ Taking into accountφ in Eq. (48), we can obtain Furthermore, solving Eq. (53) gives the ray trajectories [30]. The results for Kerr spacetime are shown in Fig. 5.

VII. RESULTS AND DISCUSSION
Here, we are in a position to solve the equations of two approaches and simulate the results for comparison. Equation (21) is solved in the Cartesian coordinates and a rectangular 2D geometry of xy space for a T E polarized wave of frequency ω = 6.3 × 10 9 (rad/s) injected from the right and the results are compared with the ray path (red line) calculated from the equations of motion. The computational domain is surrounded by a perfectly matched layer that absorbs the outward waves to ensure that there are no unwanted reflections, and the simulations are done by means of a standard software solver (COMSOL).
For Reissner-Nordström spacetime the medium parameters are given by Eq. (27) and the results are depicted in Table I and Figs. 1-4. In order to obtain dimensionless parameters, we normalize all parameters to r 0 = rs 2 in which r s = 2m is the event horizon of Schwarzschild black hole where the gravity is so strong, that light cannot escape [30].  Taking into account the plotted figures, we find that the ray path follows the center of the beam. Moreover, the ray trajectories are consistent with the behavior expected from theory. Speaking more precisely, the case of D = D c represents photon sphere [32] which is a spherical region of space where the gravitational effect is strong enough that photons are forced to travel in an orbit (in the xy plane). This area is patently obvious in Figs. 1b, 2b, 3b and 4b. The photon sphere is located farther from the center of a black hole than the event horizon, the radius of photon sphere can be determined as r c in Eq. (37) and some numerical results of r c can be found in Table I. In addition, for D < D c , we expect that the photons are drown into the event horizon, it can be seen in Figs. 1a, 2a, 3a and 4a. Moreover, in the case of D > D c , deflection of ray trajectory can be observed in Figs. 1c, 2c, 3c and 4c. Furthermore, it is notable that the charge of black hole, Q, affects the critical value of impact parameter, D c , and as a result the radius of photon sphere, r c is changed. It can be seen from Table I that when Q increases, the critical impact parameter, D c , and radius of photon sphere, r c , decrease.
It can be observed from all figures that the beam splits into a set of rays or sub-beams, one part falls into black hole due to having D < D c , the other part escapes from black hole because of D > D c and another part bends around the photon sphere of the black hole and interferes with the primary beam.
It is expected from theory that the case of Q = 0 represents Schwarzschild spacetime. The radius of photon sphere is determined r c = 3m for Q = 0 that is similar to our expectation, also Fig. 1 matches the results reported in Ref. [17].
For the Kerr spacetime the medium parameters are given by Eqs. (39) and (40), while the numerical and simulated results are shown in Table II and Fig. 5. It is notable that in order to have dimensionless parameters, they are normalized to r s , as before. As we mentioned before, computational analysis of rotating spacetime is complicated, and therefore, we just focus on the critical values of the impact parameter which lead to photon sphere. In contrast to the Reissner-Nordström black hole, the Kerr black hole is not spherical symmetry, but only enjoys axially symmetry, which has profound consequences for the photon orbits. A circular orbit can only exist in the equatorial plane [32]. Fortunately, both methods of simulations, numeric solutions to Maxwell equations and the ray trajectory arisen from the Hamiltonian formalism, are in agreement to each others. As it is depicted in Table II, by decreasing the spin parameter (a), the critical impact parameter (D c ) and the radius of photon sphere (r c ) increase. In addition, in the case of static case (a = 0), the results of the Kerr black hole match to that of driven from Schwarzschild black hole, as expected.

VIII. CONCLUDING REMARKS
Here, we compared geometrical description of gravitational systems with optical media in flat space with two approaches. We observed a good correspondence between the two methods.
We considered two supplementary factors of the Schwarzschild black hole; the electric charge and rotation parameters. More precisely, we investigated the effects of electric charge and rotation factors on the trajectory of light with various impact parameters. We found that the medium parameters given for Reissner-Nordström and Kerr spacetime can mimic the ray path near the black hole in the flat spacetime. Interestingly, some theoretical phenomena like photon sphere are observed and the effects of charge, Q, and spin parameter, a, are investigated. We showed that the critical impact parameter is a decreasing function of both the electric charge and spin parameters. This behavior is expected due to the fact that increasing the electric charge and rotation lead to decreasing the event horizon radius (weaken the gravitational effect) in the classical black hole scenario.
It will be interesting to apply these calculations to the KerrNewman black hole and other interesting nontrivial solutions of Einstein gravity.

IX. ACKNOWLEDGEMENTS
We wish to thank Shiraz University Research Council. This work has been supported financially by the Research Institute for Astronomy and Astrophysics of Maragha, Iran.