Topology-optimized source shifter for optical location camouflaging

Garuda Fujii; Garuda Fujii

doi:10.1364/OE.503183

1. Introduction

Since the proposal of transformation optics [1,2], which is arguably one of the most innovative theories in optics, the concept of illusion optics [3] has spread widely through the enthusiastic interest of many researchers. Metastructures for generating illusions are usually designed under transformation theory, not only in optics but also in other fields through analogous physical scenarios [4]. To construct these metastructures, metamaterials with artificial microstructures are frequently used because of their high tunability in realizing the required specifications in material parameters. Some innovative devices, such as invisibility cloaks, have been developed [5,6].

Using source transformations, source illusions [7] have also been proposed [8]. Several attempts have demonstrated the manipulation of radiation patterns from sources [9] and generated illusions in which the number of sources are misperceived [10]. The source shifter is one such metastructure that generates a source illusion in camouflaging location so that observers viewing at a distance misperceive the location of the light source as if the light waves emanate from a different location. Camouflaging the location of an electromagnetic source has been proposed using metamaterials [11] and ultra-transparent photonic crystals [12]. However, the approaches for generating optical source illusions remain substantially limited to the above source transformation, which is based on the composition of metamaterials or photonic crystals.

In this work, we propose a topology optimization of light-source shifters to generate a high-quality illusion in location-camouflaging of a light source. Topology optimization [13] is one of the structural optimization methods with the highest design flexibility. This flexibility allows the generation of new pores that change the structural topology during optimization; higher performances are therefore expected through optimum design processing. To formulate the topology optimization of light-source shifters, we propose a scheme, an objective function, and constraints for generating a high-quality source illusion, and present a way to solve the optimization problem, overcoming difficulties such as multimodality, interdependence among design variables, and ill-posed-ness in topology optimization.

2. Methods

In this section, we describe the topology optimization problem to be solved in this work. In the following first subsection, we illustrate a design scheme for the light-source shifter and explain the domain settings of the optimization problem and the situations assumed in the numerical simulation. In the second subsection, we formulate the topology optimization problem, including an objective function, constraints, and structural modeling. The third subsection describes the algorithm behind the covariance matrix adaptation evolution strategy (CMA-ES) used as a method to search for an optimal solution to the proposed optimization problem.

2.1 Scheme

For the design problem concerning location camouflaging (see Fig. 1), we set two sources; one is a light source that is physically present at $\boldsymbol{x}_\mathrm {src}$ [Fig. 1(c)], the other is a virtual source at $\boldsymbol{x}_\mathrm {ref}$. Although the virtual light source is not physical, the source shifter causes an observer viewing from outside to misperceive the light source to be at $\boldsymbol{x}_\mathrm {ref}$, when in fact it is at $\boldsymbol{x}_\mathrm {src}$. The $x$-axis is set on the straight line connecting the locations of the source and the virtual source, $\boldsymbol{x}_\mathrm {src}$ and $\boldsymbol{x}_\mathrm {ref}$, and the $y$-axis in the middle between the two sources. For topology optimization as applied to optical location camouflaging, a circular fixed design domain $\Omega _\mathrm {D}$ of radius $R_\mathrm {D}$ is set at the center of the analysis domain. A cavity domain of air $\Omega _\mathrm {src}$ in which to place the light source [Fig. 1(c)] at $\boldsymbol{x}_\mathrm {src}=(R_\mathrm {D}/2, 0)$ is assumed. The structures of dielectrics, denoted as $\Omega _\mathrm {d}$, are assumed to be made of SiO$_2$; they are transformed in $\Omega _\mathrm {D}$ to camouflage the actual location of the light source at $\boldsymbol{x}_\mathrm {src}$, placing it instead at a virtual location perceived from a distance to be at $\boldsymbol{x}_\mathrm {ref}=(-R_\mathrm {D}/2, 0)$. The positional relation between the domains and the two sources is important because the dielectric structure $\Omega _\mathrm {d}$ changes within the fixed design domain while optimizing the structural topology. We set perfectly matched layers $\Omega _\mathrm {PML}$ around the whole analysis domain to implement an absorbing boundary condition for simulating light scattering in an open region.

Fig. 1. (a) Schematic illustrating topology optimization for optical location camouflaging. Domain sizes are set to $L_x=L_y=1.5\times R_\mathrm {D}$, $R_\mathrm {src}=R_\mathrm {D}/8$, $L_\mathrm {grid}= R_\mathrm {D}/60$. (b) Discretized level-set functions $\phi _j$ and structural expression highlighting the clear interface between SiO$_2$ and air. (c) Domain $\Omega _\mathrm {src}$ in which the real light source is located.

Download Full Size | PDF

2.2 Optimization problem

In this simulation, we assume $E_z$-polarized waves for which the electric field oscillates along the $z$-axis. The objective function for camouflaging the light source location is defined as:

(1)$$\Psi= \frac{1}{\Psi_\mathrm{bare}} \int_\mathrm{\Omega_\mathrm{out}} \big|E_z-E_\mathrm{ref} \big|^2 d\Omega,$$

where $E_z$ represents the electric field emanating from source location $\boldsymbol{x}_\mathrm {src}$ surrounded by a source shifter, $E_\mathrm {ref}$ the reference field emanating from a virtual source at $\boldsymbol{x}_\mathrm {ref}$ with no surrounding dielectrics to be reproduced for location camouflaging, and $\Omega _\mathrm {out}$ the outer domain of $\Omega _\mathrm {D}$ [Fig. 1(a)] in which $E_\mathrm {ref}$ is reproduced. By decreasing the value of this function, $E_z$ moves closer to $E_\mathrm {ref}$, reproducing $E_\mathrm {ref}$ despite the presence of the source at $\boldsymbol{x}_\mathrm {src}$, and thus camouflages the location of the light source. The objective function is normalized by $\Psi _\mathrm {bare}$, which is given as:

(2)$$\Psi_\mathrm{bare}=\int_\mathrm{\Omega_\mathrm{out}} \big| E_\mathrm{bare}-E_\mathrm{ref} \big|^2 d\Omega,$$

