Limitations on bandwidth-integrated passive cloaking

Abstract

We present a general framework for the computation of structure-agnostic bounds on the performance of passive cloaks over a nonzero bandwidth. We apply this framework in 2D to the canonical scenario of cloaking a circular object. We find that perfect cloaking using a finite-sized isotropic cloak is impossible over any bandwidth, with the bounds scaling linearly with the bandwidth before saturating due to the finite size of the cloak and the presence of material loss. The bounds also exhibit linear scaling with material loss in the cloak and linear scaling with the inverse of the radial thickness of the design region before saturation due to finite-size effects or the presence of material loss. The formulation could readily find applications in the development of cloaking devices, setting expectations and benchmarks for optimal performance.

Introduction.— The ability to engineer optical cloaking of objects to make them undetectable by observers has wide-ranging implications for radar [1], obscurance [2, 3], and sensing [4] applications. Several techniques have been proposed to achieve (near-)invisibility at a single operating optical wavelength, including transformation optics [5], scattering-cancellation via plasmonic [6] or mantle cloaking [7], waveguide cloaking [8], transmission-line networks [9], or the use of anomalous localized resonances [10], to name a few. A fundamental and practically relevant question is whether one can achieve (near-)invisibility over a wide spectral window, e.g., the entire visible spectrum. As pointed out by Pendry et al. using an argument based on the phase and group velocities of light in the presence of dispersion [5], and more recently by others using more formal arguments for the inevitable distortion of a pulse wave by materials with local response [11], perfect cloaking of an isolated object in vacuum over a nonzero bandwidth is impossible due to causality; the finite speed of light and the presence of dispersion limits fully effective cloaking to a single frequency. However, the Pendry et al. and Miller proofs say little about imperfect cloaking (a nonzero but small scattering cross section). For instance, engineering a cloak to attain near-perfect invisibility at a single frequency does not preclude “good” cloaking over some nonzero bandwidth around that frequency. Hashemi et al. [12] presented an alternative proof that agrees with the results of Pendry et al. and Miller for the impossibility of perfect cloaking over a nonzero bandwidth for physical, causal materials, but additionally derived scaling relations with implications for design, e.g., the allowed bandwidth over which one may expect practical cloaking scales inversely with the diameter of the cloaked object.

In this article, we adapt recently developed constrained optimization approaches [13, 14, 15] to study the following question related to bounds on cloaking performance: Given a specified cloak material along with an operating frequency and bandwidth, what is the best possible performance of a passive cloaking device that can be formed out of that material within a prescribed region surrounding a cloaked object? In addition to capturing scaling behaviors, the presented formalism allows for quantitative assessments of expected performance of cloaks. Typical realizations of cloaks require anisotropic metamaterials with extreme and exotic properties, such as permittivities and permeabilities less than one and even equal to zero in some regions of the cloak [5, 16, 17]. This raises the natural question of what is the best achievable performance using the simplest materials possible. While the numerical bounds lack the transparency and intuition of analytic expressions, they can still provide insight into scaling properties and general trends, and even lead to fairly tight limits (within an order of magnitude of structures discovered via topology optimization).

Formulation.— In line with the majority of previous works on cloaking systems [12, 18, 19, 15], and as a proof of concept, we study the common scenario of a cloaked object illuminated by an incident plane wave emitted by a far-away source. We consider extinction, a nonnegative quantity for passive systems [20], as the fundamental figure of merit, for the following reasons. First, it is better suited for describing the practical goals of cloaks as opposed to scattered or absorbed power alone: minimizing extinction minimizes the scattered field over both near- (absorbed power) and far-field (scattered power) domains [15], implying that an object becomes truly undetectable regardless of the position of an observer as the figure of merit approaches zero. Second, analyticity in the upper-half of the complex frequency plane means that one can exploit the residue theorem to evaluate the spectral average over a nonzero bandwidth via a single scattering calculation [21, 12], which greatly simplifies computations.

Throughout, we use the language of scattering theory laid out in previous descriptions of electromagnetic bounds [14, 22, 23], wherein the $\mathbb{T}$ operator describing bound polarization currents in the medium is defined by the relation

$$\mathbb{I}=\mathbb{T}\left(\mathbb{V}^{-1}-\mathbb{G}_{\mathrm{bg}}\right)=\left(\mathbb{V}^{-1}-\mathbb{G}_{\mathrm{bg}}\right)\mathbb{T}.$$ (1)

The operator $\mathbb{G}_{\mathrm{bg}}$ represents the background Green’s function, which for vacuum satisfies $\left[\nabla\times\nabla\times-\frac{\omega_{0}^{2}}{c^{2}}\right]\mathbb{G}_{\mathrm{vac}}(\mathbf{r},\mathbf{r}^{\prime};\omega_{0})=\frac{\omega_{0}^{2}}{c^{2}}\delta^{(3)}(\mathbf{r}-\mathbf{r}^{\prime})$. (Note that in keeping with prior work, this definition includes an additional factor of $\frac{\omega_{0}^{2}}{c^{2}}$ compared to the standard convention [24]). The $\mathbb{V}$ operator is the scattering potential (susceptibility) relative to this background medium (whatever additional material response was not included in the definition of $\mathbb{G}_{\mathrm{bg}}$), and $\left|\mathbf{E}_{i}\right>$ and $\left|\mathbf{J}_{i}\right>$ are defined as the incident electric fields and electric currents in the background, respectively. When acting on incident fields, the $\mathbb{T}$ operator produces the generated current $|\mathbf{J}_{g}\rangle=-\frac{ik_{0}}{Z}\mathbb{T}|\mathbf{E}_{i}\rangle$, with $Z=\sqrt{\mu_{0}/\epsilon_{0}}$ denoting the impedance of free space.

Employing the $\mathbb{T}$ operator relation to a scattering problem involving a prescribed incident field [14], the extinguished power may be written as,

