1. Introduction

The design of electromagnetic devices is a challenging engineering effort due to the complex dynamics of electromagnetic waves in dielectric and conductive media. Those rich dynamics, even in linear media, give rise to a multitude of effects—from resonances and waveguiding to bandgaps, metamaterials and topological effects—that empower researchers and engineers to design progressively more complex and efficient systems and devices. Optimization plays a key role in this process and it can be as simple as tuning a few geometrical parameters of the structure, or as complex as designing the whole structure itself based on high-level functionality goals. The latter case is usually referred to as topology optimization, where the shape of the device, or its material composition and distribution, are somehow mapped to the optimization variables.

Recent decades have seen a constant progress in topology optimization in electromagnetic devices over the whole frequency spectrum between RF and UV. Since the pioneering work by Bendsøe and Kikuchi [1], considerable improvements have been made on topology optimization, giving rise to several different approaches, including the so-called Solid Isotropic Material with Penalization (SIMP) or density-based method [2], which can be considered as the most popular method employed for topology optimization.

The reviews by Jensen and Sigmund [3] and by Molesky et al. [4] offer a thorough picture of this progress, with particular interest in the field of nanophotonics. The short wavelengths combined with advanced and precise fabrication methods for large areas (with respect to the wavelength) make topology optimization in nanophotonics particularly challenging, but, at the same time, a rewarding effort. Important contributions in the field have been obtained by the density-based method, from which we highlight the pioneering works [5–10]. Time and again, non-intuitive designs obtained through the most diverse topology optimization methods have demonstrated figures of merit far beyond their intuitive counterparts, including bio-inspired methods [11–13], transformation optics [14–16], level-set methods [17–19], perturbation theory techniques [20], objective-first formulation [21] and pixel-by-pixel optimization method [22,23]. Recent proposals also deal with manufacturing uncertainties and fabrication constraints [24–29], which are important to guarantee that the optimized designs can be fabricated using current technology.

Advancing the development of efficient topology optimization methods, in this work we present a formulation for the electromagnetic scattering problem based on the topological derivative method. The concept of topological derivative was first proposed in the context of shape and topology optimization [30]. A complete description of the method can be found in [31]. This relatively new concept has applications in many different fields such as electromagnetic scatterers [32], band-gap structures [33], image processing [34], multi-scale constitutive modeling [35] and, fracture [36] and damage [37] mechanics. For an account on the theoretical development and applications of the topological derivative method, see the series of review papers [38–40] and references therein.

The main idea of the topological derivative method is to compute the variation of a predefined functional over a change in the topology of a given domain. More intuitively, if the domain is composed of two materials, it calculates the variation of the functional when those materials are swapped in a given location, similar to a gradient of the functional. The topological derivative can, therefore, be used in any gradient-based optimization algorithm. In the context of topological-derivative-based topology optimization methods, the algorithms available in the literature usually combine topological derivatives with shape derivatives or level-set methods [41–44], leading to a two-stage topology/shape optimization procedure. More precisely, new features are nucleated according to the topological derivative, while standard tools in shape optimization are used to move the new boundaries. Following the original ideas proposed by Amstutz & Andrä [45], in this paper a topology optimization algorithm based on the topological derivative together with a level-set domain representation method is devised. The basic idea consists in achieving a local optimality condition written in terms of the topological derivative and a level-set function, leading to an efficient one-stage algorithm driven by the topological derivative only. In particular, the associated topological derivative is rigorously derived and used as gradient descent direction to find local optima for the inverse electromagnetic design problem. We demonstrate that the resulting topology design algorithm is remarkably efficient and of simple computational implementation, since it features a minimal number of user-defined algorithmic parameters.

One important difference of the topological derivative formulation from more conventional approaches based on the level-set method [46], is that it does not require one to solve any auxiliary evolution problem such as the Jacobi equation and, more importantly, the resulting algorithms do not require any special features from the initial condition, such as a minimal number of holes in the domain. The implementation of the algorithm requires only the solution of the original electromagnetic problem and its adjoint, both which have the same formulation and, as such, can be solved by the same Finite Element Method (FEM) implementation (for instance, COMSOL Multiphysics’ electromagnetic solver) with just a change of the source term. From both solutions, the derivative of the objective function with respect to the device topology is directly obtained, i.e., no further calculations based on the chain rule for derivatives are necessary, which makes the algorithm straightforward to implement for a variety of objective functions.

Furthermore, whereas in density-based method a gradient vector is calculated with respect to continuous changes in the dielectric constant (or any other material parameter) over the optimization domain, the topological derivative, as the name implies, is calculated with respect to the topology itself. That means when two materials are chosen for the design problem—for example, Si and SiO₂, as is often the case in nanophotonics—the topological derivative results in designs where the dielectric constant is either that of Si or SiO₂. On the other hand, density-based methods result in dielectric constant values in the range between those of SiO₂ and Si, what is usually called a gray-scale design. It is often the case that a binarization step is then required to discretize the design, possibly modifying its calculated optimal response, although recent advances in projection techniques [47–52] are capable of both tackling the binarization issue and offering accurate control of the minimal length scales in the design, which is of utmost importance for any practical purposes in any optimization method. Fabrication constraints can also be included in the topological derivative formulation as a penalty in the level-set function used to characterize the topology [29]. As a proof-of-concept, we demonstrate the implementation of a simple hard constraint for conventional fabrication technology in integrated photonics in two of the presented examples. Further exploration of the best techniques to introduce such constraints—particularly in three-dimensional problems—are out of the scope of this work.

This paper is organized as follows: in Section 2, the electromagnetic scattering problem is defined, where the goal is to design a dielectric (possibly lossy) device based on a predefined objective function. Section 3 introduces the topological derivative method in the context of electromagnetic scattering problems. Finally, Section 4 presents several design examples obtained from different objective function definitions both in two and three spatial dimensions, before a short summary of this work. The theoretical foundations and the derivation of the topological derivative associated with the electromagnetic problem are described in the Appendix A.

2. Model problem

Let us consider an open and bounded domain $\mathcal {D} \subset \mathbb {R}^{d}$, with $d=2,3$, as illustrated in Fig. 1. A near-field domain is denoted as $\mathcal {B}_R \subset \mathcal {D}$, with boundary $\Gamma = \partial \mathcal {B}_R$. The optimization region is given by $\Omega \subset \mathcal {B}_R$, which is split into two disjoints subdomains $\Omega _a$ and $\Omega _b$, such that $\Omega = \Omega _a \cup \Omega _b$ and $\Omega _a \cap \Omega _b = \emptyset$, with $\Omega _a$ and $\Omega _b$ representing regions with refractive index $n_1$ and $n_2$ respectively (light and dark regions in Fig. 1). The weak form of the electromagnetic scattering problem (assuming linear and isotropic media) is stated in the conventional FEM formulation [53,54]: find $E \in \mathcal {V}$, such that

 

