The development of invisibility cloaks is a longstanding goal in electromagnetics [1–7]. An object that is invisible to electromagnetic waves would have widespread application in both military and commercial systems. At radio frequencies, cloaked objects have a zero radar cross section [8,9], while optical displays with low scattering and reflection provide an improved viewing experience [10,11]. It has also been proposed that cloaked sensors may provide enhanced sensitivity because of a reduced interaction with the environment [12–16].

More recently, the development of metamaterials and metasurfaces has expanded the design space to include materials with negative refractive indices and plasmonic properties [17]. Materials with negative indices can be used with transformation optics to design electromagnetic cloaks [18–20], and several implementations have been demonstrated [21–23]. Transformation optics has the considerable advantage of imposing no limitations on the size and shape of the object. However, the designs often require complex anisotropic materials with negative refractive indices, which inevitably suffer from large absorption losses and dispersion [24]. Carpet cloaks can be designed using materials with positive isotropic indices [22,25], but they do not provide invisibility.

Scattering cancellation is an alternative method to achieve invisibility [26]. Of course, antireflection coatings are used throughout industry to reduce reflection by planar surfaces [10,11]. For nonplanar structures, such as spheres and cylinders, multilayered designs involving dielectric [27,28] and plasmonic materials [29,30] have been proposed. However, these designs are either limited to simple shapes or suffer from large absorption losses [24]. Metasurfaces with tailored surface impedance or mantle cloaks have been proposed as an alternative to thicker, multilayer designs [31–34]. However, their application is limited to subwavelength objects. Therefore, a more generalized design approach applicable to objects of any size and shape that does not involve plasmonics or materials with negative indices would be a significant advance.

In recent years, inverse methods have become the predominant approach in photonic design [35–43]. Broadly described, inverse design is the optimization of an objective function with respect to the structure or material properties. Inverse design has been used to optimize photonic circuits [35–37] and nanophotonic resonant structures [38–42]. The key to developing an effective inverse design algorithm lies in the definition of the objective function. For instance, to design highly scattering structures, a suitable objective function is the determinant of the wave operator defining the poles of the scattering matrix [40]. On the other hand, nonscattering structures can be designed by minimizing the scattering cross section [44,45]. However, the need to calculate the cross section from the scattered field for all angles at each step in the optimization poses a computational challenge, and also limits the scale and type of objects that can be optimized.

In this paper, we develop a more efficient optimization method to design unidirectional invisibility cloaks. Our inverse design approach is based on minimizing the forward scattering amplitude, which by the optical theorem is equal to the total cross section [46,47]. The use of the optical theorem circumvents the need to integrate the scattered fields over an angle by evaluating the cross section from the forward scattering amplitude alone, thus providing a simpler, more computationally efficient algorithm. Our method is completely general and makes no assumptions about the size, shape, or composition of the object, though large objects still present a computational challenge. For demonstration, we apply it to design unidirectional gradient-permittivity cloaks to minimize scattering from metallic targets of different sizes and shapes. We show that orders of magnitude reductions in the cross section are achievable with subwavelength coatings. The ease and effectiveness of our approach enables the optimization of large-scale nonscattering structures and invisibility cloaks.

