Enhanced control of light and sound trajectories with three-dimensional gradient index lenses

We numerically study the focusing and bending effects of light and sound waves through heterogeneous isotropic cylindrical and spherical devices. We first point out that transformation optics and acoustics show that the control of light requires spatially varying anisotropic permittivity and permeability, while the control of sound is achieved via spatially anisotropic density and isotropic compressibility. Moreover, homogenization theory applied to electromagnetic and acoustic periodic structures leads to such artificial (although not spatially varying) anisotropic permittivity, permeability and density. We stress that homogenization is thus a natural mathematical tool for the design of structured metamaterials. To illustrate the two-step geometric transform-homogenization approach, we consider the design of cylindrical and spherical electromagnetic and acoustic lenses displaying some artificial anisotropy along their optical axis (direction of periodicity of the structural elements). Applications are sought in the design of Eaton and Luneburg lenses bending light at angles ranging from 90° to 360°, or mimicking a Schwartzchild metric, i.e. a black hole. All of these spherical metamaterials are characterized by a refractive index varying inversely with the radius which is approximated by concentric layers of homogeneous material. We finally propose some structured cylindrical metamaterials consisting of infinitely conducting or rigid toroidal channels in a homogeneous bulk material focusing light or sound waves. The functionality of these metamaterials is demonstrated via full-wave three-dimensional computations using nodal elements in the context of acoustics, and finite edge-elements in electromagnetics.


Introduction
There is currently a keen interest in transformation-based metamaterials, which offer new functionalities thanks to an enhanced control of acoustic and electromagnetic wave trajectories. The rapid growth of these emerging research areas is fueled by the electromagnetic paradigms of two invisibility cloaks independently proposed by Pendry et al [1] on the one hand and Leonhardt [2] on the other hand, in 2006. These two discoveries were an instant hit on the World Wide Web; this was also emphasized by a subsequent experiment chiefly achieved in the microwave regime by the former grouping [3]. However, the required parameters of such optical media are very difficult to achieve as the refractive index goes to zero or to infinity at some points. The need for the realization of invisible devices therefore prompted the work of Tyc and Leonhardt on the so-called transmutation of singularities [4]. After a transformation is applied, the refractive index distribution equivalent to the original one is non-singular but displays slight anisotropy, which is a weaker type of singularity where the permittivity and permeability tensors are discontinuous at one point. In order to illustrate their theory, Tyc and Leonhardt theoretically investigate in [4] the trajectories of light in invisible spheres such as the transformed Eaton lens which bends light by an angle of 180 • without any scattering, thereby making a perfect retroreflector [87]. Such a spherical lens described by the spatially varying refractive index n(r ) = √ ε r µ r = 2a r − 1 for r a and n = 1 for r > a, (1) was proposed 60 years ago by Eaton [5], where the relative permittivity ε r is spatially varying and the relative permeability µ r is assumed to be 1 (non-magnetic medium). It has been used for years in radar technology, even if its refractive index tends to infinity when r goes to 3 zero. Importantly, the lens is impedance matched to vacuum at the boundary r = a. Tyc and Leonhardt have solved this long-standing issue by noting that the above refractive index diverges like n ∼ r −1/2 near zero. Physically, this means that light slows down near the center of the lens. Tyc and Leonhardt's bold claim is to expand the space around the point r = 0 when the refractive index behaves as n(r ) = r p using the transform r = r p+1 , where r is the stretched radial coordinate. With this new coordinate system, it can be shown that the relative permittivity and permeability are now described by diagonal tensors ε = µ = Diag r 2 r 2 dr dr , dr dr , dr dr = Diag p + 1, in the neighborhood of the center of the lens. This means that we have now replaced the singular isotropic refractive index in the center of the lens by anisotropic tensors of permittivity and permeability, which in the case of the Eaton lens are associated with p = −1/2. We can see that anisotropy plays a positive role in the transmutation story. However, it is not always desirable to have anisotropy, which might be a genuine effect associated with some staircase approximation of a continuous and isotropic refractive index such as in (1). It is such a negative story that we wish to tell in this paper. Historically, such transformation techniques were first developed to simplify computational problems [6][7][8]. One area of tremendous use of such techniques is the modeling of open waveguides [9][10][11], whereby a domain of infinite extent is mapped onto a finite domain to be meshed, a computational method known as perfectly matched layers (PMLs) which was introduced by Bérenger [12] and Chew and Weeden [13] in 1994. A PML is an artificial absorbing layer for wave equations, commonly used to truncate computational regions in numerical methods to simulate problems with open boundaries, especially in the finite-element method and the finite difference time domain method [14]. The key property of a PML that distinguishes it from an ordinary absorbing material is that it is designed so that waves incident upon the PML from a non-PML medium do not reflect at the interface-this property allows the PML to strongly absorb outgoing waves from the interior of a computational region without reflecting them back into the interior. PMLs can be mathematically justified for quite general convex domains [15] and can also be extended to star domains [16]. Transformation techniques are also very useful regarding the design of complementary media [17][18][19][20][21][22][23], whereby space folding leads to cancelation of the acoustic or optical path, hence enabling the object and image planes of a negatively refracting lens to be mapped onto one another, or in invisibility cloaks, whereby a blowup of a point leads to a hole in acoustic or electromagnetic space, thus making this region invisible. One of several issues that are difficult to overcome is in particular singularities of the material parameters occurring on the invisibility cloak's inner boundary, where the eigenvalues of the rank-2 tensors of permittivity and permeability (in electromagnetics) or rank-2 tensor of density and scalar compressibility modulus (in acoustics) tend towards zero or infinity.
Structured metamaterials based upon resonant structural elements [3] have gained popularity over the last few years due to their counter-intuitive physics as well as their scope for subwavelength imaging applications (through negative refractive effective index [24]) and invisibility cloaks (via artificial magnetism and anisotropic heterogeneous features). Transformational optics is a vibrant area in metamaterial research based on newly engineered heterogeneous anisotropic media [1,2], low-index materials [25] and of course negative 4 refractive index materials leading to cloaking via anomalous resonances [26]. Importantly, the mathematical foundations behind the scene have been known from researchers working in the area of inverse conductivity problems [27]. The first experimental realization of an invisibility cloak, chiefly achieved in the microwave regime [3], suggests that cloaking will be limited to a very narrow range of frequencies. However, it will not be perfect since the cloak is necessarily dissipative and dispersive, and some of its tensor components are singular on the inner boundary. This drawback can be overcome by considering a generalized transform [28] in an upper dimensional space and then projecting the resulting metric on the physical space, and this leads to non-singular tensors of permittivity and permeability. The approach of projective optics from upper dimensional space proves to be a very fruitful one, as it opens new avenues in planar and spherical transformation optics metamaterials [29][30][31][32][33][34][35][36][37]67]. It is even possible to conceive an ultra-thin planar metamaterial using graphene [38].

