Electromagnetic cloaking in convex and concave media with surface modelled as a parameterised function

The onset of transformation optics has opened avenues for designing of a plenitude of applications related to propagation of electromagnetic waves in anisotropic media. In this paper, an algorithm is proposed using a coordinate transformation and a piecewise function for the purpose of designing a three dimensional cloak having an arbitrary geometry which could be convex or non-convex in nature. The surfaces of the cloak as well as of the body under consideration are assumed to be conformal to each other. For an arbitrary geometry, the coordinate system needed to model the surface can be a non-orthogonal system. For the purpose of verification of the algorithm, a ray tracing process is carried out for an ellipsoid as well as for a concave surface having axial symmetry. In order to solve the Hamiltonian equation for the purpose of ray tracing, the process of finding the derivatives analytically, for an arbitrary geometry as considered here, becomes very cumbersome. Here, a numerical method is described which provides a better approximation to the partial derivatives than the conventional finite difference approach based on forward differences.


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,2,3,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].

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). 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 2 , u 3 ) 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 2 , u 3 ). 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 2 , u 3 ): The unit vectorR can be expressed in terms of its components aŝ 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, where I is the Identity tensor. Relative permittivity and permeability in the transformed space can be expressed as: The expression for Jacobian in the primed coordinates is: In Eq. (4), the position vector is expressed in terms of variables (r, u 2 , u 3 ) which are considered over here as the generalised curvilinear coordinates denoted by (q 1 , q 2 , q 3 ) 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: whereê 1 ,ê 2 ,ê 3 are the standard Cartesian unit vectors. This can be expressed as: Since q 1 = 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 1 , q 2 , q 3 ) is transformed to the primed coordinate system (q ′1 , q ′2 , q ′3 ) such that (q 1 = r) → (q ′1 = r ′ ) whereas both q 2 & q 3 remain unchanged in the primed coordinate system.

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 2 . 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:

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 1 , q 2 , q 3 ) 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]: where t is the parameterisation variable and x is the position vector. Ray tracing is carried out by solving (x, k) as a function of t starting from t = 0, which corresponds to the point at which the incident wave touches the cloak.
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 1 , q 2 , q 3 ) coordinate system. That in turn implies that the ray tracing algorithm also needs to be written in terms of the generalised coordinates.
where k is still in the Cartesian system. 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 ′(1) , q ′(2) , q ′(3) ), 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) i = 1, 2, 3 is an infinitesimal variation in q ′(i) . For the purpose of illustration, we assume that i = 3 i.e. the numerical derivative with respect to q ′(3) is being calculated in the following steps.
When ∆q ′(3) 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 ′(3) ) 2 . Instead, if it is assumed that ∆q ′(3) is imaginary, and is expressed as j∆q ′(3) , the Taylor series for H(q ′(1) , q ′(2) , q ′(3) ) can now be expressed in terms of j∆q ′ (3) .

H.
Thus, The order of the first neglected term (∆q ′(3) ) 3 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 can be obtained by putting i = 1, 2 respectively and writing a Taylor series expression The differentiation of the generalised coordinates (q ′(1) , q ′(2) , q ′ (3) ) can be carried out as: where 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 1 , k 2 , k 3 ) and generalised coordinates (q ′1 , q ′2 , q ′3 ). The partial derivatives of Hamiltonian H at a given value of (q ′(1) , q ′(2) , q ′(3) ) 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 ′(1) , q ′(2) , q ′ (3) ) 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]: where k i is the incident wave and k is the wave vector inside the cloak at the point of incidence. By solving the coupled differential equations, we get expression of (q ′(1) , q ′(2) , q ′(3) ) and k as a function of t for the ray propagating in the cloak. For plotting the path of the ray, (q ′(1) , q ′(2) , q ′(3) ) are used to calculate the position vector x in the Cartesian coordinate system.

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. Since the body under consideration is axisymmetric, only one parameter is required to define the surface, and we denote it as (q ′2 = 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 the terms in Eq. (18) can now be expressed as: In order to verify the proposed algorithm, a MATLAB based simulation is carried out using R(q ′(2) ) 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.
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.
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 ′(2) ) = 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  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. 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.
The Hamiltonian H is function of the Jacobian Λ il in Eq. (15)and Eq. (19) which depend on the parametric equation of the outer surface R(q 2 , q 3 ) (generalised case) or R(q 2 ) (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.   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.