where $E_\mathrm {bare}$ denotes the electric field when a bare source is at $\boldsymbol{x}_\mathrm {src}$ without any surrounding dielectric structures. The light source is assumed to be an oscillating dipole emitting polarized dipole radiation $E_z$ governed by the differential equation:

(3)$$\nabla^2 E_z + \frac{\omega^2}{c^2}\epsilon(\boldsymbol{x})E_z ={-}\frac{\omega^2}{\epsilon_0c^2}P_\mathrm{dp}\delta(\boldsymbol{x}-\boldsymbol{x}_\mathrm{src}),$$

where $\omega$ denotes the frequency of the light wave, $c$ the speed of light in a vacuum, $\epsilon _0$ the permittivity of the vacuum, $\epsilon (\boldsymbol {x})$ the relative permittivity that depends on position $\boldsymbol{x}$, $P_\mathrm {dp}$ the $z$-component of the polarization vector of the dipole, and $\delta$ Dirac’s delta function. The total field $E_z$ is determined by the interference between the scattered and incident waves, i.e., $E_z=E_\mathrm {s}+E_\mathrm {i}$. Substituted into the above equation, we obtain:

(4)$$\nabla^2 E_\mathrm{s} + \frac{\omega^2}{c^2}\epsilon(\boldsymbol{x})E_\mathrm{s} ={-}\frac{\omega^2}{c^2}\big( \epsilon(\boldsymbol{x}) -\epsilon_\mathrm{air} \big) E_\mathrm{i},$$

where $\epsilon _\mathrm {air}$ denotes the relative permittivity of air (specifically, we assume $\epsilon _\mathrm {air}=1$), and $E_\mathrm {i}$ satisfies in the absence of dielectrics the differential equation for air:

(5)$$\nabla^2 E_ \mathrm{i}+ \frac{\omega^2}{c^2}\epsilon_ \mathrm{air}E_ \mathrm{i} ={-}\frac{\omega^2}{\epsilon_0c^2}P_\mathrm{dp}\delta(\boldsymbol{x}-\boldsymbol{x}_\mathrm{src}),$$

having solution:

(6)$$ E_{\mathrm{i}}=\frac{i \omega^2}{4 \varepsilon_0 c^2} P_{\mathrm{dp}} H_0^{(1)}\left(\frac{\omega}{c} \sqrt{\varepsilon_{\mathrm{air}}}\left\|\boldsymbol{x}-\boldsymbol{x}_{\mathrm{src}}\right\|\right) $$

where $H_0^{(1)}$ represents the Hankel function of the 1st kind, and $i=\sqrt {-1}$. For $\Omega _{\rm D}$, this design problem optimizes the distribution of $\epsilon (\boldsymbol {x})$ defined as:

(7)$$\begin{aligned}&\epsilon(\boldsymbol {x}) = \begin{cases} \epsilon_\mathrm{air}+\chi(\epsilon_\mathrm{d} -\epsilon_\mathrm{air}) & \mathrm{for}\ \boldsymbol {x} \in \Omega_{\rm D}\\ \epsilon_\mathrm{air} & \mathrm{for}\ \boldsymbol {x} \in \Omega_{\rm out}, \end{cases} \end{aligned}$$

where $\epsilon _\mathrm {d}$ is the relative permittivity of the dielectric, and $\chi$ the characteristic function defined as:

(8)$$\begin{aligned}\chi\big(\phi(\boldsymbol {x})\big) = \begin{cases} 1 & {\rm if}\quad \boldsymbol {x} \in \Omega_\mathrm{d}\\ 0 & {\rm if}\quad \boldsymbol{x}\in\Omega_{\rm D}\backslash\Omega_\mathrm{d}, \end{cases} \end{aligned}$$

where $\phi (\boldsymbol {x})$ denotes the level-set function defined on $\Omega _\mathrm {D}$ as:

(9)$$\begin{aligned} \Omega_\mathrm{D}\backslash \Omega_\mathrm{d} &= \{ \boldsymbol{x}\ |\ -1 \leq \phi(\boldsymbol {x}) < 0 \},\\& \Gamma = \{ \boldsymbol{x}\ |\ \phi(\boldsymbol {x}) = 0 \}, \\ \Omega_\mathrm{d}\backslash \Gamma &= \{ \boldsymbol{x}\ |\ 0 < \phi(\boldsymbol {x}) \leq 1 \}. \end{aligned}$$

The level-set function prescribes the distribution of the positive-valued domains of the function representing the dielectric structures, $\Omega _\mathrm {d}$:dielectric domain, the negative-valued domains representing the air-filled spaces outside the dielectric structures, $\Omega _\mathrm {D}\backslash \Omega _\mathrm {d}$:air domain, and the zero-isosurfaces representing the structural boundaries, $\Gamma$ [Fig. 1(b)]. The level-set function is discretized over grid points, and body-fitted finite elements are created along the boundaries [14] where wave scattering usually occurs. At each structure generation during an optimization computation, body-fitted finite elements are generated along the dielectric boundaries by remeshing. This clear boundary modeling is indispensable for an exact analysis of wave scattering. The grid size $L_\mathrm {grid}$ [Fig. 1(b)] relates to the size of the finite elements because the body-fitted elements are generated in each grid along the boundaries. A value of $L_\mathrm {grid}$ set to the shorter, finer finite elements is generated, and the accuracy of finite-element analyses is improved. In contrast, discretized level-set functions are arranged on the grid points; therefore, a shorter $L_\mathrm {grid}$ leads to an increase in the number of design variables.

Topology optimization is an ill-posed problem in which infinitely small structures can be generated. To regularize the ill-posed-ness, we employ a perimeter constraint that minimizes the value of the structural perimeter using the objective function and thereby the fitness which we regularize by:

(10)$$\underset{\phi}{\mathrm{inf}}\qquad F =\max\big( \Psi, \tau L_\mathrm{p} \big),$$

where $L_\mathrm {p}$ denotes the perimeter of the dielectric structures and $\tau$ the coefficient of regularization, defined as the ratio between $\Psi$ and $L_\mathrm {p}$. The value of $L_\mathrm {p}$ is computed from the total length of $\Gamma$ modeled by the iso-surfaces of level-set functions.

