1. INTRODUCTION

It was recently shown that systematic computational optimization can lead to non-intuitive photonic devices that are much more efficient and compact when compared to conventional devices [1–3]. This success, in turn, has motivated efforts to advance computational methods for photonic device optimization.

In general, for typical optimization methods, within each iteration, one starts by evaluating the performance of a structure. The structure is then modified in some fashion with the aim to improve the device performance. This process is repeated until a structure with a sufficiently good performance is found. In these optimization methods, therefore, developing a technique to determine how to modify the structure in an iteration is of central importance.

In addition to some of general considerations in optimizations, there are considerations specific to photonic device designs. First of all, most on-chip photonic devices are constructed from a discrete set of materials. In fact, many current non-intuitive devices are constructed by etching air voids into silicon [4–6]. These structures thus consist only of two materials, silicon and air for example. In this paper, we refer to a structure that consists only of two materials as a binary structure. For practical considerations, it is often desirable that the optimization procedure results in a binary structure. Second, on-chip photonic devices are typically fabricated with lithographic techniques, which requires that the minimum feature sizes of the devices must be larger than a certain threshold. Optimization methods for photonic devices needs to consider these practical considerations.

Among all the methods for photonic device optimization, the adjoint method, which is a continuous optimization method, has been widely used [2,3,7–12]. In this method, the gradient of the figure of merit (FOM) with respect to the continuous variation of the dielectric function in every pixel of the structure can be efficiently computed. One then modifies the structure based on the gradient. A key advantage of this method is that within each iteration, one simultaneously adjusts a large number of design parameters. Also, since the gradient points to the direction for the improvement of FOM, every iteration improves the device performance. However, the adjoint method typically modifies the dielectric function of the structure continuously. Thus, it may not result in an optimized structure that is binary, even if the starting structure is [3]. Neither does the method guarantee an optimized structure with a required minimum feature size.

There have been attempts seeking to modify adjoint methods to generate structures that satisfy the practical constraints. For example, one can seek to approximate the resulting optimized structure from an adjoint optimization procedure with a binary structure that satisfies the minimal feature size requirement [2,13]. Doing so, however, typically sacrifices some of the device performance. Alternatively, the adjoint method can be combined with a level set method [2,10,14], which guarantees a binary structure, but doing so limits the available search space [15,16].

Instead of using a continuous optimization method, one can use a pixel-based discrete optimization method to guarantee a binary structure with a sufficiently large minimum feature size. One could divide the structure into pixels, where each pixel has a size equal to the minimum feature size. To get a binary structure consisting of silicon and air, for example, one considers only structures where each pixel consists of either silicon or air. One of the difficulties here, however, is to design a strategy to modify the structure. Within each iteration, one could evaluate a few candidate structures before deciding how to modify the structure. This approach, known as the direct binary search [2,17], is however quite expensive computationally. Alternatively, one might consider using the adjoint method for optimization guidance. However the adjoint formulation, which mathematically provides a sensitivity to infinitesimal permittivity change in the system [8,9], is likely to make wrong predictions for binary optimization, in particular for a large refractive index contrast or large pixel size.

In this paper, we introduce a method to overcome the difficulty, as outlined above, for pixel-based discrete optimization method. We show that, once the Green’s function of an initial structure is determined, the change in the FOM, for every structure differing from the initial structure by a single pixel, can be determined efficiently. Thus, one can choose among this set of structures the best in terms of the FOM as the modified structure. Moreover, the Green’s function for the modified structure can also be updated efficiently, using Dyson equations [18]. This process can then be iterated until an optimal structure is found. This approach guarantees that every binary permittivity change improves the FOM. Also, the minimal feature size requirement is straightforwardly satisfied in this method. We will illustrate the method with numerical examples.

The rest of the paper is organized as follows: in the second section, we will demonstrate that a Green’s function formalism [18] allows for an efficient and guaranteed ascent of the FOM for pixeled-based discrete permittivity change, and we will combine it to an adjoint method for better efficiency in a general case. In the third section, we will show various device optimization examples, illustrating remarkable features of the method, including its high stability for highly resonant binary structures, and the straightforward and rigorous control of minimal feature size allowed by the method.

2. MATHEMATICAL BACKGROUND

Consider an electric field $\mathit{E}(\mathit{r})$ generated by an excitation ${\mathit{E}}_{\text{exc}}(\mathit{r})$, for a structure as described by a given permittivity $\mathit{\epsilon }(\mathit{r})$. For illustration of the main idea of the paper, suppose we would like to optimize a figure of merit:

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