$$P_{\mathrm{ext}}(\omega_{0})=\frac{1}{2}\operatorname{Re}\left[\left<\mathbf{E}_{i}|\mathbf{J}_{g}\right>\right]=\frac{1}{2Z}\operatorname{Im}\left[k_{0}\left<\mathbf{E}_{i}|\mathbb{T}|\mathbf{E}_{i}\right>\right],$$ (2)

where $k_{0}$ is left inside Im because it will be continued to the complex plane to map this quantity to a bandwidth-averaged quantity [21, 12]. Namely, since real sources emit light over a nonzero bandwidth a key figure of merit is the bandwidth average of extinct power. Considering a Lorentzian window function $L(\omega)\equiv\frac{\Delta\omega/\pi}{(\omega-\omega_{0})^{2}+\Delta\omega^{2}}$ centered at $\omega_{0}$ with a bandwidth $\Delta\omega$, the average extinction power for incident plane waves,

$$\langle P_{\mathrm{ext}}\rangle\equiv\int_{-\infty}^{\infty}P_{\mathrm{ext}}(\omega)L(\omega)\differential\omega\to P_{\mathrm{ext}}(\tilde{\omega}),$$ (3)

may be evaluated by closing the contour and picking up the residue at the pole of the Lorentzian in the upper-half complex frequency plane, yielding Eq. (2) with all vectors and operators evaluated at a complex frequency $\tilde{\omega}\equiv\omega_{0}+i\Delta\omega$ and complex wave number $\tilde{k}\equiv(\omega_{0}+i\Delta\omega)/c$, simplifying the calculation so that one no longer needs to consider power at each individual frequency within the Lorentzian window. Note also that working at a complex frequency requires the evaluation of $\chi(\tilde{\omega})$, which is mathematically equivalent to using modified materials at a real frequency $\omega_{0}$ [12].

To isolate the impact of the cloak on scattering from a fixed object, we let $\mathbb{G}_{\mathrm{bg}}$ refer to the vacuum Green’s function and shift the known scattering properties of the fixed object into the background of the scattering operator describing the cloak, with $\mathbb{G}_{\mathrm{obj}}$ denoting the Green’s function of the fixed object in isolation, so that $\mathbb{G}_{\mathrm{obj}}$ satisfies $\left[\nabla\times\nabla\times-\frac{\omega_{0}^{2}}{c^{2}}(\mathbb{V}_{\mathrm{obj}}+\mathbb{I})\right]\mathbb{G}_{\mathrm{obj}}(\mathbf{r},\mathbf{r}^{\prime};\omega_{0})=\frac{\omega_{0}^{2}}{c^{2}}\mathbb{I}\delta^{(3)}(\mathbf{r}-\mathbf{r}^{\prime})$. In particular, let $\mathbb{T}_{\mathrm{obj}}=(\mathbb{V}_{\mathrm{obj}}^{-1}-\mathbb{G}_{\mathrm{vac}})^{-1}$ be the scattering operator of the fixed object (the object to be cloaked) in isolation, and let $\mathbb{T}_{\mathrm{d}}=(\mathbb{V}_{\mathrm{clk}}^{-1}-\mathbb{G}_{\mathrm{obj}})^{-1}$ be the scattering operator of the designable object (the cloak) dressed by the fixed object. After some straightforward algebraic manipulations (see Appendix), we find

$$P_{\mathrm{ext}}(\omega_{0})=\frac{1}{2Z}\operatorname{Im}\left[k_{0}\left<\mathbf{E}_{i}\right|\mathbb{G}_{\mathrm{vac}}^{-1}\mathbb{G}_{\mathrm{obj}}\mathbb{V}_{\mathrm{obj}}\left|\mathbf{E}_{i}\right>\right]\\ +\frac{1}{2Z}\operatorname{Im}\left[k_{0}\left<\mathbf{E}_{i}\right|\mathbb{G}_{\mathrm{vac}}^{-1}\mathbb{G}_{\mathrm{obj}}\mathbb{T}_{\mathrm{d}}(\mathbb{G}_{\mathrm{vac}}\mathbb{T}_{\mathrm{obj}}+\mathbb{I})\left|\mathbf{E}_{i}\right>\right].$$ (4)

The first term is the extinguished power when no cloak is present while only the second term depends on the cloaked object and quantifies the interaction between the cloak and cloaked object.

Our derivation of bounds exploits the optimization procedure based on Lagrange duality laid out in Refs. [25, 13]. The loosest such bound only imposes that the optimal scattering operator satisfies the conservation of power (optical theorem [20]) over the entire design domain, and not the full scattering equations. Defining $|\mathbf{E}_{\mathrm{d}}\rangle\equiv(\mathbb{G}_{\mathrm{vac}}\mathbb{T}_{\mathrm{obj}}+\mathbb{I})|\mathbf{E}_{i}\rangle=\frac{iZ}{k_{0}}\mathbb{G}_{\mathrm{obj}}\left|\mathbf{J}_{i}\right>$, the total field in the presence of only the fixed object, and $|\mathbf{T}_{\mathrm{d}}\rangle\equiv\mathbb{T}_{\mathrm{d}}|\mathbf{E}_{\mathrm{d}}\rangle$ so that $-\frac{ik_{0}}{Z}|\mathbf{T}_{\mathrm{d}}\rangle$ is the induced current in the cloak, we relax the problem such that, instead of optimizing over $\mathbb{T}_{\mathrm{d}}$ with support in the design region, we optimize over $|\mathbf{T}_{\mathrm{d}}\rangle$ with support in the design region, i.e., optimize over all possible polarization currents within the design region so that the bound considered here automatically takes into account all possible distributions of vacuum and the prescribed cloak material within the design region. It should be noted, however, that this bound formalism does not yield the optimal material distribution of the cloak [13], a task that calls for more computationally demanding topology optimization methods (NP-hard problems that preclude guarantees of globally optimal solutions) [26]. Similar techniques have recently been used to derive bounds on deterministic scattering and fluctuational electrodynamic phenomena [14, 27, 23, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37]. Concretely, we wish to solve the following problem:

$$\underset{\left|\textbf{T}_{\mathrm{d}}\right>}{\text{min}}\,\,\frac{1}{2Z} \operatorname{Im}\left[k_{0}\left<\mathbf{E}_{i}\right|\mathbb{G}_{\mathrm{vac}}^{-1}\mathbb{G}_{\mathrm{obj}}\left(\mathbb{V}_{\mathrm{obj}}\left|\mathbf{E}_{i}\right>+\mathbb{P}_{\mathrm{des}}\left|\textbf{T}_{\mathrm{d}}\right>\right)\right]$$ (5a)

$$\text{s.t.}\,\forall m,\quad \langle\mathbf{E}_{\mathrm{d}}|\mathbb{P}_{m}|\mathbf{T}_{\mathrm{d}}\rangle-\langle\mathbf{T}_{\mathrm{d}}|\mathbb{U}_{m}|\mathbf{T}_{\mathrm{d}}\rangle=0,$$ (5b)

where $\mathbb{P}_{\mathrm{des}}$ is a projection operator into the design region $\Omega_{\mathrm{des}}$, $\mathbb{P}_{m}$ is a projection operator into a subregion $\Omega_{m}\subseteq\Omega_{\mathrm{des}}$ contained within the design region, and $\mathbb{U}_{m}\equiv\mathbb{P}_{\mathrm{des}}(\chi_{\mathrm{clk}}^{-1\dagger}\mathbb{P}_{m}-\mathbb{G}_{\mathrm{obj}}^{\dagger}\mathbb{P}_{m})\mathbb{P}_{\mathrm{des}}$. Taking the real and imaginary parts of the complex-valued constraints enforce the conservation of resistive and reactive power, respectively, spatially integrated over a region $\Omega_{m}\subseteq\Omega_{\mathrm{des}}$ [13]. This is a quadratically constrained quadratic program over the field $|\mathbf{T}_{\mathrm{d}}\rangle$ and we numerically compute bounds on Eq. (5a), the primal objective, by evaluating and optimizing the corresponding concave Lagrange dual function [25].

Figure 1: (a) Bounds on bandwidth-integrated extinction power as a function of the bandwidth. All curves and markers refer to performance bounds and inverse designs, respectively, expressed as a ratio of the extinguished power from incident plane waves in the presence and absence of a cloak for $R_{\mathrm{obj}}=\lambda_{0}/4$ and $R_{\mathrm{des}}=3\lambda_{0}/4$ (bottom-right schematic). Insets: Representative inverse designs corresponding to $\Delta\omega/\omega_{0}=10^{-7}$. (b) Performance of structures discovered via topology optimization for $\Delta\omega/\omega_{0}=10^{-7}$, $R_{\mathrm{obj}}=\lambda_{0}/4$, $R_{\mathrm{des}}=3\lambda_{0}/4$, $\chi_{\mathrm{clk}}(\tilde{\omega})=5+10^{-3}i$, and $\chi_{\mathrm{obj}}(\tilde{\omega})=-4+i$ (solid) or $\chi_{\mathrm{obj}}(\tilde{\omega})=4+10^{-3}i$ (dashed). Top: The cloak suppresses extinction within a frequency range around $\omega_{0}$ at the expense of worse performance in other parts of the spectrum. Bottom: Real part of the total out-of-plane electric field at $\omega_{0}$. Green circles outline the region of space where the object and cloak are contained, with an inner and outer radius $R_{\mathrm{obj}}$ and $R_{\mathrm{des}}$, respectively.

Applications.— We now exploit the framework above to numerically obtain performance bounds for passive cloaks composed of an isotropic electric susceptibility. For simplicity, we focus on 2D settings and consider the scenario of a fixed circular object of radius $R_{\mathrm{obj}}$, a common benchmark example for cloaking [38], with the cloak restricted to an annular design region with an inner radius $R_{\mathrm{obj}}$ and outer radius $R_{\mathrm{des}}$ [see Fig. 1(a) bottom-right schematic] for the case of incident TM (electric field out of plane) plane waves.

First, we consider the situation of illumination of broadband light from the far-field (plane waves). For a fixed $\chi(\omega_{0}+i\Delta\omega)$ we find that the bounds scale linearly with the bandwidth before saturating due to the finite size of the design domain and due to the presence of material loss, see Fig. 1(a). As seen in the figure, enforcing the conservation of power within the design domain leads to nontrivial bounds with interesting trends and scaling behavior that are also fairly tight (within an order of magnitude of the objective values of discovered structures). The structures found via topology optimization perform well within a bandwidth $\Delta\omega$ around the “center” frequency $\omega_{0}$, but at the expense of more scattering in nearby regions of the spectrum, see Fig. 1(b). This is in agreement with the findings in Ref. [39], which showed that any passive cloak made of a linear, nondiamagnetic material respecting causality always increases the scattering and extinction integrated over all frequencies compared to the original uncloaked object. Thus, near-invisibility in a given frequency window necessarily implies significant scattered/extinct power at other frequencies, and therefore such cloaks may be detected more easily than the original uncloaked objects when illuminated by sufficiently broadband light (short pulses). As seen in Fig. 1(b), the electric field in the lossy metallic object is expelled from and rerouted around the object with the introduction of the cloak. The dielectric object with a smaller loss, on the other hand, may try to create zero-field regions or it may allow the electric field to penetrate the object if it is advantageous as is the case for the dielectric cloak example shown in Fig. 1(b).

Next, we investigate the performance limits of 2D cloaks of isotropic susceptibility as a function of the material loss, see Fig. 2. The bounds (solid lines) follow trends seen in topology-optimized designs (circular markers) over a broad range of material loss values and both demonstrate, unsurprisingly, improved performance as the material loss in the cloak decreases, exhibiting linear scaling with decreasing material loss before eventually saturating. This scaling can be understood as follows. Writing the total $\mathbb{T}$ operator in $2\times 2$ block-form over the two regions, namely,

