Controlling electromagnetic waves has broad range‎ ‎of applications in various area of science and engineering‎, ‎for‎ ‎instance the design of perfect flat lenses‎, ‎waveguides‎, ‎invisibility cloaks or in general metamaterials‎. ‎The idea of‎ ‎controlling electromagnetic waves by introducing materials‎ ‎inhomogeneity was introduced in [1] by‎ ‎Pendry et‎. ‎al‎. ‎It stems on this fact that the Maxwell equations‎ ‎are invariant under the coordinate transformation‎. ‎Therefore‎, ‎by ‎finding an analytic mapping from the physical (undeformed) space‎ ‎to the clocked (deformed) space in which the transformation maps a‎ ‎point onto a ball with the finite reduce‎‎, ‎the mentioned ball‎ ‎will be hidden from the electromagnetic waves‎. ‎Later‎, ‎this method‎ ‎has been extensively used by various researchers for the design of‎ ‎metamaterials [2]-[4], ‎photonic band gap‎ ‎devices [5]‎, ‎active nanoplasmonic‎ ‎metamaterials [6]‎, ‎extreme-angle broadband‎ ‎metamaterial lens [7]‎, ‎all-dielectric‎, ‎polarization-independent optical cloaks‎ [8]‎.

‎    There are several efforts in the literature for the design of‎ ‎invisibility cloaks‎, ‎in particular within the frequencies spectrum‎ ‎of the visible light‎, ‎for instance the optical cloaking with‎ ‎metamaterials [9]‎, ‎full-wave invisibility of‎ ‎active devices at all frequencies [10]‎, ‎macroscopic invisibility cloak for visible light‎ [11]-[13]. ‎However‎, ‎as it is already warned by Pendry et‎. ‎al‎. ‎[1]‎, ‎the need for‎ ‎highly anisotropic inhomogeneous media is the main difficulty of‎ ‎their method‎. ‎This makes it difficult to realize these materials‎ ‎for real world applications‎. ‎To cope this limitation, some authors‎ ‎use conformal mapping to generate nearly isotropic cloaks‎, ‎for‎ ‎instance see‎ [14]-[16]. ‎The use of constrained optimization and quasi-conformal mapping [17]-[19]‎ is‎ ‎another alternative approach to overcome this limitation‎. ‎Computing‎ ‎the mapping based on the laplace equation‎ [20]-[21] was another idea to reduce‎ ‎the material anisotropy‎.

‎    The use of topology optimization method is an alternative way to‎ ‎control the propagation of electromagnetic waves in a‎ ‎heterogeneous media‎. ‎In [22]-[25]‎, ‎the topology optimization approach was used for the design of‎ ‎directional cloaks considering a plane wave with a specific ‎frequency‎. ‎In these works the optimal design problem leads to the‎ ‎solution of a distributed parameter identification problem under‎ ‎the Maxwell equation as a constraint which was computationally very ‎expensive‎. ‎In [26] a simple method for the design of‎ ‎directional cloaks has been suggested based on the laplace‎ ‎equation‎. ‎In this method‎, ‎the light trajectories are computed by‎ ‎the solution of the laplace equation and then the material‎ ‎properties‎, ‎the index of refraction‎, ‎was computed directly using‎ ‎the eikonal equation‎. ‎This method produces the isotropic material‎ ‎properties and does the perfect cloaking under the geometric‎ ‎optics approximation‎. ‎However‎, ‎there is no bound on the material‎ ‎properties‎. ‎

    In the present study‎, ‎the topology optimization method is adapted‎ ‎to formulate the original problem as a distributed parameter‎ ‎identification problem‎. ‎However‎, ‎unlike formerly mentioned works‎‎, ‎it is assumed that the light trajectories follow‎ ‎the eikonal equation‎. ‎Therefore‎, ‎we avoid to solve the Maxwell‎ ‎equation which is computationally very expensive‎. ‎Instead‎, ‎we solve the eikonal equation that can be solved very efficiently by ‎the fast marching or fast sweeping methods‎. ‎

    The remainder of this ‎paper is organized as follows‎. ‎The problem formulation is‎ ‎presented in the next section‎. ‎The sensitivity analysis and‎ ‎necessary optimality conditions are presented in section 3‎. ‎The‎ ‎solution algorithm is presented in section 4‎. ‎Section 5 is devoted‎ ‎to the space discretization and numerical algorithm‎. ‎The numerical‎ ‎results for a model problem is presented in section 6‎. ‎Finally‎, ‎section 7 summarizes this work and presents the outline of further ‎studies‎.

I.     problem formulation