Wave equations for transformed electromagnetic and acoustic heterogeneous anisotropic media: similarities and discrepancies
The equations governing the propagation of acoustic and electromagnetic waves do not respond to a geometric transform in exactly the same way. While the vector Maxwell equations are known to be transformation invariant: (where ε and µ denote the rank-2 tensors of permittivity and permeability, J is the Jacobian of the transformation and ω is the wave frequency), the elastodynamic equations do not retain their form under a geometric transform in general [40,42,47,48].
Here, we focus on the simple case of scalar pressure waves propagating in a fluid, and in this case we have where ρ denotes the rank-2 tensor of density and κ is the (spatially varying) compressibility modulus. One should note that the compressibility modulus remains a scalar under a geometric transform, unlike for the density, permittivity and permeability. The paradigm of transformation optics is the spherical invisibility cloak, which is an anisotropic heterogeneous coating in the region R 1 < r < R 2 , whose tensors can be expressed as [1] (in a spherical basis) New Journal of Physics 14 (2012) 035011 (http://www.njp.org/)

5
The paradigm of transformation acoustics is the silence cloak, whose tensor of density and scalar compressibility can be expressed as (in a spherical basis) [41][42][43][44][45] ρ −1 = Diag One can clearly see the discrepancies in the optical and acoustical cloaks' parameters. These two examples motivate the need for comparisons between optical and acoustical metamaterials. Even if the transformation rules in (3) and (4) look very similar up to a matrix inversion for the tensors of permittivity, permeability and density, this has important physical consequences, e.g. the radial component of (5) vanishes at r = R 1 , while its counterpart in (6) is infinite at r = R 1 (as is the compressibility). The interpretation of anisotropy thus depends upon the physical context.
One should note in passing that the transformation rules (3) and (4) give very different features in cylindrical coordinates (r, θ, z). Indeed, the cylindrical cloak in optics displays the following tensor of permittivity and permeability: whereas the acoustic cloak displays the following tensors of density and scalar compressibility [46]: which shows that one needs to be very careful with correspondences between electromagnetic and pressure wave equations in two and three dimensions. There seems to be a contradiction between the expressions of (7) and (8). Actually, the transformed equation associated with (3) takes the following form for cylindrical domains: with H z the longitudinal magnetic field and ε −1 HereT denotes the upper diagonal part of the transformation matrix T = J T J det(J −1 ) and T zz its third diagonal entry. This paradox can be sorted out using the duality lemma derived in [66]:  Comparisons between effective permeability and density derived from the Maxwell-Garnett formula (MG), see equation (25) and homogenization theory (HT), see theorems 1 and 2, for a cylindrical inclusion with high permittivity ε 0 = 50 (resp. density ρ 0 = 50) in a background medium of permittivity ε 1 = 1 (resp. density ρ 1 = 1) (only the contrast matters between inclusions and surrounding matrix). We note that the higher the filling fraction f , the smaller the first two diagonal entries of the effective permittivity ε 11 = ε 22 (resp. density ρ 11 = ρ 22 ) and the larger the third diagonal entry , hence the higher the effective anisotropy along the cylindrical axis. In contrast, the Maxwell-Garnett formula provides us with effective permittivity (resp. density) which is isotropic and whose value increases with the contrast. We note the strong correspondence with values of permittivity (resp. density) of the GRIN lens given by formula (24) in the fifth column and ε 11 = ε 22 (resp. ρ 11 = ρ 22 ) in the second column, unlike for ε 33 (resp. ρ 33 ) in the third column. Importantly, when the permittivity ε 0 (resp. density ρ 0 ) tends to infinity, ε 33 (resp. ρ 33 ) tends to infinity: this limiting case corresponds to the configuration of infinite conducting/rigid cylindrical inclusion. All the aforementioned results apply mutatis mutandis to the case of toroidal inclusions (in which case the effective permittivity and density tensors need to be expressed in the toroidal system of coordinates). Importantly, the expressions for effective permittivity and density are identical; see (10) and (11) in lemma 1.
This paradox between the expressions of cylindrical optical and acoustical cloaks is solved using this lemma in [49]. This lemma explains why the expressions or effective permittivity and density are the same as those in table 1.
Regarding the finite-element modeling of the equations for light and sound, one important feature is that the Maxwell equations are of course vectorial, but they also have another interesting feature. Indeed, taking the divergence in both members of (3), one sees that ∇ · (µ H) = 0, an equation also known as a constitutive or compatibility equation. If one uses nodal elements, it is customary to penalize the weak form of (3) by the positive term s D µ H · V, where s needs be chosen according to one's computational intuition, D is the computational domain (with, say, Dirichlet data on its boundary) and V is a nodal vector function. This leads to the famous problem of spurious modes, which has been solved by Bossavit [8] through the introduction of edge elements which naturally fulfill the compatibility relation: edge elements are discrete Whitney forms whose tangential components are continuous across interfaces on a tetrahedron mesh. They also behave in a nice way under geometric transforms, as these are differential one-forms [8].
The wavelength is of the same order as the fixed domain f , and hence very large compared to the small inclusions.
One may wonder what to do with the resulting transformed media as these have unusual properties not encountered in most materials: they are described by matrix-valued permittivity and permeability. Fortunately, homogenization theory shows us a route towards such metamaterials.

Homogenized wave equations for artificially anisotropic media
2.2.1. Homogenized Maxwell's system for anisotropic permittivity and permeability. Let H η be a sequence of locally square integrable (vector valued) functions over R 3 solutions of the sequence of problems where η is a small positive parameter (describing the typical size of heterogeneities, see figure 1), k 0 is the wavenumber in vacuum and ε(x, x η ) and µ(x, x η ), respectively, denote the relative permittivity and permeability of the medium equal toε(x, x η ) andμ(x, x η ) in a periodic structure f (known as metamaterial) and equal to 1 in R 3 \ f , i.e. in vacuum. We suppose that the sequence of fields H η , solutions of the sequence of problems (P H η ), has a two-scale expansion of the form where Our goal is to characterize the electromagnetic field when η tends to 0. If the functions H i do not increase too fast when η tends to 0, the limit of H η will be H 0 , a rougher approximation to 8 H η . Hence, we make the assumption that for all , so that the expansion still makes sense in the neighbourhood of 0. If the above expansion is relevant, we can state the following fundamental result [62]: Theorem 1. When η tends to zero, the H η solution of the problem (P H η ) converges weakly in L 2 ( f ) to the average of the first term of its asymptotic expansion on the basic cell Y , namely H hom (x) = Y H 0 (x, y) dy, which is the unique solution of the following problem: with the homogenized permittivity (15) and the homogenized permeability , are the unique solutions (up to an additive constant) which are square integrable with a square integrable gradient on Y of one of the following six problems (K j ) and (M j ) of electrostatic type: and Remark 1. Theorem 1 states that the homogenized permittivity and permeability derive from an electrostatic problem, which is both counter-intuitive (at first glance, one might expect that homogenized permeability derives from a magnetostatic problem) and reassuring (the transformed permittivity and permeability in (3) have exactly the same structure). This seemingly insignificant remark emphasizes that deep connections exist between the mathematics of transformation optics and homogenization theory.

Homogenized acoustic equation for artificial anisotropic density.
Let p η be a sequence of locally square integrable (scalar valued) functions over R 3 solutions of the sequence of problems where η is a small positive parameter (describing the typical size of heterogeneities), k 0 is the wavenumber in vacuum and ρ(x, x η ) and κ(x, x η ), respectively, denote the density and compressional modulus of the medium equal toρ(x, x η ) andκ(x, x η ) in a periodic structure f (known as acoustic metamaterial) and equal to 1 in R 3 \ f . We suppose that the sequence of pressure fields p η , solutions of the sequence of problems (P p η ), has a two-scale expansion of the form where Our goal is to characterize the pressure field when η tends to 0. If the coefficients p i do not increase 'too much' when η tends to 0, the limit of p η will be p 0 , a rougher approximation to p η . Hence, we make the assumption that for all , so that the expansion still makes sense in the neighbourhood of 0. If the above expansion is relevant, we can state the following fundamental result [64,65]: When η tends to zero, the p η solution of the problem (P p η ), converges weakly in L 2 ( f ) to the average of the first term of its asymptotic expansion on the basic cell Y , namely p hom (x) = Y p 0 (x, y) dy, which is the unique solution of the following problem: with are the unique solutions (up to an additive constant) which are square integrable with a square integrable gradient on Y of one of the following three problems (L j ) of acoustostatic type: Remark 2. Theorem 2 states that the homogenized anisotropic density derives from an acoustostatic problem, while the homogenized compressibility remains a scalar which is both counter-intuitive (at first glance, one might expect that the homogenized compressibility is also anisotropic) and reassuring (the transformed density and compressibility in (4) have exactly the same structure). This second seemingly insignificant remark emphasizes that deep connections exist between the mathematics of transformation acoustics and homogenization theory.
It should be noted that in theorems 1 and 2, the homogenized (matrix valued) coefficients are spatially varying, which is a feature allowed by so-called two-scale convergence [83,84]. The group of Greenleaf used this feature in order to approximate (heterogeneous anisotropic) acoustic and quantum spherical cloaks' parameters to any desired accuracy through the twoscale homogenization of concentric isotropic multilayers [85]. The same statement holds here for spherical electromagnetic and acoustic cloaks. Let us now describe the various geometric transformations considered in this paper and their connection to effective (anisotropic or not) media. We start with the design of cylindrical lenses.

Cylindrical gradient index (GRIN) lenses in optics and acoustics
In this section, we extend the design of dielectric cylindrical concentrators with toroidal air channels proposed recently in the context of electromagnetic waves [51] to the case of toroidal channels with infinite conducting boundaries (electromagnetic micro-waves) and toroidal inclusions with rigid boundaries (acoustic pressure waves). Gradient Index (GRIN) lenses are an old subject, which started with the paradigm of Maxwell's fisheye (discovered by the father of the eponym equations) more than one hundred and fifty years ago [52]. A century after Maxwell. Luneberg [53] discovered that the isotropic but spatially varying refractive index of the fisheye can be deduced from the stereographic projection of a sphere on a plane. Nearly half a century passed again before Leonhardt [54] theorised that Maxwell's fisheye contained within it all the ingredients for a perfect lens (up to an adequately placed mirror, which is reminiscent of the time-reversal experiments of Fink's group [86]). This claim led to some controversy, which renewed the interest in GRIN lenses [55-57, 68, 70, 71, 78] in conjuction with a rapid growth of metamaterials subject [80][81][82]. The ongoing debate as to whether or not negative refraction is necessary in order to achieve a perfect lens [24] is one of tremendous interest, and is still partly open. Importantly, the Maxwell fisheye is also a key ingredient in making a non-superluminal invisibility cloak [28,67], which is yet another reason to believe it is a past device promised to a great future. From now on, we adopt the notation (x, y, z) for the Cartesian coordinates, which were previously denoted by ( x 1 , x 2 , x 3 ), as there is no longer any ambiguity with rescaled variables. Importantly, all the devices in optics are assumed to be non-magnetic to ease a potential realization.
More precisely, we would like to consider the hyperbolic secant profile where R 0 is the radius of the cylindrical lens, and r = x 2 + y 2 . Moreover, n 0 and n R are, respectively, the refractive indices (invariant along the z-axis) at r = 0 (along the optical axis) and at r = R 0 (at the outer edge). Such a design can be deduced from quasi-conformal grids, as detailed for instance in the textbook [58]. We note that such a refractive index profile has already been implemented in the 2D optical and acoustical settings in [39,50]. We would like to propose a structured design of GRIN lens, which mimics the above transformed design. We consider a matrix of polymethylmethacrylate with a refractive index n 0 = 1.486 and we drill some toroidal channels of varying infinite conducting cross-section. We note that the radius of the channels r increases with the radius of lens R and thus the refractive index decreases with increasing radius in accordance with the table shown in panel (c) of figure 2.
Intuitively, the effective permittivity e (resp. density ρ e ) of the composite medium consisting of cylindrical inclusions should be reasonably well approximated by the classical Maxwell-Garnett formula for 3D composites [60,61]: where 0 (resp. ρ 0 ) and 1 (resp. ρ 1 ) are the permittivity (resp. density) of material inside the inclusions (here of high conductivity/rigidity) and background. We consider a cylindrical inclusion of radius r of length d in a cubic cell of sidelength d and hence the filling fraction f = 1 − πr 2 /d 2 . We take d = 200 nm in table 1. Importantly, the effective permeability is 1, as the structured medium is non-magnetic. Similarly, the effective acoustic medium has the same compressibility as the bulk surrounding the toroidal inclusions (assuming that the compressibility is the same in both media in order to emphasize correspondences between optics and acoustics, otherwise the effective compressibility κ e is the harmonic mean of the compressibility).
We computed the corresponding tensors of effective permittivity and density for the cylindrical structures shown in figures 3 and 4 in the following way: we first note that the Proposal for a cylindrical Luneburg lens for electromagnetic (EM) or pressure waves: (a) three-dimensional (3D) view of toroidal channels (dark blue) in a cylindrical bulk (light blue); (b) side view of (a) showing the geometric parameters R (distance between the center of the toroid to the center of the lens) and r (radius of a toroidal channel); (c) geometric parameters R and r of the lens, and corresponding effective refractive index (equal to the square root of the effective permittivity in the electromagnetic case, as it is non-magnetic); (d) numerical validation at wavelength λ = 700 nm for an EM gradient index (GRIN) lens with a thickness of 600 nm and a radius R 0 of 1000 nm, for the refractive index profile (23). Left: side view of the real part of the y-component of electric field (phase); middle: side view of the norm of the electric field; right: 3D plot of the norm of the electric field.
homogenized Maxwell system is simply described by an effective anisotropic permittivity and a constant permeability when the structured medium is a dielectric medium with infinite conducting inclusions [63], and in the same way the homogenized pressure wave equation is amply described by an effective anisotropic density for a periodic set of rigid inclusions in a homogeneous fluid [64]. The only thing to do in the annex problems K j and L j (see equations (17) and (23)) is to consider that ε and ρ take infinite values inside the infinite conducting/rigid inclusion of the periodic cell. It is therefore clear that these two annex problems will lead to the same numerical values for the entries of the tensors of permittivity and density if we choose the same geometry for the electromagnetic and acoustic problems. It is thus enough to compute the annex problem for either the electromagnetic or acoustic problem to find out the required set of geometric and (normalized) material parameters. We choose to do it for the electromagnetic problem.
We then note that the topology of a toroidal inclusion in a toroidal cell with periodic boundary conditions is the same as that of a straight cylindrical inclusion in a periodic cubic cell with the two ends of the inclusion touching the left-and right-hand sides of the cell. The topological equivalence allows us to simplify the problem of finding effective properties for a periodic set of infinite conducting/rigid toroidal inclusions to that of cylindrical inclusions in a perforated cubic cell, which is easily tractable with the homogenization theory framework discussed above.
Using the COMSOL finite-element package, we computed the anisotropic matrix of effective permittivity (resp. density) from the annex problems and compared it to what the Maxwell-Garnett formula provides as a (scalar) estimate. These results are reported in table 1 for various radii of the toroidal inclusions corresponding to given filling fractions f in the Maxwell-Garnett formula. This formula obviously breaks down as the inclusions are far from being spherical and, moreover, the contrast is large. Importantly, the effective properties are isotropic in the x z-plane shown in figures 3 and 4 (they correspond to ε 11 = ε 22 , resp. ρ 11 = ρ 22 ). The anisotropy is along the axis of the toroids (that is, in the xy-plane in figures 3 and 4) and therefore the lensing effect in the x z-plane in figures 3 and 4 cannot be attributed to some guidance phenomenon. The large anisotropy might be held responsible for the poor concentration effect of the infinite conducting toroidal channels compared to the analogous design with air channels we proposed in [51]. Interestingly, the concentration effect is weaker in the case of pressure waves, compared to electromagnetic waves, as can be quantitatively observed in panels (a) of figure 3 and 4. Moreover, the focusing effect appears further away from the lens in the acoustic case, which bears some resemblance to jet effects, discussed later on in this paper.

Spherical GRIN lenses in optics
Let us now turn our attention to spherical GRIN lenses. In that case, it goes without saying that there are still some issues in artificial anisotropy induced by the homogenization approach, the effect of which spoils to a certain extent the functionality of multilayered GRIN lenses, which are assumed to have a spatially varying, but usually isotropic, refractive index. In contrast, this artificial anisotropy is exactly what one needs in order to design spherical invisibility cloaks, concentrators or optical black holes.
Indeed, it follows from theorem 1 that in a multilayered spherical configuration, the tensors of homogenized permittivity and permeability can be written in the spherical basis (r, θ, φ) as where f Y = 1 0 f (r, θ, φ) dr , which obviously induces some artificial anisotropy.
14 Importantly, from Hölder's inequality [65], we know that for any functions h and g such that h p and g q are integrable on Y with 1/ p + 1/q = 1 Hence, taking p = q = 2, h = √ f and g = h −1 , squaring both members of the inequality, we are assured that f −1 −1 Y f Y (that is, the harmonic mean is smaller than the arithmetic mean), which shows that the smaller component of the anisotropic tensors [ε hom ] and [µ hom ] is along the radial direction. Importantly, this inequality holds even when p = 1 and q = ∞, in which case the second integral in the right-hand side of the inequality is the essential supremum norm of |g| on Y , which is inf{a ∈ R | L(x : |g(x)| > a) = 0}, where L is the Lebesgue measure on Y . This allows for a variety of inequalities between arithmetic, harmonic and geometric means, the former being always greater than the latter (which is itself greater than harmonic mean). Indeed, the rule of thumb in the homogenization theory of Maxwell's equations [62] is that the anisotropy is always stronger along the directions displaying elongated structures, which in the present case is along θ and φ. Fortunately, it is exactly along these directions that one needs anisotropy for invisibility cloaks [1,2], as well as optical black holes that mimic the Schwarzschild metric [76].
From theorem 2, we derive that similar results hold for the homogenized density and compressibility: However, we note that the harmonic and arithmetic means in [ρ hom ] are interchanged compared to (25), which is due to the fact that the curl operator rotates the homogenized matrix inside the gradient operator of (27). Such a duality result is known in the context of homogenization of checkerboards [66]. Fortunately, the same duality occurs between the transformed Maxwell and wave equations. Let us now explore the focusing properties of a few spherical GRIN lenses, keeping in mind the fact that artificial anisotropy is ubiquitous whenever some multilayered design is implemented, and can thus play an antagonistic role. From now on, we assume that the optical devices consist of non-magnetic media, i.e. the relative permeability is one.

Eaton lens
There are a number of variants of the Eaton lens worth looking at, but we shall focus in what follows on the following three cases [69]: • Case 1. Bending by 90 • .
The refractive index for this Eaton lens is the solution of the following equation: We report in figure 5 the corresponding material parameters and control of the light trajectory. The 90 • bend is indeed achieved with a multilayered Eaton lens, and we checked that this works over the finite interval of wavelengths [600, 1000] nm. We note that a Chinese group (Mei and Cui) has successfully implemented some metamaterial device allowing for a 90 • bend in a 2D setting [77].
The refractive index for this electromagnetic Eaton lens is the solution of the following equation: We report in figure 6 the corresponding material parameters and control of the light trajectories.
It should be noted that most of the light is leaking away from this device, as a U turn is an extremely sharp twist on the light trajectory. Moreover, the refractive index varies over a larger interval than for the 90 • bend, see figure 6. One way to improve this retroreflector would be to make use of Ulf and Tyc's transmutation theory, see [4] and also our introductory section.
• Case 3. Bending by 360 • (i.e. straightforward trajectory with a loop inside the lens). The refractive index for this Eaton lens is the solution of the following equation: We report in figure 7 the corresponding material parameters and control of the light trajectories, where the scattering worsens compared to 90 • and 180 • bends. Moreover, the refractive index varies over a larger interval than for the 90 • and 180 • bends, see figures 5 and 6. It is fairly hard to visualize the loop inside the multilayered Eaton lens, but one can see that some light indeed travels through this less than usual device.

Luneburg lens
Back in 1944, Luneburg [53] proposed a refractive index which creates two conjugate foci outside of a spherical lens. Each point on the surface of a Luneburg lens is the focal point for parallel radiation incident on the opposite side. Ideally, the permittivity of the material composing the lens decreases from 2 at its center to 1 at its surface (or equivalently, the refractive index n decreases from √ 2 to 1), according to where R is the radius of the sphere. The fact that the refractive index is 1 on the surface of the sphere ensures that the lens is impedance matched to vacuum and hence transparent. We report in figure 8 the corresponding material parameters and control of the light trajectories. It is clear here that a focusing effect occurs on the boundary of the spherical lens.

Optical black hole
According to [76], one way to mimic a real Schwarzschild black hole with transformation optics is to consider some anisotropic heterogeneous permittivity and permeability, which have the following form in spherical coordinates: where L is the event horizon of the approximate black hole. We note that for r = L, there is a singularity which can be avoided by replacing L by L + δ with δ being small positive parameter on the sphere r = L of the event horizon. In Cartesian coordinates these tensors take the following form: We report some finite-element simulation for such an anisotropic heterogeneous optical black hole in figure 9(b).
Such anisotropic tensors can be approximated using homogenization theory in a way similar to what was done in the previous section.
For this, we note that if the spherical lens consists of an alternation of two homogeneous isotropic dielectric layers of thicknesses d A and d B and permittivities ε A , ε B , we have [63] 1 where η = d B /d A is the ratio of thicknesses for layers A and B and d A + d B = 1.
If the layers are also magnetic and alternate permeabilities µ A and µ B , the resulting effective permeability is given by [62] 1 Interestingly, the structure of (33) and (35) and (36) is very similar, and this suggests that approximating the black hole parameters by a piecewise constant permittivity and permeability might work. Furthermore, the ray trajectories should be nearly preserved even for a constant permeability (at the cost of a slight impedance mismatch). If we throw away the magnetic parameter, we find that the staircase approximation in figure 9(c) gives indeed a reasonable estimate of the required anisotropy, and this is validated by the plot of the electromagnetic field in figure 9(d).
However, there are other works which achieve very similar light trajectories for some refractive index profile of the Coulomb potential type. More precisely, the relative permittivity for the optical black hole is given by [72,73] ε = (r 2 /r 1 ) 2 + iδ, 0 r r 1 , where r 1 , r 2 are, respectively, the radii of the internal core and the coating and δ is a small positive parameter which accounts for absorption. Such optical black holes were first introduced by the grouping of Ulf Leonhardt [74] in the context of slow light and thus were not expected to mimic the quantum effects of a black hole (no Hawking radiation). These optical black holes, nevertheless, mimic the classical properties of a gravitational black hole [74,75], making it potentially useful in studying other properties of gravitational black holes.
We report in figure 10 the results of our investigations for such a Coulomb profile when δ 1. These results are indeed similar to that of figure 9.

Spherical GRIN lenses in acoustics
The refractive index for acoustic Eaton lenses is the solution of the same equations as in the previous section in agreement with the fact that the transformed permittivity as well as homogenized permittivity and density have the same structure. The fact that pressure waves are solutions of a scalar wave equation transpires in the scalar nature of the transformed and homogenized compressibility (unlike for the permeability), but otherwise all things are very similar. In this section, we therefore adopt a more sketchy description of the material properties; see the previous section for further details.
It is actually fairly easy to deduce from the annex problem L j , see [64], that if the multilayered lens consists of an alternation of two homogeneous isotropic layers of thicknesses

Radius (µm)
Permittivity [ε] Figure 9. Multilayer approximation (11 layers) that mimics the black hole signature with r 1 = 1.2 µm, r 2 = 2.3 µm. The wavelength for the incident beam is λ = 700 nm; (a) diagrammatic view of the computational domain for an optical black hole: the sphere corresponds to its event horizon and the anisotropic space given by (34) is of infinite extent; here the box is surrounded by PMLs specially designed for the black hole medium; (b) (x-y)-plane view of the norm of the electric field for a Gaussian beam incident upon a heterogeneous anisotropic black hole given by (34); (c) staircase approximation of the profile of permittivity ε θ = ε φ given by (33); (d) (x-y)-plane view of the norm of the electric field for a Gaussian beam incident upon a multilayered black hole given by (37): the black hole effect is less visible than in (b), since only the region within the event horizon r L approximates the spatially varying anisotropic Schwarzschild metric with concentric layers of anisotropic homogeneous media (each anisotropic layer being deduced from a homogenization of thinner isotropic layers), while for r > L, the permittivity is that of vacuum. d A and d B and compressibility moduli κ A , κ B and densities ρ A and ρ B , we have where η = d B /d A is the ratio of thicknesses for layers A and B and d A + d B = 1. This formula has been rediscovered in many instances, but has already appeared in the mathematics literature in 1978 [65].

Eaton lens
We once again consider three cases. • Case 1. 90 • bend Parameters for the Eaton lens are given by where R is a parameter to match impedance (ideal case R = r 2 ) and r 1 , r 2 are, respectively, the radii of the internal core and the coating. We report in figure 11 the corresponding material parameters and control of the sound trajectories. The beam is indeed bent at an angle of 90 • . We note that the functionality of the Eaton lens is enhanced compared to its optical counterpart, as there is less unwanted scattering.
• Case 2. 180 • bend Parameters for the Eaton lens are given by   where R is a parameter to match impedance (ideal case R = r 2 ) and r 1 , r 2 are, respectively, the radii of the internal core and the coating. We report in figure 12 the corresponding material parameters and control of the sound trajectories. A bending of 180 • is achieved and is markedly enhanced compared to the optical counterpart.
• Case 3. 360 • bend Parameters for the Eaton lens are given by where R is a parameter to match impedance (ideal case R = r 2 ) and r 1 , r 2 are, respectively, the radii of the internal core and the coating. We report in figure 13 the corresponding material parameters and control of the sound trajectories for a multilayered Eaton lens. The Gaussian beam is bent to a certain extent.

Luneburg lens
Parameters for the acoustic Luneburg lens are given by   where R is a parameter to match impedance and control the size of the focus point (ideal case R = r s ) and r s is the radius of the spherical lens. We report in figure 14 the corresponding material parameters and control of the sound trajectories when R = r s . It is clear here that a focusing effect occurs on the boundary of the spherical lens.
Let us now consider the case when R > r s . In this case, there is an impedance mismatch of the refractive index on the sphere boundary, as exemplified by the jump in the profile of the density in figure 15(c). One should, however, note that the focusing effect now occurs further away from the boundary of the spherical lens: the lens now behaves as an acoustic jet, in a similar way to what Bonod et al achieved for photonics [79].

Black hole
Let us, finally, consider some refractive index profile of the Coulomb potential type for an acoustic black hole. More precisely, the density is given by ρ = (r 2 /r 1 ) 2 + iδ, 0 r r 1 , ρ = (r/r 2 ) 2 , r 1 r r 2 , ρ = 1, where r 1 , r 2 are, respectively, the radii of the internal core and the coating and δ is a small positive parameter which accounts for viscosity. The more the viscosity, the more the absorption inside the acoustic black hole. We report in figure 16 the corresponding material parameters and control of the sound trajectories. Sound is trapped to a certain extent within the acoustic black hole for δ 1.  (d) (x-y)-plane view of the norm of the pressure field.

Conclusion
In this paper, we have reviewed the bending and focusing properties of a class of cylindrical and spherical Eaton and Luneberg lenses in conjunction with an effective medium approach. All these metamaterials have been studied for electromagnetic and acoustic waves, with special emphasis on the similarities of the corresponding governing equations. We noted that while both transformed permittivity and permeability are rank-2 tensors, only the transformed density is a rank-2 tensor (the transformed compressibility is a scalar function). In the same way, both homogenized permittivity and permeability are anisotropic, whereas only the homogenized density is anisotropic. Such correspondences between the structure of the transformed and homogenized equations allowed us to design cylindrical Luneberg lenses with toroidal channels (infinite conducting for electromagnetic waves and rigid for acoustic waves) as well as spherical multilayered Luneberg and Eaton lenses for light and sound focusing and bending. We finally looked at optical and acoustic black holes which trap light and sound to a certain extent. We hope our study will foster experiments in these directions.