Fig. 1. Domain definitions for the calculation of the topological derivative. The electromagnetic scattering problem is defined in the domain $\mathcal {D}$, including sources, boundary conditions and PMLs as usual in FEM formulations. The optimization domain $\Omega$ is partitioned between $\Omega _a$ and $\Omega _b$, which are represented by light and dark regions, respectively. The topology and shape of those regions are the subject of the optimization of a shape functional, defined in terms of integrals over the boundary of a near-field domain $\mathcal {B}_R$.

Download Full Size | PDF

The goal is to maximize or minimize a shape functional given, in general form, by

3. The topological derivative method

A quite general approach for dealing with such topology optimization problem is based on the concept of topological derivative [31,56], which, in our particular case, represents the (topological) sensitivity of the shape functional $\psi (E)$ with respect to the nucleation of a small inclusion confined in $\Omega$ endowed with different material properties from the background. To be more precise, the associated topological derivative is given here by the following result:

Theorem 1 Let $\psi (E)$ be the shape functional (5), with $E$ used to denote the solution to the variational problem Eq. (1). Then, its topological derivative can be written as

Proof 1 The proof of this results is presented in the Appendix A.

Corollary 2 From the definition for the contrast $\gamma (x)$ in Theorem 1, the topological derivative can be rewritten as

It is important to notice that the solution of Eq. (9), required to find $D_T \psi (x)$ in Eq. (10), can be obtained by the same solver used to evaluate the original FEM formulation Eq. (1), except for the change in the source term on the right-hand side. Therefore, the implementation of the topological derivative-based optimization algorithm does not depend on the implementation of new algorithms. Many commercial and open source FEM solvers can be used to implement Eq. (9), including COMSOL Multiphysics and FEniCS [57].

4. Numerical results

In this section, a topology optimization algorithm based on the topological derivative combined with a level-set domain representation method is presented. It has been proposed in [45] and it consists basically in achieving a local optimality condition for the optimization problem we are dealing with, given in terms of the topological derivative and a level-set function. In particular, the domains $\Omega _a$ and $\Omega _b$ are characterized by a level-set function $\zeta$ of the form:

The process ends when the condition $\theta _j \leq \varepsilon _\theta$ is satisfied at some iteration, where $\varepsilon _\theta$ is a given small numerical tolerance. If $\alpha$ is found to be smaller then a given numerical tolerance and the local optimality condition is not yet satisfied, namely $\theta _j\;>\;\varepsilon _\theta$, then a mesh refinement of the domain $\Omega$ is carried out in iteration $j$ and the process is continued.

In practice, the objective functions in all the following examples converged to acceptable values much sooner than the optimality condition on $\theta _j$, which may require extreme mesh refinements. That means that stopping criteria based on electromagnetic design requirements, such as insertion loss above −1 dB or reflection below −20 dB, are reached faster than the absolute convergence of the problem, which can save on computation time.

In the following sections we present a few examples of optimization results for two and three-dimensional devices. Each example uses a different objective function $\Psi (\Omega )$ defined in terms of transmission or reflection coefficients (scattering parameters) or far field radiation. They are intended to demonstrate the flexibility offered by the optimization method. Information on the required computational resources is reported in Appendix B.

4.1 Reflector

The goal in this first example is to design a reflector for a conventional silicon rib waveguide at 1.55 µm. The objective of the optimization is to maximize the reflection coefficient, i.e., maximize $|S_{11}|^{2}$ for a 1-port device at a single wavelength.

To define the objective function in terms of scattering parameters $S_{\ell m}$, the input and output ports of the device must be defined over $\Gamma$. With the normalized input mode $m$ excited, the coefficient can be calculated by projecting the electromagnetic fields onto the normalized output port mode $(e_\ell , h_\ell )$ [60–62]:

For the reflector, we can define the objective function to be minimized

The optimization domain is defined as a 2.5 µm × 2.5 µm area at the end of a 450 nm-wide waveguide, illustrated in Fig. 2(a). The electric field is perpendicular to the simulation domain and the refractive indices used are of SiO₂ and that of a 220 nm Si slab waveguide in SiO₂ with the appropriate polarization. All other 2D simulations presented here use these same polarization and refractive indices.

 

Fig. 2. Optimized reflector design. (a) Simulated electric field magnitude $|E_z|$ in linear scale (PML regions excluded). The white contours indicate material boundaries. (b) Evolution of the reflection coefficient $|S_{11}|^{2}$ and the absolute convergence parameter $\theta$ during the optimization. The values above the plot indicate the refinement of the optimization mesh in nanometers. (c) Refractive index distributions at 4 iteration steps: initial guess, first iteration, 15^th iteration, and final result. Dark and light areas represent Si and SiO₂, respectively.

Download Full Size | PDF

The results of the optimization can be seen in Fig. 2. The electric field magnitude distribution along the waveguide in Fig. 2(a) distinctively presents the standing wave profile expected from a good reflector. The final design presents a reflection coefficient of 0.94 (−0.29 dB). The convergence plot of Fig. 2(b) shows that the objective function quickly converges, surpassing −1 dB even before the first mesh refinement. The absolute convergence parameter $\theta$ is more noisy and converges more slowly than the objective function. In larger and more complex problems it is a good idea to use the objective function as a stopping criterion, based on practical device requirements.

It is interesting that the shape of the final reflector resembles a conventional grating, except for a few deformations. A grating would be the intuitive way to manually design a reflector and, although the initial condition for the optimization was a homogeneous block of Si, the algorithm started approaching the grating-like design already on the first step, as shown in Fig. 2(c). It is possible to identify regions at the sides of the main radiation direction where the design seems to have an “inverted phase” and a higher fraction of Si. Towards the far end of the reflector it also loses its spherical appearance. Finally, we note that the region directly in touch with the input waveguide appears to form a rounded cap that shapes the outgoing waves to conform to the grating.

4.2 Power splitter

An unbalanced power splitter can be designed by combining the transmission coefficients from several outputs in the objective function. This example demonstrates the 2D optimization of a $2\mathbin {:}1$ power splitter through the following objective function and topological derivative:

The optimization area measures 4.0 µm × 3.0 µm and is initialized as a uniform block of Si. The resulting device is presented in Fig. 3(b) in conjunction with the connecting waveguides. The plot in Fig. 3(a) shows the wavelength response achieved by optimizing at the 2 selected wavelengths: the splitter presents a flat response in most of the 300 nm range simulated. The wideband simulation was performed by Finite Difference in Time Domain (FDTD) method using MEEP [63], and it matches the FEM simulations at both wavelengths used for optimization. The electric field magnitude distributions at those wavelengths in Figs. 3(c) and (d) does not seem to contain any resonant structures, which probably assists in the flat splitter response.

 