Consider  () as an open bounded domain with sufficiently regular boundaries occupied with two isotropic materials. Assume the characteristic length scale of  is denoted by. Consider arbitrary electromagnetic wave with wavelengths (considering non-monochromatic waves includes no technical difficulty). It is well known that under assumption of the wave (light) trajectories obeys the following eikonal equation:

                                                                                        (1)

Where denotes the first arrival time of wave at , denotes the index of refraction and  denotes the set of point source in. Note that isosurfaces of are equivalent to the wave fronts, i.e. the wave trajectory at is orthogonal to -isosurface.

The goal of optimal design problem in the present study is to determine the optimal distribution of  in the design domain  such that the waves propagate according to a-priori desired trajectories in the control domain. For this purpose, the objective functional, to be minimized, is defined as the  distance of and target wave fronts over,

                                                                                                            (2)

 It is worth mentioning that the objective function can defined in more general form to result other desires. For instance, to make hidden from the light, i.e. to make invisibility cloaks for visible light spectrum, we can maximize the arrival time of light into , and on the other hand minimize the arrival time of light in the remained areas, i.e. in . For this purpose, the objective function can be defined alternatively as follows:

                                                                                               (3)

To make the optimal design manufacturable, is constrained to lies in the admissible design space  which is defined as follows:

                                                                                 (4)

Where denotes the admissible lower and upper bound for pointwise refractive index respectively. Therefor the optimal design problem in the present study can be summarized as follows:

                                                                                            (5)

To minimize the objective functional, the constrained gradient flow of (2) or (3) is used in the present study. To do this, it is required to compute the fist variation of objective functional (2) or (3) with respect to the design variable. Because after the space discretization, the number of design parameters will be possibly very large, the adjoint approach is adapted here to compute the first variation of (2) or (3) with respect to .

A.     Adjoint sensitivity analysis

Consider  as the Lagrange multiplier corresponding to PDE constraint (1). Ignoring the bound contestants  (they are treated by the projected gradient approach later), the local solutions of (5) are subsets of stationary points of the following augmented lagrangian:

                                                                                   (6)

Taking the first variation of with respect to  and equating it to zero will be identical to the direct PDE (1). The first variation of with respect to  will be computed as follows (consider objective functional (2)):

                                               (7)

Consider as a scalar field and as a vector field, applying identity  to the second term in right hand side of (6) results in (take and):

                                                                            (8)

Applying the divergence theorem on the first term on the right hand side of (8) results in:

                                                                            (9)

Where denotes the boundaries of, i.e.  and denotes the outward unit normal vector on. Considering (7) to gather with (9) and equating to zero the first variation of with respect to  results in the following adjoint PDE:

                                                             (10)

Where denotes the characteristic function corresponding to, i.e.  in and in. Note that the boundary conditions of adjoint PDE is computed by equating to zero the first term in the right hand side of (9). If alternatively the objective functional (3) is considered, by similar derivation we reach to the following adjoint PDE:

                                                      (11)

Finally, taking the first variation of with respect to leads to:

                                                                  (12)

Therefore, equating it to zero the first variation of with respect to results in:

                                                                                                               (13)

II.     first order necessary optimality conditions

Treating the bound constraints by the projected gradient approach (cf. [27], [28]), the first order KKT necessary optimality conditions (cf. [2]) corresponding to the optimization problem (5) can be expressed as follows:

where operator  denotes the orthogonal projection of trial fields onto design space . By straightforward computation, the explicit form of operator can be computed as follows (cf. [2]):

For the case of objective functional (3) the necessary first order optimality conditions will has the following form:

III.     solution algorithm

Considering results of section III, the solution for (5) based on the constrained gradient flow could be expressed as follows:

where  denotes the characteristic function of the design space, i.e.  in  and  in . Because the gradient is computed based on the norm, the regularity of the new design  may not be better than (is a square integrable field). In fact, it can possesses discontinuity in which is not desirable from the practical and computational point of views. It can decrease the convergence rate of algorithm 1, or even leads to oscillatory behavior in the course of computations. To overcome this problem, we perform a regularization step in algorithm 1 and lift the new design variable from the function space to the more regular space in which both of function and its first order derivative are square integrable. For this purpose, the following Helmholtz equation should be solved at every step:

                                                 (14)

whereis a regularization parameter that controls the smoothness of solution and  denotes the regularized design variable. According to the above discussion, algorithm 1 could be modified as follows:

IV.     numerical algorithm

For the purpose of numerical solution, the spatial domain is assumed to be a rectangular domain  with dimension  unites and is discretized into a uniform squared grid with edge length. The spatial variables are defined at center of finite difference cells, and  is used to denote the discretized counterpart of  at . Because, direct equation (1) is a nonlinear first order partial differential equation, using conventional methods for its solution leads to a nonlinear system of equations which is computational expensive to solve in practice. However, exploiting the hyperbolic and causality nature of this equation leads to a very efficient solution algorithm called the fast marching method [28] with  computational complexity; where. Later, the fast sweeping method has been introduced in [31] by Zhao