$$\mathbb{G}^{\mathrm{vac}}=\begin{bmatrix}\mathbb{G}^{\mathrm{vac}}_{\mathrm{obj},\mathrm{obj}}&\mathbb{G}^{\mathrm{vac}}_{\mathrm{obj},\mathrm{clk}}\\ \mathbb{G}^{\mathrm{vac}}_{\mathrm{clk},\mathrm{obj}}&\mathbb{G}^{\mathrm{vac}}_{\mathrm{clk},\mathrm{clk}}\end{bmatrix},$$ (6)

$$\mathbb{V}^{-1}=\begin{bmatrix}\mathbb{V}^{-1}_{\mathrm{obj}}&0\\ 0&\mathbb{V}^{-1}_{\mathrm{clk}}\end{bmatrix},$$ (7)

$$\mathbb{T}^{-1}=\begin{bmatrix}\mathbb{T}_{\mathrm{obj}}^{-1}&-\mathbb{G}^{\mathrm{vac}}_{\mathrm{obj},\mathrm{clk}}\\ -\mathbb{G}^{\mathrm{vac}}_{\mathrm{clk},\mathrm{obj}}&\mathbb{T}_{\mathrm{clk}}^{-1}\end{bmatrix},$$ (8)

leads, after the use of the Woodbury formula for a matrix inverse, to a total $\mathbb{T}$ operator expression of the form

$$\mathbb{T} \equiv\begin{bmatrix}\mathbb{Y}_{\mathrm{obj}}&\mathbb{T}_{\mathrm{obj}}\mathbb{G}^{\mathrm{vac}}_{\mathrm{obj},\mathrm{clk}}\mathbb{Y}_{\mathrm{clk}}\\ \mathbb{T}_{\mathrm{clk}}\mathbb{G}^{\mathrm{vac}}_{\mathrm{clk},\mathrm{obj}}\mathbb{Y}_{\mathrm{obj}}&\mathbb{Y}_{\mathrm{clk}}\end{bmatrix}.$$ (9)

Here, we introduced $\mathbb{Y}_{\mathrm{clk}}$ which is the scattering operator of the designable cloak dressed by the presence of the fixed object, as seen from $\mathbb{Y}_{\mathrm{clk}}\equiv(\mathbb{T}_{\mathrm{clk}}^{-1}-\mathbb{G}^{\mathrm{vac}}_{\mathrm{clk},\mathrm{obj}}\mathbb{T}_{\mathrm{obj}}\mathbb{G}^{\mathrm{vac}}_{\mathrm{obj},\mathrm{clk}})^{-1}=(\mathbb{V}_{\mathrm{clk}}^{-1}-(\mathbb{G}^{\mathrm{vac}}_{\mathrm{clk},\mathrm{clk}}+\mathbb{G}^{\mathrm{vac}}_{\mathrm{clk},\mathrm{obj}}\mathbb{T}_{\mathrm{obj}}\mathbb{G}^{\mathrm{vac}}_{\mathrm{obj},\mathrm{clk}}))^{-1}$ and recognizing $\mathbb{G}_{\mathrm{vac}}+\mathbb{G}_{\mathrm{vac}}\mathbb{T}_{\mathrm{obj}}\mathbb{G}_{\mathrm{vac}}=\mathbb{G}_{\mathrm{obj}}$ as the Green’s function of the fixed object in isolation in vacuum. Likewise for $\mathbb{Y}_{\mathrm{obj}}$ but with the fixed and design regions switching roles. Letting $\chi_{\mathrm{clk}}=\chi_{\mathrm{clk}}^{\prime}+i\chi_{\mathrm{clk}}^{\prime\prime}$ be the separation of the susceptibility into its real and imaginary parts, we find

$$\mathbb{T}_{\mathrm{clk}}=(\chi_{\mathrm{clk}}^{\prime}\mathbb{P}_{\mathrm{clk}})[(\mathbb{I}-\mathbb{G}^{\mathrm{vac}}_{\mathrm{clk},\mathrm{clk}}\mathbb{V}_{\mathrm{clk}}^{\prime})^{-1}\\ +\sum_{n=0}^{\infty}n(\mathbb{G}^{\mathrm{vac}}_{\mathrm{clk},\mathrm{clk}})^{n}(\chi_{\mathrm{clk}}^{\prime})^{n-1}(i\chi_{\mathrm{clk}}^{\prime\prime})]\\ +(i\chi_{\mathrm{clk}}^{\prime\prime}\mathbb{P}_{\mathrm{clk}})(\mathbb{I}-\mathbb{G}^{\mathrm{vac}}_{\mathrm{clk},\mathrm{clk}}\mathbb{V}_{\mathrm{clk}}^{\prime})^{-1}+O[(\chi_{\mathrm{clk}}^{\prime\prime})^{2}]$$ (10)

and similar for $\mathbb{Y}_{\mathrm{clk}}$ but with $\mathbb{G}_{\mathrm{vac}}$ replaced with $\mathbb{G}_{\mathrm{obj}}$. Using these expansions for $\mathbb{T}_{\mathrm{clk}}$ and $\mathbb{Y}_{\mathrm{clk}}$ in the $2\times 2$ block-form for the total scattering operator, Eq. (9), leads to the prediction that the extinction at a single frequency, $P_{\mathrm{ext}}(\omega)=\frac{k_{0}}{2Z}\left<\mathbf{E}_{i}|\operatorname{Im}\left[\mathbb{T}\right]|\mathbf{E}_{i}\right>$, scales linearly with $\chi_{\mathrm{clk}}^{\prime\prime}$ towards a constant for a fixed geometric structure as $\chi_{\mathrm{clk}}^{\prime\prime}\to 0$. The numerical Lagrange dual bounds follow the same scaling behavior as the primal objective function for a fixed structure, demonstrating that faster scaling is not possible. The approach to a positive constant rather than 0 as the material loss of the cloak vanishes is in agreement with the findings of Hashemi et al. [40], which showed that if the attainable refractive index contrast is bounded from below (as is the case for a cloak of isotropic susceptibility considered here) then there is a bound on the reduction of the scattering cross section for transformation-based invisibility cloaking of an isolated object; in particular, a minimum achievable refractive index contrast greater than 0 necessarily implies a positive cross section. A Pendry cloak can achieve perfect invisibility at a single frequency with a finite device footprint but it requires a vanishing permittivity and permeability at the inner surface of the cloak [5, 16]. As shown, however, the presented formalism can provide quantitative lower bounds supporting such qualitative observations.