To solve Eq. (3), one can discretize the volume $V$ in pixels, determine first the electric field ${\mathit{E}}^{(\text{new})}({\mathit{r}}^{\prime \prime })$ for each pixel of $V$ by solving a square linear system of equations, and then determine the electric field in the whole space by application of (3) outside volume $V$. This approach is commonly known as the discrete dipole approximation (DDA) [19–21]. The DDA method becomes exact when the pixel size is far smaller than the wavelength. This procedure can be used to determine the change in the figure of merit $\mathrm{\Delta }\mathrm{FOM}$ as described in Eq. (2).

However, in the case of single-pixel modification, $\mathrm{\Delta }\epsilon$ is non-zero only in a small volume $V$ centered at ${\mathit{r}}^{\prime \prime }$. It is then in fact possible to determine $\mathrm{\Delta }\mathrm{FOM}$ in a much simpler way. To see this, by combining Eq. (3) with Eq. (2), we have

We can easily determine ${\mathit{E}}^{(\text{new})}({\mathit{r}}^{\prime \prime })$ for every new structure that differs from the old structure by a single pixel modification at ${\mathit{r}}^{\prime \prime }$. Indeed, from Eq. (3), ${\mathit{E}}^{(\text{new})}({\mathit{r}}^{\prime \prime })$ can actually be determined straightforwardly and efficiently, by solving the simple following Lippman–Schwinger equation:

Therefore, we have shown that, if we know the Green’s function of a given (“old”) structure, from Eqs. (4) and (5), one can determine the change of the FOM for all single pixel modification in an efficient and rigorous way, valid even for high-permittivity change such as binary change. Using this information, among all the structures that differ from the old structure by a single-pixel modification, we can then select the structure that offers the best improvement in the FOM as our new candidate (“new”) structure. This is a major improvement in comparison with traditional optimization techniques, such as adjoint method with standard EM methods, that perform a first-order derivative perturbation with the only available old field ${\mathit{E}}^{(\text{old})}({\mathit{r}}^{\prime \prime })$ in Eq. (4) and that therefore require a continuous permittivity change. In the following then, we use the same ${\mathit{r}}^{\prime \prime }$ to denote the pixel where we have carried out the modification of the structure.

To complete the iterative process, we will also need to compute the field in the new structure, as well as its Green’s function. The field in the new structure can be straightforwardly calculated as

We note that solving an electromagnetic scattering problem using Eqs. (7) and (8), by adding sequentially the different pixels constituting a device, was proposed in Ref. [22]. Our contribution here is to point out the effectiveness of this approach in the context of optimization.

In terms of the computational time, if our design region consists of $N$ pixels, after simulating an initial structure, at each iteration, the cost of determining and solving for the new structure is far smaller than the cost for solving the initial structure. This is in contrast with traditional approaches in discrete optimization, where in each iteration one needs to fully solve $N$ structures to determine the new structure.

In terms of the memory requirement, in the simplest implementation of Eqs. (5)–(8), one needs to store the full Green’s function of the structure. This is doable for a modest-size structure, but it can become expensive when the structure becomes larger. In this paper, in the numerical demonstration of the next section, we do store the full Green’s function. In a subsequent paper [23], we show that there exists a more efficient way of storing the Green’s function.

Above, for illustration purposes, we have focused on a FOM that depends only on the field value at a single spatial point. In practice, other forms of FOM may be used. As an example, later in this paper, we will optimize the power coupling efficiency into a particular waveguide mode, as described by the modal field pattern $|\begin{array}{c}{\mathit{E}}_{\mathit{m}}\\ {\mathit{H}}_{\mathit{m}}\end{array}$, for which case the FOM has the form

Equation (9) is of the general form,

With a single pixel variation, the change in FOM thus becomes

We end this section by commenting on the connection and the contrast of the present method with the standard adjoint variable method. When $\mathrm{\Delta }\epsilon$ is small, the change of FOM can be obtained with either Eq. (4) or Eq. (13), but replacing ${\mathit{E}}^{(\text{new})}$ with ${\mathit{E}}^{(\text{old})}$ in these equations, since the difference between ${\mathit{E}}^{(\text{new})}$ and ${\mathit{E}}^{(\text{old})}$ is of the order $O(\mathrm{\Delta }\epsilon )$. Such a replacement is equivalent to the standard adjoint variable method, where the gradient of FOM with respect to an infinitesimal change of $\mathrm{\Delta }\epsilon$ is computed in this way. Thus, in a certain sense, one can view our work as a generalization of the standard adjoint variable method to the case of a discrete and large change of dielectric function.

