1. INTRODUCTION

For on-chip photonic device optimization, there are a few important constraints that need to be satisfied. First, these devices are typically constructed out of a discrete set of materials with a corresponding discrete set of permittivities, and the final optimized design must use only this set of permittivities. In many cases, this set in fact consists of two materials, and the resulting design must be a binary design [1]. Second, the design must be amenable to nanofabrication technology; in particular, the design needs to satisfy the minimum feature size requirement.

While continuous optimization methods, in particular those based on the adjoint variable method [2–8], have been widely used for the optimization of photonic devices, it remains a considerable challenge to incorporate the constraints mentioned above, in a systematic fashion, while maintaining computational efficiency and high performance of the final design [9]. Discrete optimization methods, e.g., a pixel-by-pixel optimization strategy, can directly incorporate these constraints from the outset, and in fact can be implemented such that the search occurs only in the parameter space that satisfies these constraints; however, the computational cost of pixel-by-pixel optimization has been high [10,11]. In particular, in traditional implementation of pixel-by-pixel optimization, where one updates a single pixel in each iteration, one needs to simulate the performance of all single-pixel variations in order to accurately decide the optimal pixel to update in each iteration, leading to substantial computational costs.

In part I of this series of papers [12], we have introduced a method, based on Dyson’s equation [13] in a Green’s function formalism, to overcome the challenge associated with discrete pixel-by-pixel optimization. We show that starting from the Green’s function of a device configuration, the change in the figure of merit (FOM) for all single-pixel variations can be efficiently computed. This then allows us to update the single pixel that leads to greatest improvement of the FOM to arrive at a new device configuration. Moreover, the Green’s function of the new device configuration can also be efficiently determined. Thus, the process can continue iteratively until an optimal configuration is reached. In part I [12], we provide an implementation of this method using the discrete dipole approximation (DDA) [13–15] as the simulation tool.

There are several drawbacks in the implementation of part I [12]. First of all, the implementation requires storage of the Green’s function. Thus the memory requirement for storage could be substantial: the Green’s function in real space is fully dense. Thus, the memory requirement scales as $N^2$, where $N$ is the number of pixels in the design region. This can become prohibitive if the design region is large. Second, many standard electromagnetic (EM) solvers that are commercially available are based on methods such as finite-difference time-domain (FDTD) [16,17], or iterative finite-element methods [18,19]. These methods efficiently compute the field distribution for a given source configuration in a structure, but usually do not compute the Green’s function. To use these methods to compute the Green’s function directly would require the calculation for a large number of source configurations, which can become very expensive.

In this paper, we show that the basic idea for efficient pixel-by-pixel optimization, as outlined in part I, can alternatively be implemented with the use of standard commercially available electromagnetic solvers, with a memory requirement that scales linearly with $N$. For this purpose, we note that to determine the position of the pixel, the modification of which will lead to the best improvement in the FOM, one requires only the diagonal part of the Green’s function, the storage of which scales linearly with $N$. The diagonal part of the Green’s function, in turn, can be updated with only a few full simulations of the structures. This idea can be generalized for an efficient block-by-block optimization, and therefore enables the direct incorporation of fabrication constraints.

The rest of the paper is organized as follows: in Section 2, we discuss the algorithms. In Section 3, we illustrate the efficiency of this algorithm with a few optimization examples. We conclude in Section 4.

2. MATHEMATICAL BACKGROUND

A. Single Pixel Modification

The same as in part I [12], we consider an electric field $\mathit{E}(\mathit{r})$ generated by an excitation for a structure as described by a given permittivity $\epsilon (\mathit{r})$. For illustration purposes, we suppose we would like to optimize a simple figure of merit that consists of the intensity of the electric field at a target point ${\mathit{r}}_{\mathbf{0}}$:

Define the change in permittivity distribution $\mathrm{\Delta }\epsilon \equiv {\epsilon }^{(\text{new})}-{\epsilon }^{(\text{old})}$ between two iterations. Using the Lippmann–Schwinger equation [15], we can determine the new electric field in the whole space:

In the case where we consider a single-pixel modification in each iteration, $\mathrm{\Delta }\epsilon$ is non-zero only in a small volume V centered at ${\mathit{r}}^{\prime \prime }$, and the change in FOM consists of

