1. Introduction

The development in the fields of metamaterials and transformation optics in the past decade has enabled scientists to envision and realise devices/applications such as electromagnetic cloaks [1–4], perfect lenses [5] and spacetime cloaks [6]. A formalism of transformation optics which also encompasses gravitational curvature has been illustrated in [7]. A technique which uses transformation between spacetime manifolds with application to cloaking is shown in [8]. Cloaks with non-conformal inner and outer boundaries have been analysed at by Li and Li [9]. The concept of cloaking or invisibility which is one such application, can be explained as bending the path of a ray around the body of interest so as to give an impression that the ray is travelling in a straight line uninterrupted by any obstacle [10]. In 2006, Pendry, et al. [4] and Leonhardt [11] carried out a coordinate transformation in order to alter the material characteristics, which would in turn cause a change in the formulation of Maxwell’s equations. The idea of ray tracing in the transformed media was used by Pendry, et al. [12] to verify the concept for spherical and cylindrical cloaks using Cartesian tensors. A generalised method for designing arbitrarily shaped cloaks using the approach of coordinate transformation has been discussed by Li and Li in [13]. In [14], an algorithm is proposed for tracing the path of a ray in an arbitrarily shaped cloak in three dimensions using spherical polar coordinates (θ, ϕ), which requires that the contour of the outer boundary of the cloak should be a single valued function of (θ, ϕ). In this paper, an algorithm is propounded for cloaking of an arbitrarily shaped body in three dimensions using a parametric formulation without the above constraint of the need of a single valued function. The requirement of a numerical method to calculate the derivatives of the Hamiltonian stems from the limitation that when the geometry under consideration is complex, analytical solutions for transformation equations may not exist. The body under consideration could either have a concave or convex geometry. The surfaces of the cloak as well as of the body under consideration are assumed to be conformal to each other. This representation can either be a piecewise continuous polynomial similar to Finite Element Method [15] or be defined in terms of Non-Uniform Rational B-Spline (NURBS) [16]. For the purpose of illustrating the technique, an axisymmetric body is represented using a piecewise linear function. That enables the algorithm to cloak a concave shaped body. Higher order functions can also be used since the procedure remains unchanged. In numerical differentiation, truncation and round-off errors cause numerical instability while using the conventional finite difference approximation based on forward differences [17]. For the purpose of ray tracing, a Hamiltonian is calculated using a technique which has been adapted from that reported in [12].

2. Mathematical formulation

In order to cloak an arbitrarily shaped object, it is coated with a thin layer of a material such that any electromagnetic wave incident on the body will move around the object with zero reflection. For practical reason, the coating thickness is uniform, hence it can be considered that the cloak would have the same topology as the object to be cloaked. Then, the formulation can be stated as follows:

Consider O to be the centre of the concentric bodies as shown in (Fig. 1).

 

Fig. 1 Coordinate transformation for an arbitrary geometry.

Download Full Size | PDF

Consider three points A, B and Q on Γ. A surface Γ′ is obtained by scaling Γ by a suitable factor.

In the figure, the scaling implies:

Here τ < 1 is the measure of the thickness of the cloak. Any surface can be represented by an implicit equation f (x, y, z) = 0. A parametric representation, which need not be unique, of the same surface can be written in terms of two parameters:

The outer surface of the cloak is described using parametric coordinates as:

The magnitude of the position vector r will be R(u₂, u₃) for any point on the surface. In [14], a parametric representation has been used in terms of spherical polar coordinates (θ, ϕ). Such an approach has a limitation in modelling general surfaces that can be concave because (θ, ϕ) representation requires R to be a single valued function of (θ, ϕ). The outer surface of the cloak in (Fig. 1) cannot be modelled by the algorithm mentioned in [14] since points C, D and E lie along the same radial path, hence for the same values of (θ, ϕ), R(θ, ϕ) will have different values.

Other parametric representations like NURBS or planar polynomial parametric approximation do not have such limitations. Here we describe a formulation in terms of an arbitrary parametric representation (u₂, u₃). The following coordinate transformation is applied.

