1. Introduction

Recent advances in the development of optical metamaterials, artificial structures with the desired local electromagnetic response determined by artificial subwavelength patterning of the media, have opened a wide range of new concepts and applications, from super-resolution imaging to optoelectronics and photovoltaics [1–8]. In particular, the recently proposed [8] “optical black hole”, an omnidirectional light concentrator and absorber, could significantly enhance the performance of solar energy harvesting systems and lead to novel nonlinear-optical devices. In this work, we present a detailed theoretical description of this device, and develop appropriate numerical methods, which – in addition to the accurate analysis of the optical black hole, - could be applied to a broad range of metamaterial-based optical systems with centrally symmetric gradient index. We note that in the context of this paper, the term “black hole” simply refers to the highly efficient “capture” of the light incident on the device, and does not imply any profound analogy to General Relativity.

The remaining part of the paper is organized as follows. In Section 2, we introduce the initial assumptions, state the geometry of the cylindrical optical system with an axially symmetric graded index, and show basic examples of ray-tracing inside such devices, including the fundamental optical black-hole [8]. Section 3 deals with the ray optical description of the system followed by its full-wave description in Section 4. The analytical and numerical apparatus for approximating the ideal fundamental black-hole with concentric lamellar systems of uniform isotropic layers is then discussed in Section 5. Section 5 also covers derivations of the scattering and absorption cross-sections, which are then compared with their semiclassical approximations in Section 6. Finally, we summarize our work in Section 7.

2. The system

In its original proposal [8], the optical “black hole” (see Fig. 1(a) ) consists of the “payload” core with the radius $r_c$, supporting the desired functionality (such as e.g. a solar cell for concentrated photovoltaics applications), surrounded by the concentrating shell of radius $r_s$with an axially symmetric dielectric function, $\epsilon \left(r\right)$, changing in the radial direction as

 

Fig. 1 (a) the schematics of the system; panels (b), (c), (d), and (e) show light ray trajectories for p equal to −1, 1, 2, and 3 respectively, where core radius r _c much smaller than shell radius r _s (r _c = r _s).

Download Full Size | PDF

3. Ray-optical description

In the semi-classical limit, the radial variation of the refractive index of the shell (where $n=\sqrt{\epsilon }$) corresponds to the effective potential

Case $p=1$ corresponds to the gravitational potential (or electrostatic potential) and emulates the motion of a non-accelerated particle in central gravitational field; so that the parabolic trajectories (Fig. 1c) are expected. Case $p=2$ gives the potential of fatal attraction [9]. In fact for all $p\ge 2$ (6) corresponds to the ray trajectories “falling” into the core of the “black hole” – see Fig. 1(b). Figures 1(b)-1(e) show the traces of light rays for p equal to −1, 1, 2, and 3 respectively; the trivial case of the beam propagating in free-space ($p=0$) is omitted.

This concept can be realized in both the spherical and cylindrical geometries, with the latter suitable for optoelectronic devices where the light propagation is generally realized in a waveguide geometry. In the present paper we will focus on this cylindrical realization of the optical black hole, leaving its spherical counterpart to a forthcoming work.

4. The wave optical description

In the present section of the paper, we develop the wave-optical description of the “black hole” light concentrator. We define a given monochromatic electromagnetic vector field by the vectors of electric field $E\exp \left(-\iota \omega t\right)$ and magnetic field $H\exp \left(-\iota \omega t\right)$, and by its angular frequency ω. We start with the monochromatic Maxwell’s equations $\nabla \times E=\iota \omega B$, $\nabla \times H=-\iota \omega D$ arriving first at [10]:

Here we assume only transverse electric (TE, $E=\widehat{z}e(x,y)$) and transverse magnetic (TM, $H=\widehat{z}h(x,y)$) waves within a nonmagnetic medium ($\mu =1$), so we reduce the above vector equations and use in-plane scalars e and h, which are further on collectively denoted as f.

4.1 Scalar wave equations and separation of variables

Hence, both TE and TM wave Eqs. (9) and (10) reduce to a common scalar wave equation

The cylindrical symmetry of the device allows a straightforward separation of radial and angular variables in the cylindrical coordinates, (r, ϕ, z), linked to the Cartesian coordinates (x, y, z) through $x=r\cos \varphi$, $y=r\sin \varphi$. Exploiting the cylindrical symmetry of the dielectric function, $\epsilon =\epsilon \left(r\right)$, and decomposing the solution of (11) into the cylindrical waves, $f={\displaystyle {\sum }_{m=-\infty }^{+\infty }{f}_m\left(r\right)\exp \iota m\varphi }$, we arrive at equation for scalar modes $f_m$