Figure 3 plots limit values as a function of the radial size of the design region, for a given object radius. Unsurprisingly, the lower bounds are monotonically nonincreasing as the allowed cloak footprint increases since any polarization current with a given performance can be contained in an enlarged design domain. In the single frequency case, $\Delta\omega\to 0$, the bounds scale inversely with the thickness of the design region, $R_{\mathrm{des}}-R_{\mathrm{obj}}$, over a notable range of thicknesses. Hashemi et al. [40] argued that for bounded refractive indices, the cloak thickness must scale proportionally to the thickness of the object being cloaked, and the scaling of our numerical bounds demonstrates this behavior. For positive bandwidths, $\Delta\omega>0$, the bounds saturate to a positive value, with saturation occurring rather quickly (around $R_{\mathrm{des}}/R_{\mathrm{obj}}\approx 2$) for bandwidths $\Delta\omega/\omega_{0}\gtrsim 10^{-5}$. This has interesting implications for the design of cloaks since it is typically not clear how large a device needs to be to achieve near-optimal or even reasonable performance. For our chosen parameters, the bounds show that for nonzero bandwidths a cloak design region with a thickness equal to that of the cloaked object essentially saturates the possible performance; larger device footprints will not always lead to significant improvements in performance.

Lastly, the presented inverse designs and bounds were restricted to plane-wave sources coming from a single direction of incidence. To consider cloaking robust to the direction of incidence, the figure of merit may be modified to encompass multiple plane wave directions (including fully angle-integrated extinction), such that $\langle P_{\mathrm{ext}}\rangle=\sum_{a=1}^{N}\frac{1}{2Z}\operatorname{Im}\left[\tilde{k}\left<\mathbf{E}_{i}^{(a)}|\mathbb{T}|\mathbf{E}_{i}^{(a)}\right>\right]$ where each $|\mathbf{E}_{i}^{(a)}\rangle$ is a plane wave incident from a different angle. This leads to an optimization problem of the form

$$\underset{\{\left|\textbf{T}_{\mathrm{d}}^{(a)}\right>\}}{\text{min}}\,\,\sum_{a=1}^{N}\frac{1}{2Z} \operatorname{Im}\left[\tilde{k}\left<\mathbf{E}_{i}^{(a)}\right|\mathbb{G}_{\mathrm{vac}}^{-1}\mathbb{G}_{\mathrm{obj}}\mathbb{V}_{\mathrm{obj}}\left|\mathbf{E}_{i}^{(a)}\right>\right]$$

$$+\sum_{a=1}^{N}\frac{1}{2Z} \operatorname{Im}\left[\tilde{k}\left<\mathbf{E}_{i}^{(a)}\right|\mathbb{G}_{\mathrm{vac}}^{-1}\mathbb{G}_{\mathrm{obj}}\mathbb{P}_{\mathrm{des}}\left|\textbf{T}_{\mathrm{d}}^{(a)}\right>\right]$$ (11a)

$$\text{s.t.}\,\forall a,b,m,\quad \langle\mathbf{E}_{\mathrm{d}}^{(a)}|\mathbb{P}_{m}|\mathbf{T}_{\mathrm{d}}^{(b)}\rangle-\langle\mathbf{T}_{\mathrm{d}}^{(a)}|\mathbb{U}_{m}|\mathbf{T}_{\mathrm{d}}^{(b)}\rangle=0,$$ (11b)

where $|\mathbf{E}_{\mathrm{d}}^{(a)}\rangle\equiv(\mathbb{G}_{\mathrm{vac}}\mathbb{T}_{\mathrm{obj}}+\mathbb{I})|\mathbf{E}_{i}^{(a)}\rangle$ and $|\mathbf{T}_{\mathrm{d}}^{(a)}\rangle\equiv\mathbb{T}_{\mathrm{d}}|\mathbf{E}_{\mathrm{d}}^{(a)}\rangle$. Without the constraints between the $a\neq b$ terms, the problem reduces to a sum of decoupled optimization problems. Although each angle is a different scattering problem, the polarization currents induced in each scenario are generated by the same structured media and this fact is enforced by the $a\neq b$ “cross constraints”, leading to additional tightening of the bounds [41, 42]. Since the incident field is the same up to a rotation, if the object and design region are invariant under rotations then the bounds are the same for decoupled problems and any tightening should come from the cross-constraints. To consider a more complicated scenario involving more than a single incident source, Fig. 4 shows computed bounds and inverse design performance values for a rectangular object discovered starting from random initializations for the case where $N=2$, comprising horizontally and vertically propagating plane waves.

Figure 2: Bounds (solid lines) and inverse designs (circles) relating to net extinction power (cloaking performance) for the system of Fig. 1, at a single frequency $\omega_{0}=2\pi c/\lambda_{0}$ and as a function of the material loss in the cloak for $R_{\mathrm{obj}}=\lambda_{0}/4$ and $R_{\mathrm{des}}=3\lambda_{0}/4$, and for both dielectric and metallic objects.

Figure 3: Bounds on extinction power for the same system as Fig. 1 but as a function of the outer radius of the design region $R_{\mathrm{des}}/R_{\mathrm{obj}}$, for a fixed object radius $R_{\mathrm{obj}}=\lambda_{0}/2$, and for different values of $\Delta\omega/\omega_{0}=10^{-3}$ (dashed), $10^{-5}$ (dotted), and zero bandwidth (solid); for comparison, crosses, circles and stars, respectively, are performance values obtained from inverse designs.

