1. Introduction

It was recently found that transformation optics open new avenues in electromagnetic cloaking, either through their heterogeneous anisotropic effective material parameters (transformation optics, [1, 2]) or through low index materials [3] or negative refractive index materials [4]. Interestingly, the invisibility is preserved in the case of an intense near field [5], when the ray optics picture breaks down. The mathematics behind the scene have been known from researchers working in the area of inverse conductivity problems [6]. The first experimental realization of an invisibility cloak, chiefly achieved in the microwave regime [7], 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. The latter drawback can be overcome by considering a generalized transform [8, 9] 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. Alternatively, one can design an approximate structured cloak via homogenization [10, 11, 12], although the rapid growth of the field, fueled by a keen interest of the optics community, promises a large panel of new technological applications.

A non trivial question to ask is whether one can design cloaks of non-spherical shapes. A parameterization of the cloak’s boundaries was proposed by three of us to design cylindrical cloaks of an arbitrary cross-section [17]. Its extension to the general three-dimensional case now requires to parameterize the inner and outer boundaries of the cloak as some surfaces with varying radii

which is somewhat out of reach running the COMSOL multiphysics package on a 256 Gb RAM computer: Compared to a spherical cloak, the tetrahedral mesh needs be further refined due to the complexity of surfaces involved. However, to demonstrate the versatility of our Fourier-based approach, it is enough to analyse arbitrary cloaks built with surfaces of revolution i.e. generated by rotating a curve about an axis.

In the present paper, we discuss the design of axially invariant three-dimensional cloaks with an arbitrary cross-section described by two functions R ₁(ϕ) and R ₂(ϕ) giving an angle dependent distance from the origin. These functions correspond respectively to the interior and exterior boundary of the cloak. We shall only assume that these two boundaries can be represented by a differentiable function. Their finite Fourier expansions are thought in the form R_j(ϕ)=a ^j _0,0 +∑^p _n=1 a ^j _0,n cos(nϕ), j=1,2, where p can be a small integer, and a ^j _0,n=0 for n>p for computational easiness. Note that our approach encompasses the case of toroidal cloaks [13] but cloaks with irregular boundaries such as cubes, fall beyond the scope of this study. Nevertheless, contrarily to previous three dimensional numerical studies [13, 12], our scheme can be used in the design of cloaks with non-convex boundaries.

To illustrate our methodology, we compute the electromagnetic field diffracted by a cloak with rotational symmetry about the z-axis. We perform full-wave finite element simulations in the commercial package COMSOL, when the cloak is illuminated by an approximate plane wave (generated by a constant electric field on the upper surface of the computational domain).

2. Change of coordinates and pullbacks from optical to physical space

2.1. Material properties of the heterogeneous anisotropic cloak

The geometric transformation which maps the field within the full domain ρ≤R ₂(ϕ) onto the annular domain R ₁(ϕ)≤ρ′≤R ²(ϕ) can be expressed as:

where 0≤ρ≤R ₂(ϕ). Note that the transformation maps the field for ρ>R ₂(ϕ) onto itself through the identity transformation.

This change of co-ordinates is characterized by the transformation of the differentials through the Jacobian:

This change of coordinates amounts to replacing a homogeneous isotropic medium with scalar permittivity and permeability ε and µ, by a metamaterial described by anisotropic heterogeneous matrices of permittivity and permeability given by [5, 14]

where T=J ^T J/det(J) is a representation of the metric tensor in the so called stretched radial coordinates. Note that there is no change in the impedance of the media since the permittivity and permeability undergo the same transformation.

After some elementary algebra, we find that

where

and

for R ₁(ϕ′)≤ρ′≤R ₂(ϕ′). Elsewhere, T ^-1 reduces to the identity matrix (c ₁₁=1, c ₁₃=0 and ρ=ρ′ for ρ′>R ₂(ϕ′)).

2.2. Singularity analysis of the transformation matrix

To exemplify the symmetric nature of the coefficients of T ^-1, and due to the role played by its first entry in the following singularity analysis, we show the variation of (T ^-1)11 in Fig. 1. We note that this coefficient varies between zero and three. The fine mesh on the inner boundary of the cloak is also apparent from Fig. 1, as are the non-convex and symmetric features of the cloak.