Fig. 3. Optimized $2\mathbin {:}1$ power splitter. (a) Wavelength dependency of the transmission coefficients for the optimized device. The square and circle markers indicate the results and wavelengths used during optimization, and the lines correspond to an FDTD evaluation of the design. The cross and plus markers indicate the results for a fabrication-constrained optimization. (b) Refractive index distribution at the last optimization iteration. (c, d) Simulated electric field magnitude $|E_z|$ in linear scale at the optimization wavelengths. The white contours indicate material boundaries. (e) Refractive index distribution of the optimized device with a fabrication constraint to ensure curvatures of at least 50 nm radius.

Download Full Size | PDF

The device was also optimized with hard fabrication constraints applied at each iteration in the form of opening and closing operations [59] in the material distribution to remove small features and ensure a minimal radius of 50 nm throughout the optimization region. The resulting topology, displayed in Fig. 3(e), is similar to the previous one, except for the removal of small features. The transmission coefficients, which can be seen as the cross and plus markers in Fig. 3(a), are within 5% of the unconstrained values, thus showing that fabrication constraints can be included in the proposed algorithm.

4.3 Modal multiplexer

This 2D example demonstrates the possibility of using the same 4.0 µm × 3.0 µm area as before to optimize a more complex device: a wideband 2-mode multiplexer. The objective function has to be altered to include both input modes, and, therefore, reads:

The optimization area, initial conditions and wavelengths are the same as for the power splitter, but the resulting topology in Fig. 4(b) is quite different. The spectral response is extremely flat, even though the optimization was performed at only 2 wavelengths, covering more than 300 nm with transmission above −1.5 dB. Similar to the power splitter case, this happens because the small wavelength difference between the simulations (less than 5%) is enough to inhibit the appearance of strongly resonant structures, as can be seen in in the field plots of Figs. 4(c)–(f).

 

Fig. 4. Optimized 2-mode multiplexer. (a) Wavelength dependency of the transmission coefficients for the optimized device. The markers indicate the wavelengths used during optimization, and the lines correspond to FDTD evaluations of the design. (b) Refractive index distribution at the last optimization iteration. (c, d) Simulated electric field magnitude $|E_z|$ in linear scale at the optimization wavelengths for the first input waveguide mode. The white contours indicate material boundaries. (e, f) Same as before but for the second mode of the input waveguide.

Download Full Size | PDF

4.4 Diplexer

The optimization over several wavelengths can also be used to tailor the spectral response of the device beyond a flat curve. Specific filters can be designed using this feature. As an example, we design a compact diplexer that splits signals from the O and C bands using $\lambda \in$ {1.31 µm, 1.55 µm}. In this case the objective function is

 

Fig. 5. Optimized diplexer for 1.31 µm and 1.55 µm. (a) Wavelength dependency of the transmission coefficients for the optimized device. The square and circle markers indicate the results and wavelengths used during optimization, and the lines correspond to an FDTD evaluation of the design. The FDTD simulation was divided in several wavelength windows due to the wide frequency range desired and the computational effort required. Stitching the results naturally leads to the small discontinuities found in the plots. The cross and plus markers indicate the results for a fabrication-constrained optimization. (b) Refractive index distribution at the last optimization iteration. (c, d) Simulated electric field magnitude $|E_z|$ in linear scale at the optimization wavelengths. The white contours indicate material boundaries. (e) Refractive index distribution of the optimized device with fabrication constraints.

Download Full Size | PDF

In contrast to the previous designs, the resulting topology presents some resonant features, as can be seen in the electric field magnitude distributions in Figs. 5(c) and (d) in the form of more intense and localized spots. Such resonant features help define the spectral response of the diplexer, presented in Fig. 5(a), which directs each wavelength range towards a different output port with less than 1 dB insertion loss. The exact values for the insertion losses at 1.31 µm and 1.55 µm are, respectively, −0.88 dB and −0.77 dB, which are very similar to the results reported by Piggott et al. [64] using an objective-first, adjoint-based method. Once again, the lines are results from a FDTD simulation to evaluate the device response in a wider range of wavelengths.

This example was also re-optimized with hard fabrication constraints in place in order to demonstrate that they can also be included in 3D designs. In this case, besides the opening and closing operations in the propagation plane to remove small features (exactly as used in the power splitter example) another constraint in the perpendicular direction was enforced that eliminated material changes in this direction, i.e., the device can be fabricated in a single etch step that removes the full 220 nm silicon layer. For each position of the optimization domain in the horizontal plane, the material for all elements along the perpendicular direction (perpendicular to the figure plane) is chosen by simple majority of the signs of the level-set function for those elements. The results for the insertion losses at 1.31 µm and 1.55 µm are, respectively, −2.4 dB and −2.1 dB, significantly lower than the unconstrained results. This is a direct result of the small features found in the unconstrained run. The room for improvement of those results indicate that further exploration of alternative constraint formulations should be pursued in the future, including the allowance of pre-determined etch depths and the use of a penalty function instead of a hard constraint.

4.5 Dual-polarization fiber coupler

This example proposes the design of a coupling structure analogous to a grating coupler for a Si waveguide with cross-section 450 nm × 220 nm. However, the design targets the modes of a photonic crystal fiber (PCF) instead of a conventional single-mode fiber, and it aims to couple both $x$ and $y$ polarization modes into the waveguide’s TE and TM modes, respectively. The challenge in this example is that no intuitive design exists that produces reasonable coupling coefficients, specially considering vertical fiber coupling. In Figs. 6(a) and (b) we can see the electric fields distributions for the fiber and waveguide modes.

 

Fig. 6. Optimized dual-polarization coupler. (a) Profile for the main component of the electric field for the $x$ ($|E_x|$) and $y$ ($|E_y|$) polarization modes of the PCF. (b) Same as (a) for the TE ($|E_x|$) and TM ($|E_y|$) waveguide modes. (c) Refractive index distribution at the last optimization iteration (central cuts in the $xy$ and $zx$ planes). (d) Simulated electric field magnitude in linear scale for the TE mode (same cut planes). The white contours indicate material boundaries. (e) Same as (d) for the TM mode.

Download Full Size | PDF

The objective function is similar to Eq. (27), except that the optimization is performed at a single wavelength: 1.55 µm. The resulting topology, in Fig. 6(c), displays more details and higher complexity than previous examples. This is both because the optimization region is larger (5.0 µm × 5.0 µm × 220 nm) and the problem itself is more complex. Nonetheless, the dual-polarization coupler presents coupling efficiencies of −2.50 dB and −2.89 dB for the TE and TM modes, respectively. The cross-talks from TE to $y$ polarization and from TM to $x$ polarization are, respectively, −19.6 dB and −29.1 dB. Although fabrication constraints would still have to be imposed on the geometry, we believe these are encouraging results for a long-sought, compact, dual polarization coupler for silicon photonics.

4.6 Dielectric nanoantenna