Consider an object in free space with isotropic relative permittivity ${ϵ}_r(\mathbf{r})$. The electric field scattered by the object is obtained by solving the time-harmonic Maxwell’s equation with the Silver–Müller radiation boundary condition. The solution for the scattered field is [2,48–50]$${\mathbf{E}}_s(\mathbf{r})=\int \mathbf{G}(\mathbf{r}-{\mathbf{r}}^{\prime })F({\mathbf{r}}^{\prime })\mathbf{E}({\mathbf{r}}^{\prime }){\mathrm{d}}^3{r}^{\prime },$$(1)where $F(\mathbf{r})=k^2[{ϵ}_r(\mathbf{r})-14]/(4\pi )$ and $\mathbf{G}(\mathbf{r}-{\mathbf{r}}^{\prime })$ is the dyadic Green’s function given by $$\mathbf{G}(\mathbf{r}-{\mathbf{r}}^{\prime })=(\mathbf{I}+\frac{1}{k^2}\nabla \nabla )\frac{e^{ik|\mathbf{r}-\mathbf{r}\prime |}}{|\mathbf{r}-{\mathbf{r}}^{\prime }|}.$$(2)The total field $\mathbf{E}(\mathbf{r})$ is given by the sum of the scattered field and the incident field ${\mathbf{E}}_i(\mathbf{r})$, which we assume is a plane wave of the form ${\widehat{\mathbf{e}}}_i{e}^{i{\mathbf{k}}_i·\mathbf{r}}$. In the far-field approximation ($|\mathbf{r}|\gg |{\mathbf{r}}^{\prime }|$), the dyadic Green’s function reduces to [49,50]$$\mathbf{G}(\mathbf{r}-{\mathbf{r}}^{\prime })\simeq (\mathbf{I}-\widehat{\mathbf{k}}\widehat{\mathbf{k}})\frac{e^{\mathit{ikr}}}{r}{e}^{-i\mathbf{k}·{\mathbf{r}}^{\prime }},$$(3)where $r=|\mathbf{r}|$ and $\mathbf{k}=k\mathbf{r}/r$ is the scattered wavevector. The scattered field then takes the form $${\mathbf{E}}_s(\mathbf{r})\simeq \frac{e^{\mathit{ikr}}}{r}\mathbf{f}(\mathbf{k},{\mathbf{k}}_i){\widehat{\mathbf{e}}}_i,$$(4)where $\mathbf{f}(\mathbf{k},{\mathbf{k}}_i)$ is the scattering amplitude dyadic [49,50]. Normally, the scattering cross section is evaluated by integrating the square modulus of the scattering amplitude over a solid angle. This expression can be used as the objective function to design nonscattering structures [44,45]. However, the need to evaluate and integrate the scattering amplitude over an angle leads to a significant computational expense. A much simpler method is to apply the optical theorem [46,47,49], which relates the total cross section to the imaginary part of the scattering amplitude evaluated in the forward ($\mathbf{k}={\mathbf{k}}_i$) direction as $$\mathit{\sigma }=\frac{4\pi }{k}\mathrm{Im}[\mathbf{f}({\mathbf{k}}_i,{\mathbf{k}}_i)].$$(5)The optical theorem implies that the scattering amplitude evaluated in just one direction can be used as the objective for a minimization function. However, to evaluate the scattering amplitude in Eq. (5), we must express the total field in terms of the incident field. This can be accomplished by discretizing space into $N$ grid points of length $h$ as $r(j)=jh$, where $j$ is a positive integer ($j\in \mathbb{N}$). In discrete form, Eq. (1) can be written as [40,51,52]$$E_{s,n}(jh)=\sum _{i,\mathit{m}}{h}^3{G}_{mn}(jh-ih)F(ih)E_m(ih),i\in \mathbb{N},m=\mathrm{1,}\mathrm{2,}\mathrm{3,}$$(6)or in matrix-vector form $${\mathbf{E}}_s=GF\mathbf{E},$$(7)where $F$ and $G$ are both $3N\times 3N$ matrices, $F$ is a diagonal matrix, and $h^3$ has been absorbed into the definition of $G$. This matrix equation can be solved to obtain the total field as [40,51,52]$$\mathbf{E}={(I-GF)}^{-1}{\mathbf{E}}_i.$$(8)Substituting this expression into the discrete form of Eq. (5), we obtain $$\mathit{\sigma }=\frac{4\pi }{k}(\mathbf{I}-\widehat{\mathbf{k}}\widehat{\mathbf{k}})\mathrm{Im}[{\mathbf{E}}_i^{*}F{(I-GF)}^{-1}{\mathbf{E}}_i].$$(9)Equation (9) forms the objective function of our inverse design optimization algorithm. For a particular incident direction, we optimize nonscattering structures by minimizing $\mathit{\sigma }$ with respect to the permittivity values using a nonlinear least squares algorithm. Note that in the continuum limit $h\to 0$, Eq. (9) is an exact expression for the cross section, accounting for all multiple scattering and resonance effects.

