Open Access
Add to BibSonomyAbstract
The inverse design of a three-dimensional nanophotonic resonator is presented. The design methodology is computationally fast (10 minutes on a standard desktop workstation) and utilizes a 2.5-dimensional approximation of the full three-dimensional structure. As an example, we employ the proposed method to design a resonator which exhibits a mode volume of 0.32(λ/n)3 and a quality factor of 7063.
©2011 Optical Society of America
1. Motivation
To date, the design of nanophotonic devices has generally involved a lengthy (days, weeks) process in which one perturbs, by trial-and-error, canonical structures such as photonic crystals or waveguide-coupled rings to achieve the desired performance for the device.
Instead, an inverse design method in which the user only specifies the desired electromagnetic field, or characteristics thereof, and then leaves the computer to find a dielectric structure satisfying these requirements, would be a much more intuitive and computationally efficient design strategy. Furthermore, such a method may even be the only feasible option available for the design of more complex multi-wavelength and multi-functional nanophotonic devices needed for on-chip integration [1].
2. 2.5-dimensional approximation
Previously, we demonstrated the inverse design of various nanophotonic devices in two dimensions [2]. Unfortunately, to directly extend our previous method to full three-dimensional space requires solving for a very large (∼ 107 × 107) and ill-conditioned matrix, namely that given by the time-harmonic Maxwell’s equation in three dimensions, which for the electric (E) field is
Therefore, to make our problem more tractable, we formulate a 2.5-dimensional approximation of the electromagnetic fields of a planar three-dimensional nanophotonic device. As shown in Fig. 1(a), this involves treating a planar three-dimensional nanophotonic device as a truncated holey fiber with identical dielectric structure. We take advantage of the planar nature of most nanophotonic devices in this way to produce a tractable, computationally-fast method. The E-field inside such a holey fiber is expressed as
where β is the wave-vector along the fiber axis. Solving for Eq. (1) now requires solving for a much smaller (∼ 105 × 105) matrix, which can be accomplished using standard linear algebra packages. This simulation technique is thus very similar to a two-dimensional finite-difference frequency-domain solver.
Fig. 1 (a) For computational feasibility, the resonant fields of a planar nanophotonic device are approximated as those of a truncated holey fiber with identical dielectric structure. (b) An example of our approximation using an L3 photonic crystal resonator. Most of the characteristics of the full three-dimensional field at the center of the slab appear in the approximate solution.
As an example, in Fig. 1(b), we compare the 2.5-dimensional electric field against the three-dimensional simulated electric field (at the central plane of the slab), given by a standard finite-difference time-domain (FDTD) software package, for an L3 photonic crystal resonator. We see that with the appropriate choice of β, which we obtain by fitting a sinusoid to the out-of-plane decay in the three-dimensional FDTD solution, our approximation captures most of the characteristics of the three-dimensional field at the center plane of the slab.
In this approximation, converting from a 2.5-dimensional structure to a full three-dimensional structure only requires the correct choice of slab thickness, t, since the in-plane values of ε remain the same. To estimate t for the full three-dimensional device, we treat the resonant mode as a guided mode in a slab waveguide surrounded by a cladding of n = 1. The expression for t is then the thickness of the slab for the lowest-order TE mode of a slab waveguide with effective refractive index n eff [3]:
Here, α is the wave-vector which determines the evanescent decay of the slab waveguide mode into air, while kx is an approximated in-plane k-vector. To determine n eff, the effective refractive index, we use the following approximation, In the above formulation and for the remainder of the manuscript, the norm (|| · ||) refers to the 2-norm over the entire grid and is equivalent to the spatial integral (∫ | · |2 dr)1/2.3. Inverse design method
3.1. Original formulation
In our inverse design method, one proceeds by specifying either the resulting field or its characteristics and then computing a structure which will produce such a field. In this work, we chose to specify the desired field characteristics, namely, small mode-volume and large quality factor. The inverse design problem can then be formulated as,
where ε and E are the permittivity and electric vector field respectively, defined within the 2.5-dimensional approximation and discretized along a standard Yee grid [4] with periodic boundary conditions. ε was constrained to be isotropic by only varying εz in each cell, and then determining εx and εy via spatial interpolation. FT is the two-dimensional Fourier transform.In this formulation, the minimization objective (Eq. (8)) is the mode volume, which we desire to be as small as possible. Similarily, Eq. (11) is a constraint on the electric field to produce a large quality factor by forbidding the existence of any Fourier components within the light cone. Eqs. (9) and (10) are physical constraints that ε and E must satisfy, that is, the wave equation for the 2.5-dimensional fiber mode and the transversality condition. Lastly, Eq. (12) is a binary constraint on epsilon, denoting that we only want the structure to be composed of air or silicon.
Although our 2.5-dimensional approximation has reduced the size and complexity of the matrices found in the above formulation, the form of the problem in Eq. (8)–(12) is actually incredibly difficult to solve, if for no other reason that Eq. (8), (9), and (12) are non-convex [5] for joint minimization on both ε and E. To make matters worse, the binary constraint in Eq. (12) generally results in a problem which is NP-hard [5]. For these reasons, our strategy is to employ an alternating directions scheme [6], in which we break Eq. (8)–(12) into two separate, tractable sub-problems which we then solve iteratively.
3.2. Alternating directions: field optimization sub-problem
The first sub-problem in the alternating directions scheme involves optimizing E while holding ε constant,
which is a quadratic problem with linear equality constraints, and is easily solved using a standard factor-solve method for sparse matrices [7].In this sub-problem, the most significant modification to the original formulation in Eq. (8)–(12) is that the constraint in Eq. (9) has been relaxed and moved into the minimization objective (Eq. (13)). We will denote this term (||∇ × ∇ × E – μω 2 εE||) as the physics residual.
At the same time, the term for the mode volume from Eq. (8) has also been modified in Eq. (13). We denote this term (||εE 2||) as the design objective. Note that although the denominator in the original formulation has been removed in order to make the term convex, the present formulation is still equivalent because we fix the E-field in the center of the structure (Eq. (14)), where we want the maximum to occur. Also note that the constraint in Eq. (14) is crucial in order to avoid the trivial solution E = 0.
Lastly, the η coefficient in Eq. (13) allows us to trade-off between the physics residual and the design objective. Thus, in order to achieve a small mode volume, the initial value of η is large and is subsequently exponentially reduced once per iteration at all points in the grid in order to bring the physics residual to 0. Furthermore, η can also be given a non-constant spatial weighting in order to reduce in-plane losses; this strategy was implemented in the results presented here.
3.3. Alternating directions: structure optimization sub-problem
The second sub-problem in the alternating directions scheme involves optimizing ε while holding E constant,
which is a quadratic problem with linear inequality constraints and is solved using CVX [8], a convex optimization package for Matlab.As in the field optimization sub-problem, the physics residual has been relaxed and placed in the minimization objective (Eq. (17)). However, we have chosen not to include the design objective term (mode volume). This modification was implemented as a heurestic that seemed to improve the aesthetic quality of the resulting structure.
The second major modification from the original formulation is found in the constraint in Eq. (18). Since the combinatorial problem given by the original constraint in Eq. (11) proved to be intractable, the constraint on the values of ε is relaxed to allow for any value in the range from ε air to ε silicon. The sub-problem is now relatively straightforward to solve, although it is very unlikely that one obtains a completely binary structure. In practice various digitization schemes can be implemented [9]; however, such schemes were not needed here since a nearly binary structure was produced fortuitously.
4. Result
Using the formulations for the field optimization and structure optimization sub-problems presented above, along with our 2.5-dimensional approximation, the inverse design of a planar three-dimensional nanophotonic device is performed.
In order to do so, the frequency, ω = 0.16c/Δx, and the out-of-plane wave vector, β = 0.24(Δx)−1, are chosen by the user. Here, Δx is the grid spacing and c is the speed of light in vacuum. Lastly, an initial dielectric structure consisting entirely of silion, ε = ε silicon everywhere on a 160 × 160 grid, is used as a starting point—although other initial structures (e. g. air everywhere or random ε) yield nearly identical structures in this case.
Having already determined the desired characteristics of our resonator (small mode-volume and large quality factor) in the formulation of the two sub-problems, the inverse design proceeds by alternately solving the field optimization and then the structure optimization sub-problems. That is to say, in every iteration of the algorithm, E is updated to be the solution of Eq. (13)–(16) and is then plugged into Eq. (17)–(18), the solution of which becomes the updated value of ε. The subsequent iteration then uses the updated ε variable to re-optimize E once again. And the algorithm continues to proceed in this way until the physics residual (||∇ × ∇ × E – μω 2 εE||) is reduced to an acceptable value.
Figures 2(a) and 2(b) are the resulting values of ε and E (Ey shown only), respectively, after 75 iterations of our inverse design method. The entire inverse design process takes only 10 minutes on a standard desktop computer for the 160 × 160 grid. An animation of the values of ε and Ey at each iteration of the algorithm is included in the supplementary material.
Fig. 2 The (a) dielectric structure, ε, and (b) target field, E (Ey shown only), produced after 75 iterations (10 minutes) on a 160 × 160 grid. The resulting ε is almost completely binary, and relatively smooth (Media 1).
The values of the physics residual and the design objective (mode volume) at each iteration are shown in Fig. 3. Most importantly, we see that the physics residual seems to converge linearly, indicating that the algorithm is reasonably efficient. The primary factor in determining this rate of convergence is the exponential decrease in η (Eq. (13)), since as η decreases, increasing emphasis is placed on minimizing the physics residual in the field sub-problem.
Fig. 3 Value of the physics residual (blue) and design objective, or mode volume, (red) at each iteration. The physics residual seems to exhibit linear convergence, while the mode volume quickly saturates after roughly 25 iterations.
5. Verification
5.1. Accuracy
To evaluate the accuracy of the inverse design method, we compared the field resulting from the iterations (target field) with the actual field of the resulting structure solved by the 2.5-dimensional approximation (2.5D fiber mode) as well as the field of a full three-dimensional FDTD simulation of the resulting structure (full 3D field), as shown in Figs. 4 and 5.
Fig. 4 Comparison (Ey) of (a) the target field from the inverse design method (from Fig. 2), (b) the actual 2.5-dimensional fiber mode, and (c) the field from the full three-dimensional FDTD simulation. The target field matches well with the full three-dimensional field.
Fig. 5 Comparison of the cross sections of the Ey field along the x-axis from the target field (blue), the actual 2.5-dimensional fiber mode (green), and the field from the full three-dimensional FDTD simulation (red). There is some discrepancy between the target field and the full three-dimensional field, but even that is confined to the edges and is only on the order of ∼ 1% of the maximum field amplitude.
As detailed above, the thickness of the corresponding three-dimensional structure is determined by Eq. (4), which yielded a value of 8.16Δx. However, a small sampling of thicknesses in that vicinity resulted in a more accurate thickness of 8.5Δx, which matched the resonant frequency of the target field and 2.5-dimensional fiber mode.
Figure 4 plots the Ey component of the target field, 2.5-dimensional fiber mode and the full three-dimensional FDTD field side-by-side. Figure 5 plots the horizontal cross sections of the magnitudes of the fields on a logarithmic scale. We see that the target field and fiber mode match very well, while there is significantly more discrepancy between the target field and the full three-dimensional field. This is expected with the use of our approximation; however, the error is still relatively small (below 1% of the maximum field strength) and confined mostly to the edges of the structure.
5.2. Performance
The fields from the full three-dimensional FDTD simulation were used to evaluate the performance of the device. The resulting mode volume was 0.32(λ/n)3 and the total quality factor was 7063. The radiative (out-of-plane) quality factor was 8808 and the in-plane quality factor was 35648.
Figure 6 shows the Fourier transforms of the target, fiber and three-dimensional Ey fields taken at the center of the slab. Although there are virtually no field components in the light cone in the case of the target and fiber fields, additional components are present in the full three-dimensional case. This explains why even when leaky radiative components were strictly disallowed in the field optimization sub-problem, the error in the 2.5-dimensional approximation unavoidably introduces some leaky components in the case of the actual three-dimensional structure.
Fig. 6 Comparison of the Fourier transforms of the Ey fields of the (a) the target field from the inverse design method, (b) the actual 2.5-dimensional fiber mode, and (c) the field from the full three-dimensional FDTD simulation. The error in our approximation introduces some small additional Fourier components into the light cone.
Finally, note that the only the Ey field is shown because it is the dominant field component and contributes most heavily to the out-of-plane leakage [10]. Specifically, the Ey component dominates because the symmetry of the structure dictates that it will contain a non-zero DC component (central component in the light cone), unlike the Ex and Ez components.
6. Conclusion
In summary, by using a 2.5-dimension approximation, we demonstrate the inverse design of a three-dimensional nanophotonic resonator. Further development of our method includes applying our inverse method to design three-dimensional devices which support multiple fields at different frequencies. This includes resonant devices such as a multi-wavelength cavity, as well as waveguiding devices such as N-port couplers, multiplexers, and grating couplers.
Acknowledgments
This work has been supported in part by the AFOSR MURI for Complex and Robust On-chip Nanophotonics (Dr. Gernot Pomrenke), grant number FA9550-09-1-0704.
References and links
1. D. A. B. Miller, “Rationale and challenges for optical interconnects to electronic chips,” Proc. IEEE 88, 728–749 (2000). [CrossRef]
2. J. Lu and J. Vuckovic, “Inverse design of nanophotonic structures using complementary convex optimization,” Opt. Express 18, 3793–3804 (2010). [CrossRef] [PubMed]
3. U. Inan and A. Inan, Electromagnetic Waves (Prentice Hall, 2000), p. 296.
4. K. Yee, “Numerical solution of initial boundary value problems involving maxwell’s equations in isotropic media,” IEEE Trans. Antennas Propag. Mag. 14, 302–307 (1966). [CrossRef]
5. S. Boyd and L. Vandenberghe, Convex Optimization (Cambridge University Press, 2004).
6. S. Boyd, N. Parikh, E. Chu, B. Peleato, and J. Eckstein are preparing a manuscript to be called, “Distributed Optimization and Statistical Learning via the Alternating Direction Method of Multipliers,” http://www.stanford.edu/~boyd/papers/distr_opt_stat_learning_admm.html.
7. Y. Chen, T. A. Davis, W. W. Hager, and S. Rajamanickam, “Algorithm 887: CHOLMOD, supernodal sparse Cholesky factorization and update/downdate,” ACM Trans. Math. Software 35(3), 2 (2009).
8. M. Grant and S. Boyd, CVX: Matlab software for disciplined convex programming, version 1.21. http://cvxr.com/cvx, January 2011.
9. D. Englund, I. Fushman, and J. Vuckovic, “General recipe for designing photonic crystal cavities,” Opt. Express 13, 5961–5975 (2005). [CrossRef] [PubMed]
10. J. Vuckovic, M. Loncar, H. Mabuchi, and A. Scherer, “Optimization of Q-factor in photonic crystal microcavities,” IEEE J. Quantum Electron. 38, 850–856 (2002). [CrossRef]