Here, r is the coordinate of any point P inside the volume enclosed by Γ′. Owing to this transformation, any point P within the outer surface Γ′ is transformed to P in the annular space between Γ and Γ′. The unprimed coordinates apply to the whole domain while the primed coordinates apply to the annular region. A point transforms from P(x, y, z) in the original system to P (x′, y′, z′) in the primed system.

Due to transformation invariance, the unit vectors in the original system and the transformed system must be equal.

Expressing the equation in the tensor notation and using Eq. (2),

Here the subscripts vary from 1 to 3 and we shall use Einstein summation convention. Any point inside the cloak can be described in terms of three variables (r, u₂, u₃):

The unit vector $\widehat{\mathit{R}}$ can be expressed in terms of its components as

In terms of the primed variables, the point can denoted by:

The Jacobian for the coordinate transformation can be written as:

In the tensor notation, it is expressed as:

The constitutive properties in the transformed space are:

If the original space is assumed to be free space,

The expression for Jacobian in the primed coordinates is:

In Eq. (4), the position vector is expressed in terms of variables (r, u₂, u₃) which are considered over here as the generalised curvilinear coordinates denoted by (q¹, q², q³) where

The basis vectors a_i for (i = 1, 2, 3) are defined as follows:

The differential position vector can be expressed in terms of Cartesian components:

Since q¹ = r, we can write using Eq. (4):

In order to define the gradient in terms of generalised coordinates, a differential volume element dV is required.

The gradient of a scalar function F can be expressed as:

The Jacobian, Eq. (11) can now be expressed as

Using [8] where transformation optics is expressed as a geometric formulation on spacetime manifolds, an expression similar to Eq. (15) has been generated.

The original coordinate system (q¹, q², q³) is transformed to the primed coordinate system (q′¹, q′², q′³) such that (q¹ = r) → (q′¹ = r′) whereas both q² & q³ remain unchanged in the primed coordinate system.

3. Axisymmetric cloak

The special case of an axisymmetric object is considered.

The object is modelled in the XZ plane as a curve with a single parameter q². A three dimensional object is generated by rotating the curve about the Z axis by 2π.

The basis vectors can be written using Eq. (13):

The Jacobian now becomes:

3.1. Hamiltonian

Since there is no loss of energy while the wave propagates through the medium, the Hamiltonian can be found in order to trace the path of the wave in the cloaked medium [12].

Here k is the propagation vector in the cloak.

The first and second terms of Eq. (20) become:

The Hamiltonian equation is:

When Eq. (19) is substituted in Eq. (21), it is observed that the Hamiltonian H in Eq. (21) is a function of generalised coordinates (q¹, q², q³) and propagation vector k.

For the purpose of ray tracing, the path can be parameterised using the Hamiltonian in Eq. (21) as explained in [19]:

In the technique proposed by Pendry, et al. [12], the independent variables for the Hamiltonian are in terms of spherical polar coordinates (r, θ, ϕ), whereas the ray tracing is done in terms of Cartesian coordinates. At every step in the algorithm, an inverse transformation from Cartesian to the spherical polar coordinate system is applied, which is trivial in that case. But in the coordinate system used in this paper, there is no explicit relation for inverse transformation from (X, Y, Z) to (q¹, q², q³) coordinate system. That in turn implies that the ray tracing algorithm also needs to be written in terms of the generalised coordinates.

In order to solve the coupled differential equations, Eqs. (24) and Eq. (25) mentioned later, the derivative of the Hamiltonian is required to be evaluated with respect to the generalised coordinates (q′⁽¹⁾, q′⁽²⁾, q′⁽³⁾), see Eq. (15). The analytical expression for that would be too tedious to evaluate. We describe a numerical method which has better accuracy than the conventional finite difference approximation based on forward differences. In the method we follow, we assume that the differential term Δq′⁽ⁱ⁾ i = 1, 2, 3 is an infinitesimal variation in q′⁽ⁱ⁾. For the purpose of illustration, we assume that i = 3 i.e. the numerical derivative with respect to q′⁽³⁾ is being calculated in the following steps.

When Δq′⁽³⁾ is real, the Taylor Series of the Hamiltonian H (for the first four terms neglecting higher order terms) can be expressed as