As we show in part I [12], the new electric field at ${\mathit{r}}^{\prime \prime }$ after a pixel has been updated in ${\mathit{r}}^{\prime \prime }$ can be efficiently computed, provided that the old Green’s function at this location is known, according to the Lippmann–Schwinger equation [12,13,15] applied in ${\mathit{r}}^{\prime \prime }$:

In the DDA-type approach described in part I [12], for the “new” structure where the dielectric constant is modified at the best-improvement pixel, the location of which is denoted as ${\mathit{r}}_{\mathit{b}}$, the field was updated everywhere using the Lippmann–Schwinger equation:

In Ref. [12], when implementing the algorithms above, we stored the full Green’s functions $\mathit{G}({\mathit{r}}^{\prime \prime },{\mathit{r}}^{\prime })$ in each iteration process. When using a numerical solver based on the discrete dipole approximation, this approach is natural. On the other hand, for most other algorithms, this is costly, since $\mathit{G}({\mathit{r}}^{\prime \prime },{\mathit{r}}^{\prime })$is fully dense. Thus, it is useful to develop a scheme in which one needs not store the full Green’s function. From Eqs. (7) and (8), it is clear that only the diagonal part of the Green’s function ${\mathit{G}}^{(\text{old})}({\mathit{r}}^{\prime \prime },{\mathit{r}}^{\prime \prime })$, where ${\mathit{r}}^{\prime \prime }$ can vary over the entire design region, is needed in the iteration process in order to determine the best improvement pixel. The required storage for ${\mathit{G}}^{(\text{old})}({\mathit{r}}^{\prime \prime },{\mathit{r}}^{\prime \prime })$ scales as $N$, where $N$ is the number of pixels in the design region. Also, instead of using Eq. (9), which involves the full Green’s function $\mathit{G}(\mathit{r},{\mathit{r}}^{\prime \prime })$ to update the field, the field for the new structure can alternatively be determined by using a standard EM solver, without the need of knowing the Green’s function of the old structure. Moreover, within each iteration, the diagonal part of the Green’s function can be updated for the new structure. Indeed, from Eq. (10), we have

In practice, the number of simulations required to determine ${\mathit{G}}^{(\text{old})}(\mathit{r},{\mathit{r}}_{\mathit{b}})$ and ${\mathit{G}}^{(\text{new})}(\mathit{r},{\mathit{r}}_{\mathit{b}})$ can be significantly reduced, by taking advantage of the simulations of the fields that have been already performed.

To start, we show a more efficient way to determine ${\mathit{G}}^{(\text{old})}(\mathit{r},{\mathit{r}}_{\mathit{b}})$.

In the case of transverse electric (TE) polarization in two dimensions, ${\mathit{G}}^{(\text{old})}(\mathit{r},{\mathit{r}}_{\mathit{b}})$ is scalar, and therefore Eq. (9) is sufficient to determine it. The required information is: ${\mathit{E}}^{(\text{old})}(\mathit{r})$ before the best-improvement pixel in ${\mathit{r}}_{\mathit{b}}$ is updated, and ${\mathit{E}}^{(\text{new})}(\mathit{r})$ and ${\mathit{E}}^{(\text{new})}({\mathit{r}}_{\mathit{b}})$ after the best-improvement pixel is updated. All of this information has already been calculated. Thus, ${\mathit{G}}^{(\text{old})}(\mathit{r},{\mathit{r}}_{\mathit{b}})$ can be straightforwardly determined by solving a simple linear equation:

In the case of the transverse magnetic (TM) polarization in two dimensions, the Green’s function is a $2\times 2$ tensor. Utilizing both ${\mathit{E}}^{(\text{new})}(\mathit{r})$ as well as ${\mathit{E}}_{\mathit{adj}}^{(\text{new})}(\mathit{r})$, which will be needed in order to compute the gradient of the FOM in Eq. (7), we have