In this example, our goal is to design a compact Si nanoantenna with main radiation lobe perpendicular to the substrate. The shape functional that will be optimized is the far field of the antenna, given by

Because $\varphi : \mathcal {V} \mapsto \mathbb {C}^{d}$ is a complex-valued vector function, then $\langle \varphi ^{\ast }(E), \varphi (E) \rangle = \varphi ^{\ast }(E) \cdot \varphi (E)$. In this particular case, the adjoint problem Eq. (9) reads: find $V \in \mathcal {V}$, such that

In order to maximize the far field, we define the objective function to be minimized as

The optimization results in a design with reflection coefficient of −16.9 dB at wavelength 1.55 µm and antenna gain of 18.8 dB, as shown in Fig. 7(c). The device topology and electric field magnitude distribution at 1.55 µm can be seen in Figs. 7(a) and (b). It is particularly interesting to observe the field intensities radiating in the forward and upward directions in Fig. 7(b), which match the 2 largest radiation lobes in the gain pattern of Fig. 7(c) in the $zx$ plane, as it would be expected.

 

Fig. 7. Optimized dielectric nanonatenna. (a) Refractive index distribution at the last optimization iteration (central cuts in the $xy$ and $zx$ planes). (b) Simulated electric field magnitude in linear scale (same cut planes). The white contours indicate material boundaries. (c) Antenna gain along 2 orthogonal radiation directions. The main lobe at 0° (upward emission) is aligned with the $z$ axis and, for the continuous curve, the 90° direction is aligned with the $x$ axis (forward emission).

Download Full Size | PDF

5. Conclusions

In this work an optimization algorithm for the electromagnetic design problem based on the topological derivative method has been proposed. Because of the general formulation of the problem, the possibilities for the shape functional and objective function are unlimited. Several examples of optimized devices in two and three spatial dimensions are presented with different objective functions, including both near field and far field formulations. All examples converge to successful results starting from a block of Si or SiO₂, indicating the robustness of the algorithm to initial conditions. Other useful formulations for the objective function that can be explored in the future include the total radiated power or the scattering cross-section of the optimized domain. Those could be used to design dielectric particles with enhanced emission or specific effective areas. Besides its robustness and the flexibility in the definition the objective function, the topological derivative-based optimization method has two other important features: straightforward implementation, since the adjoint problem can be solved by the same solver used for the electromagnetic scattering problem; and binary search space by definition, which avoids the need for a binarization of the optimal solution. For these reasons, we believe the topological derivative-based method will become a valuable tool in the design of compact and complex optical components in the future.

Appendix A. Topological sensitivity analysis

In this appendix, we present the proof of the main theoretical result of the paper stated through Theorem 1. Let us start by introducing a topological perturbation confined in a small ball $\mathcal {B}_\varepsilon (x_0) \subset \Omega$ of size $\varepsilon$ and center at $x_0 \in \Omega$ of the form (see sketch in Fig. 8)

(33)$$\gamma_\varepsilon(x) = \begin{cases} \gamma(x), & x\in \mathcal{B}_\varepsilon(x_0) \\ 1, & \textrm{otherwise} \end{cases}$$

Therefore, the topologically perturbed counterpart of the shape functional Eq. (5) is given by

(34)$$\psi(E_\varepsilon) = \langle \varphi(E_\varepsilon), \varphi^{{\ast}}(E_\varepsilon) \rangle,$$

where $E_\varepsilon$ is solution of the following variational problem: Find $E_\varepsilon \in \mathcal {V}$, such that

(35)$$\int_{\mathcal{D}}\left(\nabla\times E_\varepsilon \cdot\nabla\times W - k_0^{2} n^{2} \gamma_\varepsilon^{2} E_\varepsilon \cdot W\right)\mathrm dx = S(W) \quad \forall \; W \in \mathcal{V}.$$

Now, we have all elements to derive a close form for the associated topological derivative, which represents the main theoretical result of the paper.

 

Fig. 8. The topological derivative of a shape functional $\psi (E)$ with respect to an inclusion at $x_0$ in the optimization domain is related to the variation of the functional in similar fashion to a conventional derivative. The variation of the functional appears as the result of a refractive index exchange in a small ball of size $\varepsilon$ centered at $x_0$.

Download Full Size | PDF

Appendix A.1. Existence of the topological derivative

The existence of the topological derivative associated to the problem we are dealing with is ensured by the following result:

Lemma 3 Let $E$ and $E_\varepsilon$ be the solutions of the variational problems Eq. (1) and Eq. (35), respectively. Then, the following a priori estimate holds true:

(36)$$\| \widetilde{E}_\varepsilon \|_{H_{\mathrm{curl}}(\mathcal{D};\mathbb{C}^{d})} \leq C \varepsilon^{\frac{d}{2} + \delta},$$

with constant $C$ independent of the small parameter $\varepsilon$ and $0\;<\;\delta\;<\;d/2$, where $\widetilde {E}_\varepsilon = E_\varepsilon -E$.

Proof 2 We start by rewriting Eq. (1) as follows:

(37)$$\int_{\mathcal{D}}\left(\nabla\times E \cdot\nabla\times W - k_0^{2} n^{2} \gamma_\varepsilon^{2} E \cdot W\right)\mathrm dx = S(W) - \int_{\mathcal{B}_\varepsilon}(\gamma^{2}-1)k_0^{2}n^{2} E \cdot W\,\mathrm dx \quad \forall \; W \in \mathcal{V},$$

where we have used the definition for the contrast Eq. (33). By subtracting Eq. (37) from Eq. (35) we obtain

(38)$$\int_{\mathcal{D}}\left(\nabla\times \widetilde{E}_\varepsilon \cdot\nabla\times W - k_0^{2} n^{2} \gamma_\varepsilon^{2} \widetilde{E}_\varepsilon \cdot W\right)\mathrm dx = \int_{\mathcal{B}_\varepsilon}(\gamma^{2}-1)k_0^{2}n^{2} E \cdot W\,\mathrm dx \quad \forall \; W \in \mathcal{V},$$

with $\widetilde {E}_\varepsilon = E_\varepsilon -E$. Let us consider a decomposition for $\widetilde {E}_\varepsilon$ of the form

(39)$$\widetilde{E}_\varepsilon = G_\varepsilon + H_\varepsilon,$$

where $G_\varepsilon \in \mathcal {V}$ is solution to

(40)$$\int_{\mathcal{D}} \nabla \times G_\varepsilon \cdot\nabla\times W \, \mathrm dx = \int_{\mathcal{B}_\varepsilon}(\gamma^{2}-1)k_0^{2}n^{2} E \cdot W\,\mathrm dx \quad \forall \; W \in \mathcal{V},$$

whereas $H_\varepsilon \in \mathcal {V}$ is solution to

