Metal–dielectric metamaterials for transformation-optics and gradient-index devices in the visible regime

We investigate the application of metamaterials made of dielectrics with embedded metal nanoparticles for transformation-optics and gradient-index devices. We demonstrate that the permittivity and permeability distribution required for cylindrical cloaking at optical frequencies can be obtained by a composite material made of silver spheres embedded in a PMMA host. We show that this approach can be applied also to devices with an isotropic permittivity distribution and illustrate this for the optical black hole. The key ingredient of our system is the tunability of the plasmonic interactions between the metallic nanoparticles, depending on distance and size. In addition, we discuss possible realizations of the proposed structures.


Metamaterials and transformation optics
The research on metamaterials over the last decade has unveiled some material properties that had been assumed to be inaccessible in nature: the realization of negative refractive index materials from microwave to optical frequencies [1][2][3][4][5][6][7][8][9] as well as applications, e.g. for imaging beyond the diffraction limit [10][11][12][13][14][15]. Furthermore, as a consequence of Maxwell's equations being invariant under coordinate transformations, metamaterials enable physicists to manipulate electromagnetic waves in almost any fashion by tuning the refractive index to desired values. This opened up the path for the emerging field of transformation optics [16,17]. It was shown theoretically that transformation optics, in principle, allows the implementation of many intriguing features, such as invisibility devices [16,18,19], superlenses [20,21], field rotators [22], object shrinkers [23], Luneburg [24,25] and Eaton lenses [25,26] or even spacetime-related objects such as omnidirectional light absorbers [27][28][29][30] (the so-called 'optical black hole' (OBH)) and wormholes [31]. However, owing to the highly anisotropic materials that are required, most of the objects mentioned above have only been experimentally demonstrated in the gigahertz to terahertz regime. For example, a cylindrical cloaking device [32] and an OBH [33] were demonstrated for microwave frequencies by utilizing concentric layers of split-ring resonators or electric-field-coupled resonators. Nevertheless, especially for these two devices, the terms 'invisibility' and 'optical' are usually related to the visible frequency spectrum, although these high frequencies cannot be accessed with conventional resonator arrangements, for example due to magnetic-response saturation [34].

Metal-dielectric metamaterials for the visible regime
For the realization of a cylindrical cloak for TM-polarized radiation in the visible regime, Cai et al [35] recently proposed a three-dimensional (3D) metal-dielectric metamaterial consisting of arrays of radially elongated metal wires of subwavelength dimensions that are embedded in a dielectric host material. In addition, Xie et al [36] proposed a similar cloaking device, 3 but only extended in two dimensions. In these proposals, the metal wires have the form of oblate ellipsoids and can generate the necessary permittivity and permeability distribution, provided that they are all aligned in the radial direction. The underlying theoretical model that describes the permittivity of the composite material in these two publications is based on the Maxwell-Garnett effective-medium theory (MGT) [37]. In this model, the shape of the metallic inclusions is taken into account by the depolarization factor; however, due to the point-dipole approximation used to derive the MGT formula, interactions between the metal particles and plasmon excitations are not considered [38,39]. In the work of Cai et al and Xie et al, such effects can be neglected with good reason owing to the low metal filling fractions and the elongated shape of the wires. However, this shape and the necessity to radially orientate the wires in the host medium are a major drawback with regard to an experimental realization.
In this paper, we extend the studies to arrays of spherical nanoparticles. Concerning the preparation process, the major advantage of spherically shaped inclusions in comparison with wires would be the fact that they are symmetric and consequently do not need to be specifically orientated within the medium. In the approach presented, we investigate silver spheres included in a host medium of polymethylmethacrylate (PMMA) and exploit the plasmonic excitations that emerge when the metal spheres are very close together. Garcia de Abajo and co-workers calculated theoretically that the plasmon resonances can be engineered by varying the distances of two metal particles from each other and hence the interaction strength between them [38,40]. We show that with metal nanoparticles we can tune the effective permittivity of the composite material to desired values. We demonstrate by full-wave simulations, taking all interactions and excitations into account, that the required permittivity and permeability distribution for cloaking at optical frequencies can be obtained by using spherical metallic nanoparticles with typical diameters up to 20 nm. In addition, we show that the outlined concept can be applied to realize various other transformation-optics and gradient-index devices. We demonstrate this for the OBH.