Finally, for 3D problems, where electric fields are three dimensional and hence ${\mathit{G}}^{(\text{old})}(\mathit{r},{\mathit{r}}_{\mathit{b}})$ is a $3\times 3$ tensor, one has to run a third simulation (that we will denote with a subscript Green), in addition to forward and adjoint simulations. This simulation can be done with any excitation source as long as such an excitation source is linearly independent of the sources used in the simulations for ${\mathit{E}}^{(\text{new})}(\mathit{r})$ and ${\mathit{E}}_{\mathit{adj}}^{(\text{new})}(\mathit{r})$. We can use a dipole source located in ${\mathit{r}}_{\mathit{b}}$, and with the dipole moment orthogonal to both ${\mathit{E}}^{(\text{old})}({\mathit{r}}_{\mathit{b}})$ and ${\mathit{E}}_{\mathit{adj}}^{(\text{old})}({\mathit{r}}_{\mathit{b}})$, which provides a linear system with good conditioning [15]. The Lippmann–Schwinger equations for the three new fields then read

We now switch to the determination of ${\mathit{G}}^{(\text{new})}(\mathit{r},{\mathit{r}}_{\mathit{b}})$. We note that Eqs. (14)–(16) are all in the form of the Lippmann–Schwinger equation, which has the following properties between the “old” and the “new” fields: one can simply switch the subscript “new”, and “old,” and replace $\mathrm{\Delta }\epsilon ({\mathit{r}}_{\mathit{b}})$ with $-\mathrm{\Delta }\epsilon ({\mathit{r}}_{\mathit{b}})$, to arrive at a new set of equations that involve ${\mathit{G}}^{(\text{new})}(\mathit{r},{\mathit{r}}_{\mathit{b}})$. As a result, Eq. (14) transforms to

We therefore arrive at the following algorithm: at the start of each iteration, we assume that we know the field distribution ${\mathit{E}}^{(\text{old})}(\mathit{r})$, the adjoint field distribution ${\mathit{E}}_{\mathit{adj}}^{(\text{old})}(\mathit{r})$, as well as the diagonal part of the Green’s function ${\mathit{G}}^{(\text{old})}({\mathit{r}}^{\prime },{\mathit{r}}^{\prime })$, where ${\mathit{r}}^{\prime }$ takes the value over the design region. We use Eqs. (7) and (8) to determine the best-improvement pixel at which the modification of the structure leads to the best improvement of the FOM. We then update the structure by modifying the best-improvement pixel, arriving at the new structure. Afterwards we determine the field distributions ${\mathit{E}}^{(\text{new})}(\mathit{r})$ and ${\mathit{E}}_{\mathit{adj}}^{(\text{new})}(\mathit{r})$ of the new structure using a standard electromagnetic solver, and determine the diagonal part of the Green’s function ${\mathit{G}}^{(\text{new})}({\mathit{r}}^{\prime },{\mathit{r}}^{\prime })$ using Eqs. (12)–(19), which may involve using the standard solver for one additional field excited by a source placed at the best-improvement pixel in the 3D case. This completes a single iteration. This process iterates until the FOM no longer improves, at which point the algorithm converges to an optimal solution. In this algorithm, the storage scales linearly as the number of pixels of the system. The main computational costs for each iteration are two or three full simulations of the fields in each iteration with a standard electromagnetic solver.