Fig.  1. Contour plot of first arrival time for iteration 10 (legend in time unit).

and it has been shown that it exhibits  computational complexity to achieve a moderate (reasonable) accuracy. Due to the high computational performance and ease of implementation, the fast sweeping method will be used in the present study to solve (1). The numerical method in the

present study is closely related to method introduced in [32]. This method will be briefly discussed in this paper in the appendix section.

V.     results and discussion

In this section, the performance and success of the presented algorithm is studied by a numerical experiment. For this purpose consider  as a rectangular domain with edge lengths and  units along and directions respectively. Moreover, assume that the spatial domain is decomposed into uniform square grids with grid size. The control domain, , is defined as a circle with radius located at the center of the spatial domain. The design domain, cloak in this example, is defined as a ring around the control domain with external radius. Our goal is to hide the control domain from light rays propagated from a point source which is located at the center of left edge of the spatial domain. For this purpose, the objective function is defined such that its minimizer leads to very large arrival time in and minimal arrival time in other regions, i.e.  objective functional (3) is considered here. Figs. 1-5 show the contour plot of first arrival times after 10, 20, 30, 50 and 100 iterations.

Fig. 2. Contour plot of first arrival time for iteration 20 (legend in time unit).

Fig. 3. Contour plot of first arrival time for iteration 30 (legend in time unit).      

Fig. 4. Contour plot of first arrival time for iteration 50 (legend in time unit).

Fig. 5. Contour plot of first arrival time for iteration 100 (legend in time unit).

Fig.  6. The variation of the objective function in the course of optimization.

According to plots the first arrival time is increased inside the control domain as the optimization proceeds. This observation confirms the success of the optimization algorithm to move toward the optimal solution. Fig. 6 shows the variation of objective functional as a function of iterations during the first 100 iterations. As the figure shows, the objective functional decreased monotonically by iterations. This observation confirms the convergence of the presented numerical algorithm.

VI.     conclusion

A mathematical model and its corresponding computational method is‎ ‎presented in this study to control the light propagation in a‎ ‎media by means of optimal determination of material inhomogeneity‎. ‎It‎ ‎is based on solving a distributed parameter identification problem ‎constrained to an eikonal equation‎. ‎In contrast to available ‎alternative methods‎, ‎the main benefit of the presented approach is ‎its ability to produce designs with bounded isotropic material‎ ‎properties‎. ‎The success of the presented method is shown by‎ ‎solving a model problem corresponding to the invisibility‎ ‎cloaking‎. ‎According to our numerical results the performance of‎ ‎optimized design‎, ‎the objective functional‎, ‎is improved by an‎ ‎order of magnitude by the presented numerical algorithm‎.

VII.     Appendix A

In this section the numerical solution algorithm for the eikonal equation and its corresponding adjoint equation used in the present study are briefly discussed.

Applying the Godunov numerical Hamiltonian to the eikonal equation (1), for , we have:

                                                                        (A.1)

where

and  denotes the positive part of, i.e. . The one sided finite difference method is used at the boundary points. The fast sweeping algorithm can be summarized as follows:

1. Initialize the point source condition  by assigning the exact value if  is a mesh point, or assigning to grid points near exact values which are computed by using the constant velocity at the point source. These values are fixed in later iterations. Assign larger positive values at all other grid points, and these values will be updated later.
2. Update the solution by Gauss-Seidel iterations with alternating sweeping (cf. [29]). At each grid point whose value was not fixed during the initialization, compute the candidate solution, denoted by, by (A.1) from the current values of its neighbors,  and then update  to be the smaller one between  and its former value; i.e., . We sweep the whole domain with four alternate ordering repeatedly: (i)  and ; (ii)  and ; (iii)  and ; (iv)  and . Here  and  are the running indices along  and directions respectively.
3. Test the convergence: given convergence criterion, check whether.

Following [30], the fast sweeping method is used in the present study to solve (11). Considering two dimensional case, equation (11) can be written in the following form:

                                                                                                               (A.2)

where ,  and  are function of  and , i.e. , , . The space discretized form of (A.2) can be expressed as follow:

                                        (A.3)

where  and  are determined according to the direction of light propagation. For instance, when, . Moreover,

Following [30], we introduce the following notation:

 Using this notation results in:

                                                   (A.4)

(A.4) can be written as follows:

                                        (A.5)

 (A.5) is in fact a system of linear equation that can be solved efficiently to a reasonable accuracy by the fast sweeping method.