Theoretical background
Transformation optics provides the possibility of controlling the path of light by applying a transformation to a coordinate grid and translating the rendered grid into equivalent permittivity and permeability tensors. In the case of the cylindrical cloak (see figure 1(a)), using cylindrical coordinates, this coordinate transformation is given by where a and b are the inner and outer radii of the cloak, and r, φ, z and r , φ , z are the transformed and untransformed coordinates, respectively. This transformation squeezes the circular area of radius b into a ring of thickness b − a and leads to diagonal permittivity and permeability tensors in the cloaking shell described by [16] (2)-(4), with the coordinate system having cylindrical coordinates r , φ and z. (b) Components of the permittivity and permeability tensors for a cloak with ideal (orange) and reduced parameters (blue).
Problematic for the realization of a material obeying such optical parameters are the radial dependence of every component and the divergence of φ and µ φ as r approaches the inner radius a. However, if we limit ourselves to an incident wave where the magnetic field vector is parallel to the z-axis ( H = (0, 0, H z ) t , i.e. the electric field vector oscillates in the crosssectional plane of the cloaking shell; in the following we label this polarization as 'TM') and take an impedance mismatch at the outer boundary into account, it was shown that the permeability dependence can be entirely removed [18,35]. In that case only µ z , r and φ enter Maxwell's equations and this so-called reduced-parameter set reads µ z and φ are constant throughout the reduced cloak, with r varying between 0 and 1 from the inner to the outer perimeter. Both the ideal and reduced parameters are shown in figure 1 Other possible applications are devices requiring a gradient-index distribution. In this paper, we exemplarily discuss the OBH, a broadband omnidirectional light absorber. This device possesses an isotropic permittivity, which is described by [29] (r ) = where a and b denote the inner and outer radius, s the permittivity of the surrounding and c the complex permittivity of the absorbing core. With our approach of a metal-dielectric 5 metamaterial we show that it is possible to fabricate an OBH that is suitable for rather high c , thus mimicking, e.g., a solar cell made of silicon.