(41)$$\int_{\mathcal{D}}\left(\nabla\times H_\varepsilon \cdot\nabla\times W - k_0^{2} n^{2} \gamma_\varepsilon^{2} H_\varepsilon \cdot W\right)\mathrm dx = \int_{\mathcal{D}} k_0^{2} n^{2} \gamma_\varepsilon^{2} G_\varepsilon \cdot W \, \mathrm dx \quad \forall \; W \in \mathcal{V},$$

From the well-posedness of the above variational problem [53,54], we have

(42)$$\| H_\varepsilon \|_{H_{\mathrm{curl}}(\mathcal{D};\mathbb{C}^{d})} \leq C_1 \| G_\varepsilon \|_{L^{2}(\mathcal{D};\mathbb{C}^{d})} \leq C_1 \| G_\varepsilon \|_{H_{\mathrm{curl}}(\mathcal{D};\mathbb{C}^{d})}.$$

By setting $W = \overline {G}_\varepsilon$ as test function in the variational problem Eq. (40) we obtain the equality

(43)$$\int_{\mathcal{D}} \nabla\times G_\varepsilon \cdot\nabla\times \overline{G}_\varepsilon \, \mathrm dx = \int_{\mathcal{B}_\varepsilon}(\gamma^{2}-1)k_0^{2}n^{2} E \cdot \overline{G}_\varepsilon \,\mathrm dx.$$

From the Poincaré inequality [55], the above bilinear form is bounded by below as follows

(44)$$c \| G_\varepsilon \|^{2}_{H_{\mathrm{curl}}(\mathcal{D};\mathbb{C}^{d})} \leq \int_{\mathcal{D}} \nabla\times G_\varepsilon \cdot \nabla \times \overline{G}_\varepsilon \, \mathrm dx,$$

so that we get

(45)$$\| G_\varepsilon \|^{2}_{H_{\mathrm{curl}}(\mathcal{D};\mathbb{C}^{d})} \leq C_2 \int_{\mathcal{B}_\varepsilon}(\gamma^{2}-1)k_0^{2}n^{2} E \cdot \overline{G}_\varepsilon \,\mathrm dx,$$

with $C_2 = 1/c$. The Cauchy-Schwartz inequality [55] yields

(46)$$\| G_\varepsilon \|^{2}_{H_{\mathrm{curl}}(\mathcal{D};\mathbb{C}^{d})} \leq C_3 \varepsilon^{\frac{d}{2}} \| G_\varepsilon \|_{L^{2}(\mathcal{B}_\varepsilon;\mathbb{C}^{d})}.$$

Notice that, Hölder inequality [55] and the Sobolev embedding theorem [65] can be used to derive

(47)$$\| G_\varepsilon \|_{L^{2}(\mathcal{B}_\varepsilon;\mathbb{C}^{d})} \leq C_4 \varepsilon^{\frac{d}{2q}} \| G_\varepsilon \|_{L^{2p}(\mathcal{B}_\varepsilon;\mathbb{C}^{d})} \leq C_4 \varepsilon^{\delta} \| G_\varepsilon \|_{H_{\mathrm{curl}}(\mathcal{D};\mathbb{C}^{d})},$$

for any $1\;<\;q\;<\;\infty$, with $1/p+1/q = 1$, and $\delta = \frac {d}{2q}$. Therefore

(48)$$\| G_\varepsilon \|_{H_{\mathrm{curl}}(\mathcal{D};\mathbb{C}^{d})} \leq C_5 \varepsilon^{\frac{d}{2} + \delta} \quad \textrm{and} \quad \| H_\varepsilon \|_{H_{\mathrm{curl}}(\mathcal{D};\mathbb{C}^{d})} \leq C_6 \varepsilon^{\frac{d}{2} + \delta}.$$

Finally, from the triangular inequality in Eq. (39) combined with the above estimates, we obtain the required result for any $0\;<\;\delta\;<\;d/2$.

Remark 4 The estimate Eq. (36) in Lemma 3 can be written simply as

(49)$$\| \widetilde{E}_\varepsilon \|_{H_{\mathrm{curl}}(\mathcal{D};\mathbb{C}^{d})} = o(\varepsilon^{\frac{d}{2}}),$$

provided that $\delta\;>\;0$. The notation $o(f(\varepsilon ))$ has to be understood as

(50)$$\lim_{\varepsilon \rightarrow 0}\frac{o(f(\varepsilon))}{f(\varepsilon)} = 0,$$

with $f(\varepsilon )$ used to denote a real-valued function that goes monotonically to zero when $\varepsilon \rightarrow 0$.

Appendix A.2. Proof of the main result

Let us subtract Eq. (5) from Eq. (34) to obtain

(51)$$\begin{aligned} \psi(E_\varepsilon) - \psi(E) &= \langle \varphi(E_\varepsilon), \varphi^{{\ast}}(E_\varepsilon) \rangle - \langle \varphi(E), \varphi^{{\ast}}(E) \rangle\\ &= 2\Re \{\langle \varphi^{{\ast}}(E), \varphi(\widetilde{E}_\varepsilon) \rangle \} + \mathcal{E}_1(\varepsilon), \end{aligned}$$

in which $\widetilde {E}_\varepsilon = E_\varepsilon - E$ and the remainder $\mathcal {E}_1(\varepsilon )$ is defined as

(52)$$\mathcal{E}_1(\varepsilon) := \langle \varphi(\widetilde{E}_\varepsilon), \varphi^{{\ast}}(\widetilde{E}_\varepsilon) \rangle.$$

The Cauchy-Schwartz inequality together with trace theorem [55] and Lemma 3 yields

(53)$$|\mathcal{E}_1(\varepsilon)| \leq \| \widetilde{E}_\varepsilon \|^{2}_{L^{2}(\Gamma;\mathbb{C}^{d})} \leq \| \widetilde{E}_\varepsilon \|^{2}_{L^{2}(\mathcal{D};\mathbb{C}^{d})} \leq \| \widetilde{E}_\varepsilon \|^{2}_{H_{\mathrm{curl}}(\mathcal{D};\mathbb{C}^{d})} = o(\varepsilon^{d}),$$

provided that $\varphi (\widetilde {E}_\varepsilon )$ is assumed to be linear with respect to the argument $\widetilde {E}_\varepsilon$ and written in terms of integrals concentrated on $\Gamma$.

From the definition for the contrast Eq. (33), the variational problem (35) can be written as: find $E_\varepsilon \in \mathcal {V}$, such that

(54)$$\int_{\mathcal{D}}\left(\nabla\times E_\varepsilon\cdot\nabla\times W - k_0^{2}n^{2}E_\varepsilon\cdot W\right)\mathrm dx = S(W) + \int_{\mathcal{B}_\varepsilon} (\gamma^{2} - 1) k_0^{2} n^{2} E_\varepsilon \cdot W \,\mathrm dx \quad \forall \; W \in \mathcal{V}.$$