In the algorithm above, we assume that at the initial step of the iteration, we have the information about the diagonal part of the Green’s function $\mathit{G}(\mathit{r},\mathit{r})$ over the entire design region. In some cases, e.g., if the initial structure consists of only vacuum, $\mathit{G}(\mathit{r},\mathit{r})$ is known analytically. On the other hand, there are common situations where the Green’s function for the initial structure is not known analytically. Instead of numerically determining $\mathit{G}(\mathit{r},\mathit{r})$ for the initial structure, which can be costly with many standard EM solvers, we can adopt the following procedure where we build up the knowledge about $\mathit{G}(\mathit{r},\mathit{r})$ as the iteration progresses. Starting from an initial structure labeled “0,” we can then randomly flip a pixel at $r_1$ to arrive at the structure labeled “1.” Using Eqs. (17)–(19), part of the Green’s function for structure 1, ${\mathit{G}}^{(1)}({\mathit{r}}_1,{\mathit{r}}_1)$ can be determined using the information of the fields for structures 0 and 1, without the knowledge of ${\mathit{G}}^{(0)}({\mathit{r}}_1,{\mathit{r}}_1)$. At this point, we approximate ${\mathit{G}}^{(1)}(\mathit{r},\mathit{r})\approx {\mathit{G}}^{(1)}({\mathit{r}}_1,{\mathit{r}}_1)$, and proceed to determine the best-improvement pixel ${\mathit{r}}_2$ and generate the structure labeled “2.” We note that ${\mathit{G}}^{(2)}({\mathit{r}}_2,{\mathit{r}}_2)$ can be determined rigorously in the same process as described above using Eqs. (17)–(19), whereas ${\mathit{G}}^{(2)}({\mathit{r}}_1,{\mathit{r}}_1)$ can be determined accurately from ${\mathit{G}}^{(1)}({\mathit{r}}_1,{\mathit{r}}_1)$ using Eq. (12). We then construct an approximation of ${\mathit{G}}^{(2)}(\mathit{r},\mathit{r})$ from ${\mathit{G}}^{(2)}({\mathit{r}}_1,{\mathit{r}}_1)$ and ${\mathit{G}}^{(2)}({\mathit{r}}_2,{\mathit{r}}_2)$. We repeat this process as outlined above in the iterations. As the iterations continue, the approximation to $\mathit{G}(\mathit{r},\mathit{r})$ becomes increasingly accurate, since $\mathit{G}(\mathit{r},\mathit{r})$ at the best-improvement pixels at all previous iterations are accurately determined. We also note that the accurate information of $\mathit{G}(\mathit{r},\mathit{r})$ is not crucial for the beginning part of the iteration process, but becomes increasingly important as the iteration progresses. Thus, the algorithm above allows us to progressively build upon the information about $\mathit{G}(\mathit{r},\mathit{r})$ without compromising the optimization process.

Here, for illustration purposes, we focus on a simple FOM as described in Eq. (1). The generalization to other types of FOMs has been described in Ref. [12], which involves modification of the adjoint simulation. The procedure for updating the diagonal part of the Green’s function, as described above, is unchanged.

B. Single-Block Modification

As already stated in Ref. [12], instead of updating a single pixel in each iteration, one can generalize the algorithm above to describe the case of single-block modification, where each block contains multiple pixels, and in the optimization process, a single block is modified in every step. Such single-block modification can be used to handle technological constraints, such as minimum feature size requirement, when the minimum feature size exceeds the size of one pixel (we recall that the size of the pixel needs to be sufficiently small as compared to the wavelength so that the discretization of the field is accurate). It is also useful in a 3D situation for the design of photonic circuits in a guiding layer. In this case, a block consists of multiple pixels within, and stacked in the direction perpendicular to, the guiding layer.

In this case, the derivation largely parallels the case of single-pixel modification. Equation (7) becomes

The electric field of the new structure can be determined using a standard electromagnetic solver. The remaining task to complete the iteration is to update the block-diagonal part of the Green’s function $\mathit{G}({\mathit{r}}_2,{\mathit{r}}_1)$. This is done by noting that updating a single block is equivalent to updating sequentially the pixels inside the block. Thus, we can update $\mathit{G}({\mathit{r}}_2,{\mathit{r}}_1)$ in a pixel-by-pixel fashion. Denoting the pixel being updated as ${\mathit{r}}_{\mathit{b}}$, the block-diagonal part of the Green’s function is updated using

3. DESIGN EXAMPLES

We now provide examples of implementation for the algorithm as described above. We use a commercial code RSOFT-Fullwave [20], which is based on the FDTD method. The FDTD method is one of the most popular EM methods: it is versatile, able to deal with arbitrary geometries and many different kinds of materials, and it is also very memory efficient in 3D, as its memory requirement scales linearly with the system size. We implement the optimization process as a code in MATLAB. At each iteration, the optimization code communicates with RSOFT-Fullwave to get the updated electromagnetic fields, and provides the new permittivity distribution to RSOFT-Fullwave. Communications between the two codes are done by exchanging text data files [20].

A. 2D TE Example in Homogeneous Background

We start by reproducing the results of part I of the paper. We optimize in 2D a device that maximizes light intensity at a target point excited from the source at a distance of one wavelength away [Fig. 1(a)]. The pixel size is $50\mathrm{nm}\times 50\mathrm{nm}$. We consider the TE polarization where the electric field is perpendicular to the 2D plane.

 