Simulation method for obtaining the effective material parameters
In this section, we describe how the effective material parameters of the composite of PMMA and the silver nanoparticles have been determined. We assume that the cylindrical cloaking shell is composed of unit cells consisting of one silver sphere included in a PMMA host medium (figure 2). We remark that although the shape of one unit cell is cylindrical, it can be approximated as a rectangular prism since the curvature is small (cf figure 2(b)). This is due to the small dimensions of the unit cell compared to the circumference of the cloaking shell. We performed three-dimensional (3D) full-wave simulations utilizing the commercial software CST Microwave Studio 2 to obtain the complex reflection and transmission coefficients S 11 and S 21 . To take the interaction between the silver spheres into account, we chose periodic boundary conditions in the directions orthogonal to the k vector of the incident light. The dispersion relation of silver is a linear interpolation on experimental data from Palik [41]. The permittivity of PMMA has been set to PMMA = 2.25. We used two independent simulation configurations, one with incident TM-polarized light propagating along the φ-direction (cf figure 2(c)) and one with TM-polarized light propagating along the r -direction (cf figure 2(d)). As explained later, this allows us to independently derive the φ and r components of the permittivity tensor. In each simulation configuration, the extensions of the unit cell in the r -, φand z-directions have been changed with a parametric sweep. We then analyzed the behavior of S 11 (reflection) and S 21 (transmission) with respect to the extension of the unit cell. In figure 3(a), as an example we show a density plot of the magnitude of the transmission coefficient |S 21 | dependent on the extension of the unit cell in the r -direction and frequency f when the incident wave is propagating in the φ-direction (cf figure 2(c)). The diameter d of the nanoparticle has been assumed to be d = 12 nm. We observe a pronounced minimum branch in the transmission that shifts to higher frequencies when the dimension l r , and hence the distance of the silver spheres in the radial direction, is increased. Furthermore, for small l r , a second branch, occurring at higher frequencies but with weaker intensity, can be seen (the emphasized region in figure 3(a)). As already outlined in section 1.2, the shape of both branches can be traced back to plasmonic resonances, which have been analyzed on metal spheres in previous publications. This resembles the plasmonic interaction between two metal spheres as discussed by Garcia de Abajo and co-workers. For example, a graph similar to figure 3(a) can be found in [38]. Figure 3(a) depicts that the resonance is very pronounced and the frequency is strongly distance dependent when the particles are very close together. For larger particle distances (in the example of figure 3(a) for l r > 18 nm, corresponding to an inter-particle distance in the r -direction of 6 nm), the curve is very flat, i.e. the resonance frequencies are almost constant when the distance of the particles increases. This demonstrates that the excitations along one direction are dominant for the behavior of the particle array if one dimension out of l r , l φ , l z is significantly smaller than the others. Figure 3 f 3 = 6.7 × 10 14 Hz (green curve) as l r increases. At small distances in the r -direction, e.g. in the l r = 12.5 nm case, a further resonance occurs at f 4 = 6.9 × 10 14 Hz (red arrow). It vanishes for higher values of l r , i.e. the curve with l r = 15.1 nm shows no second resonance.
To determine the effective material parameters, we made use of the parameter-retrieval method of Smith et al [42,43]. With knowledge of S 11 and S 21 we calculated the effective refractive index n ret and the effective impedance z ret of the unit cell. From these values we obtained the effective permittivity and permeability. We remark that, if the light propagates in . At low distances in the r -direction between the spheres (l r = 12.5 nm), a further resonance occurs at f 4 = 6.9 × 10 14 Hz (red arrow). It vanishes for higher values of l r (e.g. the transmission curve with l r = 15.1 nm shows no second resonance). Increasing l r shifts the resonance to higher frequencies. (c) Retrieved radial permittivity r,ret with respect to the extension of the unit cell in the r -direction for l r = 12.5 nm (red curve), l r = 15.1 nm (blue curve) and l r = 18 nm (green curve). Like the transmission curves, they also show a pronounced resonance that shifts to higher frequencies if l r increases.
the unit cell along the φ-direction (simulation configuration shown in figure 2(c)), the electric field vector points into the r -direction. Therefore, in this case, the parameter retrieval yields the radial permittivity r,ret of the permittivity tensor. If, on the contrary, we let the light enter from the r -direction, i.e. the electric field vector points in the φ-direction (simulation configuration shown in figure 2(d)), we obtain the tangential permittivity component φ,ret of the permittivity tensor. In figure 3(c), the retrieved r,ret , obtained from the curves of |S 21 | in figure 3(b), are shown. At the frequencies of the transmission minima, they show a pronounced resonance that also shifts to higher frequencies if l r increases. For all the investigated cases, the permeability has been very close to 1, as desired (not shown). We point out that the individual simulations contain only one unit cell in the propagation direction of the wave (cf figures 2(c) and (d)). To make sure that the effective permittivity is not affected by the interaction between adjacent cells in the direction of propagation, we also performed simulations with several unit cells. We find 8 that the influence of these additional unit cells is relatively small, in particular at the frequencies of interest.

