## Abstract

An optical black hole is studied using a parallel radially dependent finite-difference time-domain (FDTD) simulation technique. The device requires non-dispersive metamaterial structures and is capable of broadband operation, based on transformation optics. Excellent absorption is demonstrated for different angles of wave incidence and illumination excitation types. In addition, a practical device, which is made to be matched to free space, is proposed, and the relevant physics is explored. Finally, peculiar phase distributions of the electromagnetic waves are observed inside the radially dependent permittivity material of the devices.

© 2010 Optical Society of America

## 1. INTRODUCTION

Transformation optics provides methodologies for engineers to manipulate electromagnetic radiation at will [1, 2]. Several microwave and optical device designs have been proposed, including broadband ground-plane cloaking structures [3, 4, 5, 6, 7, 8] and exotic absorbing devices [9, 10]. More recently, an interesting application that has been proposed is the design of artificial optical black holes [11, 12]. Mimicking their celestial counterparts, these devices bend the electromagnetic waves appropriately and can almost completely absorb the incident electromagnetic radiation crossing into their region. In addition, they can operate for all angles of incidence over a broad frequency spectrum and can be constructed with non-resonant metamaterials. Potential applications of these exotic devices are perfect absorbers [13], efficient solar energy harvesting photovoltaic systems [14], thermal light emitting sources, and optoelectronic devices [15].

It should be noted that despite their light-absorbing capabilities, the name “black hole” is a misnomer in this case, because these devices do not possess the main characteristic of gravitational black holes, namely, the event horizon, the artificial boundary around the device beyond which no light can escape. As we also show here, radiation generated inside the device does indeed escape into the surrounding environment. However, the original term black hole for the device as introduced in [11] will be maintained here for consistency, despite the fact that a term like “optical attractor” might be more appropriate [12].

We investigate the performance of the spherical optical black hole embedded in a background medium using a full-wave simulation technique, the well-established finite-difference time-domain (FDTD) method [16]. Due to the large dimensions of the three-dimensional (3-D) device, a parallel version of the FDTD technique is utilized, which divides the simulation domain into sub-domains that are processed in parallel and significantly decrease the total simulation time. Different excitation types are chosen to illuminate the device, namely, plane waves and spatially Gaussian pulses. Moreover, the performance of an alternative black hole, which is not embedded in a particular material [11] but is matched directly to free space, is examined. The latter is a more practical approach to the optical black hole, in a similar fashion that the ground-plane quasi-cloaks operate in free space, as proposed in [17] and verified experimentally in [18]. The performance of the latter device is found to be similar to the performance of the embedded structure. Finally, the losses at the core of the black hole are removed, and a point source is placed inside. We observe that the cylindrical wavefronts of the source are disturbed from the radially dependent permittivity of the metamaterial structure, leading to peculiar field phase distributions.

## 2. PARALLEL RADIALLY DEPENDENT FDTD TECHNIQUE

The FDTD method is used to explore the physics of the absorbing device. The full-wave numerical technique can efficiently demonstrate the performance of the broadband device. The proposed technique can also model the diffraction effects and near fields, which are neglected with ray tracing methods used previously [11]. As a result, a more clear and complete picture of the behavior of the modeled device is achieved.

The FDTD method is based on the temporal and spatial discretization of Faraday’s and Ampere’s laws, which are

where $\overline{E}$, $\overline{H}$, $\overline{D}$, and $\overline{B}$ are the electric field, magnetic field, electric flux density, and magnetic flux density components, respectively. Note that harmonic time dependence $\text{exp}\left(\u0131\omega t\right)$ of the field components is assumed throughout this paper. Faraday’s and Ampere’s laws are discretized with a common procedure [19], and the conventional updating FDTD equations are obtained for the field components:*n*is the number of the current time steps. The temporal discretization $\Delta t$ has to satisfy the Courant stability criterion $\Delta t=\Delta x/\left(\sqrt{3}c\right)$ [19] to achieve stable FDTD simulations. $\Delta x$ is the uniform spatial resolution of the FDTD cubical Cartesian mesh. The relative permittivity $\epsilon \left(x,y,z\right)$ and permeability $\mu \left(x,y,z\right)$ can be spatially dependent according to the geometrical features of the material or metamaterial structure that is simulated. The derived FDTD equations are similar, and simpler, compared to the previously proposed radially dependent numerical technique [20, 21] used to model another interesting transformation-based device: the cloak of invisibility. The difference is that unlike most invisibility cloaks, the black hole is composed of non-dispersive materials, as it will be shown in the next section. Thus, there is no reason to discretize the constitutive equations leading to a more complicated FDTD algorithm.