In the case of a spherical cloak, c ₁₃ vanishes and T ^-1 reduces to $\mathrm{Diag}(\frac{{\rho }^2}{c_{11}{\rho \prime }^2},c_{11},c_{11})$. We note that the first eigenvalue vanishes on the inner boundary and the other two remain constant, which is consistent with the singularity analysis led in [15]. This is unlike the circular cylindrical case whereby one eigenvalue goes to zero while the other one goes to infinity on the inner boundary [5].

 

Fig. 1. 3D plot of (T ^-1)₁₁, as given by Eq. 5, Eq. 6 and Eq. 7 within the cloak with boundaries given by Eq. 10. The symmetries are noted.

Download Full Size | PDF

However, in our case, the cloak is of an arbitrary shape, and it is therefore illuminating to look at the behaviour of its permittivity and permeability tensors’s eigenvalues. The eigenvalues of (5) are found to be

We can therefore see that λ _j,j=1,2 are spatially varying functions of ρ′(ρ,ϕ), such that λ ₁=0 when ρ=0 i.e. at the inner boundary, while λ ₂>0. Importantly, λ ₃ is a strictly positive constant for a given angle ϕ′(i.e. a function independent upon ρ). Last, we checked that when there is no longer a symmetry of revolution about one axis, all three eigenvalues are also spatially varying with θ′, but λ ₃ remains independent upon ρ. This reflects the fact that the geodesics for light within the arbitrary cloak follow more and more complex trajectories when we perturb the geometry away from the spherical cloak.

3. Finite elements computations

We would now like to further investigate the electromagnetic response of the cloak to an incident plane wave from above, see Fig. 2. For this, we choose the electric field E as the unknown:

where $k=\omega \sqrt{{\mu }_0{\epsilon }_0}=\omega ⁄c$ is the wavenumber, c being the speed of light in vacuum, and ε̳′ and µ̳′ are defined by Eqs. (4). Also, E=E _i+E _d, where E _i is the incident field and E _d is the diffracted field which satisfies the usual outgoing wave conditions (to ensure existence and uniqueness of the solution). The weak formulation associated with Eq. 9 is discretised using second order finite edge elements (or Whitney forms) which behave nicely under geometric transforms (pull-back properties) [14].

 

Fig. 2. 3D plot of the magnitude $\underset{=}{\epsilon \prime }$ of the total electric field for a plane wave of wavenumber k=2π/0.3 incident from above on an non-convex invisibility cloak.

Download Full Size | PDF

For the sake of illustration, cf. Fig. 2, let us consider a cloak with inner and outer boundaries expressed as

 

Fig. 3. 2D plot of $\sqrt{E_1^2+E_2^2+E_3^2}$ generated by a slice of Fig. 2 in the xz-plane for y=0.

Download Full Size | PDF

We report the computations for the magnitude of the total electric field in Fig. 2 (3D plot), Fig. 3 (2D plot in the xz-plane) and Fig. 4 (2D plot in the yz-plane) for a plane wave incident from above at wavenumber k=2π/0.3 (units are in inverse of a length, say µm-1 for nearly visible light (UV)).

Around 4.10⁵ elements were used in this computation, which corresponds to about 2.7 10⁶ degrees of freedom. While the convergence of the numerical scheme has been checked by considering different types of meshes, the large size of the system means we were not able to further refine the mesh for the computational resources at hand.

 

Fig. 4. 2D plot of $\sqrt{E_1^2+E_2^2+E_3^2}$ generated by a slice of Fig. 2 in the yz-plane for x=0.

Download Full Size | PDF

4. Conclusion

In conclusion, we have proposed a design of an arbitrarily shaped cloak using a Fourier approach. We only assumed that the cloak displays a symmetry of revolution about one axis in order to reduce the computational complexity of the problem (extension to completely arbitrarily shaped cloaks is a straightforward matter). Cloaking has been confirmed numerically for an incident plane wave in resonance with the concealed region and an analysis of the cloak’s singularity has been carried out.

A. Diatta and S. Guenneau acknowledge funding from EPSRC grant EP/F027125/1.