Because of the design problem required, a structural symmetry about the $x$-axis is implemented through a symmetrical implementation of the level-set function. Since the structure is symmetrical about the $x$-axis, discretized level-set functions at symmetrical positions to the $x$-axis have the same value. Then, the spatially discretized level-set functions over the grid points [Fig. 1(b)] located in the half-plane $y\geq 0$ become design variables. This symmetry reduces the number of design variables optimized by approximately half, but not exactly half because there is also a level-set function on the $x$-axis. The optimal set of design variables, $\boldsymbol {\phi }=\{ \phi _1, \ldots, \phi _j, \ldots, \phi _n \}$, where $n$ denotes the number of design variables, is explored for minimizing the fitness using the CMA-ES [15] with box-constraint handling [16], $-1\leq \phi _j \leq 1$.

2.3 CMA-ES

To overcome the difficult properties emerging in topology optimization problems that manipulate waves, we employ CMA-ES [15], which is one of the most powerful stochastic search algorithms among many evolutional approaches for exploring optimal design variables $\boldsymbol{\phi }$. CMA-ES adapts its parameters for a sampling distribution to the landscape of the minimized fitness function; the method is known as the robust search algorithm to combat difficult properties [17] in optimization problems such as multimodality and the interdependence of design variables, without implementing a trial-and-error approach required in making initial guesses. The methods for box-constraint handling [16,18] help to explore optimal solutions by CMA-ES by restricting the range of design variables to intervals $-1\leq \phi _j \leq 1$. The CMA-ES has been successfully employed in topology optimizations for wave manipulations [19–21].

Among some variations of the algorithm of CMA-ES [22,23], we use the following CMA-ES algorithm. At the first generation $g=0$, the distribution parameters are initialized:

(11)$$\boldsymbol{m}^{(0)} =0.5(b_\mathrm{u} +b_\mathrm{l}) \mathbf{1}_n, \hspace{3mm} \sigma^{(0)}=0.3(b_\mathrm{u}-b_\mathrm{l}), \hspace{3mm} \boldsymbol{C}^{(0)} =\boldsymbol{I}_n, \hspace{3mm} \boldsymbol{p}_{\sigma}^{(0)}= \boldsymbol{0}_n, \hspace{3mm} \boldsymbol{p}_{\boldsymbol C}^{(0)}=\boldsymbol{0}_n,$$

where $\boldsymbol{m}^{(g)}$ denotes the mean of the top $\mu$ solutions, $b_\mathrm {u}$ and $b_\mathrm {l}$ denote the upper and lower bounds of the design variable, respectively, $\boldsymbol{1}_n$ and $\boldsymbol{0}_n$ the $n$-dimensional vectors in which all elements are $1$ and $0$, respectively, $\sigma ^{(g)}$ denotes the spread of the sampling distribution, $\boldsymbol{C}^{(g)}$ the covariance matrix of the distribution for sampling, ${\boldsymbol I}_n$ the identity matrix of size $n$, and $\boldsymbol{p}_{\sigma }^{(g)}$ and $\boldsymbol{p}_{\boldsymbol C}^{(g)}$ denote evolution paths. Generating candidate solutions for sampling is based on the distribution parameters:

(12)$$\boldsymbol{\phi}_i^{(g+1)}=\boldsymbol{m}^{(g)} + \sigma^{(g)}\boldsymbol{BD} \boldsymbol{z}_{i}^{(g+1)},$$

where $\boldsymbol{\phi }_i^{(g+1)}$ denotes a candidate solution, $\boldsymbol{z}_{i}^{(g+1)}$ the vector of random numbers based on a multivariate normal distribution, denoted $\boldsymbol{z}_i^{(g+1)} \sim \mathcal {N}(\boldsymbol {0}_n, \boldsymbol {I}_n)$, $\boldsymbol{B}$ an orthogonal matrix, the columns of which are the eigenvectors of $\boldsymbol{C}^{(g)}$, and $\boldsymbol{D}$ a diagonal matrix, the elements of which equal the square root of the eigenvalues of $\boldsymbol{C}^{(g)}$. The eigen-decomposition of $\boldsymbol{C}^{(g)}$ is used to obtain $\boldsymbol{B}$ and $\boldsymbol{D}$ as

(13)$$\boldsymbol{C}^{(g)} = \boldsymbol{B}\boldsymbol{D}^2\boldsymbol{B}^\top\enspace.$$

Only when a certain condition is satisfied, indicating that $\boldsymbol{C}$ has been sufficiently learned, do we perform an eigen-decomposition because of the high computational cost of eigen-decompositions; specifically,

(14)$$10n( c_1 + c_\mu ) \Delta g > 1,$$

where $\Delta g$ denotes the number of generations that has passed since the last eigenvalue decomposition, $c_1$ the learning rates of the rank-one update, and $c_\mu$ the learning rates of the rank-$\mu$ update. Using $F$ with an adaptive penalty function for handling the box constraint [16], $-1\leq \phi _j \leq 1$, the candidate solutions are ranked. The differences between the mean vector $\boldsymbol{m}^{(g)}$ and the $i$-th-ranked candidate out of $n_\mathrm {smp}$-solutions at the next generation $g+1$, $\boldsymbol{\phi }_{i:n_\mathrm {smp}}^{(g+1)}$, are rescaled as

(15)$$\boldsymbol{y}_{i:n_\mathrm{smp}}^{(g+1)} =\left(\boldsymbol{\phi}_{i:n_\mathrm{smp}}^{(g+1)}-\boldsymbol{m}^{(g)}\right)\Big/\sigma^{(g)}.$$

Updating the distribution parameters is performed using the procedure:

(16)$$\boldsymbol{m}^{(g+1)} = \sum_{i=1}^\mu w_i \boldsymbol{\phi}_{i:n_\mathrm{smp}}^{(g+1)},$$

(17)$$\boldsymbol{C}^{(g+1)} =\boldsymbol{C}^{(g)}+h_\sigma^{(g+1)} c_1\bigg( \boldsymbol{p}_{\boldsymbol C}^{(g+1)} \Big(\boldsymbol{p}_{\boldsymbol C}^{(g+1)} \Big)^\top \hspace{-1mm} -\boldsymbol{C}^{(g)} \bigg) +c_\mu\hspace{-1mm}\sum_{i=1}^{\mu}w_i \bigg( \boldsymbol{y}_{i:n_\mathrm{smp}}^{(g+1)} \Big(\boldsymbol{y}_{i:n_\mathrm{smp}}^{(g+1)} \Big)^\top \hspace{-1mm} -\boldsymbol{C}^{(g)}\bigg),$$