Figure 4: Bounds on bandwidth-integrated extinction as a function of the bandwidth with either vertically or horizontally incident plane waves, or an objective function with an equally weighted average of both (bottom-right schematic). The rectangular object has lengths $L_{\mathrm{obj},x}=3\lambda_{0}/4$ and $L_{\mathrm{obj},y}=\lambda_{0}/2$ while the rectangular design region has lengths $L_{\mathrm{des},x}=L_{\mathrm{des},y}=\lambda_{0}$. The susceptibility values are $\chi_{\mathrm{obj}}(\tilde{\omega})=4+10^{-3}i$ and $\chi_{\mathrm{clk}}(\tilde{\omega})=5+10^{-3}i$.

Conclusion.— In conclusion, we presented a formalism for computing bounds on passive, broadband cloaking systems that can provide benchmarks and set expectations for future cloaking devices, particularly those developed through the use of topology optimization. The many variables available for tuning means an exhaustive study of their relationship to the bounds on performance is complex. We presented a few key examples to demonstrate the capabilities of the formalism and to extract intuition for the obvious parameters one may wish to vary in experiments. Our numerical results support the claims that broadband cloaking of electrically large objects is practically impossible using passive linear cloaks [12, 38, 18] and provide a quantitative assessment of the limits on performance. We believe this represents an important result for the science of cloaking, suggesting that new concepts and designs, including opening to the field of nonlinear and active metamaterials [38], are necessary for increasing the cloaking bandwidth and compensating the unwanted scattering required by linearity and causality in the presence of absorption losses [18]. In closing, we note that while our examples focused on 2D and isotropic permittivities, the formalism is valid in 3D settings, and may be easily extended to consider anisotropic permittivities and permeabilities.

Acknowledgments.— We acknowledge the support of a Princeton SEAS Innovation Grant and by the Cornell Center for Materials Research (MRSEC) through Award DMR-1719875. S.M. acknowledges financial support from the Canada First Research Excellence Fund via the Institut de Valorisation des Données (IVADO) collaboration. The simulations presented in this article were performed on computational resources managed and supported by Princeton Research Computing, a consortium of groups including the Princeton Institute for Computational Science and Engineering (PICSciE) and the Office of Information Technology’s High Performance Computing Center and Visualization Laboratory at Princeton University.

Appendix A Objective function derivation

In this section, we derive an expression for the total (including the fixed and designable objects) scattering operator which isolates the contribution of the designable object. A useful operator for this purpose is the $\mathbb{W}$ operator defined by

$$\mathbb{I}=\mathbb{W}\left(\mathbb{I}-\mathbb{V}\mathbb{G}_{\mathrm{vac}}\right)=\left(\mathbb{I}-\mathbb{V}\mathbb{G}_{\mathrm{vac}}\right)\mathbb{W}.$$ (12)

From Eqs. (12) and (1), it is clear that $\mathbb{I}=\mathbb{W}\mathbb{V}\mathbb{T}^{-1}$ and $\mathbb{W}=\mathbb{I}+\mathbb{T}\mathbb{G}_{\mathrm{vac}}$. The $\mathbb{T}$ operator takes an incident field and gives $\propto$ the generated current, $|\mathbf{J}_{g}\rangle=-\frac{ik_{0}}{Z}\mathbb{T}|\mathbf{E}_{i}\rangle$, while the $\mathbb{W}$ operator takes an initial current and gives the total current, $|\mathbf{J}_{t}\rangle=\mathbb{T}\mathbb{V}^{-1}|\mathbf{J}_{i}\rangle=\mathbb{W}|\mathbf{J}_{i}\rangle$. We will make use of the nesting property of the $\mathbb{W}$ operator

$$\mathbb{W}_{\mathrm{tot}}=\mathbb{W}_{\mathrm{obj}}\mathbb{W}_{\mathrm{d}},$$ (13)

where $\mathbb{W}_{\mathrm{tot}}$ is the $\mathbb{W}$ operator of the total system (including all the scatterers), $\mathbb{W}_{\mathrm{obj}}$ is the $\mathbb{W}$ operator of the fixed object only, and $\mathbb{W}_{\mathrm{d}}$ is the $\mathbb{W}$ operator for the designable body with the fixed object contained in the background, i.e., $\mathbb{G}_{\mathrm{vac}}$ replaced by $\mathbb{G}_{\mathrm{vac}}\mathbb{W}_{\mathrm{obj}}=\mathbb{G}_{\mathrm{vac}}(\mathbb{I}+\mathbb{T}_{\mathrm{obj}}\mathbb{G}_{\mathrm{vac}})=\mathbb{G}_{\mathrm{vac}}+\mathbb{G}_{\mathrm{vac}}\mathbb{T}_{\mathrm{obj}}\mathbb{G}_{\mathrm{vac}}$ (which is the total Green’s function of the fixed object in the background, $\mathbb{G}_{\mathrm{obj}}$) in the defining relation for $\mathbb{W}_{\mathrm{d}}$. The nesting property holds since $\mathbb{W}_{\mathrm{d}}^{-1}\mathbb{W}_{\mathrm{obj}}^{-1}=(\mathbb{I}-\mathbb{V}_{\mathrm{clk}}\mathbb{G}_{\mathrm{obj}})\mathbb{W}_{\mathrm{obj}}^{-1}=(\mathbb{I}-\mathbb{V}_{\mathrm{clk}}\mathbb{G}_{\mathrm{vac}}\mathbb{W}_{\mathrm{obj}})\mathbb{W}_{\mathrm{obj}}^{-1}=(\mathbb{W}_{\mathrm{obj}}^{-1}-\mathbb{V}_{\mathrm{clk}}\mathbb{G}_{\mathrm{vac}})=(\mathbb{I}-\mathbb{V}_{\mathrm{obj}}\mathbb{G}_{\mathrm{vac}}-\mathbb{V}_{\mathrm{clk}}\mathbb{G}_{\mathrm{vac}})=\mathbb{W}_{\mathrm{tot}}^{-1}$. Using $\mathbb{W}_{\mathrm{tot}}\mathbb{V}_{\mathrm{tot}}=\mathbb{T}_{\mathrm{tot}}$ and the nesting relation Eq. (13) yields