Fig. 1. 2D TE example in homogeneous background. (a) Schematic of the setup of the optimization, which aims to maximize light intensity at a target point, separated from an emitting source by a distance of one wavelength ($\lambda =1.55\mathrm{\mu m}$). We consider a binary system with two permittivities of 2.25 and 6.25, respectively. The design space is a $2\lambda$-sized square, divided in 50-nm-sized square pixels, with the source at the center. (b) Optimized structure obtained in part I [12] with the DDA approach. (c) FOM as a function of iteration number. (d) Optimized structure. (c)–(d) are obtained using the approach described in this paper. (e) Electric field of the optimized structure, as determined using FDTD simulation.

Download Full Size | PDF

In the optimization process, we start from a silica homogeneous background, for which the Green’s function can be analytically obtained [15]. At each iteration, we follow the algorithms as discussed in Section 2.A, with the required fields determined from FDTD simulations.

We reproduce in Fig. 1(b) the optimized structure as obtained using the DDA approach and presented in part I. In Fig. 1(c), we show the FOM evolution as a function of iteration number, and in Figs. 1(d) and 1(e), we show the optimized structure and the field distribution of the optimized structure. The FOM of the optimized structure, as well as the number of iterations required to reach convergence, is very similar to those from the DDA approach in Ref. [12]. The final optimized structure [Fig. 1(d)] is also very similar to the structure obtained using the DDA approach [Fig. 1(b)]. However, the memory required by the optimization code and the RSOFT code [$90*90$ 50-nm-sized square pixels including the perfectly matched layers (PML)s] together is only 257 kB, compared to 236 MB required for the DDA approach in Ref. [12], which represents a significant decrease in memory usage. In the total time spent, the optimization code accounts for only 0.38% and is thus negligible. Almost all simulation times are spent on FDTD simulations.

B. 2D TE Example in Complex Background

We then reproduce the results obtained in part I for the 2D TE device design problem with complex background [12]. The device operates as a two-wavelength demultiplexer, which separates two incident waves with wavelengths ${\lambda }_1=1.31\mathrm{\mu m}$ and ${\lambda }_2=1.55\mathrm{\mu m}$ from the fundamental mode of a 1-μm-wide input multimode waveguide, to two 0.3-μm-wide single-mode output guides [Fig. 2(a)]. The FOM is the average of power coupling efficiencies to the two output waveguides [12].

 

Fig. 2. 2D TE example in complex background. (a) Schematic of the setup of the optimization: light incident as the fundamental mode of a 1-μm-wide multimode input waveguide is routed to two 0.3-μm-wide single-mode output waveguides in a wavelength-dependent fashion. Light at wavelength ${\lambda }_1=1.31\mathrm{\mu m}$ is routed to the left, and light at the wavelength ${\lambda }_2=1.55\mathrm{\mu m}$ is routed to the right. (b) Optimized structure obtained in part I [12] with the DDA approach. (c) FOM as a function of iteration number, starting from a uniform silica background in the design region. (d) Optimized structure. (c)–(d) are obtained using the approach described in this paper. (e) Corresponding electric fields at wavelengths ${\lambda }_1=1.31\mathrm{\mu m}$ (left panel) and ${\lambda }_2=1.55\mathrm{\mu m}$ (right panel), as obtained using FDTD simulations of the optimized structure.

Download Full Size | PDF

In this structure, the initial structure consists of the input and output waveguides, with the design region being uniform. The Green’s function of this structure is not known analytically. In contrast to the DDA-based method in part I, for which the Green’s function needed to be initialized, requiring time-costly pre-calculations or matrix inversion [12], we use the method discussed towards the end of Section 2.A, where we approximate the diagonal part of the Green’s function and gradually improve the approximation as the iteration progresses. The final structure produced in this process [Fig. 2(d)] does operate as intended [Fig. 2(e)], and is very similar to the structure produced using the DDA approach [Fig. 2(b)], with a comparable FOM of 78% [Fig. 2(c)]. The small discrepancies between the structures in Figs. 2(b) and 2(d) are expected, since the Green’s function used in the two methods are not identical. The results here indicate that the approximate Green’s function used is sufficient to reach high-performing structures. The memory required by the optimization code and the RSOFT code ($110*110$ 50-nm-sized square pixels including PMLs) together is only 316 kB, compared to 41 MB required for the DDA approach.