3. DESIGN EXAMPLES

To demonstrate the power of the method, we consider first a 2D problem, consisting of a structure and a point source. The aim of the optimization is to maximize the light intensity at a target point separated from the source by a distance of one wavelength [Fig. 1(a)]. The figure of merit is therefore simply ${|\mathit{E}({\mathit{r}}_0)|}^2$, where ${\mathit{r}}_0$ is the position of the detector. The wavelength considered is the telecom wavelength $\lambda =1.55\mathrm{\mu m}$, and the design space has a size $2\lambda \times 2\lambda$ and is divided into $62\times 62$ pixels. Each pixel has a size of $50\mathrm{nm}\times 50\mathrm{nm}$. We consider the polarization where the electric field is perpendicular to the 2D plane.

 

Fig. 1. Concentration of light intensity at a target point. The optimization aims at maximizing light intensity at a target point (the black dot), separated from the source (the red dot) by a distance of one wavelength. The wavelength is $\lambda =1.55\mathrm{\mu m}$. We consider a scalar electric field perpendicular to the plane. The permittivities allowed for the binary system are respectively 2.25 and 6.25. The space for optimization has a square size with a dimension of $2\lambda \times 2\lambda$, divided into $50\mathrm{nm}\times 50\mathrm{nm}$ square pixels.

Download Full Size | PDF

We perform the optimization for a binary system, consisting of a background material with permittivity ${\epsilon }_{\mathrm{bck}}=2.25$, representative of silica, and a second material, with a permittivity differing from that of the background material by $\mathrm{\Delta }{\epsilon }_{\max }=4$, which is representative of the effective index of a Silicon slab. In the optimization process, the starting structure consists of a uniform background.

We first solve this problem with a few established approaches. All these approaches utilize the finite difference frequency domain method (FDFD) [24] as the method to solve Maxwell’s equations. These approaches all utilize the standard adjoint variable method, where the gradient of FOM to an infinitesimal change of $\mathrm{\Delta }\epsilon$ is computed. In Fig. 2, we show the performance of the FOM as a function of iteration steps for a few of these methods.

 

Fig. 2. Performance of various continuous optimization techniques for the light concentration setup in Fig. 1. (a) FOM in a log10 scale as a function of the iteration number; black curve: pixel-by-pixel binary-adjoint optimization guided by gradient information, starting from silica background; green curve: level-set optimization of a binary structure from a random binary configuration; blue and red curves: continuous optimization followed by level-set optimization starting either from silica background or a random starting point, respectively. The constant maximal permittivity change in the domain in each iteration is constrained to be 0.25. The crosses on the curve show the end of continuous optimization. (b) Corresponding to each case in (a), the left column shows the initial structure. The middle column shows the final structure with the corresponding FOM in a log10 scale. The right column shows the electric field intensity distribution for the final structure, in a log10 scale. (c) FOM in a log10 scale as a function of the iteration number. Both curves correspond to a continuous adaptive optimization process, starting from the same initial conditions as continuous optimization of (a), but with a maximal permittivity change between iterations continuously reduced to guarantee a monotonic ascent of the FOM during the optimization process. (d) Corresponding to the two cases in (c), the left column shows the structure at the end of the continuous optimization process. The middle column shows the structure after a thresholding process to create a binary structure. In both columns, the corresponding FOM in a log10 scale for each structure is also shown. The right column shows the electric field intensity distributions in a log10 scale for the final binary structures. Each iteration takes approximately 0.1 s on an Intel core i5-7200u CPU at 2.50 Ghz, RAM 8 GB.

Download Full Size | PDF

In Fig. 2(a), the blackline shows the performance the “binary-adjoint” method, starting from silica background. In this method, at each iteration, one performs a binary change for the single pixel where there is the highest gradient of FOM as calculated by the adjoint method. With this method, before reaching the final value, the FOM increases at every iteration. The FOM saturates and no longer improves after approximately 100 iterations. While not visible in the graph due to the use of a log scale in the figure, the FOM in fact oscillates in this regime: the iterations consist of adding and removing the same pixel. The behavior here indicates that the gradient as calculated in the adjoint variable method, which is for an infinitesimal change, does not accurately predict the change of FOM for a discrete large change of the dielectric function. The final value of the FOM is relatively low ($\mathrm{FOM}={10}^{1.97}$) compared with the other methods. The final structure features a one-dimensional array of a high-dielectric region, reminiscent of a waveguide microcavity [25], that increases light intensity at the target point [Fig. 2(b)].