The cloak
We used the software COMSOL Multiphysics 3 to model the cloak with (i) continuous-ideal, (ii) continuous-reduced and (iii) discrete-reduced parameters (figures 4(a)-(c), respectively). In detail, we performed a 2D simulation in the x-y plane, modeling the cloak with perfectly matched layer (PML) boundaries around the simulation area at a frequency of f = 6.65 × 10 14 Hz, corresponding to a wavelength in the visible of λ = 451 nm. The outer radius of the cloaking shell has been assumed to be b = 2 µm and the inner radius a = 1 µm. The permittivity of the surrounding area has been set to vacuum with s = 1. The material in the cloaked region was assumed to be a perfect conducting (PEC) cylinder. Figure 4(a) shows the E y field induced by a TM-polarized plane wave, which excites the ideal cloak with its continuous material parameters (equations (2)-(4)) as well as the power flow, which is guided smoothly around the cloaked area. Figure 4(b) presents the cloak excited by a TM-polarized plane wave when applying reduced continuous material parameters to the cloaking shell (equations (5)- (7)). The power flow behind the cloak is slightly less than in the case of ideal material parameters, which results from the impedance mismatch of the reduced cloak. To come closer to a possible realization, we assumed the cloak to be composed of n cylindrical layers of equal thickness (b − a)/n (cf figure 2(a)). We substituted the continuous-reduced parameters with discretereduced parameters within each layer. For this purpose, we calculated the permittivity values in the middle of each layer and set it constant within the layer. We varied the number of layers and compared the cloaking effect to that of the reduced-continuous case. We found that a good cloaking effect is already obtained when using n = 8 layers. In figure 4(c), the reduced cloak with discretized material parameters (eight layers) is displayed. As in the previous case, a TMpolarized plane wave excites the discretized, reduced cloak. The E y field is disturbed slightly more than in the reduced-continuous case with the power flow damped a bit more. In table 1, columns r and φ show the discrete values that have been applied to the individual layers in this simulation.
To realize these constant permittivity values within each layer, we simulated in a fourth step a unit cell with PMMA host medium and an inclusion of a silver sphere as outlined in section 3. We kept the length in the z-direction (l z = 20 nm) and φ-direction (l φ = 28 nm) as well as the diameter of the silver sphere (d = 12 nm) constant. Using the simulation configuration with incident plane waves in the φ-direction (figure 2) and changing the length of the unit cell in the r -direction from l r = 12.5 to l r = 18 nm leads to a change of the radial permittivity r,ret , as shown in figure 3(c). For example, at the frequency f = 6.65 × 10 14 Hz we find a change in the radial permittivity from 0 to 1 for values of l r between 12.56 and 13.37 nm (table 1, column l r ). We chose eight values of the radial permittivity, which coincide best with the values that we calculated as the constant values within each layer in the discrete-reduced case. The simulation configuration with the incident plane wave in the r -direction gave us the change in tangential permittivity φ at these values of l r . It turned out that the tangential permittivity varies only slightly, with an average value of φ = 4.09. We remark that from equation (6)  a and b of the cloaking shell can be adjusted to desired values of φ . For example, if we fix the inner radius of the cloaking shell at a = 1 µm, we can calculate for the average φ of φ = 4.09 an outer radius of b = 1.98 µm. Substituting the retrieved values of the radial and tangential permittivity into COMSOL, we get a good cloaking effect, which is shown in figure 4(d). In particular, inside the cloaking shell the field is still guided smoothly around the inner core. Table 1. Calculated radial and tangential permittivity values r and φ for the discrete, reduced cloak composed of n = 8 cylindrical layers. The radial permittivity values r have been assumed to be constant within each layer and were calculated using equation (5) for each layers middle. The retrieved radial permittivity values r,ret have been calculated at a frequency of f = 6.65 × 10 14 Hz (λ = 451 nm). The l r values are the corresponding extensions of the unit cells in the r -direction. l z = 20 nm and l φ = 28 nm have been set constant. The retrieved tangential permittivity φ,ret changes only slightly, with an average of φ = 4.09. The metal filling factor is denoted by t. The corresponding COMSOL simulation is shown in figure 4(d). However, the small interval of change in l r , and therefore the almost identical metal filling factors t, is a problem in realizing the cloak. Therefore, we investigated a unit cell whose extension is changed in the l r and in the l z extension simultaneously. In table 2, we provide an example at a frequency of f = 6.34 × 10 14 Hz (corresponding to a wavelength of λ = 473 nm) where the permittivity distribution required for cloaking is satisfied when the r -extension ranges from 12.5 to 14 nm and the z-extension from 12.5 to 45.83 nm. The extension in the φ-direction and the diameter of the silver sphere are fixed at l φ = 28 nm and d = 12 nm, respectively. In this case with a variation of two parameters the change of l r is larger (≈0.5 nm), compared to the change of l r in the case of the discrete-reduced cloak with only one variable parameter (≈0.1 nm). Furthermore, the metal filling factor varies over a broader range from 0.05 to 0.21. The disadvantage of introducing a second variable parameter is that φ,ret varies over a broader range of values (from 2.35 to 3.85). This is due to the fact that changing the z-extension of the unit cell affects r,ret as well as φ,ret , while we found that changing the r -extension has an effect on r,ret alone (compare with table 1). In the configuration presented, the average tangential permittivity is 2.91. Therefore, if a cylinder of radius a = 1 µm should be cloaked, equation (6) would lead to an outer radius of the cloaking shell of b = 2.42 µm. The corresponding COMSOL simulation is depicted in figure 4(e). It can be observed that although the spread in φ,ret is 3.85 − 2.35 = 1.50, the cloaking performance is not significantly altered.
In a third approach, we kept the extension in the φand r -direction constant at l φ = 28 nm and l r = 12.5 nm, varying only the z-extension of the unit cell. We found that at the frequency of f = 6.34 × 10 14 Hz the radial permittivity can also be changed from 0 to 1, comprising changes in l z from 12.5 to 28.79 nm (table 3). The corresponding filling factors vary between 0.21 and 0.09 from the inner to the outer layer. φ,ret varies from 2.94 to 3.83 with an average value of 3.55. The outer radius of the cloaking shell with inner radius a = 1 µm can be calculated Table 2. Simulations and calculations following the same way as described in table 1. In this table, we vary two parameters of the unit cell (l r and l z ). l φ was set constant at l φ = 28 nm. The operating frequency is f = 6.34 × 10 14 Hz (λ = 473 nm). This data set corresponds to the COMSOL simulation in figure 4( figure 4(f). The cloaking performance is similar to that in figures 4(d) and (e). This third approach combines the advantages of varying only one structure parameter (l z ) with a relatively large spread in the metal filling fractions. Therefore we conclude that this would be the easiest way for a possible experimental realization. The imaginary part of the composite material is about 0.1 throughout the cloaking shell, a value similar to that in [35]. As mentioned in section 3 the interaction between adjacent unit cells in the direction of propagation has only a minor influence on the retrieved parameters at the frequencies of interest. For example, for the cloak in table 1, the retrieved permittivities r,ret for the first and eighth layers change to −0.01 and 0.95 (0.06 and 1.10), respectively, when two (three) adjacent unit cells in the direction of propagation are used.
Comparing the far-field scattering patterns of the discussed cloaks to that of a bare PEC cylinder, we find that all cloaks reduce the scattering significantly. It turns out that the cloak where only the dimension l r is changed (table 1) has advantages at large scattering angles, whereas the other cloaks (tables 2 and 3) are advantageous for smaller scattering angles.

The optical black hole
The OBH has been modeled in the first and second steps with COMSOL, i.e. for continuous permittivity (figure 5(a)) and discretized permittivity ( figure 5(b)). The outer and inner radii have been set to b = 20 µm and a = 10 µm, respectively. The permittivity of the surrounding material is s = 2.1, chosen as that of silica. The absorbing core within the middle of the OBH has been assumed to consist of n-doped silicon (n ≈ 2.7 × 10 20 cm −3 ) and has a complex permittivity of c = 12 + 0.7i as outlined in [29]. In figure 5(a), the electric field intensity |E| and the power flow of a Gaussian beam incident on the continuous OBH are shown. The Gaussian beam is directed to the middle of the OBH and gets absorbed by the core. Figure 5(b) shows |E| and the power flow of a Gaussian beam passing a discretized OBH with six layers. As described above for the discretized cloaking device, we calculated the permittivity in the middle of each layer and set it constant within the layer. We obtained good agreement of the beam path with that of the ideal case already for six layers. Due to impedance mismatch, minor reflections occur at the boundaries of the layers. To obtain unit cell dimensions with effective permittivity components matching the discretized permittivity values, we simulated with CST Microwave Studio in a third step a symmetric unit cell (l r = l φ = l z ), since in the case of an OBH the permittivity is Table 4. Calculated permittivity values (r ) for the discrete OBH with equal change in l r , l φ and l z . The retrieved permittivity values (r ) ret have been obtained from the CST simulations of the unit cell of PMMA host medium with an inclusion of a silver sphere. The l r,φ,z values are the extension of the unit cells in the r, φ, z-direction, which belongs to the retrieved permittivity values. The metal filling factor is denoted by t.

Layer
(r ) (r ) ret l r,φ,z (nm) t  (8)). We like to mention that due to the isotropic permittivity tensor an OBH would be easier to realize than a cloaking device, which requires anisotropy. A symmetric unit cell fulfills the requirement that the permittivity value in every direction is constant. Varying the length of the unit cell equally in every direction, we obtained retrieved permittivity values for the layers of the OBH with values of l r = l φ = l z between 12.5 and 18.72 nm (table 4) at a frequency of f = 3 × 10 14 Hz (corresponding to a wavelength of λ = 1000 nm). The obtained values have been substituted into COMSOL. In figure 5(c), the intensity |E| and the power flow of the discretized OBH with the retrieved permittivity values are shown. It can be seen that the ray path is almost the same as that in figure 5(c). The reflections at the interfaces are slightly more pronounced, because the impedance mismatch at the inner layers has increased slightly. The major challenge is the permittivity values of the inner layers, which have to be rather high. In that case, one needs to be near the resonance frequency of the composite medium, but this imposes a high imaginary part, at which the OBH would lose its performance. However, at a frequency of f = 3 × 10 14 Hz the complex permittivity of the inner layer is r = 8.49 + 0.24i, so that the imaginary part is still reasonably low for a metal-dielectric system with high metal filling factors at the inner layers (see in table 4).
We note that in our example we adopted material values of the inner and outer materials similarly to [29]. In that work a semiconductor with a rather high permittivity (Re( c ) = 12) was chosen for the inner core to mimic a solar cell that can potentially be integrated into the structure, exploiting the radiation-attraction capability of the OBH. If one assumes a lower inner-core permittivity of Re( c ) = 4, e.g. as used in [44,45], the problem of unfavorable imaginary parts can be avoided. Since all permittivities have been chosen to be larger than one, the realization of an OBH is not limited to metal inclusions, but can also be achieved by a dielectric-dielectric composite (not shown), which would not cause absorption due to a metal compound. However, as mentioned above, such a dielectric-dielectric composite would be restricted to rather low inner-core permittivities. Large effective permittivities are a unique property of metal-dielectric composites since only there do plasmonic resonances occur.