In the case of the computational intensive modeling of the 3-D spherical black hole, a parallel version of the FDTD algorithm is used. The FDTD updating equations are the same as before. The only difference is that the computational domain is divided into smaller sub-domains based on the space decomposition technique [22], and every domain is then assigned to one processor. The tangential field components are communicated between the adjacent interfaces at each time step with an appropriate synchronization procedure, which is provided by the message passing interface library. After the calculations are completed and the solver has finished, the resulting fields from each sub-domain are retrieved and are combined together to obtain the results in the whole initial domain. The parallel version of the algorithm is ideal for handling complicated and computational intensive problems, comprised of huge computational domains and complex electromagnetic parameters.

## 3. PARAMETERS OF SPHERICAL/CYLINDRICAL OPTICAL BLACK HOLE

The spherical 3-D and the cylindrical two-dimensional (2-D) optical black holes will be investigated. Each device is divided into two regions: the core (usually absorbing) and the shell. The radially dependent permittivity distributions of the spherical and cylindrical black holes are given by the formula [11]

The objective is to achieve a matched device with the surrounding space, reducing the reflection of the impinging electromagnetic waves to a minimum. To obtain this effect at the surface of the device, the radius of the device’s core has to vary according to the parameters of the surrounding medium, for a given core material with permittivity ${\epsilon}_{c}$. It is given by

## 4. NUMERICAL RESULTS OF THE SPHERICAL/CYLINDRICAL OPTICAL BLACK HOLE

#### 4A. Spherical Optical Black Hole Embedded in Dielectric Material

In this subsection we investigate a spherical black hole operating as an electromagnetic concentrator for improving the light capturing capabilities of solar cells. The spherical black hole is embedded in silica glass $\left({\text{SiO}}_{2}\right)$ with a relative permittivity of ${\epsilon}_{0}=2.1$. The core of the device is composed of *n*-doped silicon with relative permittivity ${\epsilon}_{c}+\u0131\gamma =12+\u01310.7$, which is a typical material of a thin film solar cell. The device is designed to operate at the infrared section of the spectrum with a central frequency of $f=200\text{\hspace{0.17em} THz}$. The range of permittivity values, required to construct the structure, is between 2.1 and 12, as indicated by Eq. (5). The device can have broadband operation, as any frequency dispersion is introduced only by whatever material is chosen to build it. Finally, all the six field components ${E}_{x}$, ${E}_{y}$, ${E}_{z}$, ${H}_{x}$, ${H}_{y}$, ${H}_{z}$ exist at the 3-D FDTD simulations.

The parameters of the device have isotropic values, where $\epsilon =\epsilon \left(r\right)$ is given in Eq. (5), and the magnetic parameters have the permeability of free space. A uniform spatial resolution is chosen for the Cartesian FDTD mesh given by $\Delta x=\Delta y=\Delta z=\lambda /20$, where $\lambda =1.5\text{\hspace{0.17em}}\mu \text{m}$ is the wavelength of the excitation signal in free space. The Courant stability condition is satisfied, and the computational domain is terminated with Berenger’s perfectly matched layers (PML) [23], modified to be embedded in the silica glass surrounding material. The excitation field is chosen to be a temporally infinite spatially Gaussian pulse with a similar wave trajectory to the ray tracing results presented in [11]. The FDTD computational domain, required for accurate modeling of the spherical black hole, is equal to $30\lambda \times 30\lambda \times 30\lambda $. It is divided along the *z*-axis to 60 sub-domains, where individual processors are solving the parallel FDTD updating equations, as explained earlier. Large amount of memory is required for every simulation, roughly 50 Gbytes of random access memory, which is also distributed between the computing nodes. Each simulation lasts approximately 18 h (30,000 time steps), until the steady state is reached and the field values no longer evolve significantly with time.