The most computationally intensive calculation in Eq. (9) is the matrix inversion, and the desire is to minimize the number of matrix inversions required to minimize the scattering cross section. Unfortunately, integral equations give rise to dense matrices, so sparse matrix techniques such as a conjugate gradient and its derivatives cannot be directly applied without introducing sparsifying approximations [51,52]. Proceeding with a conventional gradient descent, we must evaluate the derivative of the objective with respect to $F(\mathbf{r})$, which involves a numerical evaluation of the Jacobian and an additional matrix inversion. Fortunately, since the objective function has an analytical form, we can obtain a closed-form expression for the Jacobian as $$\frac{\delta }{\delta F}F{(I-GF)}^{-1}{\mathbf{E}}_i=[I+F{(I-GF)}^{-1}G]{(I-GF)}^{-1}{\mathbf{E}}_i.$$(10)This expression provides a rapid determination of the Jacobian matrix at each step in the gradient descent and, more importantly, only requires the same matrix inversion already performed to calculate the cross section. This removes extraneous matrix inversions at each step in the gradient descent and greatly speeds up the nonlinear least squares algorithm.

We tested this methodology using MATLAB on a standard PC, and the memory limitations provided a practical limitation on the size of the matrix that could be inverted, thus limiting the size of the scattering object compared to the wavelength. To relax this size limitation, we tested the methodology using a 2D object and transverse electric field. In this case, Eq. (1) reduces to a scalar equation for the transverse component, where the integral is over two dimensions and Green’s function is the scalar 2D Helmholtz kernel [3,51]$$G_{2D}(\mathbf{r}-{\mathbf{r}}^{\prime })=\frac{i}{4}{H}_0^{(1)}(k|\mathbf{r}-{\mathbf{r}}^{\prime }|),$$(11)where $H_0^{(1)}(x)$ is a Hankel function of the first kind. While this form of Green’s function is used in Eq. (8), the far-field form replacing Eq. (3) is $$G_{2D}(\mathbf{r}-{\mathbf{r}}^{\prime })\simeq \sqrt{\frac{i}{8\pi kr}}{e}^{ikr}{e}^{-i\mathbf{k}·{\mathbf{r}}^{\prime }}.$$(12)Rederiving the optical theorem in 2D results in an expression for the cross section similar to Eq. (9), absent the dyadic prefactor and with volume elements replaced by area elements [46,47].

To test this methodology to design unidirectional nonscattering objects, we considered various cases of an optimized cloak applied to a metal rod. The model was parameterized with the dimensions of the metal rod, cloak, and grid elements all defined relative to a wavelength $\lambda$. The permittivity of the metal rod was chosen to be $1+i{10}^{10}$, corresponding to a highly conductive metal. Although we quickly found that neither the resulting cloak design nor the scattering cross-section reduction was highly dependent on the exact permittivity of the rod as long as it was sufficiently conductive, we did find that using realistic metal permittivity values at optical wavelengths resulted in significant absorption cross section values that put a lower bound on the total cross-section reduction. Using a highly conductive rod minimizes absorption in the rod so the scattering cross-section reduction is equivalent to the total cross-section reduction as computed using the optical theorem. Figure 1 shows the model setup, which in this case included the highly conductive metal cylinder with a $0.5\lambda$ radius, a $0.5\lambda$ thick cloak, and a grid size of $0.0417\lambda$ ($\lambda /24$), resulting in 1804 grid elements across the rod and cloak. The plane wave was incident in the $+\widehat{\mathbf{x}}$ direction with a polarization parallel to the length of the rod. The top row of plots [Fig. 1(a)] shows the permittivity distribution with no cloak as well as the real part of the scattered and total field inside and outside the rod calculated using the discrete integral form of Eqs. (7) and (8). The total cross section of the uncloaked rod calculated using Eq. (9) is $2.4\lambda$. As expected, the rod strongly scatters the incident wave, creating interference patterns throughout the region around the rod.