C. 2D TM Example in Complex Background

Here we perform the same 2D problem as example 3.B, but for TM-polarized light where the electric field is an in-plane (Fig. 3). In this case, in each iteration, we use both forward and adjoint simulations to determine the diagonal part of the Green’s function at the best-improvement pixel, as described by Eq. (18) with $\mathit{r}={\mathit{r}}_{\mathit{b}}$. Then, the update of the diagonal part of the Green’s functions for all former best-improvement pixels positions is done with Eq. (12), using Eqs. (15) and (18). An approximate form of the diagonal part of the Green’s function is built in the process as described in Section 2.A. We confirm, here again, this process produces an optimal structure with a high FOM. The memory required by the optimization code and the RSOFT code ($110*110$ 50-nm-sized square pixels including PMLs) together is only 392 kB, compared to 164 MB required for the DDA approach.

 

Fig. 3. 2D TM example in complex background. Same as in Fig. 2, but for TM polarization. (a) FOM as a function of iteration number, starting from a uniform silica background in the design region. (b) Optimized structure. (c) Corresponding magnetic fields at wavelengths ${\lambda }_1=1.31\mathrm{\mu m}$ (left panel) and ${\lambda }_2=1.55\mathrm{\mu m}$ (right panel), as obtained using FDTD simulations of the optimized structure.

Download Full Size | PDF

D. 3D Example for Optimization of a Slab Structure

We then perform a similar problem as example 3.A, but in a 3D configuration [Fig. 4(a)]. We consider a silicon [$\epsilon =12.25$) slab with 0.3 μm thickness, surrounded by silica ($\epsilon =2.25)$]. The aim of the optimization is to maximize light intensity emitted from a $y$-polarized source dipole, toward a target point at half-wavelength distance away. Both the source and target are located in the middle of the slab along the $y$ direction [Fig. 4(a)]. In the optimization process, we consider structures where various patterns are introduced into the slab. Since these patterns are generally fabricated by lithography and etching, we assume that the patterning results in structures that are uniform in $y$ within the slab. Therefore, we use the algorithms as discussed in Section 2.B, where each block contains a line of pixels aligned along the $y$ direction.

 

Fig. 4. 3D example for optimization of a slab structure. (a) Schematic of the setup of the optimization, which aims to maximize light intensity at a target point, separated from a $y$-polarized emitting source dipole by a distance of a half-wavelength ($\lambda =1.55\mathrm{\mu m}$), with the source at the center. Source and target point are located in the middle of the slab along the $y$ direction. We consider a binary system with two permittivities of 2.25 (silica) and 12.25 (silicon), respectively. The design space is a $1.55\mathrm{\mu m}\times 0.3\mathrm{\mu m}\times 1.55\mathrm{\mu m}$ slab, divided in 50-nm-sized cubic pixels. (b) FOM as a function of iteration number, starting from a uniform silica background. (c) Optimized structure, using the approach described in this paper. (d) Ey electric field of the optimized structure in the middle of the slab, as determined using FDTD simulation.

Download Full Size | PDF

For this optimization, we choose the initial structure consisting of uniform silica regions. As the initial background is homogenous, made of silica, one can initialize analytically the Green’s function. However, this initialization should be consistent with the way the Green’s function is updated with 3D-FDTD. In 3D-FDTD, as other 3D finite-difference methods, the electric field components are calculated at edges of the Yee cells [16,17]. Therefore, one has to analytically initialize the Green’s function correspondingly, in particular, the diagonal part, in a way different from the standard procedure defined in literature where the spatial coordinates of the Green’s function is assumed to be at the center of the cell [15,21].