If only the first order term is retained, the standard forward difference approximation for the derivative is obtained. The first neglected term is the 2^nd order term (Δq′⁽³⁾)². Instead, if it is assumed that Δq′⁽³⁾ is imaginary, and is expressed as jΔq′⁽³⁾, the Taylor series for H(q′⁽¹⁾,q′⁽²⁾,q′⁽³⁾) can now be expressed in terms of jΔq′⁽³⁾.

Thus,

The order of the first neglected term (Δq′⁽³⁾)³ is three instead of two in the earlier case. This technique described by Eq. (23) is used for evaluating the partial derivative numerically instead of the conventional finite difference method. This method does not suffer from round-off errors if the derivative is small. Further, as mentioned above, the lowest nonzero derivative is of 3^rd order in the residue compared to 2^nd order in the conventional finite difference method. Similar expressions for $\frac{\partial H}{\partial {q^{\prime }}^{(1)}}$ & $\frac{\partial H}{\partial {q^{\prime }}^{(2)}}$ can be obtained by putting i = 1, 2 respectively and writing a Taylor series expression The differentiation of the generalised coordinates (q′⁽¹⁾,q′⁽²⁾,q′⁽³⁾) can be carried out as:

Let ${\left[\tilde{\mathbf{\Lambda }}{\tilde{\mathbf{\Lambda }}}^T\right]}_{jl}{k}_l=V_j$ v is a vector whose j^th component in Cartesian coordinates is V_j, and k_l are the Cartesian components of vector k.

The differential of the Hamiltonian can be written as

To determine the path of the ray moving through the medium, the coupled differential equations, Eqs. (24) and Eq. (25), have to be solved. In each of these equations, the right hand side is a function of Cartesian components of k(k₁, k₂, k₃) and generalised coordinates (q′¹, q′², q′³). The partial derivatives of Hamiltonian H at a given value of (q′⁽¹⁾, q′⁽²⁾, q′⁽³⁾) are evaluated using Eq. (23) and substituted in Eq. (25). A program has been written in MATLAB and the differential equation toolbox has been used to solve the equations. The initial value of (q′⁽¹⁾,q′⁽²⁾,q′⁽³⁾) is taken as the coordinate of the point where the incident wave touches the outer cloak boundary. The initial values for components of k are obtained by solving the Hamiltonian Eq. (21) at the point of incidence using the following boundary condition on the surface [12]:

3.2. Piecewise linear model

An example of a piecewise linear model is shown in (Fig. 2). It is modelled as a union of piecewise linear functions between adjacent nodes (similar to that used in FEM). In the figure, the line joining nodes 1 and N (here N = 7) is the Z-axis. By rotating the figure about the Z-axis, a three dimensional object is generated.

 

Fig. 2 Piecewise linear model.

Download Full Size | PDF

Since the body under consideration is axisymmetric, only one parameter is required to define the surface, and we denote it as (q′² = u). Between k^th and (k + 1)^th nodes, any point in the unprimed system is modelled in terms of a parameter u, which takes values as (k − 1) ≤ u ≤ k. By assuming the piecewise function to be

In order to verify the proposed algorithm, a MATLAB based simulation is carried out using R(q′⁽²⁾) as a piecewise function defined by Eq. (26) to interpolate between successive points. In the first case, the nodal points are chosen to be on the surface of an ellipsoid. The coordinates of the position vector r are plotted inside the cloak for two of three different sets of nodes as shown in (Fig. 3) and (Fig. 4). In (Fig. 3), the solid coloured lines depict the paths of the rays to traverse half the cloak: the path in green is of the ray when only 5 nodes are used, the path in black when 25 nodes are used, while red denotes the path when 65 nodes are used. The dashed lines in the figure, which are in red, indicate the inner and outer boundaries of the cloak for 65 nodes while the dashed lines in green indicate the boundaries for 5 nodes. It can be observed from the figure, that not only there is a prominent difference in the paths traced by the rays for 5 nodes (solid green curve) and 65 nodes (solid red curve), the geometry of the cloak has also been affected (observe the change in geometry from a polygon for 5 nodes to an ellipsoid for 65 nodes) causing a change in the angle of incidence. In (Fig. 4), the solid coloured lines depict the paths of the rays to traverse half the cloak: the path in blue is of the ray when only 10 nodes are used, the path in black when 25 nodes are used, while red denotes the path when 65 nodes are used. The dashed lines in the figure, which are in red indicate the inner and outer boundaries of the cloak for 65 nodes while the dashed lines in blue indicate the boundaries for 10 nodes. It can be observed in this case also, that there is a difference in the paths traced by the rays for 10 nodes (solid blue curve) and 65 nodes (solid red curve), and the geometry of the cloak has also been affected, but that change is less conspicuous than in the previous case for 5 nodes.

 