After subtracting Eq. (1) from Eq. (54) we obtain

(55)$$\int_{\mathcal{D}}\left(\nabla\times\widetilde{E}_\varepsilon\cdot\nabla\times W - k_0^{2}n^{2}\widetilde{E}_\varepsilon\cdot W\right)\mathrm dx = \int_{\mathcal{B}_\varepsilon} (\gamma^{2} - 1) k_0^{2} n^{2} E_\varepsilon\cdot W \,\mathrm dx.$$

By setting $W = \widetilde {E}_\varepsilon$ in Eq. (9) and $W = V$ in Eq. (55), we obtain the following equalities, respectively:

(56)$$\int_{\mathcal{D}}\left(\nabla\times V\cdot\nabla\times\widetilde{E}_\varepsilon - k_0^{2}n^{2}V\cdot\widetilde{E}_\varepsilon\right)\mathrm dx = \langle \varphi^{{\ast}}(E), \varphi(\widetilde{E}_\varepsilon) \rangle, $$

(57)$$\int_{\mathcal{D}}\left(\nabla\times\widetilde{E}_\varepsilon\cdot\nabla\times V - k_0^{2}n^{2}\widetilde{E}_\varepsilon\cdot V\right)\mathrm dx = \int_{\mathcal{B}_\varepsilon}(\gamma^{2} - 1)k_0^{2}n^{2}E_\varepsilon\cdot V\,\mathrm dx. $$

From the symmetry of both bilinear forms and after taking the real part of the above equalities, the following important result holds true:

(58)$$\begin{aligned} 2 \Re\{\langle \varphi^{{\ast}}(E), \varphi(\widetilde{E}_\varepsilon)\rangle \} &= 2 \Re\{\int_{\mathcal{B}_\varepsilon}(\gamma^{2} - 1)k_0^{2}n^{2}E_\varepsilon\cdot V\,\mathrm dx \}\\ &= 2 \Re\{\int_{\mathcal{B}_\varepsilon}(\gamma^{2} - 1)k_0^{2}n^{2}E\cdot V\,\mathrm dx \} + \mathcal{E}_2(\varepsilon), \end{aligned}$$

where we have considered (51). The remainder $\mathcal {E}_2(\varepsilon )$ is defined as

(59)$$\mathcal{E}_2(\varepsilon) := 2 \Re\{\int_{\mathcal{B}_\varepsilon}(\gamma^{2} - 1)k_0^{2}n^{2}\widetilde{E}_\varepsilon\cdot V\,\mathrm dx\},$$

which is bounded as follows:

(60)$$|\mathcal{E}_2(\varepsilon)| \leq C \varepsilon^{\frac{d}{2}} \| \widetilde{E}_\varepsilon \|_{L^{2}(\mathcal{B}_\varepsilon;\mathbb{C}^{d})} \leq C \varepsilon^{\frac{d}{2}} \| \widetilde{E}_\varepsilon \|_{H_{\mathrm{curl}}(\mathcal{D};\mathbb{C}^{d})} = o(\varepsilon^{d}),$$

where we have used the Cauchy-Schwartz inequality and Lemma 3. From the interior elliptic regularity [65] of functions $E$ and $V$, we obtain

(61)$$2 \Re\{\langle \varphi^{{\ast}}(E), \varphi(\widetilde{E}_\varepsilon)\rangle \} = 2 |\mathcal{B}_\varepsilon(x_0)| \, \Re\{(\gamma^{2}(x_0) - 1) k_0^{2} n^{2}(x_0) E(x_0) \cdot V(x_0)\} + \sum_{i=2}^{3} \mathcal{E}_i(\varepsilon).$$

The remainder $\mathcal {E}_3(\varepsilon )$ is defined as follows:

(62)$$\mathcal{E}_3(\varepsilon):= 2 \Re \{\int_{\mathcal{B}_\varepsilon}[(\gamma^{2} - 1)k_0^{2} n^{2}E \cdot V - (\gamma^{2}(x_0) - 1) k_0^{2} n^{2}(x_0) E(x_0) \cdot V(x_0)] \mathrm dx \},$$

which is trivially bounded as

(63)$$|\mathcal{E}_3(\varepsilon)| \leq C \int_{\mathcal{B}_\varepsilon} \| x-x_0 \| = O(\varepsilon^{d+1}) = o(\varepsilon^{d}),$$

since all arguments in the integrand are Lipschitz continuous almost everywhere in $\Omega$. After collecting all terms from Eqs. (51), (58), and (61), we obtain

(64)$$\psi(E_\varepsilon) - \psi(E) = 2 |\mathcal{B}_\varepsilon(x_0)| \, \Re\{(\gamma^{2}(x_0) - 1)k_0^{2} n^{2}(x_0)E(x_0) \cdot V(x_0)\} + \mathcal{E}(\varepsilon),$$

where, according to Eq. (53), Eq. (60) and Eq. (63), $\mathcal {E}(\varepsilon ) = \sum _{i=1}^{3} \mathcal {E}_i(\varepsilon ) = o(\varepsilon ^{d})$.

Finally, the leading term of the expansion Eq. (64) can be identified as the topological derivative of the shape functional $\psi (E)$, which concludes the proof of Theorem 1.

Appendix B. Computational resources

As an indication of the computational resources required to perform the optimizations presented in the previous examples, Table 1 presents the total computational time spent in each of them as well as the number of iterations performed. The stopping criteria for all examples was not defined in terms of their objective function; the optimizations were allowed to continue until no further improvements could be found and the optimization mesh could not be further refined due to the memory limits in the computers. It is important to notice that, as a consequence, the final iterations are much slower than the initial ones.

Table 1. Computational resources for the optimizations. LS and FEM are the number of elements in the optimization regions for the level-set function and the simulation mesh, respectively.

View Table

Funding

Fundação de Amparo à Pesquisa do Estado de São Paulo (2015/24517-8, 2016/19270-6, 2018/25339-4); Fundação Carlos Chagas Filho de Amparo à Pesquisa do Estado do Rio de Janeiro (E-26/203.041/2017); Conselho Nacional de Desenvolvimento Científico e Tecnológico (302036/2018-0, 310512/2017-4, 408274/2018-2, 438272/2018-8); Coordenação de Aperfeiçoamento de Pessoal de Nível Superior (88881.311020/2018).

Disclosures

The authors declare no conflicts of interest.

References

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

2. M. P. Bendsøe and O. Sigmund, Topology optimization. Theory, methods and applications (Springer-Verlag, Berlin, 2003).

3. J. S. Jensen and O. Sigmund, “Topology Optimization for Nano-Photonics,” Laser Photonics Rev. 5(2), 308–321 (2011). [CrossRef]  

4. S. Molesky, Z. Lin, A. Y. Piggott, W. Jin, J. Vucković, and A. W. Rodriguez, “Inverse design in nanophotonics,” Nat. Photonics 12(11), 659–670 (2018). [CrossRef]  