The green line in Fig. 2(a) shows the performance of the “binary level-set” method. In this method, one starts from a random binary configuration. The gradient of the FOM with respect to the shape of the level set that defines the boundaries between the binary regions is then calculated with the adjoint method. Based on such gradient, the shape of the binary regions is updated based on the Hamilton–Jacobi evolution [10,12,14]. To satisfy the minimum feature size requirement, however, in each iteration, one updates the regions in a pixelated fashion. Compared with the binary-adjoint method, in this method, in each iteration, one simultaneously changes multiple pixels, which enables exploration of a larger phase space. The final FOM is higher ($\mathrm{FOM}={10}^{2.65}$), corresponding to a structure where the entire design region is used to enhance the field strength at the detector point. The convergence behavior, in terms of FOM as a function of the iteration number, is not monotonic. There are iteration steps where the FOM in fact decreases significantly, which again points to the fact that the adjoint variable method, which is designed to describe the change of FOM for the infinitesimal structure change, does not accurately predict the change of FOM for a discrete large change of the geometrical shape.

In Fig. 2(a), the blue and red lines show the performance of the method of continuous permittivity optimization, from two different starting points, of either a uniform silica background, or a random starting point, both of which are shown in Fig. 2(b). In this method, at each iteration, all pixels undergo a continuous permittivity change proportional to the calculated gradient of FOM. In gradient-based optimization methods, the magnitude of the change along the gradient direction still needs to be specified. Here, we set a constant maximal permittivity change in the domain in each iteration to be 0.25. With such method, the FOM improves monotonically until it saturates, at which point oscillation starts to occur. Once the oscillation occurs, we use a binary thresholding and a level-set optimization to generate a binary structure and to further improve the performance. The resulting $\mathrm{FOM}={10}^3$ is significantly better than the previous method, since the gradient information computed with the adjoint variable method is more applicable to guide permittivity change in each iteration, when such permittivity change is small.

Finally, in Fig. 2(c), we show the results from a continuous adaptive optimization method. This method is similar to the method of continuous permittivity optimization as discussed above, with the only difference being that the maximal permittivity change between iterations is continuously reduced in the optimization process. For this particular optimization problem, the optimal structures are highly resonant and can be very sensitive to small dielectric permittivity change. Reducing the maximal permittivity change adaptively thus ensures that the gradient information as obtained by the adjoint variable method remains relevant for structure update. In that case, FOM close to ${10}^4$ can be obtained, but the number of iterations is significantly larger. Also, in this method, we obtain the final binary structure through permittivity thresholding. But it is no longer possible to further improve the FOM on such binary structure through level-set optimization. As a result, the binary structure has a FOM that is quite a bit lowered as compared to the structure obtained from continuous optimization [Fig. 2(d)]. The final optimized structure is strongly non-intuitive. The field pattern corresponding to such a structure indicates both an increase of local density of states at the source location and a focusing of light toward the detector position.

We now optimize the system with our approach as described in Section 2. Section 2 shows that it is possible to efficiently compute the change of FOM with respect to all single-pixel variations, and to update the Green’s function of the system for each single-pixel variation. Based on this algorithm, in each iteration, we vary the single pixel to achieve the largest improvement of the FOM. There are a few attractive aspects of this algorithm. (1) In this method, there is no empirical parameter in the optimization process. This is in contrast with the methods discussed above, which use the parameters such as maximal permittivity change steps in continuous optimization, or velocity in level-sets optimization. In practice, there is substantial subtlety in choosing the value of these parameters to achieve good performance of the optimization algorithms. Removing the use of such empirical parameters is therefore an advantage of our proposed method. (2) In this method, the FOM improves at each iteration by design. (3) The optimization is guaranteed to end within a finite number of steps, and it ends when no single-pixel variation can improve the FOM.

We illustrate the performance of our method in Fig. 3, where we show the FOM as a function of iteration steps for a few random starting points. Since our method requires the knowledge of the Green’s function of such starting points, we start with a silica initial background, where the Green’s function is known, and add random pixels within the first iterations (250, 500, 750, and 1000 iterations), before starting the optimization process. For all starting points, our method results in final structures with a FOM above ${10}^4$, which significantly exceeds the best FOM obtained from the adaptive continuous optimization method. The optimized structures are automatically binary. The number of iterations required to reach the optimal structure is also smaller when compared to the adaptive continuous optimization method.

 