Fig. 3 Two-dimensional plot of an ellipsoidal cloak in the X–Z plane: the solid lines are the paths of the rays for different number of nodes; green colour for 5 nodes, black colour for 25 nodes and red for 65 nodes. The dashed lines in red indicate the inner and outer boundaries of the cloak for 65 nodes while the dashed lines in green indicate the boundaries for 5 nodes.

Download Full Size | PDF

 

Fig. 4 Two-dimensional plot of an ellipsoidal cloak in the X–Z plane: the solid lines are the paths of the rays for different number of nodes; blue colour for 10 nodes, black colour for 25 nodes and red for 65 nodes. The dashed lines in red indicate the inner and outer boundaries of the cloak for 65 nodes while the dashed lines in blue indicate the boundaries for 10 nodes.

Download Full Size | PDF

It can be concluded that as the number of nodes is changing, the geometry gets modified because of the piecewise linear function involved, causing a change in the point of incidence. The modelling with different sets of points is done to bring out the dependence of the path of a ray on the number of points used in the model. To elucidate the importance of nodes required for plotting a curve, an ellipsoid as shown (Fig. 5) is plotted for two sets of nodes: 5 and 20. For the smaller number of nodes, the curve generated is a polygon and the point of incidence gets changed because of the distorted geometry. It can be observed from (Fig. 3) and (Fig. 4), that the smoothness and accuracy of the path traced improves as the number of nodes is increased.

 

Fig. 5 Plotting of an ellipsoid for 2 sets of nodes: blue coloured curve is for 5 nodes while maroon coloured curve is for 20 nodes.

Download Full Size | PDF

In order to verify the validity of the algorithm for concave surfaces, the nodal points are chosen to be on a concave surface and the MATLAB based simulation is carried out using R(q′⁽²⁾) = R(u) as the piecewise function defined by Eq. (26). The path of the ray is plotted inside the concave surfaced cloak with moderate curvature for three different sets of nodes as shown in (Fig. 6). It can be observed that due to the usage of only 10 nodes for the blue curve, the algorithm is not able to capture the cusp in the geometry. As the number of nodes is increased, the variations in geometry are faithfully reproduced and the path of the ray is in accordance with the actual geometry of the cloak. In order to test the robustness of the algorithm, the sharpness of the concave curve is increased further and ray tracing is carried out as seen in (Fig. 7). The cyan coloured curve depicts the path of the position vector r turning just next to the sharp edge of the concave curve, and it is also within the outer boundary of the cloak as desired. As the sharpness of the curve increases further, piecewise linear model might not be sufficient to handle those sudden variations. (Fig. 7) shows the limiting case of the usage of linear function for parameterising, where the ray to be traced just manages to turn sufficiently to avoid entering the body to be cloaked. A quadratic or even higher order function would be required to model the surface with greater sharpness. Further, if NURBS is used, it would be possible to model any surface with complicated geometry. In (Fig. 7), in order to compare the two numerical differentiation techniques, one using imaginary method Eq. (23) and the other, conventional finite difference method,the paths of the ray are tested on linear fitting and a parameter du is varied such that there is a difference between the paths using both the numerical differentiation techniques. The results are plotted in region between (−0.5) and (0.5) of the Z axis to highlight the region around the cusp. It can be observed in (Fig. 8), that the cyan coloured ray which is plotted using conventional finite difference method hits the outer boundary of the cloak while the black dotted ray plotted using imaginary method differential technique manages to circumvent around the outer edge of the cloak and continue its path further.

 