(18)$$\sigma^{(g+1)} = \sigma^{ (g)}\ \exp\Biggl(\frac{c_{\mathrm \sigma}}{d_\sigma} \biggl(\frac{||\boldsymbol{p}_{\sigma}^{(g+1)}||} {\mathbb{E}[|| \mathcal{N}(\boldsymbol{0}_n,{\boldsymbol I}_n)||]}-1\biggr)\Biggr),$$

where $\mu$ denotes the number of higher-ranked solutions, $w_i$ recombination weights, $c_\sigma$ the learning rate for the cumulations of the step-size adaptation, $d_\sigma$ the damping parameter for the step-size adaptation, and $\mathbb {E}[|| \mathcal {N} (\boldsymbol {0}_n,{\boldsymbol I}_n)||]$ the expected value of the norm of an $n$-variate normal distribution; the expectation is approximated using $\mathbb {E}[|| \mathcal {N} (\boldsymbol {0}_n,{\boldsymbol I}_n)||]=\sqrt {n} \left (1 - \frac {1}{4n} + \frac {1}{ 21n^2}\right )$. The Heaviside step function $h_\sigma ^{(g+1)}$ is defined as

(19)$$\begin{aligned}h_\sigma^{(g+1)} = \begin{cases} 1 & \ || \boldsymbol{p}_{\sigma}^{(g+1)} || < \left(1.5 + \frac{1}{n-0.5}\right)\mathbb{E}[|| \mathcal{N}(\boldsymbol{0}_n,{\boldsymbol I}_n)||] \\ 0 & \ \rm{otherwise}, \end{cases} \end{aligned}$$

and the evolution paths are updated using the scheme

(20)$$\begin{aligned}\boldsymbol{p}_{\sigma}^{(g+1)}=(1-c_{\mathrm \sigma})\boldsymbol{p}_{\sigma}^{(g)} + \sqrt{c_{\mathrm \sigma}(2-c_{\mathrm \sigma})\mu_\mathrm{eff}} \sqrt{\boldsymbol{C}^{(g)}}^{{-}1}\ \frac{\boldsymbol{m}^{(g+1)}-\boldsymbol{m}^{(g)}}{\sigma^{(g)}}, \end{aligned}$$

(21)$$\begin{aligned}\boldsymbol{p}_{\boldsymbol C}^{(g+1)} = (1 - c_{\boldsymbol C})\boldsymbol{p}_{\boldsymbol C}^{(g)} +h_\sigma^{(g+1)}\sqrt{c_{\boldsymbol C}(2-c_{\boldsymbol C})\mu_\mathrm{eff}} \ \frac{\boldsymbol{m}^{(g+1)} -\boldsymbol{m}^{(g)}}{\sigma^{(g)}}, \end{aligned}$$

where $c_{\boldsymbol {C}}$ denotes the learning rate for the cumulations of the rank-one update, $\mu _\mathrm {eff}$ the variance effective selection mass for the mean, and $\sqrt {\boldsymbol {C}^{(g)}}^{-1}$ is defined as

(22)$$\sqrt{\boldsymbol{C}^{(g)}}^{{-}1}=\boldsymbol{B}\boldsymbol{D}^{{-}1}\boldsymbol{B}^\top.$$

We list the default values of the strategy parameters in CMA-ES [24] in Table 1.

Table 1. Default strategy parameters of the CMA-ES

View Table | View all tables in this article

The optimization continues until the following convergence criterion:

(23)$$\Delta m=\frac{||\boldsymbol{m}^{(g+1)}-\boldsymbol{m}^{(g)}||}{||\boldsymbol{m}^{(g+1)}||} \leq 0.01,$$

obtained using the error of the mean vector $\boldsymbol{m}^{(g)}$ is satisfied.

3. Results

In this section, we show the result of topology optimizations. The outline of this section is as follows: The first subsection introduces the reference electric field of a light wave emanating from a bare source located at $\boldsymbol{x}_\mathrm {ref}$. This reference field is a target field to be reproduced outside the light-source shifter under the source location at $\boldsymbol{x}_\mathrm {src}$ through the topology optimization presented above. Also, the electric field of a light wave emanating from a bare source at $\boldsymbol{x}_\mathrm {src}$ is evaluated. The field is used for the normalization of the objective function that evaluates the performance of the light-source shifter. The proposed topology optimization is demonstrated in the second subsection under various settings of the regularization coefficient $\tau$ and distance between the two source locations. Dependencies of the camouflaging performance on frequency of light and relative permittivity of the structural material are investigated. The third subsection presents attempts at optimum designs of light-source shifters operating at multiple frequencies.

3.1 Reference field and the electric field emitted by a bare source

We show results for reference fields in Fig. 2(a)–(c). If a light source is located at the virtual location $\boldsymbol{x}_\mathrm {ref}$ [Fig. 2(a)], the electric field emanating from $\boldsymbol{x}_\mathrm {ref}$ becomes the reference field to be reproduced for location camouflaging [Fig. 2(b)]. The evaluated difference in the objective function then vanishes completely everywhere [Fig. 2(c)]. If the source location changes from $\boldsymbol{x}_\mathrm {ref}$ to $\boldsymbol{x}_\mathrm {src}$ [Fig. 2(d)], the electric field emitted by the source also changes from $E_\mathrm {ref}$ [Fig. 2(b)] to $E_\mathrm {bare}$ [Fig. 2(e)]. The difference between the two fields is then evident [Fig. 2(f)]. Because of the prescribed normalization of the objective function, the value of the objective function for the bare source at $\boldsymbol{x}_\mathrm {src}$ is that shown in Fig. 2(f). We note that the normalized frequency is set to $\omega R_\mathrm {D}/(2\pi c)= 2.5$ corresponding to the wavelength of $R_\mathrm {D}/2.5$; then, one wavelength of light in air is numerically expressed by approximately 24 finite elements under discretization $L_\mathrm {grid}=R_\mathrm {D}/60$ [Fig. 1(b)].

