In part I of this series of papers, we introduced an efficient pixel-by-pixel optimization technique, utilizing the Dyson’s equation in a Green’s function formalism. This technique directly incorporates material and minimum feature size constraints in practical photonic devices. The implementation in part I, however, requires full storage of the Green’s function, which is costly. In this paper, we further develop this optimization technique by showing that only the diagonal part of the Green’s function needs to be stored, and by showing that the optimization process can be implemented using standard commercially available electromagnetic solvers. The development here should enable a wider use of this optimization technique.
© 2019 Optical Society of America under the terms of the OSA Open Access Publishing Agreement
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 . 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 . 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 , we have introduced a method, based on Dyson’s equation  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 , 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 . 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 , where 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 . 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 . 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 , we consider an electric field generated by an excitation for a structure as described by a given permittivity . 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 :
Define the change in permittivity distribution between two iterations. Using the Lippmann–Schwinger equation , we can determine the new electric field in the whole space:
In the case where we consider a single-pixel modification in each iteration, is non-zero only in a small volume V centered at , and the change in FOM consists of12] 7) is justified as long as the size of the volume V is much smaller than the wavelength.
As we show in part I , the new electric field at after a pixel has been updated in 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 :7) and (8), the best pixel to update can be rigorously and efficiently determined even for large and discrete permittivity changes, as required in binary systems . We refer to this pixel as the “best-improvement pixel” below.
In the DDA-type approach described in part I , for the “new” structure where the dielectric constant is modified at the best-improvement pixel, the location of which is denoted as , the field was updated everywhere using the Lippmann–Schwinger equation:12]: 10) first for : 10).
In Ref. , when implementing the algorithms above, we stored the full Green’s functions 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 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 , where 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 scales as , where is the number of pixels in the design region. Also, instead of using Eq. (9), which involves the full Green’s function 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 have15].
In practice, the number of simulations required to determine and 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 .
In the case of transverse electric (TE) polarization in two dimensions, is scalar, and therefore Eq. (9) is sufficient to determine it. The required information is: before the best-improvement pixel in is updated, and and after the best-improvement pixel is updated. All of this information has already been calculated. Thus, 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 tensor. Utilizing both as well as , 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 is a 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 and . We can use a dipole source located in , and with the dipole moment orthogonal to both and , which provides a linear system with good conditioning . The Lippmann–Schwinger equations for the three new fields then read
We now switch to the determination of . 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 with , to arrive at a new set of equations that involve . As a result, Eq. (14) transforms to15) and (16) transform to 17)–(19).
We therefore arrive at the following algorithm: at the start of each iteration, we assume that we know the field distribution , the adjoint field distribution , as well as the diagonal part of the Green’s function , where 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 and of the new structure using a standard electromagnetic solver, and determine the diagonal part of the Green’s function 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 over the entire design region. In some cases, e.g., if the initial structure consists of only vacuum, 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 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 as the iteration progresses. Starting from an initial structure labeled “0,” we can then randomly flip a pixel at to arrive at the structure labeled “1.” Using Eqs. (17)–(19), part of the Green’s function for structure 1, can be determined using the information of the fields for structures 0 and 1, without the knowledge of . At this point, we approximate , and proceed to determine the best-improvement pixel and generate the structure labeled “2.” We note that can be determined rigorously in the same process as described above using Eqs. (17)–(19), whereas can be determined accurately from using Eq. (12). We then construct an approximation of from and . We repeat this process as outlined above in the iterations. As the iterations continue, the approximation to becomes increasingly accurate, since at the best-improvement pixels at all previous iterations are accurately determined. We also note that the accurate information of 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 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. , 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. , 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) becomes8) becomes 20) involves solving Eq. (21), which is a small linear system with the dimension linearly proportional to the number of pixels () inside a block. With this information, we then modify the block that leads to the best improvement of FOM, to arrive at the new structure; below, we refer to this block as the “best-improvement block.”
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 . This is done by noting that updating a single block is equivalent to updating sequentially the pixels inside the block. Thus, we can update in a pixel-by-pixel fashion. Denoting the pixel being updated as , the block-diagonal part of the Green’s function is updated using13)–(19).
3. DESIGN EXAMPLES
We now provide examples of implementation for the algorithm as described above. We use a commercial code RSOFT-Fullwave , 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 .
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 . We consider the TE polarization where the electric field is perpendicular to the 2D plane.
In the optimization process, we start from a silica homogeneous background, for which the Green’s function can be analytically obtained . 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. . 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 [ 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. , 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 . The device operates as a two-wavelength demultiplexer, which separates two incident waves with wavelengths and 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 .
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 , 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 ( 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 . 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 ( 50-nm-sized square pixels including PMLs) together is only 392 kB, compared to 164 MB required for the DDA approach.
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 [) slab with 0.3 μm thickness, surrounded by silica (]. The aim of the optimization is to maximize light intensity emitted from a -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 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 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 direction.
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 of a given block as the “new” structure. Equation (16) can then be used to deduce for all within the same block as . Since the initial structure is uniform, from this we can then deduce for and 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 , is here again large. The memory required by the optimization code and the RSOFT code ( 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 ().
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 mode into 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-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.
In part I  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.
Air Force Office of Scientific Research (AFOSR) (FA9550-17-1-0002); FP7 People: Marie-Curie Actions (PEOPLE) (600382).
1. D. J. Lockwood and L. Pavesi, “Silicon fundamentals for photonics applications,” in Silicon Photonics (Springer, 2004), pp. 1–50.
2. N. K. Georgieva, S. Glavic, M. H. Bakr, and J. W. Bandler, “Feasible adjoint sensitivity technique for EM design optimization,” IEEE Trans. Microwave Theory Tech. 50, 2751–2758 (2002). [CrossRef]
3. G. Veronis, R. W. Dutton, and S. Fan, “Method for sensitivity analysis of photonic crystal devices,” Opt. Lett. 29, 2288–2290(2004). [CrossRef]
4. C. M. Lalau-Keraly, S. Bhargava, O. D. Miller, and E. Yablonovitch, “Adjoint shape optimization applied to electromagnetic design,” Opt. Express 21, 21693–21701 (2013). [CrossRef]
5. O. D. Miller, “Photonic design: from fundamental solar cell physics to computational inverse design,” arXiv:1308.0212 (2013).
6. J. Lu and J. Vučković, “Objective-first design of high-efficiency, small-footprint couplers between arbitrary nanophotonic waveguide modes,” Opt. Express 20, 7221–7236 (2012). [CrossRef]
7. A. Y. Piggott, J. Lu, K. G. Lagoudakis, J. Petykiewicz, T. M. Babinec, and J. Vučković, “Inverse design and demonstration of a compact and broadband on-chip wavelength demultiplexer,” Nat. Photonics 9, 374–377 (2015). [CrossRef]
8. V. Liu, Y. Jiao, D. A. Miller, and S. Fan, “Design methodology for compact photonic-crystal-based wavelength division multiplexers,” Opt. Lett. 36, 591–593 (2011). [CrossRef]
9. A. Y. Piggott, J. Petykiewicz, L. Su, and J. Vučković, “Fabrication-constrained nanophotonic inverse design,” Sci. Rep. 7, 1786 (2017). [CrossRef]
10. M. A. Seldowitz, J. P. Allebach, and D. W. Sweeney, “Synthesis of digital holograms by direct binary search,” Appl. Opt. 26, 2788–2798 (1987). [CrossRef]
11. B. Shen, P. Wang, R. Polson, and R. Menon, “An integrated-nanophotonics polarization beamsplitter with 2.4× 2.4 μm2 footprint,” Nat. Photonics 9, 378–382 (2015). [CrossRef]
12. S. Boutami and S. Fan, “Efficient pixel-by-pixel optimization of photonic devices utilizing the Dyson’s equation in a Green function formalism: Part I: implementation with the method of discrete dipole approximation,” J. Opt. Soc. Am. B 36, 2378–2386 (2019). [CrossRef]
13. O. J. Martin, A. Dereux, and C. Girard, “Iterative scheme for computing exactly the total field propagating in dielectric structures of arbitrary shape,” J. Opt. Soc. Am. A 11, 1073–1080 (1994). [CrossRef]
14. B. T. Draine, “The discrete-dipole approximation and its application to interstellar graphite grains,” Astrophys. J. 333, 848–872 (1988). [CrossRef]
15. O. J. Martin and N. B. Piller, “Electromagnetic scattering in polarizable backgrounds,” Phys. Rev. E 58, 3909 (1998). [CrossRef]
16. K. Yee, “Numerical solution of initial boundary value problems involving Maxwell’s equations in isotropic media,” IEEE Trans. Antennas Propag. 14, 302–307 (1966). [CrossRef]
17. A. Taflove, “Application of the finite-difference time-domain method to sinusoidal steady-state electromagnetic-penetration problems,” IEEE Trans. Electromagn. Compat. EMC-22, 191–202 (1980). [CrossRef]
18. J. M. Jin, The Finite Element Method in Electromagnetics (Wiley, 2015).
19. W. Shin and S. Fan, “Choice of the perfectly matched layer boundary condition for frequency-domain Maxwell’s equations solvers,” J. Comput. Phys. 231, 3406–3431 (2012). [CrossRef]
20. Fullwave® and Synopsis®, “Rsoft Design,”.
21. A. D. Yaghjian, “Electric dyadic Green’s functions in the source region,” Proc. IEEE 68, 248–263 (1980). [CrossRef]