Fig. 3. Performance of our Green-function-based optimization techniques for the light concentration setup in Fig. 1. (a) FOM in a log10 scale as a function of the iteration number. The solid lines are from our optimization techniques, starting from a uniform silica design region. Each iteration takes 0.15 s on an Intel(r) core(tm) i5-7200u CPU at 2.50 ghz, RAM 8 GB. Solid black line: the optimization process starts from the first iteration. Other solid color lines: the first 250, 500, 750, or 1000 iterations involve randomly adding high index pixels. The optimization process starts afterwards. The dashed dotted lines are from adaptive continuous optimization included for comparison purposes. (b) Evolution of structure and the corresponding electric field intensity profiles during the optimization process. The structures are taken at the iterations numbers 336, 500, 1250, 2308 along the black line in (a). The FOM in a log10 scale of the final structure is also indicated. (c) The electric field intensity distribution of the final structure shown in (b), as calculated using the FDFD method. (d) The top and bottom rows show the initial and final structures, and the corresponding electric field intensity distribution shown in a log10 scale for the four cases as presented in Fig. 1(a), where the initial structures consist of randomly added pixels. The FOM in a log10 scale of the final structure is also shown.

Download Full Size | PDF

In Fig. 3(b), we show the evolution of the structure in the optimization process. The high-index regions first extend along the source–target axis, and then, as they reach the upper and bottom borders of the optimization space, it progressively expands inside the whole domain. This progressive pixel-by-pixel shape evolution is reminiscent of the cell-by-cell morphogenesis evolution in nature [26]. The high FOM of the solution obtained by our algorithm has been confirmed with FDFD [Fig. 3(c)].

We consider in Fig. 4 a second example for the design of an integrated demultiplexer. The input port of the device consists of a 1-μm-wide waveguide, and the output ports of the device consists of two 0.3-μm-wide monomodal waveguides. We inject light at two wavelengths ${\lambda }_1=1.31\mathrm{\mu m}$, and ${\lambda }_2=1.55$, into the fundamental mode of the input waveguide. The input waveguide is multi-moded at both wavelengths. The aim of the design is to separate the two wavelengths to two output waveguides. This device thus operates both as a mode converter and a wavelength sorter. To achieve this aim, we modify the structure in a $2\mathrm{\mu m}\times 2\mathrm{\mu m}$ design region located between the input and output waveguides. Like the example above, we consider only binary structures with permittivities of 2.25 and 6.25, respectively. The design region is divided into $40\times 40$ pixels. Each pixel has a size of $50\mathrm{nm}\times 50\mathrm{nm}$. We consider the polarization where the electric field is perpendicular to the 2D plane.

 

Fig. 4. Mode and frequency splitter. Light coming from the fundamental mode of a 1-μm-wide multimode waveguide is sent into either of the 0.3-μm-wide single-mode waveguides depending on its wavelength. Light with the wavelength ${\lambda }_1=1.31\mathrm{\mu m}$ is sent to the left waveguide, and light with the wavelength ${\lambda }_2=1.55\mathrm{\mu m}$ is sent to the right waveguide.

Download Full Size | PDF

In this design, the FOM is the average of power-coupling efficiencies, from the input mode $|\begin{array}{c}{\mathit{E}}_{\text{in}-1}\\ {\mathit{H}}_{\text{in}-1}\end{array}$ toward the output mode $|\begin{array}{c}{\mathit{E}}_{\text{out}-1}\\ {\mathit{H}}_{\text{out}-1}\end{array}$ at wavelengths ${\lambda }_1$, and the input mode $|\begin{array}{c}{\mathit{E}}_{\text{in}-2}\\ {\mathit{H}}_{\text{in}-2}\end{array}$ toward the output mode $|\begin{array}{c}{\mathit{E}}_{\text{out}-2}\\ {\mathit{H}}_{\text{out}-2}\end{array}$ at wavelengths ${\lambda }_2$, respectively:

As a comparison, we first optimize this system using a standard method, where we perform a continuous adjoint-based optimization, followed by binarization and binary-level-set optimization. We obtain a structure with a $\mathrm{FOM}=79\%$ for an empty initial structure [Figs. 5(a) and 5(b)], and a $\mathrm{FOM}=74\%$ for a filled initial structure [Figs. 5(c) and 5(d)]. Our method reaches a comparable high performance of $\mathrm{FOM}=78\%$ and $\mathrm{FOM}=80\%$, respectively, for both initial configurations (Fig. 5). In the previous example, where the structure is highly resonant, we see that in the optimization based on continuous optimization, it is difficult to reach a high-performing structure that is binary and that satisfies minimal feature constraint, since the structure is highly sensitive to detailed changes. In contrast, this system is not strongly resonant, and the FOM of the final structures obtained by either the continuous optimization methods or by our method are more comparable.

 