Fig. 2. Results for the reference and normalization. (a) a light source at the virtual location $\boldsymbol{x}_\mathrm {ref}$ covered by $\Omega _{\text {D}}$ filled with air, (b) the reference field $E_\mathrm {ref}$, (c) the absolute difference evaluated for $E_z$ corresponding to $E_\mathrm {ref}$, (d) a bare light source at actual location $\boldsymbol{x}_\mathrm {src}$, (e) the electric field emanating from a bare source $E_\mathrm {bare}$, and (f) the absolute difference evaluated for normalizing the objective function. The normalized frequency is set to $\omega R_\mathrm {D}/(2\pi c)=2.5$. Note that the source location $\boldsymbol{x}_\mathrm {src}$ in Eq. (6) changes to $\boldsymbol{x}_\mathrm {ref}$ in the results for reference fields (a)–(c).

Download Full Size | PDF

3.2 Topology optimization for light-source shifters

We show the histories of the values of $F$, $\Delta m$, $\Psi$, and $L_\mathrm {p}$ in Fig. 3. Also, intermediate topologies during optimizations, optimal configurations, and their camouflaging performances are shown in Fig. 4. Figure 3(a) shows the best fitness value $F$ out of $n_\mathrm {smp}$-candidates at each generation $g$. Although the $F$-value exhibits a slight fluctuation caused by sampling generations based on random numbers, the value decreases almost monotonically. The relative movement of the mean vector, $\Delta m$ [Fig. 3(b)] used for the termination criteria in Eq. (23) ultimately reaches $\Delta m \leq 0.01$; the topology optimizations demonstrated are then judged to have converged. The objective function for location camouflage also improved as $g$ proceeds [Fig. 3(c)] under the perimeter constraint reducing the structural perimeter [Fig. 3(d)].

Fig. 3. History of (a) the fitness, (b) the convergence error for termination criteria, (c) the objective function for camouflaging performance, and (d) the perimeter of dielectric structures for regularization.

Download Full Size | PDF

Fig. 4. Results of topology optimizations: (a, e, i) transition of the configuration during topology optimization, (b, f, j) topology-optimized source shifters composed of SiO$_2$, (c, g, k) electric field around the light-source shifters, (d, h, l) the evaluated difference between electric field manipulated by a topology-optimized source shifter, and the reproduced reference field. We assume SiO$_2$ as the dielectric material with permittivity $\epsilon _\mathrm {d}=3.90$. The regularization coefficient is set to (a–d) $\tau =0.001$, (e–h) $\tau =0.0001$, and (i–l) $\tau =0.00001$. The number of samplings in the CMA-ES is $n_\mathrm {smp}=450$; the number of design variables is then $n=5642$.

Download Full Size | PDF

By tuning the coefficient $\tau$ in the perimeter constraint, priority between simplifying optimum configurations and improving camouflaging performance can be adjusted. Larger $\tau$ prioritizes reducing $L_\mathrm {p}$ [$\tau =0.001$ in Fig. 3(d)] rather than reducing $\Psi$ [$\tau =0.001$ in Fig. 3(c)] and provides a simpler configuration with a smaller perimeter [Fig. 4(b)]. Smaller $\tau$ in contrast improves substantially $\Psi$ [$\tau =0.00001$ in Fig. 3(c)], with less priority given to reducing the perimeter [$\tau =0.00001$ in Fig. 3(d)]; it also provides a complex configuration with a larger perimeter but achieves superior performance by prioritizing performance enhancements [Fig. 4(j)–(l)]. All optimal configurations succeed in reproducing $E_\mathrm {ref}$ in the outer domain $\Omega _\mathrm {out}$ [Fig. 4(c)(g)(k)], and observers viewing at a distance from the source shifters perceive a virtual light source at $\boldsymbol{x}_\mathrm {ref}$ as the light source. The minimized difference between the two fields nears zero in $\Omega _\mathrm {out}$ [Fig. 4(d)(h)(l)] and, in the best case [Fig. 4(j)–(l)], the value of the objective function falls below $0.05{\%}$ of that when a bare source is present at $\boldsymbol{x}_\mathrm {src}$ [Fig. 2(f)].

How far can the virtual source location be separated from the actual source location in the generation of source illusions by topology optimization? We demonstrate topology optimizations for different distances between virtual and actual locations, $||\boldsymbol {x}_\mathrm {src}-\boldsymbol {x}_\mathrm {ref}||$, and show the relationship between distance and camouflaging performance from optimal configurations obtained [Fig. 5]. Despite the increase in distance between virtual and actual locations, the objective function value remains less than 2% of that without the source shifter and does not degrade performance.

Fig. 5. Distance between actual and virtual sources, normalized by the wavelength in air $\lambda$, versus camouflaging performances. The normalized frequency is set to $\omega R_\mathrm {D}/(2\pi c)=R_\mathrm {D}/\lambda =2.5$. The coefficient of regularization is set to $\tau =0.001$. The number of samplings for CMA-ES is $n_\mathrm {smp}=450$; the number of design variables is then $n=5642$. Note that the value for normalizing the objective function, $\Psi _\mathrm {bare}$, is computed for each $||\boldsymbol {x}_\mathrm {src}- \boldsymbol {x}_\mathrm {ref}||$.

Download Full Size | PDF

For the topology-optimized source shifters [Fig. 4], Fig. 6 displays the responses in camouflaging performance for frequency $\omega$ and permittivity $\epsilon _\mathrm {d}$. Only a single dip in the objective function is found at the surface of the response; this implies that the topology-optimized source shifters work only at the frequency and only for the permittivity of the material assumed in the topology optimizations demonstrated.

Fig. 6. Response of the camouflaging performance of the optimum configuration to frequency and relative permittivity: (a) $\tau =0.001$, (b) $0.0001$, and (c) $0.00001$.

Download Full Size | PDF

3.3 Topology-optimized light-source shifters operating for multifrequency

Here, we widen the camouflaging frequency by minimizing the fitness for improving the performance at multiple frequencies,