A spatially confined *y*-polarized $\left({E}_{y}\right)$ Gaussian pulse is chosen to impinge on the spherical black hole from two different angles of incidence to evaluate its omnidirectional absorbing performance. Both pulses are $\lambda /2$ wide in space [full width at half-maximum (FWHM)], and the frequency of operation is chosen as $f=200\text{\hspace{0.17em} THz}$. The excitation pulses are infinite in time. After steady state is reached, the 3-D electric field amplitude distribution ${E}_{y}$ is retrieved. In order to visualize the results, we plot the 2-D amplitude distribution on three different planes, each one parallel to each of the Cartesian planes of the domain. The results for two different simulation scenarios are shown in Figs. 1, 2 . The locations of the pulse entrance for each figure are $\left(15\lambda ,15\lambda ,1\lambda \right)$ and $\left(15\lambda ,7.5\lambda ,4\lambda \right)$, respectively. It can be clearly seen that the incoming field power is totally absorbed inside the core of the device for both cases. The field trajectories rapidly bend toward the core of the device (see Fig. 2), as expected. The device is matched to the surrounding material and, as a result, the reflections of the device are almost zero. Furthermore, the absorption of the radiation is excellent, reaching approximately 95%. Finally, note that the 3-D figures were constructed after downsampling the computational domain to half its original size, due to memory constraints in the post-processing.

#### 4B. Cylindrical Optical Black Hole Embedded in Dielectric Material

In the next subsections of the paper, the study will be focused on the 2-D black hole design due to its simplicity in FDTD modeling compared to the parallel FDTD simulations. As it shall be shown, the results of the cylindrical device simulations are very similar to the case of the 3-D spherical black hole studied in the previous subsection. The cylindrical black hole is embedded in ${\text{SiO}}_{2}$, and the core of the device is composed of *n*-doped silicon. The frequency of interest is $f=200\text{\hspace{0.17em} THz}$, and the device can have broadband operation. The cylindrical coating inner and outer radii are ${R}_{c}=5.6\lambda $ and ${R}_{sh}=13.33\lambda $, respectively, the same with those of the spherical embedded black hole. Transverse magnetic polarization is used throughout the 2-D simulations of the optical black hole, without loss of generality. Only three field components exist ${E}_{x}$, ${E}_{y}$, ${H}_{z}$, and the parameters are isotropic and equal to ${\epsilon}_{x}={\epsilon}_{y}=\epsilon \left(r\right)$ and $\mu ={\mu}_{0}$, given in Eq. (5).

The spatial resolution is chosen uniform and equal to $\Delta x=\Delta y=\lambda /30$, where *λ* is the wavelength of the excitation signal at free space. The temporal resolution of the 2-D simulation is chosen as $\Delta t=\Delta x/\sqrt{2}c$, according to Courant stability condition [19], where *c* is the speed of light in free space. The computational domain is terminated again with PMLs [23]. During all the 2-D FDTD simulations presented here, the steady state is reached after approximately 8000 time steps or 1.5 h. A temporally infinite *z*-polarized $\left({H}_{z}\right)$ spatially Gaussian pulse is used, $2\lambda /3$ wide in space (FWHM), to illuminate the black hole, in order to mitigate the ray tracing simulations [11]. Finally, the FDTD domain is noticeably large $28\lambda \times 28\lambda $, because the device operates correctly only when its dimensions are much larger than the wavelength. The real part and the amplitude of the magnetic field distribution when a Gaussian beam is impinging on the proposed black hole, after steady state is reached, are reported in Figs. 3a, 3b (Media 1) respectively. The performance of the device is excellent, achieving almost full absorption of the beam, similar to the performance of the spherical black hole shown in Fig. 2.

Next, a *z*-polarized $\left({H}_{z}\right)$ plane wave illuminates the device. The plane wave is constructed after replacing the top and bottom PMLs with periodic boundary conditions [19]. The results are shown in Figs. 4a, 4b , where the real part and the amplitude of the magnetic field distribution are shown. A large shadow is cast at the back of the structure, and the reflections are negligible. The largest part of the plane wave’s energy is absorbed at the core, as can be clearly seen in Fig. 4b. The behavior of the device is that of an efficient electromagnetic absorber, similar to the one proposed in [10].