Alternatively, it is also possible to initialize numerically the Green’s function for a homogeneous background. Here, we use Eq. (16) and we refer to the initial uniform silicon structure as the “old” structure, and the structure where we add a first pixel at ${\mathit{r}}_{\mathit{b}}$ of a given block as the “new” structure. Equation (16) can then be used to deduce ${\mathit{G}}^{(\text{old})}({\mathit{r}}_2,{\mathit{r}}_{\mathit{b}})$ for all ${\mathit{r}}_2$ within the same block as ${\mathit{r}}_{\mathit{b}}$. Since the initial structure is uniform, from this we can then deduce ${\mathit{G}}^{(\text{old})}({\mathit{r}}_2,{\mathit{r}}_1)$ for ${\mathit{r}}_1$ and ${\mathit{r}}_2$ belonging to the same block, which concludes the initialization process for the block-diagonal part of the Green’s function. This process is directly compatible with 3D-FDTD simulations and results in a rigorous initialization of the block-diagonal part of the Green’s function. After the initialization process, the iteration process can then proceed according to the description in Section 2.

For this optimization problem, the FOM as a function of iteration number is shown in Fig. 4(b). As an important characteristic of this class of methods we propose here, the FOM improves for each interaction until convergence. The cross section of the resulting structure is shown in Fig. 4(c). The structure can be fabricated by etching such a cross-sectional pattern into the slab, as intended for this optimization tool. The optimal structure shows strong field concentration at the target point [Fig. 4(d)].

The reduction in memory use obtained with our approach, as compared to the DDA stand-alone method of part I [12], is here again large. The memory required by the optimization code and the RSOFT code ($50*38*50$ 50-nm-sized cubic pixels including PMLs) together is only 9.5 MB, whereas the DDA approach would require 4.8 GB, for final results with similar performance. The additional time for the optimization algorithm accounts for 35.6% of the simulation time, due mainly to the additional Green’s simulation ($\sim 33\%$).

E. 3D Example of a Slab Structure in Complex Background with Minimum Feature Size Constraints

We finally end with a 3D slab structure example in a complex background, with minimum feature size constraints. The design domain is 1 μm wide and 0.3 μm thick, made of silicon and silica, optimized to convert ${\mathrm{TE}}_{00}$ mode into ${\mathrm{TE}}_{01}$ mode, for identical 1-μm-wide and 0.3-μm-thick input and output Si waveguides [Fig. 5(a)]. We consider three minimum feature sizes: 50 nm, 100 nm, and 200 nm, leading, respectively, to vertical blocks of six pixels, 24 pixels, and 96 pixels. The Green’s function is not analytical for this problem in a complex background, and is thus progressively determined according to the procedure described earlier. As shown in Fig. 5, as the minimum feature size increases, the device performance naturally decreases, but the method finds the optimal device, within a given minimum feature size constraint. The memory required by the optimization code and the RSOFT code ($50*38*50$ 50-nm-sized cubic pixels including PMLs) together is 6.6 MB, 12.9 MB, and 37.6 MB, respectively, for the 50 nm, 100 nm, and 200 nm minimum feature size cases, whereas it would be 829 MB for the DDA approach.

 

Fig. 5. 3D example of a slab structure in complex background with minimum feature size technological constraints. (a) Schematic of the setup of the optimization, which aims to maximize conversion of ${\mathrm{TE}}_{00}$ mode into ${\mathrm{TE}}_{01}$ mode, for 1-μm-wide 0.3-μm-thick input and output waveguides ($\lambda =1.55\mathrm{\mu m}$). We consider a binary system with two permittivities of 2.25 (silica) and 12.25 (silicon), respectively. The design space is a $1\mathrm{\mu m}\times 0.3\mathrm{\mu m}\times 1\mathrm{\mu m}$ slab, divided in 50-nm-sized cubic pixels. (b) Optimized structures, using the approach described in this paper, for minimum feature sizes of (i) 50 nm, (ii) 100 nm, and (iii) 200 nm. (c) Ex electric field of the optimized structures in the middle of the slab, as determined using FDTD simulations, with indicated conversion efficiencies.

Download Full Size | PDF

4. CONCLUSION

In part I [12] of the paper, we introduced an efficient pixel-by-pixel optimization technique, based on Dyson’s equation in a Green’s function formalism. This technique naturally incorporates material and minimum feature size constraints in practical devices. In this paper, we have taken a step further, by showing that this technique can be implemented with standard commercially available electromagnetic solvers, with a memory requirement that scales linearly with the number of pixels in the design region. Our work should lead to further development of discrete optimization techniques for design and implementation of nanophotonic devices.