We present a technique for large-scale optimization of optical microcavities based on the frequency-averaged local density of states (LDOS), which circumvents computational difficulties posed by previous eigenproblem-based formulations and allows us to perform full topology optimization of three-dimensional (3d) leaky cavity modes. We present theoretical results for both 2d and fully 3d computations in which every pixel of the design pattern is a degree of freedom (“topology optimization”), e.g. for lithographic patterning of dielectric slabs in 3d. More importantly, we argue that such optimization techniques can be applied to design cavities for which (unlike silicon-slab single-mode cavities) hand designs are difficult or unavailable, and in particular we design minimal-volume multi-mode cavities (e.g. for nonlinear frequency-conversion applications).
© 2013 Optical Society of America
In this paper, we present a new technique for large-scale optimization of optical microcavities based on the frequency-averaged local density of states (LDOS), which circumvents computational difficulties posed by previous eigenproblem-based formulations and allows us to perform full topology optimization of three-dimensional (3d) leaky cavity modes. Essentially, this technique allows us to minimize modal volume V for a given cavity quality factor Q, or equivalently a given bandwidth, by solving a single complex-frequency scattering problem at each optimization step, which is both computationally easier than solving an eigenproblem and avoids the difficulty of selecting which eigenvalue to optimize. We present proof-of-concept results in both 2d and 3d computations in which every pixel of the design pattern is a degree of freedom (“topology optimization”), e.g. for lithographic patterning of dielectric slabs in 3d. For a design Q of ≈ 104 in 3d for silicon slabs (index n = 3.52), we obtain a modal volume V = 0.06(λ/n)3 (relative to the vacuum wavelength λ), which is 4× smaller than the best comparable-Q volume found in the literature . Here, the focus is on illustrating the technique rather than on designing practical cavities (for which many good designs are already available in the single-mode silicon case), so we do not incorporate regularizations [2–8] to force the optimization to find easily fabricated structures. However, our results are still useful in establishing theoretical bounds, showing that significant room for improvement remains even for silicon cavities. More importantly, we show that such optimization techniques can be applied to design cavities for which (unlike silicon-slab single-mode cavities) hand designs are difficult or unavailable, and in particular we design proof-of-concept minimal-volume multi-mode cavities (e.g. for nonlinear frequency-conversion applications [9, 10]) in Secs. 8.2.3 and 8.3. We find that the optimum doubly degenerate cavity appears to be three-fold symmetric, while two-frequency cavities have more complex shapes depending on the polarization, although further work is needed to obtain manufacturable structures. Theoretical optimization is also useful to investigate the tradeoff between the size of the design region and the maximum attainable Q in slab-like situations (Sec. 8.4), and find it to be roughly exponential. All of these results are enabled by a sequence of mathematical transformations of the original cavity-design problem, depicted in Fig. 1: from the Purcell factor Q/V to the more well-posed problem of minimizing V for a given Q (Sec. 2), from the eigenproblem to a scattering problem via the LDOS (Sec. 3), from a frequency-averaged LDOS to a single complex-frequency scattering problem  via contour integration (Secs. 4–5), and finally from maximizing LDOS to minimizing 1/LDOS in order to circumvent optimization problems that arise for sharply peaked objectives (Sec. 7). Even though the optimization problem is nonlinear and probably non-convex , we obtain similar optima for many initial structures (including vacuum or random pixels), suggesting that the result may be close to a global optimum.
Microcavity design, which seeks to confine a cavity “mode” for a long time (dimensionless lifetime Q) in a small “modal volume” V, has a long history [13, 14], and until recently has been dominated by hand designs, from ring resonators [15, 16] to photonic-crystal slab cavities [17, 18], in which the parameters are manually tweaked to obtain the desired results (or in some cases computer optimization is performed over a handful of parameters [19–21]). These manual designs have been tremendously successful, attaining experimentally verified Q of 106 or more in modal volumes only a few wavelengths in diameter [19–21]. However, fundamental questions remain: how close are the existing designs to the theoretical optima (e.g. the minimal V for any given Q), and do the optimal designs resemble the hand designs or are they entirely different? (As discussed in Sec. 2, the question is not what is the maximum Q or Q/V, because those questions have a trivial theoretical answer: ∞.) Furthermore, while design of single-mode cavities in silicon slabs has been heavily studied and many design heuristics are known, very different designs and grueling effort may be required in radically different circumstances, from new materials to new design goals such as multi-mode cavities for nonlinear optics . Large-scale “topology” optimization, in which the design pattern is completely determined by computational search (often with thousands of free parameters) with little or no a priori information, offers a different route to addressing these questions. Although several authors have developed topology-optimization techniques for microcavity design [22–25], most of that previous work was limited to 2d calculations [22, 23] and none of the previous work fully addressed the design tradeoff between Q and V. Some topology optimization dealt with lossless systems and avoided Q entirely , or maximized a 2d Q without controlling V  (which we argue in Sec. 2 leads to an ill-posed optimization problem). Other work used a “2.5d” heuristic for the cavity radiation loss , which permitted 2d calculations and led to intriguing designs, but limits the attainable Q (to < 104) because there appears to be no systematic way to improve this heuristic loss estimate. One group did perform fully 3d calculations with absorbing boundaries to capture radiation loss , but limited their degrees of freedom to a small region inside a photonic crystal and included an unphysical absorbing material within the cavity itself, leading to a heuristic objective whose calculation does not appear to have a rigorous quantitative relationship to key cavity properties. As we review in Sec. 3, maximizing LDOS corresponds mathematically to maximizing the power radiated by a dipole current source inside the cavity. Naively, this may seem backwards: isn’t the purpose of a cavity to minimize radiated power? However, that intuition stems from the situation in which the energy density inside the cavity is held constant: there is an equation Q = ωU/P relating the cavity Q and frequency ω with the radiated power P and the confined energy U , so that P ∼ 1/Q for fixed U. In our case, however, we are holding the current amplitude in the cavity fixed, in which case the radiated power P ∼ Q (as reviewed in Sec. 3) while one can show from a coupled-mode framework  that U ∼ Q2. Conversely, minimizing radiated power for a fixed current would correspond to minimizing LDOS, which would result in the absence of resonances. (Interestingly, Frei et al.  actually minimized radiated power from a fixed-current dipole, but their introduction of a heuristic absorbing region inside the cavity apparently compensated for this inverted objective and resulted in a resonant mode, albeit a mode optimizing an unclear objective.) On the other hand, our work adapts two crucial ideas from : solving a scattering problem rather than an eigenproblem, and optimizing over a finite bandwidth. Nor did any previous work, to our knowledge, address topology optimization of multi-mode cavities.
2. Eigenproblem formulation
There are two key figures of merit for a resonant mode En(x) of a cavity: quality factor Q and modal volume V. The quality factor Q is a dimensionless lifetime, and 1/Q is a dimensionless decay rate . Mathematically, Q is related to the frequency-domain Maxwell eigenvalue problem:14] is 26]; however, we will circumvent all of these difficulties in this paper by transition to a Green’s-function LDOS approach in Sec. 3 that does not deal explicitly with eigenvalue problems. The modal volume V , defined as Eq. (3) is a standard expression, it should really be viewed as the modal volume in the limit of a lossless cavity. For a lossy cavity, ∫ε(x)|En(x)|2dx does not converge in an infinite open system; a more rigorous treatment is to define the LDOS in terms of the Green’s function as in Sec. 3 and obtain the Purcell factor from the ratio with the LDOS of the homogeneous medium.) The Purcell factor , which describes the enhancement of spontaneous emission rate, can be written as [29,30]
For applications with light-matter interactions (such as lasers, sensors, and nonlinear frequency converters), maximal lifetime Q and minimal modal volume V are desirable . It is therefore tempting to use the Purcell factor in Eq. (4) or Q/V as the figure of merit for cavity optimization. Unfortunately, maximizing Q/V leads to an ill-posed problem, because the maximum of Q/V is ∞; for example, Q/V grows exponentially with radius for a ring resonator . In practice, any optimization in a finite computation cell will obtain a finite Q and V , but the values are then just an artifact of the finite computational domain; in this sense, maximizing Q/V is not well-posed because the solution does not converge as one increases the size of the computational domain.
In practice, however, there is an upper bound on the useful Q for two reasons. First, besides the intrinsic radiation loss (Qrad) in a cavity, there are also radiation losses due to surface roughness (Qroughness) and material absorption (Qabsorption). The total loss rate 1/Qtotal is the sum of these three effects :33], although microdisk and microsphere resonators can achieve higher Q by delocalizing the mode away from the surface [34, 35]. Second, there is another quality factor in the system. For any cavity-based device, the cavity is always intentionally coupled to some channels (e.g., waveguides) to get light in and out. That coupling process will be described by its own lifetime Qcoupling. It turns out that the losses in such a coupled device are proportional to Qcoupling/Qtotal . Once these losses are decreased below the desired loss budget, it often does not matter in practice if one decreases them further.
A better optimization problem might be, instead, to maximize Q/V subject to Q = Q̃, or alternatively:Eq. (1), one could obtain Q and V from eigenvectors En(x) and eigenvalues ωn through Eq. (2) and Eq. (3). Then a natural question to ask is which eigenvalue one should optimize. In practice, one has some design frequency ω̃ given by the application, so one could optimize the eigenvalue closest to ω̃. (Equivalently, thanks to the scale-invariance of Maxwell’s equations , we choose units so that ω̃ = 2π, in which case the main computational choice is the resolution, i.e. the number of pixels per wavelength λ, and in some cases a slab thickness/λ.) However, asking for the mode closest to ω̃ leads to discontinuities: as the structure changes during optimization, the mode closest to ω̃ will tend to hop discontinuously. Although there are some ways to deal with this , the problem becomes worse when one simulates the radiation loss in a patterned dielectric slab, because in this case the finite cell is approximating a continuum of radiation modes above and below the slab. As a result, there are more and more closely spaced modes as the cell size increases. Hence, we want to circumvent this difficulty by adopting a new approach: turning the eigenproblem into a linear scattering problem. This also reduces the computational expense, because we will now require only a single linear solve per structure, whereas finding eigenvalues in the interior of a spectrum (e.g. by the shift-and-invert algorithm) requires many linear solves .
3. LDOS formulation
The well-known Purcell factor (Q/V and variations thereof) is, in fact, only an approximation of a more fundamental quantity, the local density of states (LDOS), which is defined in terms of the Green’s function of the system (and can be related to a “density of modes” per unit frequency per unit volume) . Paradoxically, the exact LDOS is easier to compute than Q/V, because the former involves only a scattering problem (a linear system of equations) whereas the latter involves a non-Hermitian eigenproblem (for a leaky mode in the interior of the spectrum). In this section, we briefly review the definition of the LDOS and how it relates to Q/V, and in the next section we describe how the LDOS can be used as the objective for a well-posed cavity optimization problem.
In particular, we begin with the per-polarization LDOS (partial LDOS), denoted by LDOSj(ω, x′) for a polarization in direction j, which is proportional to the power radiated by a dipole in direction (unit vector) êj at position x′ with a frequency ω, i.e. a current J ∼ êje−iωtδ (x − x′). (The total LDOS is simply ∑j LDOSj, proportional to power radiated by a randomly oriented dipole.) This is a key physical quantity because quantum fluctuations, e.g. spontaneous emission or thermal fluctuations, can be viewed semiclassically as dipole currents, and the LDOS yields the enhancement of these fluctuations by any given geometry. For example, the spontaneous emission rate is proportional to the LDOS [29, 38], and therefore the ratio of LDOS between two systems indicates enhancement or suppression of spontaneous emission. The Purcell factor Q/V is an approximation to this LDOS enhancement for a high-Q microcavity, assuming that the emitter (e.g. atom or quantum dot) is positioned and oriented at the location of peak LDOSj in the cavity [29, 30, 38, 39].29, 37, 38, 40, 41] definition of LDOSj: 29, 38]. (The normalization is irrelevant for optimization or for evaluating LDOS ratios in different systems to obtain enhancement factors.)
In the limit of a low-loss cavity, in which the LDOS is dominated by the contribution of a single pole in Maxwell’s equations (a single “resonant mode”), one can derive the approximation 29] obtains the traditional Purcell factor Eq. (4). In particular, the LDOS in a microcavity resonating a frequency ω̃ (1 + i/2Q̃) is approximately in the form of a Lorentzian peak centered at ω̃ with a bandwidth ω̃/2Q̃ , and the Purcell factor is the enhancement at the peak. This relationship between Q and LDOS bandwidth is the key to obtaining a tractable well-posed optimization problem in the next section.
We propose that LDOSj or its variants can be used as figure of merit for the characterization and optimization of a microcavity. Note that the precise figure of merit should really depend on the application. For example, if we are interested in the spontaneous emission rate for the dipole at a specific position x′ with a specific polarization êj, then LDOSj(ω, x′) is the most relevant figure of merit. On the other hand, if we have a dipole at a specific point x′ with a randomly distributed polarization, then the figure of merit would be ∑j LDOSj(ω, x′) (for the average case) or minj LDOSj(ω, x′) (for the worst case). In nonlinear devices for frequency conversion (in particular, second harmonic generation), a more relevant figure of merit might be minn LDOSj(ωn, x′), where ωn are different frequencies of interest. Instead of at a single point x′, if the dipoles of interest are distributed [with probability density function s(x)] in a region V with polarization êj, the most relevant figure of merit in this case is∫V LDOSj(ω, x)s(x)dx. Depending on the applications for enhancement or inhibition, we should maximize or minimize the figure of merit correspondingly. Many other variations could be devised.
For all the above mentioned applications, LDOSj is the basic building block. We therefore focus on this case: we maximize the spontaneous emission rate for a dipole at a point x′ with given polarization êj. For simplicity, from now on, we shall omit the explicitj and x′ dependence from LDOSj(ω, x′), and simply write LDOS(ω) to denote the figure of merit given in Eq. (9). In Sec. 8.2, however, we will also consider a case where the dipole polarization is randomly distributed, and in Secs. 8.2.3 and 8.3 we consider multi-mode cavities.
4. Frequency-averaged LDOS
In previous section, we proposed that one way of attacking the problem of microcavity design is to maximize the LDOS. However, simply maximizing LDOS(ω) in Eq. (9) is still ill-posed as in Sec 2, being equivalent to Q/V. As explained in Sec. 2, a well-posed formulation could be obtained by specifying a desired cavity lifetime Q̃. In terms of LDOS, this is equivalent to specifying an upper bound ω̃/2Q̃ on the bandwidth Γ̃. For computational convenience, we instead maximize the average LDOS over this finite bandwidth Γ̃:38], so that the integral Eq. (11) is determined mainly by the 1/V factor as soon as the resonance is narrower than Γ̃ (Q ≫ Q̃). The key point is that maximizing the bandwidth-averaged LDOS regularizes the optimization problem to eliminate the possibility of diverging-Q solutions as the computational cell size increases.
The main remaining question is how to efficiently compute the average LDOS of Eq. (11). The most straightforward approach would be to apply a standard numerical-integration procedure , which would involve evaluating the LDOS (i.e., solving a scattering problem) for many separate frequencies over the bandwidth Γ̃. (This is similar in spirit to the multi-frequency optimization approach of .) However, as described in the next two sections, we can exploit techniques from complex analysis to evaluate the exact LDOS integral by solving only a single scattering problem at a complex frequency ω̃ + iΓ̃. The key to this technique is that the LDOS derives from a causal Green’s function, as reviewed in Sec. 4.1. This allows us to perform a contour integration, with an appropriate choice of W(ω), in order to obtain the integral as described in Sec. 4.2. As explained in Sec. 4.3, we can reinterpret any complex-frequency scattering problem as a real-frequency scattering problem with complex materials (i.e., absorption). Finally, in Sec 5 we discuss convenient choices of the window function W(ω) that reduce the contour integral to a single complex-ω scattering problem. (We previously applied a very similar application of causality and complex analysis to analysis of electromagnetic cloaking bandwidth , and related ideas can be found in quantum field theory [44, 45].)
4.1. Causality and analyticity
Before we proceed, we first define a function f (ω), which is a complex version of LDOS(ω):Eq. (9), it is clear that LDOS(ω) = Re[f(ω)]. [From now on, we will omit the x′ dependency of f (ω, x′).] Note that the operator ℳ (ε, ω) given in Eq. (8) is a linear operator relating the electric field E(x, ω) to the (time-harmonic) input electric current J(x) at a given frequency ω. Causality [the electric field E comes after (not before) the current J] implies that E(x, ω) is analytic in the upper-half complex-ω plane . Therefore, f (ω) is also analytic in the upper-half complex-ω plane.
4.2. Contour integration
In this section, we are going to compute the mean LDOS by exploiting the analyticity of f (ω), via:40,42,47,48]. Now we want to complete our integration contour (Fig. 2) in the upper-half plane and evaluate L by residue theorem  49]. In Sec. 5, we will choose W(ω) so that these poles and residues are easy to evaluate. Furthermore, we can choose the weight function W(ω) so that it decays faster than 1/|ω|3 for large ω, in which case the contribution from arc A3 will be zero since LDOS(ω) is proportional to ω2 in 3d (or ω in 2d) free space and f (ω)W(ω) will decay faster than 1/|ω| on arc A3. We can obtain the contribution from arc A2 by evaluating the residue due to the simple pole of f (ω) at ω = 0 Eq. (13), L is just the path integral along the path A1 . Therefore, we have
4.3. From complex frequency to material absorption
To compute the residue at the complex poles, we need to solve the scattering problem at complex frequencies. More precisely, the scattering problem Eq. (8) at complex frequency ω + iΓ can be written as50].) In particular, adding a positive imaginary part to ω (for ω̄k in the upper-half plane ) corresponds to a positive imaginary part in ε̃(x) and μ̃(x), which corresponds (with our e−iωt convention) to an absorption loss.
5. Possible window functions
In this section, we discuss two convenient window functions: a simple Lorentzian and the square of a Lorentzian. (A third possibility, the difference of two Lorentzians is discussed in .)
5.1. A simple Lorentzian
The simplest window function, with only a simple pole in the upper-half plane, is a Lorentzian centered at ω̃ with half-width Γ̃. The frequency-average LDOS against this weight isEq. (17) at a single complex-frequency ω̃ + iΓ̃. L1 is a perfectly finite, well-defined quantity in a discretized simulation with a finite spatial resolution (finite grid).
In combination with Sec. 4.3, the objective Eq. (19) has a simple interpretation that coincides with the discussion in Sec. 2. Our frequency-averaged LDOS objective is equivalent to maximizing the LDOS at a single frequency, i.e. the Purcell factor Q/V of a cavity, but a cavity in which an absorption loss has been added to otherwise lossless materials. In particular, the absorption loss that arises from multiplying ε(x) and μ(x) by 1 + i/2Q̃ can be seen from perturbation theory  (for Q̃ ≫ 1) to yield an absorption lifetime of exactly Qabsorption = Q̃. From Eq. (5), this means that Qtotal ≤ Q̃, and that increasing Qrad ≫ Q̃ will have little effect on the Purcell factor Qtot/V. Hence, optimizing the frequency-averaged objective will tend to increase Qrad until it is > Q̃, and after that will mainly try to decrease V, similar in spirit to Eq. (6).
However, a careful examination reveals that this simple average does not converge as the resolution increases. There are two equivalent ways to understand this. First, in a continuous medium, the integral does not converge because the window function decays like 1/|ω|2 while LDOS(ω) behaves like |ω| (in 2d free space) or |ω|2 (in 3d free space) for large |ω|. (For finite spatial resolution, there is an upper frequency cutoff that eliminates this divergence.) Second, from the relationship between the complex-frequency scattering and lossy material discussed in the previous section, we know that the residue Re[f (ω̃ + iΓ̃)] is actually the power emitted by a dipole in lossy material, which is the sum of the power radiating to the outside of the cavity and the power absorbed by the lossy material in the cavity . It is known that this absorbed power is infinite because E(x) diverges as 1/r3 in the neighborhood of the dipole (in 3d free space) [40, 47, 52, 53]. (In a lossless medium, only Im[E(x)] diverges as 1/r3, so LDOS ∼ Re[E(x)] is finite.) In discretized space, the Green’s function is finite and diverges as (resolution)3 in 3d. To avoid this singularity, we need to choose window functions which decay faster than |ω|3 at large |ω|. Two natural candidates are the difference of two Lorentzians  and the square of a Lorentzian.
5.2. Square of a Lorentzian
To ensure that the W(ω) decays faster than 1/|ω|3, we propose the window function,Eq. (16) and Eq. (20) Eq. (12) and Eq. (22), it is clear that both f (ω̃ +iΓ̃) and f′(ω̃ +iΓ̃) can be obtained from a single scattering solution E(x, ω̃ + iΓ̃) (see appendix 10) and Eq. (17) at a complex frequency ω̃ + iΓ̃.
We know that Eq. (20) gives a finite average LDOS because it decays fast enough with ω, but it is interesting to also consider how it fixes the divergence from the second viewpoint in Sec. 5.1 (that of the infinite power absorption from a dipole in a lossy medium). The explanation is essentially that the second term in Eq. (23) is roughly a subtraction of the divergent absorbed power from the first term: Γ̃ε is ω Im(ε̃) from Sec. 5.1, and ω Im(ε̃)|E|2 is absorbed power . (A subtlety arises from having +ETE, rather than −|E|2, but the 1/r3 divergence at r → 0 should be dominant in Im(E) for small Γ̃ so one should have ETE ≈ −Im(ET)Im(E) ≈ −|E|2 as r → 0.)
Since the role of the second term in Eq. (23) is essentially to subtract off the divergent absorbed power in lossy ε̃ medium, and this divergence comes from the 1/r3 field divergence that is independent of geometry (the scattered field from the surrounding geometry is finite at r = 0), one might expect that the second term in Eq. (23) plays little role in geometry optimization at a fixed resolution. Indeed, we find in numerical experiments that the optimizations with and without the second terms in Eq. (23) for the 2d TE case (discussed in Sec. 8.2) discover similar structures. Therefore, in Sec. 8 we optimize the simpler single-Lorentzian objective of section 5.1, although Eq. (23) is computationally feasible if it becomes necessary.
6. A preliminary formulation
Now we have a preliminary formulation for our cavity optimization in terms of the frequency-averaged local density of states:Eq. (17) once. If we choose the window function W(ω) from Eq. (20), then the problem can be reformulated as Eq. (25), we need
To speed up the optimization, we must also have the gradient of the objective with respect to the design parameter εk (the dielectric constant at each “pixel” k). Applying the standard adjoint technique , only one more linear system with the same operator ℳ̃(ε, ω̃) but a different source term need be solved to obtain the gradient. We provide the detailed calculation in appendix 10 and summarize the procedure here:
- Take the real part of ∂ℒ/∂εk to obtain ∂L/∂εk.
Note that the scattering operator Eq. (26) in the sensitivity analysis is the same as the operator in the objective evaluation. We can take advantage of this by reusing the information (e.g., the preconditioner or LU factorization) from the solution of Eq. (17). We will discuss this in detail in Sec. 7.1.
7. Numerical scheme for cavity optimization
In this section, we discuss the numerical implementation for our frequency-averaged LDOS formulation given in Sec. 6. In order to solve this PDE-constrained optimization problem computationally, we need fast and efficient implementations for objective evaluation, gradient evaluation, and optimization.
7.1. Objective and gradient evaluation
As we discussed in Sec. 6, evaluating the objective LDOS involves solving the scattering problem Eq. (17). We can apply any standard frequency-domain solver technique to this problem (e.g., finite-difference, finite-element, or boundary-element methods). Here, we simply adopt the finite-difference approach [55–57]. If we impose mirror symmetry planes in the system, we can obtain an 8-fold reduction in the number of unknowns (see Sec 8.5.1).
For the finite-difference frequency-domain (FDFD) method, the most robust solution technique is a sparse-direct solver, which is excellent in 2d, but expensive (in both memory and time) in 3d . In contrast to direct solvers, iterative solvers (e.g., GMRES or BiCGStab) work quite well if one has a good preconditioner . Here, we combine both of these techniques. During the optimization, we re-solve many times for slightly different structures. Therefore, we can use sparse-direct factorization from one step as a preconditioner for iterative solvers in many subsequent steps . We implemented the FDFD solver with the parallel sparse-matrix library PETSc [59–61] and the parallel sparse-direct solver PaStiX .
7.2. Optimization scheme
We use a free-software implementation  of standard gradient-based optimization algorithms, and we find that low-storage BFGS  and conservative convex separable approximations (CCSA)  work equally well. However, we find that one additional mathematical transformation is required in order to optimize high-Q cavities: instead of maximizing L, we minimize 1/L.
If we attempt to maximize L directly as in Eq. (24), we typically find that the convergence of any standard optimization algorithm slows to a halt for Q ≳ 1000. The reason for this is that, for high-Q cavities, the objective function becomes a “narrow ridge” (a sharply-peaked function along some low-dimensional manifold in the parameter space), with a large second derivative (∼ Q3, as shown in ) in the direction perpendicular to the ridge, and it is well known that optimization along narrow ridges is problematic . Most algorithms tend to “zigzag” slowly along the ridge, and even quasi-Newton methods like BFGS break down when the second derivative is so large that the Hessian matrix becomes ill-conditioned. However, because the LDOS is strictly positive, there is a simple solution: maximizing L is equivalent to minimizing 1/L, and the reciprocal of a sharp peak is a shallow valley, so we find that the 1/L objective avoids most of the problems of slow convergence. (For example, in the 2d TM optimization to be discussed in Sec. 8.1, maximization of L makes no progress if the initial structure is a high-Q photonic-crystal cavity, whereas minimization of 1/L converges to a significantly improved structure as described below.)
In practice, there is another useful technique: a successive-refinement strategy (somewhat analogous to [67,68]). We found that gradually increasing the specified Q̃ (decreasing the bandwidth 1/Q̃) tends to reduce the likelihood that the optimization becomes trapped in a poor local optimum close to the starting structure, and usually produces a much better result at the highest Q̃. That is, we optimize for Q̃ = 10, then use that result as a starting point for optimization at Q̃ = 100, and so on.
8. Results for cavity optimization
In this section, we present some 2d and 3d cavity-optimization results, summarized as follows. We start with high-resolution 2d cases, and run simulations with different initial guesses (vacuum, photonic crystal with a defect, and random structures) and different dipole polarizations (TM, TE, and random). In the region to be optimized, we allow the dielectric constant ε ∈ [1, 12.4] at each pixel to be one degree of freedom [Fig. 3(a)]. For the 2d TE case, the optimization discovers similar structures for maximizing the spontaneous emission rate of a specific dipole polarization and a randomly polarized dipole. However, optimizing the worst case of a randomly polarized dipole finds a three-fold symmetric structure (Sec. 8.2.3), while two-frequency cavity optimizations yield more complex patterns depending on the polarization (Sec. 8.3). In another 2d scenario, to obtain the Q versus V tradeoff analogous to 3d slabs, we limit the degrees of freedom in one direction and choose a thin strip, instead of a square, as the region for optimization [Fig. 3(b)]. As the degrees of freedom increase, the radiation Q first increases and then becomes saturated, limited by the numerical precision in the computation. Finally, we ran a full 3d optimization on a supercomputer and obtained a structure with extremely small mode volume V = 0.06(λ/n)3, ( ). In the following, we are minimizing 1/L (the inverse of the frequency-averaged LDOS) and variations thereof, as described in the previous sections, but for simplicity we describe this below as “maximizing LDOS.”
Here, the focus is on illustrating the formulation and on establishing theoretical bounds. Further regularization techniques are generally required to obtain easily manufacturable structures, as discussed in Sec. 9.
8.1. 2D TM case
In this section, we maximize the LDOS for a dipole with TM polarization (E out of plane) in a 2d setting [Fig. 3(a)]. One possible starting point is a photonic crystal with a defect , like the one shown in Fig. 4(a). This is a periodic arrangement (periodicity a) of dielectric silicon rods (radius 0.2a and permittivity ε = 12.4) with one defect rod at the center (radius 0.1a). The defect TM mode is at frequency 0.32(2π/a), with quality factor Q = 1.41×108 and mode volume V = 0.097(λ/n)2. With this structure as an initial guess, we run the optimization and obtain an entirely different nested-ring structure [Fig. 4(b)] with quality factor Q = 1.01 × 1010 and mode volume V = 0.075(λ/n)2. Clearly, the optimization itself discovers a “radially periodic” structure, reminiscent of a Bragg onion . We also run the optimization with vacuum as initial guess and obtain similar structure (Fig. 5) with Q = 1.30 × 109 and V = 0.075(λ/n)2.
In these two optimizations, we gradually increase the specified Q̃ (we decrease the bandwidth 1/Q̃) from 10 to 105. The optimization at Q̃ = 10 actually gives a high Q cavity (almost the same radiation Q) with the resonance at about 1.003ω̃. The optimizations at higher Q̃ simply tune this structure so that the resonant frequency becomes much closer to ω̃. In a two-dimensional situation such as this, there is no intrinsic Q versus V tradeoff, unlike in 3d slabs , because Q → ∞ for a finite modal volume V as the number of layers in a Bragg onion increases . So, the optimization in such a 2d case is mainly minimizing V, while the Q is bounded only by the computational-cell size.
Note that we allow the dielectric permittivity of each pixel to vary continuously from εmin = 1.0 to εmax = 12.4, but almost all the pixels (except a few at the interfaces) are at either εmin or εmax in the optimized structure. This phenomenon (reminiscent of “bang-bang” solutions in control theory ) had also been observed empirically in other cavity-related optimization work [3,71,72]. There has been some recent progress in proving theoretically that this is the expected solution: Osting and Weinstein  recently analyzed optimization problems for scalar waves, and showed that maximizing an energy confinement time over the permittivity at every point in space generally leads to a solution in which the permittivity is either the maximum or the minimum allowed value at every point, excepting a set of measure zero (e.g., at the interfaces between regions) in the limit as the resolution goes to infinity. However, we also obtain counterexamples for other types of cavity optimization, e.g. in the doubly-degenerate cavity design of Sec. 8.2.3, where substantial regions converge to intermediate values. Fortunately, in cases where this occurs, there are a variety of techniques to obtain binary solutions as discussed in Sec. 9.
8.2. 2D TE case
In this section, we consider the 2d TE polarization, for which (unlike the TM polarization of Sec. 8.1) the objective function breaks rotational symmetry and we do not expect similar Bragg-onion solutions. We will start with maximizing LDOS for a fixed êx polarization in 8.2.1. For a randomly polarized dipole, we consider two cases: maximizing the average over all polarizations in 8.2.2 and maximizing the minimum (worst case) over all polarizations in 8.2.3.
8.2.1. Fixed dipole polarization: max LDOS(ω̃; êx)
For the 2d TE polarization, let us first look at the case where the dipole is polarized in the êx direction. In other words, we want to maximize LDOS(ω̃; êx). From a vacuum initial guess, the optimization discovers the structure shown in Fig. 6. This structure has quality factor Q=5.16×108 and mode volume V = 0.092(λ/n)2. [Again, the Q̃ = 10 gives an equally high-Q cavity with resonant frequency at 1.0007ω̃. The optimizations at higher Q̃ = 102 to 105 simply tune the resonant frequency to ω̃.] For a random initial guess (uniform random pixels), we also obtain the same structure as the one from a vacuum initial guess, which suggests that the result may be close to a global optimum.
8.2.2. Randomly polarized dipole: max meanj LDOS(ω̃; êj)
Now we want to study the case in which the dipole is randomly polarized in the plane and the objective is the mean LDOS. One might hope that the optimization will find a structure that resonates for both polarizations at the same frequency, and hence by linearity will resonate for all in-plane polarizations—this corresponds to the requirement that the microcavity be doubly degenerate, and many symmetry groups besides circular symmetry can support double degeneracies (such as three-fold, four-fold, or six-fold symmetrical structures ). However, we find that this is not the case: the optimization finds a singly resonant cavity that enhances LDOS for one polarization (chosen “randomly” depending on the initial structure) at the expense of the other polarization. As explained below, this suggests that the LDOS of the best single-mode TE cavity is more than twice the LDOS of the best doubly degenerate TE cavity. We performed the mean-LDOS optimization as follows. It is easy to show [29, problem 8.6] that maximizing the LDOS for a random polarization by averaging all polarizations is equivalent to maximizing the mean from the êx and êy polarizations, namely [LDOS(ω; êx) + LDOS(ω; êy)]/2. For this new objective, we ran many different simulations with different pseudo random initial guesses (each pixel is randomly chosen between εmin and εmax). About half gave the structure optimizing the êx polarization (Fig. 6), while the other half gave the 90°-rotated structure optimizing the êy polarization. From these results, it seems that the optimization, instead of favoring both êx and êy polarization simultaneously, simply randomly picks one direction and optimizes it. That is, there is a spontaneous symmetry breaking: it is better to optimize one polarization at the expense of the other than to try to obtain a doubly degenerate cavity that resonates for both polarizations. By selecting only one polarization to enhance, these structures sacrifice a factor of 2 in the mean LDOS, which is why we conclude that the best single-mode LDOS is at least twice as big as the best two-mode LDOS.
8.2.3. Optimization for a randomly polarized dipole: max minj LDOS(ω̃; êj)
In the previous section, we showed that maximizing the mean LDOS over all in-plane polarizations was equivalent to maximizing the LDOS for a single polarization at the expense of the orthogonal polarization. In order to obtain a structure that enhances all polarizations equally (via a doubly degenerate resonance), we consider a different objective: we maximize the minimum LDOS over all in-plane polarizations (rather than the mean). The result of this optimization is a three-fold symmetric microcavity shown in Fig. 7(a), which is the smallest symmetry group that supports a (non-accidental) doubly degenerate mode  shown in Fig. 7(b–c). (The rectangular FDFD grid breaks the three-fold symmetry, but the structure converges to true three-fold symmetry with an exact degeneracy as the resolution is increased.) Intuitively, larger symmetry groups have fewer degrees of freedom for optimization (the Cmv symmetry group  with m-fold rotational symmetry plus reflections leaves only a wedge of angle π/m in the degrees of freedom), so the optimization is picking the smallest possible symmetry group in order to maximize the number of degrees of freedom available for maximizing LDOS.
As predicted in Sec. 8.2.2, the effective modal volume V of this doubly degenerate cavity is 0.34(λ/n)2 > 2 times the volume of the single-mode cavity from Sec. 8.2.1.
Technically, maximizing the minimum LDOS over all polarizations is significantly harder than maximizing the mean LDOS. Unlike the mean LDOS, it is not sufficient to consider only the êx and êy polarizations, and in principle we must solve for the LDOS at all angles. The simplest approach is to merely sample the set of polarizations to compute the LDOS at many discrete angles, and then to maximize the minimum LDOS) over this set. Although we tried sampling up to 20 discrete angles, we found that it was sufficient (obtaining the same structure) to sample only three angles (0, 60, and 120). Furthermore, the minimum LDOS over several angles is no longer a differentiable function, but we circumvent that problem by the standard transformation  of introducing a “dummy” variable t and solving maxt subject to the non-linear constraints t ≤ LDOSj at each angle j. (The resulting nonlinear-programming problem was solved via the CCSA algorithm .)
8.3. 2D optimization for different frequencies: max minnLDOS(ω̃n; êj)
In nonlinear devices, e.g. for nonlinear frequency conversion, it is often desirable to design a cavity which resonates at multiple distinct frequencies [9, 10, 76–79], leading to a challenging multi-frequency cavity-design problem if one desires a small-volume cavity [80, 81]. The simplest approach to designing such a multi-frequency cavity with our optimization techniques is to maximize the minimum LDOS for sources at two (or more) frequencies, ω̃1 and ω̃2. For example, we begin by considering ω̃2 = 2ω̃1 (coupling TM to TM or TE to TE), the desired relationship for intra-cavity second-harmonic generation (SHG) [79, 81], with results shown in Figs. 8(a) and 8(c). These results exhibit Bragg-onion-like structures similar to the single-frequency optimization in the previous sections, which is not surprising considering that Bragg mirrors tend to have gaps at integer frequency multiples  and so the same mirror can confine both a fundamental and harmonic frequency. A more challenging case for optimization, therefore, is to design frequencies that are not integer multiples, and for illustration purposes we considered ω̃2 = 1.5ω̃1, which results in the more complicated structure shown in Figs. 8(b) and 8(d).
8.4. 2D TE thin strip case
In previous 2d TM and TE polarization cases, we expect and obtain almost no Q versus V trade off since the cavity can be surrounded by complete photonic bandgap or a Bragg onion, with Q only limited by the computational-cell size. In a 2d setting, to get a Q versus V trade off analogous to 3d slabs, we need to limit the degrees of freedom in one direction in order to force the possibility of radiation loss. In this section, we choose the region for degrees of freedom to be a thin strip [Fig. 3(b)].
Although we expect cavities in a finite-thickness strip to have intrinsic radiation loss in the direction transverse to the strip, we also expect that it should be possible to make this radiation loss arbitrarily small at the expense of modal volume. For example, if one starts with a periodic waveguide structure  and introduces a defect in the periodicity, a resonant mode can be trapped, and the Q can be made arbitrarily large by gradually tapering from the defect to the periodic structure [14,82–84]. To explore this tradeoff, we limit the modal volume by considering a finite-length d × 1λ strip of degrees of freedom. For each value of d (d = 1, 2, 3, 5λ), we plot the obtained radiation lifetime Qrad as a function of the requested bandwidth Q̃, and the result is shown in Fig. 9(a). For any given d, as Q̃ is increased then Qrad is increased as well (as explained in Sec. 4), until it saturates at a maximum determined by two factors. First, the maximum Qrad is limited by d: a longer d gives more degrees of freedom and increases the maximum possible Qrad, as expected. Second, around Qrad ≈ 107 the improvement becomes limited by the numerical precision, which prevents the optimization from making progress even though higher Qrad should theoretically be possible. For d = 5λ, the resulting structure is shown in Fig. 9(b). [Note that the data in Fig. 9(a) are from the optimization result for the objective ε(x′)LDOS(ω, x′).]
8.5. 3D case
In this section, we present the optimization results for the TE-like polarization in the 3d slab setting (Fig. 10). We also briefly discuss some variations on the optimized structure and its comparison with an air-slot cavity.
8.5.1. 3D slab optimization
With the optimization tools developed in the previous sections, we run large-scale simulations on 3d slab case (with an in-plane dipole source, coupling to even-symmetry “TE-like”  resonances). Here we choose the dimensions of the slab to be 3λ-3λ-0.19λ, where the thickness is 0.19λ = 0.67(λ/n). A sketch of the physical model is shown in Fig. 10(a), and the real computation domain (with mirror-symmetry reductions) is illustrated in Fig. 10(b). The optimization discovers a structure (Fig. 11) with quality factor Q = 30000 and extremely small mode volume V = 0.06(λ/n)3. [This result is obtained from optimizations with absorption Q̃ gradually increasing from 10 to 104. The optimization discovers structures with radiation Q = 1.18×104 at Q̃ = 100, with radiation Q = 2.55 × 104 at Q̃ = 1000, and with radiation Q = 2.98 × 104 at Q̃ = 104.]
A comparison with other large- or small-scale optimization work, such as 2.5d optimization , L3-type cavity  and H0-type cavity  optimization are given in table 1. Clearly, the optimization was able to achieve four times smaller mode volume than the smallest mode volume (for Q within one order of magnitude of ours) we found in the literature . (Note that L3-type  and H0-type  cavities are designed from small-scale optimization and can be fabricated; while both 2.5d optimization  and our results are purely theoretical and computational investigation from a large-scale optimization perspective.)
8.5.2. Sensitivity to small features
In the optimized structure, there are some filament-like structures which would be difficult to fabricate at infrared scales. As a first assessment of the importance of these tiny features, we manually remove the filaments and obtain a structure (Fig. 12) with Q = 10000 and roughly same V. Instead of post-processing the structure in this way, which has the disadvantage of no longer being a local optimum, future work should consider suppress these tiny features by imposing some explicit constraints during optimization or by some regularization and projection  as further discussed in Sec. 9.
8.5.3. Comparison with air-slot cavity
All these cavities listed in table 1 are dielectric cavities. In other words, the centers of these cavities are high-dielectric materials (Si and GaAs) and these cavities are useful for dipoles/emitters lying in these materials. It is also reflected in our unit of mode volume. For example, the mode volume of the cavity we obtain is V = 0.06(λ/n)3 = 0.06(λ/nSi)3.
It is known that air-slot cavities [85–88] can have extremely small volumes. For example, Nomura  reported an air-slot cavity with Q = 4.8 × 106 and V = 0.015(λ/nair)3. Although 0.015 is smaller than 0.06, these two kinds of cavities are not comparable in two ways. First, the two mode volumes are in different units (λ/nair)3 versus (λ/nSi)3. In our units, their modal volume is 0.65(λ/nSi)3. Second, these two types of cavities are for different applications: air-slot cavities are useful for emitters lying in air, while the semiconductor-based cavities are designed for emitters lying in Si and GaAs.
If the application is for emitters lying in air, in theory, we can also introduce an infinitesimal air-slot at the center, oriented perpendicular to the electric field, into our structure. As discussed in , after the introduction of an air-slot, the mode volume decreases by a factor of (nSi/nair)2 without changing Q. In our case, this factor is about 12.4, and the new mode volume is 4.8 × 10−3(λ/nSi)3 = 1.1×10−4(λ/nair)3. Because the resolution we used (46-pixel per wavelength in air) is not that high, the optimization discovers a dielectric cavity, instead of one with air-slot type. In future work, one could run the optimizations with high resolutions (at least in 2d cases) to investigate whether air-slot structures can be discovered; at such high resolution, the regularized figure of merit of Eq. (21) discussed in Sec. 5.2 becomes essential.
9. Manufacturability: Projection and Regularization
The primary purpose of this paper is to improve the mathematical formulation of the microcavity problem to make it amenable to large-scale topology optimization. It is well known, however, that such topology optimization can sometimes lead to non-manufacturable designs, due to three problems: regions of intermediate ε values, fine features (such as the “hairs” in Fig. 11), and extreme sensitivity to variations in the design parameters. Fortunately, there are a variety of techniques to correct these difficulties, which could easily be combined with our LDOS formulation, and we briefly review these projection, regularization and robustification techniques here for the benefit of future work.
In cases where the optimization does not lead to “bang-bang” solutions, i.e. where there are large regions of intermediate values that do not correspond to available materials, one can use a variety of penalization techniques that add a penalty (increased as needed) to the objective for intermediate ε values [3, 89]. Another possibility is to use a level-set method [4, 90, 91], which guarantee binary solutions (except possibly on a set of measure zero, since intermediate values are typically used at the level set boundary in order to ensure continuity). A third possibility is the SIMP (Solid Isotropic Material with Penalization) , a re-parameterization of the problem that is used in conjunction with the filtering techniques described below in order to project the solution towards a bang-bang result.
To eliminate small features, a common solution is to apply a smoothing filter to the degrees of freedom (combined with a smoothed Heaviside projection to transform back towards bang-bang designs) [2, 5–8, 92–95]. This is often viewed as a regularization to ensure that the problem is well-posed in the sense that the optimum should converge with resolution (rather than yielding finer and finer features) . Another approach that can eliminate small features as well as designs with extreme parameter sensitivity is “robust” optimization, in which one typically optimizes a worst-case design over a set of parameter uncertainties (e.g. a small “noise” added to each pixel) [67,96,97], which has been shown in some cases to also eliminate small features in topology optimization for photonics [68, 98, 99].
In this paper, we presented a novel formulation for large-scale optimization of optical cavities via frequency-averaged LDOS. With this formulation, we obtained various 2d and 3d cavity-optimization results for different applications. Our results show that several times smaller modal volume is theoretically possible for silicon cavities compared to previous work, without sacrificing Q. Many other possibilities present themselves for future work. First, the slab thickness in the 3d optimization (Sec. 8.5) is fixed, but one can further extend our approach by allowing the slab thickness to be one additional degree of freedom to improve the figure of merit. Second, although we focused here on silicon-based cavities, one can easily apply our approach to other materials (such as lower-index diamond for visible frequencies [100,101]). Since silicon cavity design has been so heavily studied and many adequate designs are already known, it is especially in the context of new materials systems that optimization can be valuable. Third, instead of maximizing LDOS, one can maximize the frequency-averaged total DOS over the whole computational cell, which is known to be related to the bounds of the light trapping in photovoltaic problems [102, 103]. Fourth, one can extend our work on multi-frequency optimization (Sec. 8.3) by defining a more specialized figure of merit for nonlinear frequency conversion that includes an overlap factor [79, 81] for the two modes.
Appendix A: Computation of f′(ω, x′)
In this section, we compute f′(ω, x′), the differentiation ∂/∂ω of f (ω, x′) from Eq. (12), to show that the residue of Eq. (21) can be evaluated without any additional matrix solves. Differentiating on both sides of Eq. (8), we haveEq. (12), we have
Appendix B: Computation of the objective and its gradient
In this section, we derive an efficient expression for the objective L = Re f (ω̃ + iΓ̃) − iΓ̃ f′(ω̃ + iΓ̃) defined in Sec. 5.2 as well as its gradient. The gradient of the objective with respect to the design parameters is calculated with standard adjoint methods  in order to minimize the number of required matrix solves.Eq. (17) with respect to εk Eq. (12), Eq. (22) and Eq. (31), we have Eq. (30), Eq. (32) and Eq. (33), we have 54], the gradient of ℒ with respect to design parameters ε can be evaluated with only a single additional matrix solve Eq. (35).
X.L. thanks J. Joannopoulos and B. Zhen for helpful discussions. This work was supported in part by the AFOSR MURI for Complex and Robust On-chip Nanophotonics (Dr. Gernot Pomrenke) grant number FA9550-09-1-0704 and by the Institute for Soldier Nanotechnologies under Contract W911NF-07-D0004. This work used the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by National Science Foundation grant number OCI-1053575.
References and links
2. B. Bourdin, “Filters in topology optimization,” Int. J. Numer. Methods Eng. 50, 2143–2158 (2001). [CrossRef]
3. M. P. Bendsoe and O. Sigmund, Topology Optimization: Theory, Methods and Applications, 2nd ed. (Springer, 2003).
4. M. Y. Wang, X. Wang, and D. Guo, “A level set method for structural topology optimization,” Comput. Methods Appl. Mech. Eng. 192, 227–246 (2003). [CrossRef]
5. G. Allaire, F. Jouve, and A.-M. Toader, “Structural optimization using sensitivity analysis and a level-set method,” J. Comput. Phys. 194, 363–393 (2004). [CrossRef]
6. O. Sigmund, “Manufacturing tolerant topology optimization,” Acta Mech. Sinica 25, 227–239 (2009). [CrossRef]
7. E. Andreassen, A. Clausen, M. Schevenels, B. Lazarov, and O. Sigmund, “Efficient topology optimization in MATLAB using 88 lines of code,” Struct. Multidiscip. Optim. 43, 1–16 (2011). [CrossRef]
8. J. S. Jensen and O. Sigmund, “Topology optimization for nano-photonics,” Laser Photonics Rev. 5, 308–321 (2011). [CrossRef]
9. J. Bravo-Abad, A. Rodriguez, P. Bermel, S. G. Johnson, J. D. Joannopoulos, and M. Soljacic, “Enhanced nonlinear optics in photonic-crystal microcavities,” Opt. Express 15, 16161–16176 (2007). [CrossRef] [PubMed]
10. R. W. Boyd, Nonlinear Optics, 3rd ed. (Academic, 2008).
11. H. Hashemi, C. W. Qiu, A. P. McCauley, J. D. Joannopoulos, and S. G. Johnson, “A diameter–bandwidth product limitation of isolated-object cloaking,” Phys. Rev. A 86,013804 (2012). [CrossRef]
12. D. P. Bertsekas, Nonlinear Programming, 2nd ed. (Athena Scientific, 1999).
13. B. E. A. Saleh and M. C. Teich, Fundamentals of Photonics (Wiley Series in Pure and Applied Optics), 2nd ed. (Wiley-Interscience, 2007).
14. J. D. Joannopoulos, S. G. Johnson, J. N. Winn, and R. D. Meade, Photonic Crystals: Molding the Flow of Light, 2nd ed. (Princeton University, 2008).
15. D. G. Rabus, Integrated Ring Resonators, Vol. 127 of Springer Series in Optical Sciences(Springer, 2007).
17. J. Vuckovic, M. Loncar, H. Mabuchi, and A. Scherer, “Optimization of the Q factor in photonic crystal microcavities,” IEEE J. Quantum Electron. 38, 850–856 (2002). [CrossRef]
18. C. W. Wong, P. T. Rakich, S. G. Johnson, M. Qi, H. I. Smith, E. P. Ippen, L. C. Kimerling, Y. Jeon, G. Barbastathis, and S. G. Kim, “Strain-tunable silicon photonic band gap microcavities in optical waveguides,” Appl. Phys. Lett. 84, 1242–1244 (2004). [CrossRef]
21. B.-S. Song, S. Noda, T. Asano, and Y. Akahane, “Ultra-high-Q photonic double-heterostructure nanocavity,” Nat. Mater. 4, 207–210 (2005). [CrossRef]
22. D. C. Dobson and F. Santosa, “Optimal Localization of Eigenfunctions in an Inhomogeneous Medium,” SIAM J. Appl. Math. 64, 762–774 (2004). [CrossRef]
23. C.-Y. Kao and F. Santosa, “Maximization of the quality factor of an optical resonator,” Wave Motion 45, 412–427 (2008). [CrossRef]
24. W. R. Frei, H. T. Johnson, and K. D. Choquette, “Optimization of a single defect photonic crystal laser cavity,” J. Appl. Phys. 103,033102 (2008). [CrossRef]
26. A. W. Snyder and J. Love, Optical Waveguide Theory, Science Paperbacks, 190 (Springer, 1983).
27. R. Coccioli, M. Boroditsky, K. W. Kim, Y. Rahmat-Samii, and E. Yablonovitch, “Smallest possible electromagnetic mode volume in a dielectric cavity,” IEE Proceedings - Optoelectronics 145, 391–397 (1998). [CrossRef]
28. E. M. Purcell, “Spontaneous emission probabilities at radio frequencies,” Phys. Rev. 69,674 (1946).
29. L. Novotny and B. Hecht, Principles of Nano-Optics (Cambridge University, 2006). [CrossRef]
31. B. Zhen, S.-L. Chua, J. Lee, A. W. Rodriguez, X. Liang, S. G. Johnson, J. D. Joannopoulos, M. Soljačić, and O. Shapira, “Enabling enhanced emission and low-threshold lasing of organic molecules using special Fano resonances of macroscopic photonic crystals,” Proc. Natl. Acad. Sci. U. S. A. 110, 13711–13716 (2013). [CrossRef] [PubMed]
32. E. A. J. Marcatili, “Bends in optical dielectric guides,” Bell Syst. Tech. J. 48, 2103–2132 (1969). [CrossRef]
33. F. Ladouceur, “Roughness, inhomogeneity, and integrated optics,” J. Lightwave Technol. 15, 1020–1025 (1997). [CrossRef]
34. V. S. Ilchenko, P. S. Volikov, V. L. Velichansky, F. Treussart, V. Lefèvre-Seguin, J. M. Raimond, and S. Haroche, “Strain-tunable high-Q optical microsphere resonator,” Opt. Commun. 145, 86–90 (1998). [CrossRef]
36. L. N. Trefethen and D. Bau, Numerical Linear Algebra (SIAM, 1997). [CrossRef]
37. K. Inoue and K. Ohtaka, Photonic Crystals: Physics, Fabrication and Applications, Springer Series in Optical Sciences (Springer, 2010).
38. A. Oskooi and S. G. Johnson, “Electromagnetic wave source conditions,” in Advances in FDTD Computational Electrodynamics: Photonics and Nanotechnology, A. Taflove, A. Oskooi, and S. G. Johnson, eds. (Artech, 2013), Chap. 4, pp. 65–100.
39. J.-M. Gerard and B. Gayral, “Strong Purcell effect for InAs quantum boxes in three-dimensional solid-state microcavities,” J. Lightwave Technol. 17, 2089–2095 (1999). [CrossRef]
40. O. J. F. Martin and N. B. Piller, “Electromagnetic scattering in polarizable backgrounds,” Phys. Rev. E 58, 3909–3915 (1998). [CrossRef]
41. G. D’Aguanno, N. Mattiucci, M. Centini, M. Scalora, and M. J. Bloemer, “Electromagnetic density of modes for a finite-size three-dimensional structure,” Phys. Rev. E 69,057601 (2004). [CrossRef]
42. J. D. Jackson, Classical Electrodynamics, 2nd ed. (John Wiley, 1975).
43. W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, Numerical Recipes: The Art of Scientific Computing, 3rd ed. (Cambridge University, 2007).
44. M. E. Peskin and D. V. Schroeder, An Introduction To Quantum Field Theory (Frontiers in Physics) (Westview, 1995).
45. S. G. Johnson, “Numerical methods for computing Casimir interactions,” in Casimir Physics, Vol. 834 of Lecture Notes in Physics, D. Dalvit, P. Milonni, D. Roberts, and F. da Rosa, eds. (Springer, 2011), Chap. 6. [CrossRef]
46. L. D. Landau, L. P. Pitaevskii, and E. M. Lifshitz, Electrodynamics of Continuous Media, 2nd ed. (Butterworth-Heinemann, 1984).
47. W. C. Chew, Waves and Fields in Inhomogeneous Media (IEEE, 1995).
48. N. A. P. Nicorovici, R. C. McPhedran, and L. C. Botten, “Relative local density of states for homogeneous lossy materials,” Physica B 405, 2915–2919 (2010). [CrossRef]
49. L. V. Ahlfors, Complex Analysis, 3rd ed. (McGraw-Hill, 1978).
50. A. W. Rodriguez, A. P. McCauley, J. D. Joannopoulos, and S. G. Johnson, “Theoretical ingredients of a Casimir analog computer,” Proc. Natl. Acad. Sci. U. S. A. 107, 9531–9536 (2010). [CrossRef] [PubMed]
51. X. Liang, Ph.D. thesis, Massachusetts Institute of Technology, 2013.
52. S. Scheel, L. Knöll, and D. G. Welsch, “Spontaneous decay of an excited atom in an absorbing dielectric,” Phys. Rev. A 60, 4094–4104 (1999). [CrossRef]
53. C. Van Vlack and S. Hughes, “Finite-difference time-domain technique as an efficient tool for calculating the regularized Green function: applications to the local-field problem in quantum optics for inhomogeneous lossy materials,” Opt. Lett. 37, 2880–2882 (2012). [CrossRef] [PubMed]
54. G. Strang, Computational Science and Engineering (Wellesley-Cambridge, 2007).
55. A. Christ and H. L. Hartnagel, “Three-dimensional finite-difference method for the analysis of microwave-device embedding,” IEEE Trans. Microwave Theory 35, 688–696 (1987). [CrossRef]
56. K. Yasumoto, Electromagnetic Theory and Applications for Photonic Crystals, Optical Science and Engineering (CRC, 2005). [CrossRef]
57. 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]
58. T. A. Davis, Direct Methods for Sparse Linear Systems (Fundamentals of Algorithms) (SIAM, 2006). [CrossRef]
59. S. Balay, W. D. Gropp, L. C. McInnes, and B. F. Smith, “Efficient Management of Parallelism in Object Oriented Numerical Software Libraries,” In Modern Software Tools in Scientific Computing, E. Arge, A. M. Bruaset, and H. P. Langtangen, eds. (Birkhäuser, 1997). [CrossRef]
60. S. Balay, J. Brown, K. Buschelman, V. Eijkhout, W. D. Gropp, D. Kaushik, M. G. Knepley, L. C. McInnes, B. F. Smith, and H. Zhang, “PETSc Users Manual,” Technical Report No. ANL-95/11 - Revision 3.3, Argonne National Laboratory (2012).
61. S. Balay, J. Brown, K. Buschelman, W. D. Gropp, D. Kaushik, M. G. Knepley, L. C. McInnes, B. F. Smith, and H. Zhang, “PETSc Web page,” http://www.mcs.anl.gov/petsc, 2012.
62. P. Hénon, P. Ramet, and J. Roman, “PaStiX: a high-performance parallel direct solver for sparse symmetric positive definite systems,” Parallel Computing 28, 301–321 (2002). [CrossRef]
63. S. G. Johnson, “The NLopt nonlinear-optimization package,” http://ab-initio.mit.edu/nlopt.
64. D. Liu and J. Nocedal, “On the limited memory BFGS method for large scale optimization,” Math. Program. 45, 503–528 (1989). [CrossRef]
65. K. Svanberg, “A class of globally convergent optimization methods based on conservative convex separable approximations,” SIAM J. Optimiz. 12, 555–573 (2002). [CrossRef]
66. J. Nocedal and S. J. Wright, Numerical Optimization (Springer, 2000).
67. A. Mutapcic, S. Boyd, A. Farjadpour, S. G. Johnson, and Y. Avniel, “Robust design of slow-light tapers in periodic waveguides,” Eng. Optimiz. 41, 365–384 (2009). [CrossRef]
68. A. Oskooi, A. Mutapcic, S. Noda, J. D. Joannopoulos, S. P. Boyd, and S. G. Johnson, “Robust optimization of adiabatic tapers for coupling to slow-light photonic-crystal waveguides,” Opt. Express 20, 21558–21575 (2012). [CrossRef] [PubMed]
70. Z. Artstein, “Discrete and continuous bang-bang and facial spaces or: Look for the extreme points,” SIAM Rev. 22, 172–185 (1980). [CrossRef]
71. A. F. Oskooi, Ph.D. thesis, Massachusetts Institute of Technology, 2010.
73. B. Osting and M. I. Weinstein, “Long-lived scattering resonances and Bragg structures,” SIAM J. Appl. Math. 73, 827–852 (2013). [CrossRef]
74. T. Inui, Y. Tanabe, and Y. Onodera, Group Theory and Its Applications in Physics (Springer, 1996).
75. S. Boyd and L. Vandenberghe, Convex Optimization (Cambridge University, 2004). [CrossRef]
76. P. D. Drummond, K. J. McNeil, and D. F. Walls, “Non-equilibrium transitions in sub/second harmonic generation,” Opt. Acta 27, 321–335 (1980). [CrossRef]
77. L.-A. Wu, M. Xiao, and H. J. Kimble, “Squeezed states of light from an optical parametric oscillator,” J. Opt. Soc. Am. B 4, 1465–1475 (1987). [CrossRef]
79. A. Rodriguez, M. Soljacic, J. D. Joannopoulos, and S. G. Johnson, “χ(2)and χ(3)harmonic generation at a critical power in inhomogeneous doubly resonant cavities,” Opt. Express 15, 7303–7318 (2007). [CrossRef] [PubMed]
80. I. B. Burgess, Y. Zhang, M. W. McCutcheon, A. W. Rodriguez, J. Bravo-Abad, S. G. Johnson, and M. Loncar, “Design of an efficient terahertz source using triply resonant nonlinear photonic crystal cavities,” Opt. Express 17, 20099–20108 (2009). [CrossRef] [PubMed]
81. Z.-F. Bi, A. W. Rodriguez, H. Hashemi, D. Duchesne, M. Loncar, K.-M. Wang, and S. G. Johnson, “High-efficiency second-harmonic generation in doubly-resonant χ(2)microring resonators,” Opt. Express 20, 7526–7543 (2012). [CrossRef] [PubMed]
82. M. W. McCutcheon and M. Loncar, “Design of a silicon nitride photonic crystal nanocavity with a Quality factor of one million for coupling to a diamond nanocrystal,” Opt. Express 16, 19136–19145 (2008). [CrossRef]
83. A. R. Md Zain, N. P. Johnson, M. Sorel, and R. M. De La Rue, “Ultra high quality factor one dimensional photonic crystal/photonic wire micro-cavities in silicon-on-insulator (SOI),” Opt. Express 16, 12084–12089 (2008). [CrossRef]
84. Z. M. Meng, F. Qin, Y. Liu, and Z. Y. Li, “High-Q microcavities in low-index one-dimensional photonic crystal slabs based on modal gap confinement,” J. Appl. Phys. 109,043107 (2011). [CrossRef]
85. J. T. Robinson, C. Manolatou, L. Chen, and M. Lipson, “Ultrasmall mode volumes in dielectric optical micro-cavities,” Phys. Rev. Lett. 95,143901 (2005). [CrossRef]
87. S. Kita, K. Nozaki, S. Hachuda, H. Watanabe, Y. Saito, S. Otsuka, T. Nakada, Y. Arita, and T. Baba, “Photonic Crystal Point-Shift Nanolasers With and Without Nanoslots Design, Fabrication, Lasing, and Sensing Characteristics,” IEEE J. Sel. Top. Quantum Electron. 17, 1632–1647 (2011). [CrossRef]
88. Y. Li, J. Zheng, J. Gao, J. Shu, M. S. Aras, and C. W. Wong, “Design of dispersive optomechanical coupling and cooling in ultrahigh-Q/V slot-type photonic crystal cavities,” Opt. Express 18, 23844–23856 (2010). [CrossRef] [PubMed]
89. O. Sigmund and J. Petersson, “Numerical instabilities in topology optimization: A survey on procedures dealing with checkerboards, mesh-dependencies and local minima,” Struct. Multidiscip. Optim. 16, 68–75 (1998).
90. S. J. Osher and R. P. Fedkiw, Level Set Methods and Dynamic Implicit Surfaces (Springer, 2002).
91. C. Y. Kao, S. Osher, and E. Yablonovitch, “Maximizing band gaps in two-dimensional photonic crystals by using level set methods, Appl. Phys. B 81, 235–244 (2005). [CrossRef]
92. J. K. Guest, J. H. Prévost, and T. Belytschko, “Achieving minimum length scale in topology optimization using nodal design variables and projection functions,” Int. J. Numer. Methods Eng. 61, 238–254 (2004). [CrossRef]
93. O. Sigmund, “Morphology-based black and white filters for topology optimization,” Struct. Multidiscip. Optim. 33, 401–424 (2007). [CrossRef]
94. S. Xu, Y. Cai, and G. Cheng, “Volume preserving nonlinear density filter based on heaviside functions,” Struct. Multidiscip. Optim. 41, 495–505 (2010). [CrossRef]
95. F. Wang, B. Lazarov, and O. Sigmund, “On projection methods, convergence and robust formulations in topology optimization,” Struct. Multidiscip. Optim. 43, 767–784 (2011). [CrossRef]
96. A. Ben-Tal, L. El Ghaoui, and A. S. Nemirovski, Robust Optimization (Princeton University, 2009).
97. D. Bertsimas, O. Nohadani, and K. M. Teo, “Robust optimization for unconstrained simulation-based problems,” Oper. Res. 58, 161–178 (2010). [CrossRef]
98. F. Wang, J. S. Jensen, and O. Sigmund, “Robust topology optimization of photonic crystal waveguides with tailored dispersion properties,” J. Opt. Soc. Am. B 28, 387–397 (2011). [CrossRef]
99. H. Men, R. M. Freund, N. C. Nguyen, J. Saa-Seoane, and J. Peraire, “Fabrication-adaptive optimization with an application to photonic crystal design,” arXiv:1307.5571.
100. A. Faraon, C. Santori, Z. Huang, V. M. Acosta, and R. G. Beausoleil, “Coupling of nitrogen-vacancy centers to photonic crystal cavities in monocrystalline diamond,” Phys. Rev. Lett. 109,033604 (2012). [CrossRef] [PubMed]
101. L. Li, M. Trusheim, O. Gaathon, K. Kisslinger, C.-J. Cheng, M. Lu, D. Su, X. Yao, H.-C. Huang, I. Bayn, A. Wolcott, R. M. Osgood Jr., and D. Englund, “Reactive ion etching: Optimized diamond membrane fabrication for transmission electron microscopy,” J. Vac. Sci. Technol. B 31,06FF01 (2013). [CrossRef]