5. O. Sigmund and J. Jensen, “Systematic design of phononic band–gap materials and structures by topology optimization,” Philos. Transactions Royal Soc. London. Ser. A: Math. Phys. Eng. Sci. 361(1806), 1001–1019 (2003). [CrossRef]  

6. P. I. Borel, A. Harpøth, L. H. Frandsen, M. Kristensen, P. Shi, J. S. Jensen, and O. Sigmund, “Topology optimization and fabrication of photonic crystal structures,” Opt. Express 12(9), 1996–2001 (2004). [CrossRef]  

7. L. Frandsen, A. Harpøth, P. Borel, M. Kristensen, J. Jensen, and O. Sigmund, “Broadband photonic crystal waveguide 60° bend obtained utilizing topology optimization,” Opt. Express 12(24), 5916–5921 (2004). [CrossRef]  

8. J. S. Jensen and O. Sigmund, “Systematic design of photonic crystal structures using topology optimization: Low-loss waveguide bends,” Appl. Phys. Lett. 84(12), 2022–2024 (2004). [CrossRef]  

9. J. S. Jensen and O. Sigmund, “Topology optimization of photonic crystal structures: A high-bandwidth low-loss T-junction waveguide,” J. Opt. Soc. Am. B 22(6), 1191–1198 (2005). [CrossRef]  

10. A. Têtu, M. Kristensen, L. Frandsen, A. Harpøth, P. Borel, J. Jensen, and O. Sigmund, “Broadband topology-optimized photonic crystal components for both te and tm polarizations,” Opt. Express 13(21), 8606–8611 (2005). [CrossRef]  

11. Z. Yu, H. Cui, and X. Sun, “Genetic-algorithm-optimized wideband on-chip polarization rotator with an ultrasmall footprint,” Opt. Lett. 42(16), 3093–3096 (2017). [CrossRef]  

12. J. C. Mak, C. Sideris, J. Jeong, A. Hajimiri, and J. K. Poon, “Binary particle swarm optimized 2×2 power splitters in a standard foundry silicon photonic platform,” Opt. Lett. 41(16), 3868–3871 (2016). [CrossRef]  

13. Y. Zhang, S. Yang, A. E.-J. Lim, G.-Q. Lo, C. Galland, T. Baehr-Jones, and M. Hochberg, “A cmos-compatible, low-loss, and low-crosstalk silicon waveguide crossing,” IEEE Photonics Technol. Lett. 25(5), 422–425 (2013). [CrossRef]  

14. D. Liu, L. H. Gabrielli, M. Lipson, and S. G. Johnson, “Transformation Inverse Design,” Opt. Express 21(12), 14223 (2013). [CrossRef]  

15. Y. Kim, S.-Y. Lee, J.-W. Ryu, I. Kim, J.-H. Han, H.-S. Tae, M. Choi, and B. Min, “Designing whispering gallery modes via transformation optics,” Nat. Photonics 10(10), 647–652 (2016). [CrossRef]  

16. L. Xu and H. Chen, “Conformal transformation optics,” Nat. Photonics 9(1), 15–23 (2015). [CrossRef]  

17. C. Y. Kao, S. Osher, and E. Yablonovitch, “Maximizing band gaps in two-dimensional photonic crystals by using level set methods,” Appl. Phys. B: Lasers Opt. 81(2-3), 235–244 (2005). [CrossRef]  

18. M. Burger and S. J. Osher, “A survey on level set methods for inverse problems and optimal design,” Eur. J. Appl. Math. 16(2), 263–301 (2005). [CrossRef]  

19. N. Lebbe, A. Glière, K. Hassan, C. Dapogny, and E. Oudet, “Shape optimization for the design of passive mid-infrared photonic components,” Opt. Quantum Electron. 51(5), 166 (2019). [CrossRef]  

20. G. Veronis, R. W. Dutton, and S. Fan, “Method for sensitivity analysis of photonic crystal devices,” Opt. Lett. 29(19), 2288–2290 (2004). [CrossRef]  

21. J. Lu and J. Vučković, “Nanophotonic computational design,” Opt. Express 21(11), 13351–13367 (2013). [CrossRef]  

22. S. Boutami and S. Fan, “Efficient pixel-by-pixel optimization of photonic devices utilizing the Dyson’s equation in a Green’s function formalism: Part I. Implementation with the method of discrete dipole approximation,” J. Opt. Soc. Am. B 36(9), 2378–2386 (2019). [CrossRef]  

23. S. Boutami and S. Fan, “Efficient pixel-by-pixel optimization of photonic devices utilizing the Dyson’s equation in a Green’s function formalism: Part II. Implementation using standard electromagnetic solvers,” J. Opt. Soc. Am. B 36(9), 2387–2394 (2019). [CrossRef]  

24. O. Sigmund, “Manufacturing tolerant topology optimization,” Acta Mech. Sinica 25(2), 227–239 (2009). [CrossRef]  

25. F. Wang, J. S. Jensen, and O. Sigmund, “Robust Topology Optimization of Photonic Crystal Waveguides with Tailored Dispersion Properties,” J. Opt. Soc. Am. B 28(3), 387 (2011). [CrossRef]  

26. M. Zhou, B. S. Lazarov, and O. Sigmund, “Topology Optimization for Optical Projection Lithography with Manufacturing Uncertainties,” Appl. Opt. 53(12), 2720 (2014). [CrossRef]  

27. A. Y. Piggott, J. Petykiewicz, L. Su, and J. Vuckovic, “Fabrication-constrained nanophotonic inverse design,” Sci. Rep. 7(1), 1786 (2017). [CrossRef]  

28. A. Michaels, M. C. Wu, and E. Yablonovitch, “Hierarchical Design and Optimization of Silicon Photonics,” IEEE J. Select. Topics Quantum Electron. 26(2), 1–12 (2020). [CrossRef]  

29. D. Vercruysse, N. V. Sapra, L. Su, R. Trivedi, and J. Vučković, “Analytical level set fabrication constraints for inverse design,” Sci. Rep. 9(1), 8999 (2019). [CrossRef]  

30. J. Sokołowski and A. Żochowski, “On the topological derivative in shape optimization,” SIAM J. Control Optim. 37(4), 1251–1272 (1999). [CrossRef]  

31. A. A. Novotny and J. Sokołowski, Topological derivatives in shape optimization, Interaction of Mechanics and Mathematics (Springer-Verlag, Berlin, Heidelberg, 2013).

32. F. L. Louër and M.-L. Rapún, “Topological sensitivity for solving inverse multiple scattering problems in three-dimensional electromagnetism. part i: One step method,” SIAM J. on Imaging Sci. 10(3), 1291–1321 (2017). [CrossRef]  