$$\mathbb{T}_{\mathrm{tot}} =\mathbb{W}_{\mathrm{tot}}\mathbb{V}_{\mathrm{tot}}=\mathbb{W}_{\mathrm{obj}}\mathbb{W}_{\mathrm{d}}\mathbb{V}_{\mathrm{tot}}$$ (14)

$$=\mathbb{W}_{\mathrm{obj}}\left(\mathbb{W}_{\mathrm{d}}\mathbb{V}_{\mathrm{obj}}+\mathbb{W}_{\mathrm{d}}\mathbb{V}_{\mathrm{clk}}\right)$$ (15)

$$=\mathbb{W}_{\mathrm{obj}}\left(\mathbb{W}_{\mathrm{d}}\mathbb{V}_{\mathrm{obj}}+\mathbb{T}_{\mathrm{d}}\right)$$ (16)

$$=\mathbb{W}_{\mathrm{obj}}\left(\mathbb{V}_{\mathrm{obj}}+\mathbb{T}_{\mathrm{d}}\left(\mathbb{G}_{\mathrm{vac}}\mathbb{T}_{\mathrm{obj}}+\mathbb{I}\right)\right)$$ (17)

$$=\mathbb{T}_{\mathrm{obj}}+\mathbb{G}_{\mathrm{vac}}^{-1}\mathbb{G}_{\mathrm{obj}}\mathbb{T}_{\mathrm{d}}\mathbb{G}_{\mathrm{obj}}\mathbb{G}_{\mathrm{vac}}^{-1},$$ (18)

where we made use of $\mathbb{W}_{\mathrm{d}}=\mathbb{I}+\mathbb{T}_{\mathrm{d}}\mathbb{G}_{\mathrm{obj}}=\mathbb{I}+\mathbb{T}_{\mathrm{d}}\mathbb{G}_{\mathrm{vac}}\mathbb{W}_{\mathrm{obj}}=\mathbb{I}+\mathbb{T}_{\mathrm{d}}\mathbb{G}_{\mathrm{vac}}\mathbb{T}_{\mathrm{obj}}\mathbb{V}_{\mathrm{obj}}^{-1}$. Define $\left|\textbf{E}_{\mathrm{d}}\right>\equiv\left(\mathbb{G}_{\mathrm{vac}}\mathbb{T}_{\mathrm{obj}}+\mathbb{I}\right)\left|\mathbf{E}_{i}\right>=\frac{iZ}{k_{0}}\mathbb{G}_{\mathrm{obj}}\left|\mathbf{J}_{i}\right>$ and $\left|\textbf{T}_{\mathrm{d}}\right>\equiv\mathbb{T}_{\mathrm{d}}\left|\textbf{E}_{\mathrm{d}}\right>$. Use $\mathbb{T}_{\mathrm{tot}}=\mathbb{W}_{\mathrm{obj}}(\mathbb{V}_{\mathrm{obj}}+\mathbb{T}_{\mathrm{d}}(\mathbb{G}_{\mathrm{vac}}\mathbb{T}_{\mathrm{obj}}+\mathbb{I}))$ to find that the total extinct power by all the scatterers, relative to a vacuum background, is given by

$$P_{\mathrm{ext}} =\frac{1}{2Z}\operatorname{Im}\left[k_{0}\left<\mathbf{E}_{i}\right|\mathbb{T}_{\mathrm{tot}}\left|\mathbf{E}_{i}\right>\right]$$ (19)

$$=\frac{1}{2Z}\operatorname{Im}\left[k_{0}\left<\mathbf{E}_{i}\right|\mathbb{W}_{\mathrm{obj}}\mathbb{V}_{\mathrm{obj}}\left|\mathbf{E}_{i}\right>\right]$$

$$~{}~{}+\frac{1}{2Z}\operatorname{Im}\left[k_{0}\left<\mathbf{E}_{i}\right|\mathbb{W}_{\mathrm{obj}}\mathbb{T}_{\mathrm{d}}(\mathbb{G}_{\mathrm{vac}}\mathbb{T}_{\mathrm{obj}}+\mathbb{I})\left|\mathbf{E}_{i}\right>\right]$$ (20)

$$=\frac{1}{2Z}\operatorname{Im}\left[k_{0}\left<\mathbf{E}_{i}\right|\mathbb{W}_{\mathrm{obj}}\mathbb{V}_{\mathrm{obj}}\left|\mathbf{E}_{i}\right>\right]$$

$$~{}~{}+\frac{1}{2Z}\operatorname{Im}\left[k_{0}\left<\mathbf{E}_{i}\right|\mathbb{W}_{\mathrm{obj}}\left|\textbf{T}_{\mathrm{d}}\right>\right]$$ (21)

with $\mathbb{T}_{\mathrm{d}}$ subject to the fundamental relation

$$\mathbb{I}_{\mathrm{des}} =\mathbb{I}_{\mathrm{des}}\left(\chi_{\mathrm{obj}}^{-1}\mathbb{I}-\mathbb{G}_{\mathrm{obj}}\right)\mathbb{I}_{\mathrm{des}}\mathbb{T}_{\mathrm{d}},$$ (22)

where $\chi_{\mathrm{clk}}$ is the susceptibility used for the object in the design domain (the cloak). This constraint is similar to that in previous works [14, 28], but where the relevant Green’s function is $\mathbb{G}_{\mathrm{obj}}$ and not the vacuum Green’s function. Use of $\mathbb{W}_{\mathrm{obj}}=\mathbb{G}_{\mathrm{vac}}^{-1}\mathbb{G}_{\mathrm{obj}}$ leads to the primal objective function presented in the main text.

