Open Access
Add to BibSonomyAbstract
A computationally-fast inverse design method for nanophotonic structures is presented. The method is based on two complementary convex optimization problems which modify the dielectric structure and resonant field respectively. The design of one- and two-dimensional nanophotonic resonators is demonstrated and is shown to require minimal computational resources.
©2010 Optical Society of America
1. Introduction
Numerous numerical methods have been devised to solve Maxwell’s equations in both time [1] and frequency [2,3] domains. We refer to these schemes as direct solvers, since they compute the electric and magnetic fields based on current sources, charge distributions and surrounding dielectric and/or metallic structures. While extremely useful in simulating optical components, using direct methods to design optical components, especially in two or three dimensions, typically requires an extremely time-consuming direct search in a large parameter space [4–8].
On the other hand, an inverse solver would be much more adept in such design and optimization problems [9, 10]. In this work, we define the inverse problem as that where the electromagnetic field is known, but the surrounding structure is not known. The goal in the inverse problem is, then, to find a dielectric structure that will produce that specific electromagnetic field profile.
We show that one can design nanophotonic resonators by specifying the electromagnetic field and its desirable characteristics (such as cavity quality (Q) factor and/or mode volume) and then using an inverse solver to find the corresponding dielectric structure. We show that the inverse method used is not only computationally-fast, but is also able to optimize for multiple device characteristics and produce multiple resonances, both of which are very difficult using direct methods.
2. Numerical setup
We start from the time-harmonic eigenvalue equation
where H, ε, ω and c are the magnetic field, relative permittivity, resonance frequency and speed of light respectively. To solve the problem numerically, H and ε are discretized in space using the standard Yee cell used in finite difference methods [1]. Also, the curl operators, since they are linear, are represented by the matrix A. Equation (1) can now be written as
where
A is the discretized curl operator,
Y = diag(ε -1) is the diagonal matrix representing the dielectric structure,
x is a vector representing H, and
In this form, given Y, we can solve the direct problem by computing x using an eigenvalue solver [13]. However, we note that Eq. (2) is also linear in Y, which allows us, if x is held constant, to solve the inverse problem by expressing Eq. (1) as
where
B = A·diag(Ax),
d = ξx, and
the variable for which we solve.
Here, diag(Ax) is the matrix with the values of Ax along the main diagonal and zeros elsewhere.
3. Least-squares method in 1D
3.1. Least-Squares
Fig. 1. Inverse design of a one-dimensional structure using the unmodified least-squares method. The target field is a sinusoid within a Gaussian envelope. The computed dielectric structure (green area) supports a field (red circles) that exactly matches the target field (blue line). The entire design process is also extremely fast and takes less than 1 second to complete on a generic desktop computer [12]. The periodic singularities in the dielectric structure are non-physical and will be addressed later in the article.
The result of applying this method to a simple one-dimensional problem is shown in Fig. 1. A generic least-squares solver [11] was used to find the dielectric structure, y (green region), that exactly produces the target field, x (blue line), using Eq. (2). Using a generic desktop computer, the solution was obtained in less than a second. Then a finite-difference time-domain (FDTD) solver was used to obtain the actual field (red circles) produced by the structure and to verify the accuracy of y.
As expected, Fig. 1 shows that the target field is reproduced exactly by the dielectric structure. However, the resulting structure is full of undesireable singularities. The rest of the section focuses on producing a well-behaved dielectric structure that still reproduces the target field accurately.
3.2. Regularized least-squares
The simplest way to produce a well-behaved dielectric structure is to add a regularization term to our least-squares problem, which is equivalent to solving the following optimization problem
Here y 0 represents some initial guess for the dielectric structure, which we want the values of y to stay close to, and η > 0 is a parameter used to trade off fit, i.e. ∥By - d∥2, and deviation from y 0, i.e. ∥y - y 0∥2.
We chose to constrain ε around a constant value of ε 0 = 10 and solved the least-squares system for η = 10-8, 10-6, and 10-4. The results, each still obtained in under a second, are shown in Fig. 2 and illustrate the trade-off between constraining ε and accurately reproducing H.
Fig. 2. Inverse design of one-dimensional structures using the regularized least-squares method. The same target field is used as in Fig. 1, and the computation time remains below 1 second. As the regularization parameter, η, is increased, ε is increasingly constrained to a chosen constant value of 10. At the same time, the mismatch between target and actual fields increases markedly. This illustrates the apparent trade-off between producing reasonable structures and accurately reproducing a fixed target field.
4. Complementary optimization in 1D
4.1. Motivation for a Complementary Optimization Strategy
The fundamental problem in the previous examples is actually not in the methods themselves, but in the improper selection of a target field. In fact, it is very difficult to select a suitable target resonant field because not every resonant mode even has a corresponding dielectric structure that is able to reproduce it. Furthermore, it is nearly impossible to select a multidimensional field which corresponds to a well-behaved, isotropic and discretely-valued ε, as would be needed for practical structures.
For this reason, a successful method must be allowed to either modify the target field, or specify it completely, in which case the user would only determine certain characteristics (e.g. mode-volume, Q-factor) that the target field should have. The former strategy is developed in both one and two dimensions, while the latter strategy is implemented in Section 6 in order to design two-dimensional resonators with discrete values of ε.
4.2. Complementary optimization
We start with the same target field as in the previous examples but we now formulate a method that allows for it to be modified during the design process. The formulation chosen is a complementary optimization routine, where we continually alternate between modifying ε to better fit the field, and then modifying the field to better fit ε. Here, we use the term “fit” to mean that either ε or H is solved so that the residual error from Eq. (1) is minimized. Additionally, both iterations are regularized in order to stably approach a solution. This algorithm can be summarized as follows,
choose x 0 and y 0
for i = 1,2,…
where Yi = diag(yi), Bi = A · diag(Axi) and di = ξxi. ∥AYi Axi - ξx i-1∥2 is used instead of ∥AYi Axi - ξxi∥2 to avoid the trivial xi = 0 solution and does not affect the overall accuracy since x changes very slowly.
Fig. 3. Inverse design of a one-dimensional structure using the complementary optimization method. The target field in Figs. 1 and 2 is used as the initial target field. The rates of change for both ε and H are controlled by regularization parameters η 1 = 10-4 and η 2 = 10-3 respectively. The 400 iterations used to achieve this result took 60 seconds to compute. This method results in a well-behaved ε that actually produces a field very similar to the original target field. Interestingly, the formation of a “steady-state” periodic structure toward the sides of the structure has emerged.
Figure 3 shows that the complementary optimization algorithm, after 400 iterations and with the correct choice of regularization parameters η 1 and η 2, results in a well-behaved structure that is able to closely reproduce the modified target field. Numerically, the least-squares problem must now be solved numerous times, which increases the computational time needed to around 60 seconds.
4.3. Complementary optimization with bounded ε
In order to achieve a more practical, discretely-valued dielectric structure, we can impose strict upper- and lower-bounds on ε. To this end, we modify our algorithm as such,
choose x 0 and y 0
for i = 1,2,…
In this algorithm, Eq. (7) is a convex optimization problem [14]. This allows us to impose hard constraints on ε, which in turn allows us to remove the regularization term present in Eq. (5). The CVX package [15], a Matlab-based modeling system for convex optimization, is used to solve Eq. (7), with each iteration of the algorithm now requiring roughly 1 second of computation time.
Fig. 4. Inverse design of a one-dimensional structure using the complementary optimization method with bounded ε. The parameters are identical to those used to produce Fig. 3 with the exception that only one regularization term is now needed (η 2 = 10-3). The algorithm was run for 100 iterations, which took 100 seconds. The structure turns out to be almost completely binary-valued and looks like a periodic structure with tapered duty cycle. It produces an actual field which very closely matches the final target field.
A nearly binary-valued dielectric structure is obtained in Fig. 4, which accurately produces the final target field. This is very useful for the design of practical structures, since they usually consist of two or three different materials at most. Interestingly, although the directly discreteness of ε was not enforced (since that would make the problem non-convex), a discrete, binary-valued structure has still arisen.
Fig. 5. Inverse design of an “S” resonator using the complementary optimization method without bounds on ε. The design was initialized by specifying an initial dielectric structure (ε= 1 everywhere) and a resonant field in the shape of an “S”. The final dielectric structure was produced after 50 iterations which took 90 seconds to complete in total. The grid size was 80×120. The final dielectric structure is quite unintuive, and yet reproduces the target field surprisingly well. This example demonstrates the versatility of the complementary optimization method in producing designs, from very simple specifications, which otherwise could be attained only with considerable difficulty.
5. Complementary optimization in 2D
5.1. “S” Resonator
We now demonstrate that the complementary optimization method is versatile and can be scaled to multiple dimensions. To ensure that ε is well-behaved we use a point-spread function which does not allow ε to change at a certain point in space without affecting the values surrounding it.
In order to show that our method can produce complex designs, we choose an S-shaped target field which is non-trivial to reproduce. The optimization results, using the complementary optimization method from Section 4.2, are shown in Fig. 5. The resulting dielectric structure is continuous, unbounded and contains some singularities (white dots), but the final target and actual fields match up well. Also, the computational cost remains quite reasonable; the 50 iterations needed required only 5 minutes of computation time. The resulting structure is completely unintuitive, and illustrates the kind of new capabilities offered by the inverse design strategy. Specifically, that a complex, intricate structure can be designed just by specifying the shape and frequency of a rather simple electromagnetic mode.
5.2. Multi-mode inverse design
The complementary optimization method can also be extended to produce dielectric structures with multiple resonances. To do so, multiple initial target fields are specified. The dielectric structure is first modified to simultaneously fit all target fields using a multi-objective least-squares method. Then each target field is individually modified to fit the structure; and we continue alternating between optimizing ε and H (j) in this way. A benefit of this scheme is that only the ε optimization increases in size, so the design process remains computationally tractable, even for several resonant fields. This algorithm can be summarized as follows,
Fig. 6. Inverse design of a doubly-resonant, degenerate “X” resonator produced by the complementary optimization method with unbounded ε. As in Fig. 5 the initial value of ε was 1 everywhere. The two initial target fields used are only slightly perturbed and are very similar to the two final actual target fields. Computationally, this design took 5 minutes to complete on a 120 × 120 grid and required 40 iterations. This example shows that the complementary optimization strategy can be extended to produce dielectric structures with multiple resonances. Such an “X” resonator is useful for polarization-entangled single-photon sources for example [16].
choose y 0, x 0 (1), x 0 (2), …
for i = 1,2,…
for i = 1,2,…
The design of an “X” resonator with two perpendicular, degenerate modes is shown in Fig. 6. The added complexity increases the total computation time to 5 minutes for 40 iterations of the algorithm.
5.3. Design of waveguiding devices
The multi-mode inverse design method can also be applied to the design of waveguiding devices such as multiplexers/demultiplexers, waveguide couplers, crossbars and dispersion-tailored waveguides. This can be accomplished by treating a waveguiding device as a doubly-degenerate resonator at its operational frequencies and then enforcing opposite symmetries (even/odd or cosine/sine) in the degenerate modes. A single-to-dual beam waveguide coupler designed based on Eq. (9) and (10) is shown in Fig. 7, and motivates how one might design other waveguiding devices such as channel-drop filters. In addition to macroscopic waveguiding devices, one should also be able to create novel periodic waveguiding structures. For instance, by controlling the frequency of a particular waveguiding mode at each of its k-vectors, one could create a waveguide with a customized dispersion characteristic, which would be useful in the design of slow-light waveguides.
Fig. 7. Inverse design of a single-to-dual beam waveguide coupler using the complementary optimization method without bounds on ε. Two degenerate modes with opposite symmetry (sine and cosine) are used as target fields (only one is shown). Only 4 iterations (14 seconds) are needed to achieve this solution on a 240×55 grid. This is a simple demonstration showing that the complementary optimization method can also be extended to design waveguiding devices.
6. Complementary optimization with bounded ε in 2D
6.1. Numerical method
We now use the complementary optimization method to design resonators with discrete, binary ε in two dimensions. At the same time, the algorithm is modified so that we no longer specify an initial target field; instead, only an initial dielectric structure and the maximum desired mode volume (i.e. mode area in 2D) are specified. In addition, the optimization process attempts to maximize the Q-factor. Such an algorithm now consists of iteratively solving two convex optimization problems,
choose y 0
for i = 1,2,…
In this algorithm, the field optimization, Eq. (11), differs significantly from Eq. (8). First, the eigenvalue equation term ∥AY i-1 Axi - ξxi∥2 no longer utilizes x i-1. Also, a new “Fourier-minimizing” term ∥Fxi∥2 has been added. Here, the row vectors of F consist of the field Fourier components that we do not want incorporated in the final solution. The motivation for adding this term is to design high-Q resonators using planar structures, where the Q-factor is limited by out-of-plane losses, and can thus be improved by eliminating field Fourier components that are not localized by total internal reflection (components inside the light cone) [10]. Lastly, the (Axi)T Y i-1(Axi) term allows us to specify the mode area (mode volume in three dimensions), that we desire for our resonant field. This works because the two minimization objectives in Eq. (11) will generally cause the mode area to be as large as possible, which means we always end up with (A 1 xi)T Y i-1(A 1 xi) = A mode.
The addition of the two new terms in Eq. (11) signifies that the field iteration in our algorithm attempts to do more than to just satisfy Maxwell’s equations. Rather than only trying to fit the dielectric structure, the resonant field now also minimizes some of its Fourier components while working with a limited mode area.
In contrast with the field optimization given by Eq. (11), the structure optimization given by Eq. (12) remains identical to the equation in the one-dimensional case, Eq. (7), and contains no terms related to notions of Fourier components or mode area. In fact, the only objective in the ε optimization is to better fit the field generated from the field optimization. This means that in this scheme, the field iteration is leading the structure iteration. In other words, Eq. (11) “looks ahead” and gradually adapts itself to become a more desirable field, while Eq. (12) just “follows along”. For this reason, this algorithm is unique from the other complementary optimization algorithms previously presented in the article, since in those algorithms both iterations follow one another and eventually “meet in the middle”.
Finally, in this scheme, we chose to limit the degrees of freedom of both xi and yi. Some components of xi were fixed at non-zero values and were not allowed to be modified simply to avoid the xi = 0 situation. To enhance the aesthetics of the resulting dielectric structure, some components of yi were frozen as well, which is useful in cases where one would only like to modify the dielectric structure within a waveguide only, for example.
6.2. Circular grating resonator
Figure 8 shows the design of a circular grating resonator using the method from Section 6.1. The dielectric structure emerged from a very simple choice of initial structure, namely a constant ε = 12.25 everywhere. The range of ε was limited to be from 1 to 12.25, the mode area was set to 5 and η = 10-3. Central components of xi were held constant to ensure that an electric dipole in the y-direction was produced in the center of the structure. Additionally, the components of ε outside a specified circle were held at a constant ε = 12.25 for the duration of the design process. The entire algorithm was run for 40 iterations and took 7 minutes.
After 40 iterations, we see that a strong binary-valued structure has formed. Interestingly, the central bowtie-like structure has emerged from previous genetic optimization methods as well [7]. Note also the extremely small amount of information needed to produce this result: only the frequency and mode area of the resonant field desired, and a trivial intial ε. The ability to produce a rather advanced design from such a simple problem specification highlights the potential usefulness of inverse methods in the design of novel nanophotonic devices.
6.3. Beam resonator
The same approach was used to design a beam resonator as shown in Fig. 9. The parameters used in this design are identical to those used for the circular resonator, with the exception that the initial dielectric structure now consists of an unbroken waveguide, and the region where ε can be modified is now confined to the center of the waveguide.
Fig. 8. Inverse design of a two-dimensional resonator using the complementary optimization method with strict bounds on ε. The initial specification is very simple and consists of an initial dielectric structure (ε = 12.25 everywhere), the frequency and mode-volume of the resonance field as well as a weighting factor, η, to avoid leaky field Fourier components. Additionally, the values of ε are only allowed to be modified within a central circular region and must be kept between 1 and 12.25. After 40 iterations on a 160 ×160 grid, which took 7 minutes to complete, a discrete structure emerged with excellent match between the predicted (x 40) and actual fields. The structure resembles a circular grating with a bowtie-like central structure for focusing the resonant energy to a single point.
Fig. 9. Inverse design of a beam resonator in two dimensions using the complementary optimization method with bounded ε. The initial conditions are identical to those for Fig. 8, except that the initial dielectric structure is an unbroken waveguide, and ε can only be modified within that waveguide. The structure emerged after 40 iterations on a 320 × 40 grid, which took 5 minutes of computation. The bowtie-like structure has reappeared in the center. Interestingly, the effect of the Fourier term in Eq. (11), is seen in the outward tapering of the holes.
These two, simple modifications produce a vastly different result as is plainly seen from Fig. 9. Some characteristics remain, such as the two closely spaced holes in the center of the structures which focus the electromagnetic energy to a central point. Other features are unique to the beam resonator, such as a gradual outward tapering of hole diameters and size ending in a periodic “steady-state” structure, in similar fashion to Fig. 3. This tapering is a direct result of the Fourier term in Eq. (11) as it leads to a smooth field profile variation and minimization of the components in the light cone that lead to out-of-plane losses [10, 17].
7. Conclusions and outlook
In summary, we have introduced and demonstrated a complementary optimization method for the inverse design of nanophotonic resonators as well as waveguiding structures. We have demonstrated that, in two dimensions, this method can be used to design non-trivial as well as multi-resonant structures. Furthermore, we have shown that by constraining the range of available dielectric constants, the same underlying method can produce binary-valued, discrete dielectric structures.
The methods presented here can readily be extended into three dimensions. In essence, the only hurdle is the size of the convex optimization problem for the field optimization. Since a full three-dimensional grid often consists of tens of millions of variables, solving the field optimization problem, Eq. (11), using methods such as LU factorization is no longer computationally tractable. Instead, iterative methods such as truncated Newton methods must be employed. As it turns out, such iterative methods are especially well suited to this complementary optimization scheme, since they can take advantage of the information in x i-1 in order to solve xi very quickly. This is very applicable to our algorithm since the field changes very little from one iteration to the next.
Although solving the field optimization problem is numerically challenging in three dimensions, solving the structure optimization problem, Eq. (12), in three dimensions is much simpler. This is because only structural variations in-plane need to be accounted for in the design of planar nanophotonic devices. The number of variables in Eq. (12) for both two- and three-dimensional problems is then effectively the same, and the same numerical techniques can be employed to solve both.
A full three-dimensional inverse design method based on complementary optimization would indeed enable the design of novel nanophotonic devices. For instance, a resonator for efficient wavelength-conversion could be designed by creating a structure with multiple well-placed resonant frequencies [18]. Furthermore, even effects such as mode overlap and the use of frequency-dependent refractive index materials could be taken into account with such a method.
Acknowledgements
This work has been supported in part by the AFOSR MURI for Compex and Robust On-chip Nanophotonics (Dr. Gernot Pomrenke). We thank Prof. Stephen Boyd for many fruitful discussions regarding convex optimization.
References and links
1. K. Yee, “Numerical solution of initial boundary value problems involving Maxwells equations in isotropic media,” IEEE Trans. Antennas Propag. Mag. 14, 302–307 (1966). [CrossRef]
2. M. Albani and P. Bernardi, “A Numerical Method Based on the Discretization of Maxwell Equations in Integral Form,” IEEE Trans. Microwave Theory Tech. 22, 446–450 (1974). [CrossRef]
3. J. M. Gerardy and M. Ausloos, “Absorption spectrum of clusters of spheres from the general solution of Maxwell’s equations. The long-wavelength limit,” Phys. Rev. B 22, 4950–4959 (1979). [CrossRef]
4. P. Deotare, M. McCutcheon, I. Frank, M. Khan, and M. Loncar, “High quality factor photonic crystal nanobeam cavities,” Appl. Phys. Lett. 94, 121106 (2009). [CrossRef]
5. J. Vuckovic, M. Loncar, H. Mabuchi, and A. Scherer, “Design of photonic crystal microcavities for cavity QED,” Phys. Rev. E 65, 1–11 (2002).
6. Y. Akahane, T. Asano, B. Song, and S. Noda, “Fine-tuned high-Q photonic-crystal nanocavity,” Opt. Express 13, 1202–1214 (2005). [CrossRef] [PubMed]
7. A. Gondarenko and M. Lipson, “Low modal volume dipole-like dielectric slab resonator,” Opt. Express 16, 17689–17694 (2008). [CrossRef] [PubMed]
8. A. Hakansson and J. Sanchez-Dehesa, “Inverse designed photonic crystal de-multiplex waveguide coupler,” Opt. Express 13, 5440–5449 (2005). [CrossRef] [PubMed]
9. P. Borel, A. Harpth, L. Frandsen, M. Kristensen, P. Shi, J. Jensen, and O. Sigmund, “Topology optimization and fabrication of photonic crystal structures,” Opt. Express 12, 1996–2001 (2004). [CrossRef] [PubMed]
10. D. Englund, I. Fushman, and J. Vuckovic. “General Recipe for Designing Photonic Crystal Cavities,” Opt. Express 12, 5961–5975 (2005). [CrossRef]
11. CHOLMOD software package, accessed via Matlab.
12. Intel Core 2 Quad 2.5GHz, 8Gb RAM.
13. S. G. Johnson and J. D. Joannopoulos, “Block-iterative frequency-domain methods for Maxwells equations in a planewave basis,” Opt. Express 8, 967–970 (1999).
14. S. Boyd and L. Vandenberghe, Convex Optimization (Cambridge University Press, 2004).
15. M. Grant and S. Boyd, CVX: Matlab software for disciplined convex programming, http://stanford.edu/∼boyd/cvx, June 2009.
16. K. Hennessy, C. Hogerle, E. Hu, A. Badolato, and A. Imamoglu, “Tuning photonic nanocavities by atomic force microscope nano-oxidation,” Appl. Phys. Lett. 89, 041118 (2006). [CrossRef]
17. B. -S. Song, S. Noda, T. Asano, and Y. Akahane, “Ultra-high-Q photonic double-heterostructure nanocavity,” Nat. Mater. 4, 207–210 (2005). [CrossRef]
18. K. Rivoire, Z. Lin, F. Hatami, W. Ted Masselink, and J. Vuckovic, “Second harmonic generation in gallium phosphide photonic crystal nanocavities with ultralow continuous wave pump power,” Opt. Express 17, 22609–22615 (2009). [CrossRef]