(24)$$\underset{\phi}{\mathrm{inf}}\quad F_\mathrm{mf} =\max\big( \Psi_{\omega_1^\mathrm{n}}, \ldots, \Psi_{\omega_N^\mathrm{n}}, \tau L_\mathrm{p} \big),$$

where $\Psi _{\omega _k^\mathrm {n}}$ denotes the value of the objective function at normalized frequency $\omega _k^\mathrm {n}=\omega _k R_\mathrm {D}/(2\pi c)$ with $k$ indexing frequencies, and $N$ the number of frequencies. We show the results obtained by minimizing $F_\mathrm {mf}$ in Fig. 7, for three frequencies $\omega R_\mathrm {D}/(2\pi c)=2.45, 2.50, 2.55$ in Fig. 7(a)–(e), and for five frequencies $\omega R_\mathrm {D}/(2\pi c)=2.40, 2.45, 2.50, 2.55, 2.60$ in Fig. 7(f)–(l). The banding of the camouflaging frequency is evident in Fig. 7(b,g). For all analyzed frequencies, the values of the objective function fall below 3% in Fig. 7(c)–(e) and 5% in Fig. 7(h)–(l). These values are not particularly outstanding, but the distribution of the electric field shows that the reference field is well reproduced, and therefore the performance is considered satisfactory. Despite the coefficient of regularization being set larger ($\tau =0.001$) to prioritize reducing the perimeter, the optimal configurations obtained in Fig. 7(a,f) have fine structures and remain complex structures exhibiting a large perimeter; this implies that realizing location camouflaging for multiple frequencies is difficult. We note that when three frequencies $\omega R_\mathrm {D}/(2\pi c)=2.40, 2.50, 2.60$ in Eq. (24) are set, three distinct dips emerge discontinuously in the performance response, meaning that generating a frequency band for location camouflaging fails.

Fig. 7. Topology optimizations for light-source shifters operating multifrequency. Normalized frequencies are set to (a)–(e) $\omega R_\mathrm {D}/(2\pi c)=2.45, 2.50$, and $2.55$ and (f)–(l) $\omega R_\mathrm {D}/(2\pi c)=2.40, 2.45, 2.50, 2.55$, and $2.60$. The coefficient of regularization is set to $\tau =0.001$. The number of samplings for the CMA-ES is $n_\mathrm {smp}=450$; the number of design variables then becomes $n=5642$.

Download Full Size | PDF

4. Discussion

We proposed a topology optimization for designing a light-source shifter for optical location camouflage and presented the scheme, the objective function, the constraints, and the structural model with a method to solve the problem. To reproduce the reference field emanating from different source locations than the actual one, the difference between the electric field controlled by a light-source shifter and the reference field outside the light-source shifter is defined as the minimized objective function. The use of CMA-ES for exploring optimal topology enables designers to avoid local minimum with poor camouflaging performance and succeeds in reaching optimal configurations with high performance for location camouflaging. Each optimal configuration obtained may possibly be one of many local minima (see Appendix A) but each exhibits a sufficiently high performance. Robust performances for location camouflage over a broad band were attempted and actually realized numerically through topology optimization incorporating multiple frequencies. Our topology optimization solved the two-dimensional structural design problem assuming that the dielectric structures are infinitely long, here for the $z$-direction. Sufficient height of the dielectric structure in the $z$-direction and an accurate measurement of the electric field at the propagation plane in the $xy$-direction are necessary to achieve the location camouflaging in experiments described in this paper. Also, with respect to the light source, the line current is ideal in realizing the two-dimensional set-up. We provided stereolithography (STL) data of the optimal configurations presented (see Appendix B). This STL data can be utilized not only for cross-verifying the camouflaging performance of the optimal configurations numerically but also for fabricating samples to demonstrate experiments. Light sources that can emit light over a broad spectral band have been developed [25,26], which may help in experimental demonstrations of location camouflaging over a wide bandwidth.

Several issues remain in topology optimization of this source shifter from both numerical and functional aspects. The governing equation solved was derived assuming dipole radiation comprising $E_z$-polarized harmonically oscillating waves; the camouflaging performance for the $H_z$-polarized waves or that in the transient state was not improved in the optimizations presented. Furthermore, although the optimal configurations become more complex and have fine structures, the robustness of performance was improved over optical frequencies, but in this topology optimization the robustness of performance to fluctuations in permittivity of the structural material was not considered. Using CMA-ES for topology optimization frees the designer from trial-and-error guessing of appropriate initial configurations and tuning strategy parameters, but requires many numerical analyses equal to the number of samples multiplied by the number of generations to converge. Therefore, developing a method to accelerate convergence is desirable.

The topology optimization presented is not limited solely to the design of the light-source shifter. By changing the scheme and the objective function, it can be applied to various two-dimensional design problems such as designing the cross section of optical fiber [27,28]. Challenges for analogous source illusions [29,30] and illusion multiphysics [20,31] are expected to be addressed in the future.

5. Conclusion

In summary, a topology optimization of light-source shifters was proposed for optical location camouflaging. We used the level-set method of generating body-fitted finite elements for modeling dielectric structures, the CMA-ES for overcoming the difficult properties encountered in optimization problems, and the perimeter constraint for regularizing the ill-posed-ness in the topology optimizations presented. Designing light-source shifters with high performance was successful and reproducing numerically the reference field outside the optimal configurations became possible. The camouflaging performance of the optimum configuration did not deteriorate as the distance between actual and virtual locations increased, as long as both locations were set within a fixed design domain. Simultaneous minimization of the objective function at several frequencies also led to the formation of a band of camouflaging frequencies.

Appendix A

Multimodal properties of topology optimization for light-source shifters

To demonstrate the multimodal property for the topology optimization problems presented, we presented several numerical results of topology optimizations demonstrated under different patterns of sampling generation in CMA-ES. In the results demonstrated in Fig. 8, topology optimizations using a different set of random numbers $\boldsymbol{z}_{i}^{(g+1)}$ converged to different optimal configurations. The fact that different optimal configurations were obtained with a different random number-generation means that this topology optimization problem has strong multimodality, and the obtained optimal configurations may be among several local optimal solutions. However, as evident in Fig. 4(d) and Fig. 8(d,h,l), all optimal configurations obtained with different random number-generation achieved sufficiently good performances to realize optical location camouflage without guessing an initial solution and adjusting strategy parameters. Thus, the effectiveness of our method is established.