Appendix B Topology optimization procedure

To calculate inverse designs, we made use of the NLopt package [43] and followed standard topology optimization algorithms [26, 44] based on the method of moving asymptotes [45] where the electric susceptibility value of each pixel in the specified design region is considered as an independent design parameter. Each pixel in the design region is allowed to explore a continuous range of susceptibility values varying linearly between the vacuum value $0$ and the prescribed material value $\chi_{\mathrm{clk}}$. The continuous range $[0,1]\cdot\chi_{\mathrm{clk}}$ can be interpreted, after normalization by $\chi_{\mathrm{clk}}$, as a filling fraction of the pixel.

We implemented a 2D Maxwell solver which could be used for arbitrary structuring (allow the pixels within any arbitrary 2D design region to vary). Doing a thousand function evaluations (iterations of optimization over the susceptibility profile for fixed material parameters, design region, bandwidth, etc.) in 2D takes less than 1 hour for a wavelength-scale design region. The primary challenge in such optimization problems lies in the computational cost of evaluating the objective function and gradient, which is evaluated several hundred if not thousand times during the course of any one optimization. Let $\mathbb{G}$ denote the total Green’s function and let $\mathbb{M}^{-1}=\mathbb{G}$ where $\mathbb{M}=\frac{c^{2}}{\omega_{0}^{2}}\nabla\times\nabla-\epsilon(\mathbf{r})$. After the discretization of the computational grid, then in the design region $\epsilon_{b}=1+\chi_{\mathrm{clk}}\bar{\chi}_{b}$ at pixel $b$. The gradient of the objective function $P_{\mathrm{ext}}(\{\bar{\chi}_{b}\};\omega_{0})=\frac{1}{2}\operatorname{Re}\left[\left<\mathbf{E}_{i}|\mathbf{J}_{g}\right>\right]=\frac{1}{2}\operatorname{Re}\left[\left<\mathbf{E}_{i}|\mathbb{V}\mathbb{G}|\mathbf{J}_{i}\right>\right]$ with respect to the degree of freedom $\bar{\chi}_{a}$ is given by

$$\frac{\partial P_{\mathrm{ext}}}{\partial\bar{\chi}_{a}} =\frac{1}{2}\operatorname{Re}\left[\left<\mathbf{E}_{i}|\frac{\partial(\mathbb{V}\mathbb{G})}{\partial\bar{\chi}_{a}}|\mathbf{J}_{i}\right>\right]$$ (23)

$$=\frac{1}{2}\operatorname{Re}\left[\left<\mathbf{E}_{i}|(\chi_{\mathrm{clk}}\mathbb{P}_{a}\mathbb{G}+\mathbb{V}\frac{\partial\mathbb{M}^{-1}}{\partial\bar{\chi}_{a}})|\mathbf{J}_{i}\right>\right]$$ (24)

$$=\frac{1}{2}\operatorname{Re}\left[\left<\mathbf{E}_{i}|(\chi_{\mathrm{clk}}\mathbb{P}_{a}\mathbb{G}-\mathbb{V}\mathbb{M}^{-1}\frac{\partial\mathbb{M}}{\partial\bar{\chi}_{a}}\mathbb{M}^{-1})|\mathbf{J}_{i}\right>\right]$$ (25)

$$=\frac{1}{2}\operatorname{Re}\left[\left<\mathbf{E}_{i}|(\chi_{\mathrm{clk}}\mathbb{P}_{a}\mathbb{G}+\chi_{\mathrm{clk}}\mathbb{V}\mathbb{G}\mathbb{P}_{a}\mathbb{G})|\mathbf{J}_{i}\right>\right]$$ (26)

$$=\frac{1}{2}\operatorname{Re}\left[\frac{k_{0}\chi_{\mathrm{clk}}}{iZ}\left<\mathbf{E}_{i}|(\mathbb{P}_{a}+\mathbb{V}\mathbb{G}\mathbb{P}_{a})|\mathbf{E}_{t}\right>\right]$$ (27)

$$\equiv\frac{1}{2}\operatorname{Re}\left[\frac{k_{0}\chi_{\mathrm{clk}}}{iZ}\langle\mathbf{E}_{i}+\tilde{\mathbf{E}}|\mathbb{P}_{a}|\mathbf{E}_{t}\rangle\right]$$ (28)

where $\bar{\chi}_{a}$ is a topology optimization degree of freedom, normalized to take values in $[0,1]$, so that the susceptibility at the $a$-th pixel in the design region is given by $\chi_{\mathrm{clk}}\bar{\chi}_{a}$, $\mathbb{P}_{a}$ is projection onto the pixel indexed by $a$, $|\mathbf{E}_{t}\rangle$ is the total electric field, and $|\tilde{\mathbf{E}}\rangle\equiv\mathbb{G}^{\dagger}\mathbb{V}^{\dagger}|\mathbf{E}_{i}\rangle$ which for reciprocal systems becomes $|\tilde{\mathbf{E}}\rangle\equiv(\mathbb{G}\mathbb{V}|\mathbf{E}_{i}^{*}\rangle)^{*}$. Hence, for a given incident field $|\mathbf{E}_{i}\rangle$ one only needs to perform two additional Maxwell solves to calculate $|\mathbf{E}_{t}\rangle$ and $|\tilde{\mathbf{E}}\rangle$ at each iteration of the optimization to obtain the required gradient Eq. (28) and the objective function $P_{\mathrm{ext}}(\{\bar{\chi}_{b}\};\omega_{0})=\frac{1}{2}\operatorname{Re}\left[\left<\mathbf{E}_{i}|\mathbf{J}_{g}\right>\right]=\frac{1}{2Z}\operatorname{Im}\left[k_{0}\left<\mathbf{E}_{i}|\mathbb{V}|\mathbf{E}_{t}\right>\right]$. Details on the computational complexities encountered in 3D and potential solution methods can be found in Appendix C of Ref. [14] and the therein cited references.