33. L. He, C.-Y. Kao, and S. Osher, “Incorporating topological derivatives into shape derivatives based level set methods,” J. Comput. Phys. 225(1), 891–909 (2007). [CrossRef]  

34. M. Hintermüller and A. Laurain, “Multiphase image segmentation and modulation recovery based on shape and topological sensitivity,” J. Math. Imaging Vis. 35(1), 1–22 (2009). [CrossRef]  

35. S. Amstutz, S. M. Giusti, A. A. Novotny, and E. A. de Souza Neto, “Topological derivative for multi-scale linear elasticity models applied to the synthesis of microstructures,” Int. J. for Numer. Methods Eng. 84, 733–756 (2010). [CrossRef]  

36. N. Van Goethem and A. A. Novotny, “Crack nucleation sensitivity analysis,” Math. Methods Appl. Sci. 33, 1978–1994 (2010). [CrossRef]  

37. G. Allaire, F. Jouve, and N. Van Goethem, “Damage and fracture evolution in brittle materials by shape optimization methods,” J. Comput. Phys. 230(12), 5010–5044 (2011). [CrossRef]  

38. A. A. Novotny, J. Sokołowski, and A. Żochowski, “Topological derivatives of shape functionals. Part I: Theory in singularly perturbed geometrical domains,” J. Optim. Theory Appl. 181(1), 1–22 (2019). [CrossRef]  

39. A. A. Novotny, J. Sokołowski, and A. Żochowski, “Topological derivatives of shape functionals. Part II: First order method and applications,” J. Optim. Theory Appl. 181(1), 1–22 (2019). [CrossRef]  

40. A. A. Novotny, J. Sokołowski, and A. Żochowski, “Topological derivatives of shape functionals. Part III: Second order method and applications,” J. Optim. Theory Appl. 181(1), 1–22 (2019). [CrossRef]  

41. G. Allaire, F. de Gournay, F. Jouve, and A. M. Toader, “Structural optimization using topological and shape sensitivity via a level set method,” Control and Cybernetics 34(1), 59–80 (2005).

42. M. Burger, B. Hackl, and W. Ring, “Incorporating topological derivatives into level set methods,” J. Comput. Phys. 194(1), 344–362 (2004). [CrossRef]  

43. H. Eschenauer, V. Kobelev, and A. Schumacher, “Bubble method for topology and shape optmization of structures,” Struct. Optim. 8(1), 42–51 (1994). [CrossRef]  

44. M. Otomori, T. Yamada, K. Izui, and S. Nishiwaki, “Matlab code for a level set-based topology optimization method using a reaction diffusion equation,” Struct. Multidiscip. Optim. 51(5), 1159–1172 (2015). [CrossRef]  

45. S. Amstutz and H. Andrä, “A new algorithm for topology optimization using a level-set method,” J. Comput. Phys. 216(2), 573–588 (2006). [CrossRef]  

46. S. Osher and J. A. Sethian, “Front propagating with curvature dependent speed: algorithms based on Hamilton-Jacobi formulations,” J. Comput. Phys. 79(1), 12–49 (1988). [CrossRef]  

47. J. K. Guest, J. H. Prévost, and T. Belytschko, “Achieving minimum length scale in topology optimization using nodal design variables and projection functions,” Int. J. for Numer. Methods Eng. 61, 238–254 (2004). [CrossRef]  

48. O. Sigmund, “Morphology-based black and white filters for topology optimization,” Struct. Multidiscip. Optim. 33(4-5), 401–424 (2007). [CrossRef]  

49. B. S. Lazarov and O. Sigmund, “Filters in topology optimization based on Helmholtz-type differential equations,” Int. J. for Numer. Methods Eng. 86, 765–781 (2011). [CrossRef]  

50. F. Wang, B. S. Lazarov, and O. Sigmund, “On projection methods, convergence and robust formulations in topology optimization,” Struct. Multidiscip. Optim. 43(6), 767–784 (2011). [CrossRef]  

51. J. D. Deaton and R. V. Grandhi, “A survey of structural and multidisciplinary continuum topology optimization: post 2000,” Struct. Multidiscip. Optim. 49(1), 1–38 (2014). [CrossRef]  

52. M. Zhou, B. S. Lazarov, F. Wang, and O. Sigmund, “Minimum length scale in topology optimization by geometric constraints,” Comput. Methods Appl. Mech. Eng. 293, 266–282 (2015). [CrossRef]  

53. A. Bossavit, Computational Electromagnetism, Electromagnetism (Academic Press, San Diego, 1998).

54. P. Monk, Finite Element Methods for Maxwell’s Equations, Numerical Mathematics and Scientific Computation (Oxford University Press, Oxford, UK, 2003).

55. A. Quarteroni and A. Valli, Numerical Approximation of Partial Differential Equations (Springer-Verlag, Berlin, Heidelberg, 1997).

56. A. A. Novotny, J. Sokołowski, and A. Żochowski, Applications of the topological derivative method, Studies in Systems, Decision and Control (Springer, Nature Switzerland, 2019).

57. M. S. Alnæs, J. Blechta, J. Hake, A. Johansson, B. Kehlet, A. Logg, C. Richardson, J. Ring, M. E. Rognes, and G. N. Wells, “The fenics project version 1.5,” Arch. Numer. Softw. 3(100), 9–23 (2015). [CrossRef]  

58. S. Amstutz, “Analysis of a level set method for topology optimization,” Optim. Methods Softw. 26(4-5), 555–573 (2011). [CrossRef]  

59. E. R. Dougherty and R. A. Lotufo, Hands-on Morphological Image Processing (SPIE, 1000 20th Street, Bellingham, WA 98227-0010 USA, 2003).

60. M. Mrozowski, Guided Electromagnetic Waves: Properties and Analysis, no. 3 in Electronic and Electrical Engineering Research Studies: Computer Methods in Electromagnetics (Research Studies Press, 1997).

61. D. M. Pozar, Microwave Engineering (Wiley, 2011), 4th ed.

62. C. A. Balanis, Advanced Engineering Electromagnetics (Wiley, 2012), 2nd ed.

63. A. F. Oskooi, D. Roundy, M. Ibanescu, P. Bermel, J. Joannopoulos, and S. G. Johnson, “Meep: A Flexible Free-Software Package for Electromagnetic Simulations by the FDTD Method,” Comput. Phys. Commun. 181(3), 687–702 (2010). [CrossRef]  

64. A. Y. Piggott, J. Lu, K. G. Lagoudakis, J. Petykiewicz, T. M. Babinec, and J. Vuckovic, “Inverse Design and Demonstration of a Compact and Broadband On-Chip Wavelength Demultiplexer,” Nat. Photonics 9(6), 374–377 (2015). [CrossRef]  

65. R. Dautray and J. L. Lions, Mathematical analysis and numerical methods for science and technology. Volume 2 – Functional and Variational Methods (Springer, Berlin, Heidelberg, New York, 1988).