The second and third rows [Figs. 1(b) and 1(c)] show the results of optimizing the cloak using the fmincon function in Matlab. The fmincon function was used to vary the permittivity at each grid element in the cloak to minimize the cross section as calculated using Eq. (9). The optimization was run with an allowed permittivity range of 1–16, the upper bound corresponding to germanium, and the gradient of the cross section was supplied to minimize the computational time required; otherwise, fmincon used the default algorithm and convergence criteria. Because optimization algorithms such as fmincon require an initial seed value for the optimization, we defined the initial shell permittivity to be a uniform value throughout the cloak. The results in the second row were for a uniform seed value of 3, and the optimized cloak permittivity distribution ranges from 1 to 8.7 with an expected mirror symmetry around the $x$ axis that was not imposed in the optimization algorithm. This design shows extremely small, scattered field intensities around the cloak, resulting in a total field nearly identical to the incident plane wave. Calculating the scattering cross section with the cloak using Eq. (9) gives a value of $6.04\times {10}^{-4}\lambda$, a reduction factor of approximately 4,000, or 36.0 dB. This strong reduction in the scattering amplitude is shown in the final column of Fig. 1(b). These results confirm the effectiveness of this methodology and were generated in only a few minutes of computational time. The results in the third row [Fig. 1(c)] used the same algorithm but with an initial seed permittivity of 5. With this higher seed value, the maximum permittivity in the cloak increases to 13.5 and the permittivity has sharper gradients within the cloak. Despite the very different design, the cross section reduction is just as effective as the previous design, with a total cross section of $4.86\times {10}^{-4}\lambda$, a reduction factor of 37.0 dB. Comparing these results to those in Fig. 1(b) shows the importance of the initial seed in the optimization and points to the need for global optimization algorithms in future work.

Having generated these promising results for a relatively thick cloak, we then performed similar optimizations with half the cloak thickness, or $0.25\lambda$. While the structure size was reduced, we maintained a grid size of $0.0417\lambda$ to ensure an accurate comparison with the results in Fig. 1. In this case, shown in Fig. 2, we observed a much greater difference in the cross-section reduction between the two optimization seeds. The cross-section reduction for the ${ϵ}_{\text{shell}}=3$ seed was only 11.5 dB, whereas the reduction with the ${ϵ}_{\text{shell}}=5$ seed was 36.6 dB, which is comparable to the results with the thicker cloak.

Next, we performed an additional cloak design for a metal rod with an elliptical cross section to show the applicability to general shapes aside from an ideal cylinder. In this case, the metal rod was an ellipse with a $0.8\lambda$ major axis length and a $0.4\lambda$ minor axis length, with the major axis rotated 45° off the incident field. The cloak was $0.5\lambda$ thick around the elliptical rod, and the grid size remained $0.0417\lambda$. The optimization results in this case, shown in Figs. 3(b) and 3(c), gave cross-section reduction factors of 35.8 dB for the ${ϵ}_{\text{shell}}=3$ seed and 28.0 dB for the ${ϵ}_{\text{shell}}=5$ seeds.

Finally, we evaluated the sensitivity of an optimized cloak to the wavelength and the angle of incidence. Since the optimization was performed for a single angle of incidence and wavelength, the design is expected to be sensitive to variations in these parameters. Figures 4(a) and 4(b), respectively, show the wavelength and angle of incidence dependence of the uncloaked cylinder in Fig. 1(a) and the cloaked design in Fig. 1(b). As expected, the cross section of the optimized cloak is a minimum at the design wavelength (${\lambda }_0$) and an angle of incidence of 0°. The scattering cross section of the cloaked structure is reduced over a 12% spectral bandwidth and is an order of magnitude lower over a 3.5% spectral bandwidth. With respect to the angle of incidence, the cross section of the cloak design is reduced over a 54° angular range and reduced by an order of magnitude over 14°.

Overall, these cases all show that the optical theorem methodology to design unidirectional invisibility cloaks is highly effective. In all cases examined here, the optimization of the cloak resulted in cross sections that were a small fraction of the initial metal rod, and the optimization was completed in 5 min or less using a standard PC. Areas for future work include using a global optimizer to get the best design possible for a given set of permittivity constraints. This should be possible using the global optimization functions in MATLAB. We also plan to explore omnidirectional and broadband cloak designs using a multiobjective optimization to simultaneously minimize the cross section at several wavelengths and angles.

We developed and applied an inverse optimization method to design unidirectional invisibility cloaks. The method is based on minimizing an objective function equal to the forward scattering amplitude of the cloaked object, which by the optical theorem is equivalent to the total cross section. The use of the optical theorem greatly simplifies the optimization since only the forward scattering amplitude must be calculated at each iteration, in contrast to the traditional approach based on integrating the amplitude over an angle. Using a nonlinear least squares gradient descent, we applied the method to design several unidirectional gradient-permittivity cloaks to reduce scattering from metallic circular and elliptical cylinders. We believe our generalized approach provides a simple, effective design tool that enables the optimization of large-scale nonscattering structures and invisibility cloaks.