Fig. 8. Various optimal configurations obtained under different patterns of random number-generation, (a–d) case 1, (e–h) case 2, and (i–l) case 3 in CMA-ES. The coefficient of regularization is set to $\tau =0.001$. The number of samplings in the CMA-ES is $n_\mathrm {smp}=450$; the number of design variables is then $n = 5642$.

Download Full Size | PDF

Appendix B

STL data of the presented optimal configurations

We provided STL data of the optimal configurations of light-source shifters in Ref. [32]. As shown in Fig. 9, we set radius and height of the fixed design domain as $R_\mathrm {D}=1.5\ \mu$m and $H_\mathrm {D}=0.5\ \mu$m, respectively. This size of the optimal configurations help set the operating frequency $\omega R_\mathrm {D}/(2\pi c)=2.5$ as 600 nm wavelength of light. A list of optimal configurations, the name of the corresponding STL files, and the operating normalized frequency and wavelength of light is given in Table 2.

Fig. 9. STL data of the presented optimal configurations.

Download Full Size | PDF

Table 2. Optimum structures, STL data, and operating wavelength and frequency

View Table | View all tables in this article

Acknowledgements

We thank Richard Haase, Ph.D., from Edanz Group for editing a draft of this manuscript.

Disclosures

The authors declare no conflicts of interest.

Data availability

Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.

References

1. U. Leonhardt, “Optical conformal mapping,” Science 312(5781), 1777–1780 (2006). [CrossRef]

2. J. B. Pendry, D. Schurig, and D. R. Smith, “Controlling electromagnetic fields,” Science 312(5781), 1780–1782 (2006). [CrossRef]

3. Y. Lai, H. Chen, Z.-Q. Zhang, and C. T. Chan, “Complementary media invisibility cloak that cloaks objects at a distance outside the cloaking shell,” Phys. Rev. Lett. 102(9), 093901 (2009). [CrossRef]

4. M. Wegener, “Metamaterials beyond optics,” Science 342(6161), 939–940 (2013). [CrossRef]

5. A. Alù and N. Engheta, “Achieving transparency with plasmonic and metamaterial coatings,” Phys. Rev. E 72(1), 016623 (2005). [CrossRef]

6. D. Schurig, J. J. Mock, B. J. Justice, et al., “Metamaterial electromagnetic cloak at microwave frequencies,” Science 314(5801), 977–980 (2006). [CrossRef]

7. Y. Luo, J. Zhang, L. Ran, et al., “Controlling the emission of electromagnetic source,” PIERS Online 4(7), 795–800 (2008). [CrossRef]

8. N. Kundtz, D. A. Roberts, J. Allen, et al., “Optical source transformations,” Opt. Express 16(26), 21215–21222 (2008). [CrossRef]

9. J. Allen, N. Kundtz, D. A. Roberts, et al., “Electromagnetic source transformations using superellipse equations,” Appl. Phys. Lett. 94(19), 194101 (2009). [CrossRef]

10. H. Chen, Y. Xu, H. Li, and T. Tyc, “Playing the tricks of numbers of light sources,” New J. Phys. 15(9), 093034 (2013). [CrossRef]

11. J. Yi, P.-H. Tichit, S. N. Burokur, and A. de Lustrac, “Illusion optics: Optically transforming the nature and the location of electromagnetic emissions,” J. Appl. Phys. 117(8), 084903 (2015). [CrossRef]

12. Z. Yao, J. Luo, and Y. Lai, “Illusion optics via one-dimensional ultratransparent photonic crystals with shifted spatial dispersions,” Opt. Express 25(25), 30931–30938 (2017). [CrossRef]

13. M. P. Bendsøe and N. Kikuchi, “Generating optimal topologies in structural design using a homogenization method,” Comput. Methods Appl. Mech. Eng. 71(2), 197–224 (1988). [CrossRef]

14. G. Fujii, H. Watanabe, T. Yamada, et al., “Level set based topology optimization for optical cloaks,” Appl. Phys. Lett. 102(25), 251106 (2013). [CrossRef]

15. N. Hansen and A. Ostermeier, “Completely derandomized self-adaptation in evolution strategies,” Evol. Comput. 9(2), 159–195 (2001). [CrossRef]

16. G. Fujii, M. Takahashi, and Y. Akimoto, “CMA-ES-based structural topology optimization using a level set boundary expression—Application to optical and carpet cloaks,” Comput. Methods Appl. Mech. Eng. 332, 624–643 (2018). [CrossRef]

17. L. M. Rios and N. V. Sahinidis, “Derivative-free optimization: a review of algorithms and comparison of software implementations,” J. Glob. Optim. 56(3), 1247–1293 (2013). [CrossRef]

18. N. Hansen, A. S. P. Niederberger, L. Guzzella, and P. Koumoutsakos, “A method for handling uncertainty in evolutionary optimization with an application to feedback control of combustion,” IEEE Trans. Evol. Computat. 13(1), 180–197 (2009). [CrossRef]

19. G. Fujii, M. Takahashi, and Y. Akimoto, “Acoustic cloak designed by topology optimization for acoustic-elastic coupled systems,” Appl. Phys. Lett. 118(10), 101102 (2021). [CrossRef]

20. G. Fujii and Y. Akimoto, “Electromagnetic-acoustic biphysical cloak designed through topology optimization,” Opt. Express 30(4), 6090–6106 (2022). [CrossRef]

21. G. Fujii, “Camouflaging the location of a sound source via topology-optimized source-shifter,” Journal of Sound and Vibration 559, 117768 (2023). [CrossRef]

22. R. Ros and N. Hansen, “A simple modification in CMA-ES achieving linear time and space complexity,” in Parallel Problem Solving from Nature – PPSN X, (Springer Berlin Heidelberg, 2008), pp. 296–305.

23. D. V. Arnold and N. Hansen, “Active covariance matrix adaptation for the (1+1)-CMA-ES,” in Proceedings of the 12th Annual Conference on Genetic and Evolutionary Computation," (Association for Computing Machinery, 2010), GECCO ’10, pp. 385–392.