4.2 TE and TM solutions for $\epsilon \left(r\right)=\text{const}$

The solvable cases that we are using here include the trivial instance of $\epsilon \left(r\right)=\text{const}$ for both TE and TM modes. In this case, (12) reduces to the Bessel equation for $\overline{r}=kr$, and the general solution is a linear combination of two linearly independent Bessel functions [11].

In most of the cases it is more convenient to use the following combination

For simulations we use the normalized variant of (13), which improves the condition number of the linear system matrix (Section 5.1):

4.3 TE and TM solutions for the ‘quadratic decay’ of permittivity, $r^2\epsilon \left(r\right)=\text{const}$

For the cylindrical realization of the black hole, which in its simplest form corresponds to an infinite cylinder with a radial distribution of the dielectric permittivity given by$r^2\epsilon \left(r\right)=\text{const}$or $k^2{r}^2={\nu }^2\text{=const}$, all optical modes of the system can also be separated into the TE and TM categories. For the TM wave$q(r)=-r^{2}l{n}^{\prime }\epsilon =2r$, while for TE we always have $q\left(r\right)=0$.

For a TM mode, (12) can be therefore reduced as

4.4 Boundary conditions and the mode matching

In previous sections the general analytical solutions were written in the domains, where the dielectric function $\epsilon (r)$ is defined by either the Bessel equation or the Euler equation. At each interface between such domains, the electromagnetic field should satisfy the standard boundary conditions for the continuity of the tangential components of the magnetic and electric fields. Furthermore, the physical requirement for the behavior at the infinity and at the origin, introduce additional constrains to the solution, as described earlier in Section 4.2

In particular for a nonmagnetic, cylindrically symmetric medium without any free charges the condition of continuity of electromagnetic fields at a given interface $r=r_i>0$ reads

TE:

If the function f describes the scattered field (i.e. does not include the incident wave), it should obey the Sommerfeld radiation condition, $\underset{r\to \infty }{\lim }\sqrt{r}(f^{\prime }-\iota kf)=0$. The asymptotic behavior of the Hankel functions$H_m^{\pm }(kr)\sim r^{-1/2}{e}^{\pm \iota kr}$ for $r\to \infty$ forces us to choose only one term in (13); hence, for the scattered field we have

4.5 Incident light: the Gaussian beams and plane wave

While the angular momentum formalism introduced above, represents the “natural” choice for a system with rotational symmetry, to complete the description of the scattering problem, we also need to develop the angular momentum expansion of the incident field. Here, we will consider the most common examples of the illumination field: the plane wave and the Gaussian beam.

Starting from the classical generating function for the Bessel functions [12]:

While plane wave illumination represents perhaps the most common setup of a scattering problem, it does not convey the case when the target is illuminated by a focused beam. The latter is generally adequately described by the Gaussian beam [13], which yields an accurate description of the beam focused by linear optical elements. However, one needs to keep in mind that the Gaussian beam in its standard form is not an exact solution of Maxwell’s equations in the free space, but instead corresponds to the paraxial approximation [13]. To avoid the approximate nature of this approach and resulting inaccuracy, we consider that the field, while following the Gaussian profile at the waist cross section, is in fact the exact solution of the Helmholtz wave equation with the given direction of propagation [14]:

Adding the arbitrary shift of the focus location by straightforward translation, $x\to x-x_0,y\to y-y_0$, the final expression of a focused beam in Cartesian coordinates is

5. Layered Systems

Practical realizations of concentric cylinder optical and optoelectronic devices often involve a design that includes a number of (homogeneous) layers. In particular, external layers may be incorporated as a protective cover, while an internal layer between the absorber and the shell of the optical concentrator provides additional mechanical support or works as a channel for a liquid absorber. In this Section, we develop the theoretical description of light scattering and absorption by such layered systems.

5.1 Cascading the cylindrical layers

The boundary value problem of the scattered cylindrical modes inside a set of concentric cylindrical layers leads to a system of linear equations. Utilizing a direct linear algebra method for obtaining the scattered cylindrical modes inside layers with losses or axially symmetric distribution of the dielectric function could be problematic due to typically large condition numbers of the system matrices. The algorithm that we use includes a normalization step, making it possible to improve poor conditioning. However, if the algorithm were used for getting analytical formulae the normalization would not be necessary.