#### 4C. Cylindrical Optical Black Hole Embedded in Free Space

If we replace the background material of silica glass (used before) with free space, the range of the black hole’s relative permittivity values will be wider and between 1 and 12. The wider material values are leading to a more complicated design of the metamaterial structure and a smaller radius of its core, given by Eq. (6). Hence, the cylindrical coating inner and outer radii become ${R}_{c}=3.85\lambda $ and ${R}_{sh}=13.33\lambda $, respectively. However, a device that is matched directly to free space is always more desirable for practical applications.

The free space cylindrical black hole is modeled with the same FDTD method as before and at the same simulation scenario shown in Fig. 3. The fields are obtained in Figs. 5a, 5b , and similar results are derived with the spherical and cylindrical embedded black holes in silica glass. There is slightly increased scattering, because the electromagnetic beam is traveling optically longer distance inside the radially dependent permittivity material. Nevertheless, the absorption of the structure is again excellent. Finally, when the beam is incident from a different angle in Figs. 6a, 6b , the performance is similar to the earlier results (particularly see Fig. 1). Thus, the optical black hole can be regarded as an omnidirectional absorber.

#### 4D. Phase Distribution of Source Placed Inside the Black Hole

The core of the free space optical black hole was initially composed of a lossy material with permittivity ${\epsilon}_{c}+\u0131\gamma =12+\u01310.7$. The imaginary part is now removed from the device, and a soft *z*-polarized $\left({H}_{z}\right)$ point source, again at a frequency of $f=200\text{\hspace{0.17em} THz}$, is placed inside the core asymmetrically with a slight deviation from the center of the *y*-axis at the point $\left(x,y\right)=\left(15\lambda ,12.33\lambda \right)$. This is conducted in order to further explore the wave interactions inside the radially dependent metamaterial shell. The real part and the phase of the ${H}_{z}$ field component are shown in Figs. 7a, 7b . A peculiar phase distribution of the cylindrical waves is observed, especially in Fig. 7b.

The radially dependent permittivity—given by Eq. (5)—is substituted in the formulas of phase and group velocity to explain the unusual phase distribution of the metamaterial structure. The derived formulas are

Note that if the source was placed at the center of the structure, only perfectly outgoing cylindrical wavefronts would exist. By placing the source away from the center, the device converts some of the cylindrical wavefronts into flat ones on the opposite side, exhibiting a lens-like behavior. This effect is caused by the longer optical path, where the wavefronts have to travel passing through the whole length of the dielectric core of the structure. Note that the radiation generated inside the device’s core always escapes to the surrounding environment, even if the imaginary part of the device is not removed (not shown here).

## 5. CONCLUSION

To conclude, the spherical (3-D) and the cylindrical (2-D) optical black holes were studied numerically with the FDTD method. Full-wave simulations were performed to study the performance of the device under different excitations and angles of incidence. The physics of this exotic metamaterial structure was explored. It behaves as an excellent omnidirectional absorber achieving absorption values of approximately 95%. An alternative black hole, matched to the surrounding free space, was proposed, and its excellent performance was demonstrated. Although the devices were tested with infinite in time excitation pulses, similar results are expected if the devices are radiated with non-monochromatic pulses. The devices do not have broadband limitations other than the natural dispersion of the conventional materials chosen to build them. Finally, the response of the radially dependent permittivity material to cylindrical waves was studied in detail. The device can have several potential applications, ranging from excellent absorbers to state-of-the-art solar cells and optoelectronics.

**1. **J. B. Pendry, D. Schurig, and D. R. Smith, “Controlling electromagnetic fields,” Science **312**, 1780–1782 (2006). [CrossRef] [PubMed]

**2. **U. Leonhardt, “Optical conformal mapping,” Science **312**, 1777–1780 (2006). [CrossRef] [PubMed]

**3. **J. Li and J. B. Pendry, “Hiding under the carpet: A new strategy for cloaking,” Phys. Rev. Lett. **101**, 203901 (2008). [CrossRef] [PubMed]

**4. **R. Liu, C. Ji, J. J. Mock, J. Y. Chin, T. J. Cui, and D. R. Smith, “Broadband ground-plane cloak,” Science **323**, 366–369 (2009). [CrossRef] [PubMed]