24. N. Hansen, “The CMA evolution strategy: a tutorial,” arXiv, arXiv:1604.00772 (2016). [CrossRef]

25. K. Xu, L. Huang, Z. Zhang, et al., “Light emission from a poly-silicon device with carrier injection engineering,” Mater. Sci. Eng., B 231, 28–31 (2018). [CrossRef]

26. K. Xu, “Silicon electro-optic micro-modulator fabricated in standard cmos technology as components for all silicon monolithic integrated optoelectronic systems,” J. Micromech. Microeng. 31(5), 054001 (2021). [CrossRef]

27. M. F. Domingues, C. A. Rodriguez, J. Martins, et al., “Cost-effective optical fiber pressure sensor based on intrinsic fabry-perot interferometric micro-cavities,” Opt. Fiber Technol. 42, 56–62 (2018). [CrossRef]

28. A. Leal-Junior, L. Avellar, V. Biazi, et al., “Multifunctional flexible optical waveguide sensor: on the bioinspiration for ultrasensitive sensors development,” Opto-Electron. Adv. 5(10), 210098 (2022). [CrossRef]

29. Y. Liu, Z. Liang, F. Liu, et al., “Source illusion devices for flexural lamb waves using elastic metasurfaces,” Phys. Rev. Lett. 119(3), 034301 (2017). [CrossRef]

30. X.-D. Fan, B. Liang, J. Yang, and J.-C. Cheng, “Illusion for airborne sound source by a closed layer with subwavelength thickness,” Sci. Rep. 9(1), 1750 (2019). [CrossRef]

31. G. Fujii and Y. Akimoto, “Optimizing the structural topology of bifunctional invisible cloak manipulating heat flux and direct current,” Appl. Phys. Lett. 115(17), 174101 (2019). [CrossRef]

32. G. Fujii, “STL data: Topology-optimized source shifter for optical location camouflaging,” figshare2023https://doi.org/10.6084/m9.figshare.24239320.v2..

Selection and Recombination	$μ = ⌊ n_{s m p} / 2 ⌋$ , $w_{i} = \frac{\ln ((n_{s m p} + 1) / 2) - \ln i}{\sum_{j = 1}^{μ} (\ln ((n_{s m p} + 1) / 2) - \ln j)}$ , $μ_{e f f} = \frac{{(\sum_{i = 1}^{μ} w_{i})}^{2}}{\sum_{i = 1}^{μ} w_{i}^{2}}$ ,
Step Size Adaptation	$c_{σ} = \frac{μ_{e f f} + 2}{n + μ_{e f f} + 5}$ , $d_{σ} = 1 + c_{σ} + 2 m a x (0, \sqrt{\frac{μ_{e f f} - 1}{n + 1}} - 1)$ ,
Covariance Matrix Adaptation	$c_{C} = \frac{4 + μ_{e f f} / n}{n + 4 + 2 μ_{e f f} / n}$ , $c_{1} = \frac{2}{(n + 1.3)^{2} + μ_{e f f}}$ ,
	$c_{μ} = min (1 - c_{1}, 2 \frac{μ_{e f f} - 2 + 1 / μ_{e f f}}{(n + 2)^{2} + μ_{e f f}})$ .

Selection and Recombination	$μ = ⌊ n_{s m p} / 2 ⌋$ , $w_{i} = \frac{\ln ((n_{s m p} + 1) / 2) - \ln i}{\sum_{j = 1}^{μ} (\ln ((n_{s m p} + 1) / 2) - \ln j)}$ , $μ_{e f f} = \frac{{(\sum_{i = 1}^{μ} w_{i})}^{2}}{\sum_{i = 1}^{μ} w_{i}^{2}}$ ,
Step Size Adaptation	$c_{σ} = \frac{μ_{e f f} + 2}{n + μ_{e f f} + 5}$ , $d_{σ} = 1 + c_{σ} + 2 m a x (0, \sqrt{\frac{μ_{e f f} - 1}{n + 1}} - 1)$ ,
Covariance Matrix Adaptation	$c_{C} = \frac{4 + μ_{e f f} / n}{n + 4 + 2 μ_{e f f} / n}$ , $c_{1} = \frac{2}{(n + 1.3)^{2} + μ_{e f f}}$ ,
	$c_{μ} = min (1 - c_{1}, 2 \frac{μ_{e f f} - 2 + 1 / μ_{e f f}}{(n + 2)^{2} + μ_{e f f}})$ .

Topology-optimized source shifter for optical location camouflaging

Abstract

1. Introduction

2. Methods

2.1 Scheme

2.2 Optimization problem

2.3 CMA-ES

3. Results

3.1 Reference field and the electric field emitted by a bare source

3.2 Topology optimization for light-source shifters

3.3 Topology-optimized light-source shifters operating for multifrequency

4. Discussion

5. Conclusion

Appendix A

Appendix B

Acknowledgements

Disclosures

Data availability

References

Data availability

Cited By

Figures (9)

Tables (2)

Equations (24)

Optics Express

Opt. Struct.	STL	Assumed normalized frequency in optimization	Operating wavelength [nm]
Fig. 4(b)	Fig4b.stl	2.50	600
Fig. 4(f)	Fig4f.stl	2.50	600
Fig. 4(j)	Fig4j.stl	2.50	600
Fig. 5(a)	Fig5a.stl	2.50	600
Fig. 5(b)	Fig5b.stl	2.50	600
Fig. 5(c)	Fig5c.stl	2.50	600
Fig. 5(d)	Fig5d.stl	2.50	600
Fig. 5(e)	Fig5e.stl	2.50	600
Fig. 5(f)	Fig5f.stl	2.50	600
Fig. 7(a)	Fig7a.stl	2.45, 2.50, 2.55	612, 600, 588
Fig. 7(f)	Fig7f.stl	2.40, 2.45, 2.50, 2.55, 2.60	625, 612, 600, 588, 577
Fig. 8(b)	Fig8b.stl	2.50	600
Fig. 8(f)	Fig8f.stl	2.50	600
Fig. 8(j)	Fig8j.stl	2.50	600