Fig. 5. Performance of our optimization techniques for the splitter setup in Fig. 4. (a), FOM as a function of iteration number. The red curves are from our approach. The blue curves are from the use of a continuous optimization followed by a level set optimization (FDFD). The design region in the starting structure is empty. (b) The top row shows the initial structure, the final structure, as well as the electric field intensity distribution at the two wavelengths of ${\lambda }_1=1.31\mathrm{\mu m}$ and ${\lambda }_2=1.55\mathrm{\mu m}$, respectively, with the corresponding power-coupling efficiencies at each wavelength. The bottom row is the same as the top row, except with our method. (c) Same as (a), for a starting structure with a filled design region. (d) Same as b, for starting structure with a filled design region. Green’s function initialization (FDFD Matrix inversion) requires 130 s per $\lambda$, then each iteration time is 0.1 s per $\lambda$ with our method, as for FDFD (the total domain is similar as previous example, but the design space where Green’s function is updated in our method is reduced), on an Intel core i5-7200u CPU @2.50 Ghz, RAM 8 GB.

Download Full Size | PDF

As seen in Figs. 5(a) and 5(c), in comparison to the standard continuous optimization method, the number of iterations required to reach the optimal structure is larger with our method. On the other hand, even in this case of non-resonant devices, where our method is less competitive in terms of the number of iterations, our method can still result in high-efficiency solutions that are generally not accessible with standard methods for achieving a binary structure such as the level-set methods.

A key advantage of our method is that it can straight-forwardly consider the minimum feature size requirement. We show in Fig. 6 the solutions obtained for the first device optimization example of Fig. 1(a), starting from silica background, for minimal feature sizes of 50 nm [Fig. 6(a)], 100 nm [Fig. 6(b)], and 150 nm [Fig. 6(c)], respectively. In all cases, in the simulation, the pixels as used in Eq. (3) are still 50 nm in size. When we update the structure in each iteration, however, we update a block of pixels simultaneously, with a block containing 1, 4, and 9 pixels, respectively, to satisfy the minimum feature size requirement. The resulting structures from such an update always satisfy the minimum feature size requirements.

 

Fig. 6. Optimization under the constraint of minimum feature sizes, using the example shown in Fig. 1. (a) Optimized structure and the electric field intensity distribution, for a minimum feature size of 50 nm. (b) Same as (a), for a minimum feature size of 100 nm. (c) Same as (a), for a minimum feature size of 150 nm. (d) FOM as a function of iteration number, using our method, for the three cases shown in (a)–(c). The times taken per iteration are 0.15 s, 0.75 s, 1.25 s for the minimum feature sizes of 50 nm, 100 nm, and 150 nm, respectively.

Download Full Size | PDF

In general, the FOM of the optimal structures, as well as the number of iteration steps, decreases as the minimum feature size increases. This is expected since the size of the phase space decreases as the minimum feature size increases. Examining the three optimal structures shown in Figs. 6(a)–6(c), we observe that, switching from 50 nm pixels to 100 nm pixels, the two optimal structures for 50 nm and 100 nm minimum feature sizes are quite similar to each other, whereas for 150 nm minimum feature sizes, the optimal structure is completely different. The result here clearly indicates the importance of incorporating the minimum size in the optimization process: an optimization process that disregards the minimum feature size requirement may not produce a structure that is implementable for a given fabrication technology. In addition, incorporating the minimum feature size requirement restricts the phase space and hence can result in the use of less iteration steps. The concept of simultaneous update of a group of pixels can also be used to implement other optimization constraints, for example, for imposing symmetry in a given structure.

4. CONCLUSION

In this article, we have introduced an efficient pixel-by-pixel discrete optimization technique, which rigorously guarantees a monotonic ascent of the figure of merit in photonic optimization, even for large discrete binary permittivity changes. This technique is based on the Dyson’s equation in a Green’s function formalism. Our technique guarantees an optimized structure that is binary, and it straightforwardly incorporates constraints on minimum feature size requirements in nanofabrication. Compared with the standard optimization method based on combinations of adjoint variable and level-set methods, we present numerical examples showing that our method produces higher-efficiency structures, especially for resonant systems. Our method could be important for systematical optimization of on-chip photonic devices.