Summary and discussion
In this paper, we propose and demonstrate the operation of the proposed transformation-optics and gradient-index devices for exemplary sets of parameters of the operating frequency, the diameter of the nanoparticle and the unit cell dimensions. These quantities can be the subject of a parameter sweep to find realizations at, e.g. a different operating frequency if this is desired. Also, to obtain a better coincidence with the continuous devices, the number of layers in the discrete cases could be raised. However, this would also increase the complexity.
In this context we would like to mention that the transformation-optics approach that is utilized in this paper is not the only strategy to describe a cylindrical invisibility cloak. For example, Chen and co-workers proposed an optimization theory [46] that is based on a fullwave scattering model of the cloak [47,48]. This ansatz additionally aims to minimize the scattering that is induced by the interfaces of the discretized cloak, similar to cloaks using plasmonic shells [49][50][51]. In that sense the cloak proposed by Chen and co-workers utilizes the advantages of both the scattering-cancellation and transformation-optics methods, leading to a good cloaking performance using only four discrete layers. However, this approach would require us to tailor three parameters (in the TM case r , φ and µ z ) instead of only one ( r ). In the case of the metal-dielectric composite that is studied in this paper it is essential that two parameters be constant throughout the cloak.
The changes of the unit-cell dimensions l r , l φ and l z are in the range of a few nanometers. This makes the realization of the discussed devices challenging. One possible approach for a realization is the self-assembly of metal nanoparticles. Using such a technique, Fan et al fabricated gold/silica composites where the gold nanoparticles of diameters below 10 nm have been arranged in a periodic face-center cubic lattice within a dense silica matrix. They describe that the cell dimensions are adjustable through control of the nanoparticle diameter and the tailorable properties of the materials used during synthesis [52]. Another promising concept was described in a recent patent where a vacuum deposition technique is introduced which allows the fabrication of nanoparticle arrays embedded in a polymer. If this deposition process is repeated it is possible to form 3D arrays with nanoparticle distances controllable on the nanometer scale [53]. Also, in the case of an isotropic permittivity tensor, the different values in the layers can be conveniently adjusted by the metal filling factor t, which, in the case of the OBH, varies over a broad range from t = 0.46 to t = 0.14.

Conclusion
In conclusion, we have proposed a way of realizing transformation-optics and gradient-index devices by choosing a metamaterial consisting of a host medium of PMMA and an inclusion of spherical metal nanoparticles. The variation of the extension of the unit cell increases or decreases the interaction strength of the particles, resulting in a shift of the permittivity at a given frequency of interest. We demonstrated possible realizations of a cylindrical cloaking device comprising anisotropic permittivity in the visible wavelength spectrum at wavelengths of λ = 451 nm and λ = 473 nm. Structure parameters for the realization of devices with isotropic permittivity have been suggested for an OBH at a wavelength of λ = 1000 nm.