Fig. 6 Two-dimensional plot of a concave surface cloak with moderate curvature in the X–Z plane: the blue coloured curve shows the path of the ray for 10 nodes, black coloured curve shows the path for 30 nodes, and cyan colour shows the path for 100 nodes.

Download Full Size | PDF

 

Fig. 7 Two-dimensional plot of a concave surface cloak in the X–Z plane where the cyan coloured wave is the path of the ray in a cloak with a sharp inward curve.

Download Full Size | PDF

 

Fig. 8 Comparison of paths of the ray plotted in a concave surface cloak in the X–Z plane using the imaginary method (Eq. (23)) and the conventional finite difference method where the cyan coloured wave is the path of the ray using the conventional finite difference method. while the black coloured wave is plotted using the imaginary method.

Download Full Size | PDF

The Hamiltonian H is function of the Jacobian ${\tilde{\mathbf{\Lambda }}}_{il}$ in Eq. (15) and Eq. (19) which depend on the parametric equation of the outer surface R(q², q³) (generalised case) or R(q²) (axisymmetric case) and the derivatives of R with respect to the generalised coordinates. Now R in turn is a function of nodes, Eq. (26), which implies that the Hamiltonian is a function of the nodes selected.

3.3. General parametric representation

The algorithm described in Sec. (3.2) uses a piecewise linear representation of an axisymmetric surface. However, the formulation is not limited to its parametric equation. As seen from Eq. (15), the Hamiltonian can be calculated if the parametric representation R(q′⁽²⁾, q′⁽³⁾) is given. A general parametric representation can be described by NURBS [16]. However, the derivatives $\frac{\partial \mathit{R}}{\partial {q^{\prime }}^2}$ and $\frac{\partial \mathit{R}}{\partial {q^{\prime }}^3}$ cannot be found analytically. Special routines need to be written to evaluate the derivatives, though they are not difficult to implement. The other equations required for ray tracing, Eq. (23), Eqs. (24) and Eq. (25), can all be generalised for that case. Thus any arbitrary parameterisation of a surface can be implemented.

In case of a cloak whose outer surface is not conformal with the body to be cloaked, Eq. (2) of [9] can be adapted to use the technique described here, but that is beyond the scope of the present work.

3.4. Numerical stability

The numerical method used for solving the Hamiltonian equation Eq. (21) which is an ordinary differential equation linearises it. In Eq. (11), when the input wave is close to the inner surface of the cloak, the denominator of the equation tends to zero since (r′ → τR). The relative permittivity and permeability Eq. (10) become very large. For this case, the eigenvalues of the linearized system differ by several orders of magnitude, or they also change during integration [18] which leads to instability in the solution. In the method used in [12], the independent variables are the Cartesian coordinates used to trace the path of the ray. Usage of this method leads to numerical instability for the case when the input wave is close to the inner surface of the cloak. In the method used in this paper, the independent variables are the generalised coordinates (q′⁽¹⁾, q′⁽²⁾, q′⁽³⁾) where q′⁽¹⁾ = r′. The method described here circumvents the problem because the generalised coordinate r′ continuously decreases and q′⁽²⁾ lies on a smaller curve compared to the outer boundary and which is conformal with it. The points r_i and r_i+1 are closer to each other as compared to the outer region i.e. the distance between the nodes reduces as the wave travels towarsd the inner surface. Hence the accuracy of the calculation improves.

4. Conclusion

An algorithm based on piecewise linear function has been proposed for designing an arbitrarily shaped cloak in three dimensions, which circumvents the problem of numerical instability when the input wave is close to the inner surface. A procedure for calculating the path of a ray through it is described using ray tracing. In order to demonstrate the validity of the algorithm, it has been applied to an axisymmetric ellipsoid cloak as well as an axisymmetric concave surfaced cloak. For the purpose of differentiation of the Hamiltonian, a numerical technique is used instead of a standard analytical approach which becomes too cumbersome in this case. The versatility of this algorithm stems from the usage of a piecewise function for plotting from one point to another in the code, enabling it to effectively cloak a body having non-convex geometry or having sharp edges. Also, since the algorithm does not depend on the linear nature of the function being used, higher order functions such as quadratic ones can be used to improve the efficiency further.