**5. **J. Valentine, J. Li, T. Zentgraf, G. Bartal, and X. Zhang, “An optical cloak made of dielectrics,” Nature Mater. **8**, 568–571 (2009). [CrossRef]

**6. **L. H. Gabrielli, J. Cardenas, C. B. Poitras, and M. Lipson, “Silicon nanostructure cloak operating at optical frequencies,” Nat. Photonics **3**, 461–463 (2009). [CrossRef]

**7. **J. H. Lee, J. Blair, V. A. Tamma, Q. Wu, S. J. Rhee, C. J. Summers, and W. Park, “Direct visualization of optical frequency invisibility cloak based on silicon nanorod array,” Opt. Express **17**, 12922–12928 (2009). [CrossRef] [PubMed]

**8. **T. Ergin, N. Stenger, P. Brenner, J. B. Pendry, and M. Wegener, “Three-dimensional invisibility cloak at optical wavelengths,” Science **328**, 337–339 (2010). [CrossRef] [PubMed]

**9. **J. Ng, H. Chen, and C. T. Chan, “Metamaterial frequency-selective superabsorber,” Opt. Lett. **34**, 644–646 (2009). [CrossRef] [PubMed]

**10. **C. Argyropoulos, E. Kallos, Y. Zhao, and Y. Hao, “Manipulating the loss in electromagnetic cloaks for perfect wave absorption,” Opt. Express **17**, 8467–8475 (2009). [CrossRef] [PubMed]

**11. **E. E. Narimanov and A. V. Kildishev, “Optical black hole: Broadband omnidirectional light absorber,” Appl. Phys. Lett. **95**, 041106 (2009). [CrossRef]

**12. **D. A. Genov, S. Zhang, and X. Zhang, “Mimicking celestial mechanics in metamaterials,” Nat. Phys. **5**, 687–692 (2009). [CrossRef]

**13. **N. I. Landy, S. Sajuyigbe, J. J. Mock, D. R. Smith, and W. J. Padilla, “Perfect metamaterial absorber,” Phys. Rev. Lett. **100**, 207402 (2008). [CrossRef] [PubMed]

**14. **H. A. Atwater and A. Polman, “Plasmonics for improved photovoltaic devices,” Nature Mater. **9**, 205–213 (2010). [CrossRef]

**15. **J. A. Schuller, E. S. Barnard, W. Cai, Y. C. Jun, J. S. White, and M. L. Brongersma, “Plasmonics for extreme light concentration and manipulation,” Nature Mater. **9**, 193–204 (2010). [CrossRef]

**16. **K. Yee, “Numerical solution of initial boundary value problems involving Maxwell’s equations in isotropic media,” IEEE Trans. Antennas Propag. **14**, 302–307 (1966). [CrossRef]

**17. **E. Kallos, C. Argyropoulos, and Y. Hao, “Ground-plane quasicloaking for free space,” Phys. Rev. A **79**, 063825 (2009). [CrossRef]

**18. **H. F. Ma, W. X. Jiang, X. M. Yang, X. Y. Zhou, and T. J. Cui, “Compact-sized and broadband carpet cloak and free-space cloak,” Opt. Express **17**, 19947–19959 (2009). [CrossRef] [PubMed]

**19. **A. Taflove and S. C. Hagness, *Computational Electrodynamics: The Finite-Difference Time-Domain Method*, 3rd ed. (Artech House, 2005).

**20. **C. Argyropoulos, Y. Zhao, and Y. Hao, “A radially-dependent dispersive finite-difference time-domain method for the evaluation of electromagnetic cloaks,” IEEE Trans. Antennas Propag. **57**, 1432–1441 (2009). [CrossRef]

**21. **C. Argyropoulos, E. Kallos, and Y. Hao, “Dispersive cylindrical cloaks under nonmonochromatic illumination,” Phys. Rev. E **81**, 016611 (2010). [CrossRef]

**22. **K. C. Chew and V. F. Fusco, “A parallel implementation of the finite difference time-domain algorithm,” Int. J. Numer. Model. **8**, 293–299 (1995). [CrossRef]

**23. **J. -P. Bérenger, “A perfectly matched layer for the absorption of electromagnetic waves,” J. Comput. Phys. **114**, 185–200 (1994). [CrossRef]