We consider a concentric cylindrical device having $l-1$ layers, $r_i<r<r_{i+1}$, $i=\overline{1,l-1}$, with the “outer” ($r>r_s$) layer filled with a host media, as shown in Fig. 2 .

 

Fig. 2 Schematic geometry of the layered structure.

Download Full Size | PDF

The general solution of the wave equation in i ^th layer can be represented as a sum of angular harmonics, $f^i(r,\varphi )={\displaystyle {\sum }_{m=-\infty }^{+\infty }{f}_m^i(r)\exp \iota m\varphi }$ where the radial part $f_m^i(r)$consists of two partial solutions $R_m^i(r)$ and $T_m^i(r)$with yet undefined coefficients. These coefficients can be obtained from appropriate $2l$ boundary conditions. All cylindrical modes are linearly independent, thus each of the boundary conditions can be separated into the individual equations for each m ^th mode, similar to that shown in (14) presented in a more general form,

Furthermore, the total field in the “outer” layer can be separated into the “incident” and “scattered” parts. The Sommerfeld radiation condition $R_{(m)}^l(r)=H_m^{+}(k_{0}r)$ is imposed on the scattered field (cf. Eq. (23)), while for the incident wave ${\tau }_{(m)}^l=a_m{J}_m({\overline{r}}_s)$, with the values of the coefficients $a_m$ taken from the expansion of the incident wave (25), (31).

To obtain the rest of $2l-2$ equations we use the standard boundary condition (20), (21) for TE polarization at each interface $r=r_{i+1}$,

For the TM polarization the only change is that $L_X^{i,0}\to \left({\epsilon }_{i+1}/{\epsilon }_i\right)\left(r_{i+1}\right)L_X^{i,0}$. Inverting the matrix at the left-hand-side of (33) and evaluating its determinant as $d^{i+1}=L_R^{i,1}-L_T^{i,1}=\frac{W[T^{i+1},R^{i+1}]}{T^{i+1}{R}^{i+1}}\left(r_{i+1}\right)$ we obtain

The resulting system obtained after stating all 2l boundary conditions has the form:

The system (37) can be solved in two sweeps, similar to the standard Thomas algorithm [17]. In this approach, the auxiliary normalized coefficients ${\overline{\tau }}^i={\tau }^{i+1}/{\tau }^i$, and ${\overline{\rho }}^i={\rho }^i/{\tau }^i$, $\left(i\in \overline{1,l-1}\right)$ are calculated by forward sweep using the recurrence (38)

The above algorithm includes the normalization of the partial solution in each layer by its value on the inner interface (the field in the core layer is an exception and is normalized by its value on the outer interface $r=r_2$, the latter is achieved by substitution $r_1=r_2$ in Eq. (32)). The normalization significantly improves the accuracy of the mathematical operations, performed numerically with a limitation of floating-point representations and arithmetic. Finally, the procedure yields the m ^th mode of the field by using (32); the solution to the entire problem in the i ^th layer is given as a sum of the angular momentum modes, $f^i\left(r,\varphi \right)={\displaystyle {\sum }_{m=-\infty }^{\infty }{f}_m^i\left(r\right)e^{\iota m\varphi }}$.

5.2 Scattering and absorption efficiencies

For any optical system, scattering and absorption efficiencies have the physical meaning of the relative size of the equivalent “black body” target which results in the same amount of scattered and absorbed power [18].

First, we recapitulate that the total field $E=\widehat{z}f$ (TE) or $H=\widehat{z}f$ (TM) outside the device at the external boundary of the shell, $r=r_s$ is defined as ($e^{-\iota \omega t}$ is omitted)

Equations (40)-(43) are used for the time-averaged Poynting vector $S^{abs}=\frac{1}{2}\mathrm{Re}\left(E\times H^{*}\right)$, and its normalized radial component $S_r^{abs}=\widehat{r}\cdot S/W$, where W is the total power transmitted by the incident wave through $2r_s$–wide strip per its unit length:

After applying (41) to (46), we note that $\mathrm{Im}\left(f_m^{*}{f^{\prime }}_m\right)$ includes a vanishing sum of a complex number with its conjugate, $\mathrm{Im}\left[{\left(c_m^{+}{H}_m^{+}{H}_m^{+}{}^{\prime }\right)}^{*}+c_m^{+}{H}_m^{+}{H}_m^{+}{}^{\prime }\right]=0$, and another term, $\mathrm{Im}\left[{\left(H_m^{-}{H}_m^{+}{}^{\prime }\right)}^{*}+{\left|c_m^{+}\right|}^2{H}_m^{-}{H}_m^{+}{}^{\prime }\right]=-\left(1-{\left|c_m^{+}\right|}^2\right)\mathrm{Im}\left(H_m^{-}{H}_m^{+}{}^{\prime }\right)$, so that

Similarly the scattering efficiency per unit length is calculated by taking the flux of normalized power of the scattered field (44) moving away from the device, $Q^{sca}={\displaystyle {\int }_0^{2\pi }{S}_r^{sca}{r}_{s}d\varphi }$, where likewise (45), $S_r^{sca}=\frac{1}{8r_s}\mathrm{Im}\left(f_{sca}^{*}\frac{\partial f_{sca}}{\partial {\overline{r}}_s}\right)$; and we obtain

5.3 The ideal electromagnetic black hole

We consider the ideal “fundamental” (p = 2) electromagnetic black hole, see (1), (2), (3). It is a three layered system with a gradient-index shell and absorbing core, for which the formalism of Section 5.1 can be applied to calculate the exact full-wave solution. Fig. 3 shows an example simulation of an ideal black hole with ${\epsilon }_s=\text{2}.\text{1}$, $r_s=\text{2}0\mu \text{m}$, and$r_c=\text{8}.\text{367}\mu \text{m}$. The device is illuminated with a Gaussian beam (free-space wavelength $\lambda =\text{1}.\text{5}\mu \text{m}$and full width $w=\text{2}\lambda$), which is focused at fixed $x_0$, and different $y_0$ (where $x_0=0$, and $y_0/r_s=\text{1}.\text{5},\text{1},0.\text{75}$, and 0).

 

Fig. 3 (Media 1) Ideal black hole with ε_s = 2.1, ε_c = 12, γ_c = 0.7, r _s = 20 µm, and r _c = r _s(ε_s/ε_c)^½ = 8.367 µm. The Gaussian beam with free-space wavelength λ = 1.5µm and full width w = 2λ is focused at x₀ = 0, and (a) y₀ = 1.5r _s; (b) y₀ = r _s; (c) y₀ = 0.75r _s, and (d) y₀ = 0. Range of modes used in series from (a) to (d): [-50;50], [-140;-40], [-170;-70], [-230;-130]

Download Full Size | PDF

Furthermore the approach is also used to determine the final form of the absorption efficiency. For this we use Eqs. (47), (48) and derive the scattering coefficients. We write the total field outside the device as$f_m^3\left(\overline{r}\right)=H_m^{-}\left(\overline{r}\right)+c_m^{+}{H}_m^{+}\left(\overline{r}\right)={\tau }_m^3{H}_m^{-}\left(\overline{r}\right)/H_m^{-}+{\rho }_m^3{H}_m^{+}\left(\overline{r}\right)/H_m^{+}$, where ${\tau }_m^3=H_m^{-}$, ${\rho }_m^3=c_m^{+}{H}_m^{+}$. Then the exact formula for the scattering coefficients is given as

The coefficients $w_{ij}^k$ are readily available from (35), then writing the explicit expressions for the matrix elements to (51) we obtain:

For reshaping ${\eta }_m$ we use the identity $-\iota \left(1-X\right)/\left(1+X\right)=\tan \left(\frac{\iota }{2}\ln X\right)$. Similarly in the expression for X, we use the identity $\left(1+\iota Y\right)/\left(1-\iota Y\right)=\exp \left(\iota 2{\tan }^{-1}Y\right)$ so that coefficients ${\eta }_m$ can be written for both polarizations as

5.4 The lamellar electromagnetic black hole

The formalism of a layered system for example can be used to study a non-ideal lamellar “black hole” optical concentrator and absorber, which approximates an ideal device (${\epsilon }_s=\text{2}.\text{1}$, $r_s=\text{2}0\mu \text{m}$, and$r_c=\text{8}.\text{367}\mu \text{m}$) with a different number of homogeneous layers, l. The device is illuminated by a plane wave with free-space wavelength $\lambda =\text{1}.\text{5}\mu \text{m}$; for l equal to 3, 5, 9, and 17. Scattering and absorption efficiencies are given in (49)-(50), the latter, $Q^{abs}$, yields 72%, 84%, 90%, and 94% respectively. Fig. 4(e) depicts the field map of the ideal black hole with smooth gradient, where 99% absorption efficiency is achieved.

 

Fig. 4 Simulated field maps for a lamellar “black hole” optical concentrator and absorber (ε_s = 2.1, r _s = 20 µm, and r _c = 8.367 µm) vs. the total number of layers, l. The device is illuminated by a plane wave with free-space wavelength λ = 1.5µm. (a) l = 3, and 72% absorption efficiency, (b) l = 5, and 84% absorption efficiency, (c) l = 9, and 90% absorption efficiency, (d) l = 17, and 94% absorption efficiency. Panel (e) depicts the reference case of the ideal black hole with smooth gradient, and 99% absorption efficiency.

Download Full Size | PDF

6. Semiclassical description

The primary objective of the present section is the derivation of the semiclassical expansion (${\overline{r}}_s\gg 1$) of the scattering and absorption cross sections of the “fundamental” (p = 2) electromagnetic black hole, introduced in Section 2.

Intuitively, from the ray dynamics of Section 2, we expect the absorption cross section to be close to the “geometrical” limit, so that $Q^{abs}\to 1$. The ray-optical treatment of Section 5 however gives no quantitative measure of how closely this limit is actually approached in a practical device. On the other hand, the exact description of Section 5, while giving a precise quantitative account, provides no compact lucid description of the system behavior. In contrast to this behavior, the semiclassical approach, mathematically corresponding to taking the formal limit ${\overline{r}}_s\gg 1$ of the exact solution, combines the clear physical picture of the actual dynamics in terms of the underlying ray trajectories, with the quantitative description of the interference and diffraction in the system.

The semiclassical approximation for the reflection coefficients and the absorption cross-section can be obtained by taking the limit ${\overline{r}}_s\gg 1$ in (52) and (54). Mathematically, this corresponds to the famous Debye expansion [12] of the Bessel functions

Physically, the angular momenta summation in absorption efficiency (47) can be separated into two parts – the sum over $\left|m\right|<{\overline{r}}_s$ and the sum over $\left|m\right|>{\overline{r}}_s$. As the dimensionless angular momentum m has physical meaning of the product of the (fixed) wavenumber and the impact parameter – see Fig. 2 – the sum over $\left|m\right|<{\overline{r}}_s$ corresponds to the rays that actually hit the target and as a result experience strong scattering. On the other hand, the sum for $\left|m\right|>{\overline{r}}_s$ takes into account the rays that miss the object. In the geometrical-optics limit, those rays experience no scattering and thus do not contribute to the cross-section. Beyond ray optics, the waves with $\left|m\right|>{\overline{r}}_s$, do experience some amount of scattering normally associated with diffraction, however in the limit ${\overline{r}}_s\gg 1$ the resulting contribution is exponentially small. Therefore, in the semiclassical limit ${\overline{r}}_s\gg 1$the sum over large angular momenta ($\left|m\right|>{\overline{r}}_s$) can be neglected, and we obtain

Furthermore, for ${\overline{r}}_s\gg 1$ the relative change in $\left|c_m^{+}\right|$, when one goes from m to $m+1$, is reasonably small – leading to the relative error $\sim 1/{\overline{r}}_s$– and thus the summation in Eq. (56) can be approximated by the integral

Note that in this approximation the values of the parameter ${\beta }_m$ that are different for the TE and TM polarization (cf. Eq. (53)), should be considered equal, i.e. $\sqrt{{\overline{r}}_s^2-m^2-1}\approx \sqrt{{\overline{r}}_s^2-m^2}$, then, applying Debye approximation to (52) and using identity $\left(1+\iota Y\right)/\left(1-\iota Y\right)=\exp \left(\iota 2{\tan }^{-1}Y\right)$ we obtain,

Then, dropping the small term $\frac{1}{2}{\overline{r}}_c{}^2/\left[{\beta }_m({\overline{r}}_c^2-m^2)\right]$ and taking $A\sqrt{{\overline{r}}_c^2-m^2}/{\beta }_m\approx 1$ (due to small losses) in (61), we are back to (58) and obtain,