Formulation for scalable optimization of microcavities via the frequency-averaged local density of states

: We present a technique for large-scale optimization of optical microcavities based on the frequency-averaged local density of states (LDOS), which circumvents computational difﬁculties 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 difﬁcult or unavailable, and in particular we design minimal-volume multi-mode cavities (e.g. for nonlinear frequency-conversion applications).


Introduction
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 ≈ 10 4 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 [1]. 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][3][4][5][6][7][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 twofrequency 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 [11] 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 [12], 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][20][21]). These manual designs have been tremendously successful, attaining experimentally verified Q of 10 6 or more in modal volumes only a few wavelengths in diameter [19][20][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 [10]. Largescale "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][23][24][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 [22], or maximized a 2d Q without controlling V [23] (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 [25], which permitted 2d calculations and led to intriguing designs, but limits the attainable Q (to < 10 4 ) 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 [24], 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 [14], 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 [14] that U ∼ Q 2 . 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. [24] 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 [24]: 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.

Eigenproblem formulation
There are two key figures of merit for a resonant mode E n (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 [14]. Mathematically, Q is related to the frequency-domain Maxwell eigenvalue problem: with radiation boundary conditions. Because of the lossy boundary conditions, the eigenproblem is non-Hermitian and the eigenvalues are complex. The Q for the mode E n (x) with eigen-frequency ω n [14] is Technically, if one considers an infinite open system (as opposed to a finite system with some absorbing boundary layer), a number of mathematical subtleties arise because the leaky "modes" are not true normalizable eigenfunctions, but are rather related to the residues of poles in the Green's function [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 [27], defined as is a measure of the volume within which the mode is confined. (Although 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)|E n (x)| 2 dx 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 [28], which describes the enhancement of spontaneous emission rate, can be written as [29,30] 3Q Here λ is the vacuum wavelength and n is the index of refraction. For applications with light-matter interactions (such as lasers, sensors, and nonlinear frequency converters), maximal lifetime Q and minimal modal volume V are desirable [31]. 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 [32]. In practice, any optimization in a finite computation cell will obtain a finite Q and V [23], 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 (Q rad ) in a cavity, there are also radiation losses due to surface roughness (Q roughness ) and material absorption (Q absorption ). The total loss rate 1/Q total is the sum of these three effects [14]: In real applications, therefore, the Q total cannot be arbitrarily large. For example, in integrated optics it is difficult to get Q total more than a few million due to surface roughness [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 cavitybased 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 Q coupling . It turns out that the losses in such a coupled device are proportional to Q coupling /Q total [14]. 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: whereQ is determined by the bandwidth and loss tolerance of applications. By solving the non-Hermitian Maxwell eigenproblem Eq. (1), one could obtain Q and V from eigenvectors E n (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 [14], 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 [22], 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 [36].

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) [37]. 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 LDOS j (ω, 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 ∼ê j e −iωt δ (x − x ′ ). (The total LDOS is simply ∑ j LDOS j , 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 LDOS j in the cavity [29,30,38,39].
These quantities can be defined more precisely as follows [29,37,38,40,41]. Poynting's theorem [42] implies that the power radiated by a dipole J(x) is where E(x) is the total electric field solving the frequency-domain scattering problem This yields the following well known [29,37,38,40,41] definition of LDOS j : where the 12/π factor is a conventional normalization that arises from the dual interpretation of the LDOS as a local "density" of eigenstates [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 [38] LDOS j (ω, and by taking the ratio of this quantity with the LDOS of a homogeneous medium with index n = √ ε one [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 [38], 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 LDOS j 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 LDOS j (ω, 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 LDOS j (ω, x ′ ) (for the average case) or min j LDOS j (ω, x ′ ) (for the worst case). In nonlinear devices for frequency conversion (in particular, second harmonic generation), a more relevant figure of merit might be min n LDOS j (ω 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 LDOS j (ω, 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, LDOS j 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 explicit j and x ′ dependence from LDOS j (ω, 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.

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 lifetimeQ. 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Γ: Here, W (ω) is some weight function or window function we choose, which is peaked around the design frequencyω and decays rapidly (with a finite integral) outside of a bandwidthΓ aroundω. As we will explain in Sec. 4.3, it will turn out that performing this average is mathematically equivalent to computing the LDOS at a single frequency with an absorption loss added into the system, which effectively limits Q to the desired bandwidth (more precisely, for Q ≫Q the figure of merit is dominated by the modal volume). Equivalently, as Q increases the LDOS of a resonance approaches a delta function with amplitude ∼ 1/V [38], so that the integral Eq. (11) is determined mainly by the 1/V factor as soon as the resonance is narrower thañ 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 [43], 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 [24].) 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 [11], and related ideas can be found in quantum field theory [44,45].)

Causality and analyticity
Before we proceed, we first define a function f (ω), which is a complex version of LDOS(ω): Comparing with Eq. (9), it is clear that [From now on, we will omit the x ′ dependency of f (ω, x ′ ).] Note that the operator M (ε, ω) 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 [46]. Therefore, f (ω) is also analytic in the upper-half complex-ω plane. Fig. 2. Contour integration path. The frequency-averaged LDOS is the path integral along arc A 1 in the limit of an infinite-radius arc. By choosing the proper window/weight function W (ω) for optimizing LDOS in a desired bandwidth, the contribution along arcs A 2 and A 3 can be made negligible compared to A 1 . Therefore, the residues at polesω k enclosed by this contour can be used to approximate the averaged LDOS.

Contour integration
In this section, we are going to compute the mean LDOS by exploiting the analyticity of f (ω), via: Here, p.v. denotes the Cauchy principal value, which we use because the imaginary part of f (ω) may have a singularity at ω = 0, as in the case of f (ω) in free space [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] A 1 Here,ω k denotes a pole of W (ω) in the upper-half plane and Res is its residue [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 A 3 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 A 3 . We can obtain the contribution from arc A 2 by evaluating the residue due to the simple pole of f (ω) at ω = 0 The factor −1/2 comes from the fact that the integration is along a clockwise semicircle. Since the weight function is peaked around design frequencyω with some narrow bandwidthΓ, we can require that W (0) is small, and therefore the contribution from A 2 is negligible comparing to the residues atω k . From Eq. (13), L is just the path integral along the path A 1 [46]. Therefore, we have

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 as We denote this complex scattering operator byM (ε, ω), namelỹ Clearly, this is equivalent to solving a scattering problem at real frequencyω with materials (In fact, any change to the frequency can be converted into a change of materials [50].) In particular, adding a positive imaginary part to ω (forω k in the upper-half plane [11]) corresponds to a positive imaginary part inε(x) and µ(x), which corresponds (with our e −iωt convention) to an absorption loss.

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 [51].)

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 is which only requires solving the scattering problem Eq. (17) at a single complex-frequencỹ ω + iΓ. L 1 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 [14] (forQ ≫ 1) to yield an absorption lifetime of exactly Q absorption =Q. From Eq. (5), this means that Q total ≤Q, and that increasing Q rad ≫Q will have little effect on the Purcell factor Q tot /V . Hence, optimizing the frequency-averaged objective will tend to increase Q rad 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 [41]. It is known that this absorbed power is infinite because E(x) diverges as 1/r 3 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/r 3 , 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 [51] and the square of a Lorentzian.

Square of a Lorentzian
To ensure that the W (ω) decays faster than 1/|ω| 3 , we propose the window function, which is a normalized square of a Lorentzian function. This window function has a double pole at ω =ω + iΓ in the upper-half plane. Applying the residue theorem, we have from Eq. (16) and Eq. (20) where f ′ (·) denotes differentiation with respect to ω. In appendix 10, we show that From 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 In summary, we can still obtain the entire frequency-averaged LDOS by solving a single scattering problem 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 [41]. (A subtlety arises from having +E T E, rather than −|E| 2 , but the 1/r 3 divergence at r → 0 should be dominant in Im(E) for smallΓ so one should have 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/r 3 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.

A preliminary formulation
Now we have a preliminary formulation for our cavity optimization in terms of the frequencyaveraged local density of states: We can evaluate the objective L by contour integration, which only requires us to solve the complex scattering problem Eq. (17) once. If we choose the window function W (ω) from Eq. (20), then the problem can be reformulated as Mathematically, to compute the objective in Eq. (25), we need 1. For given ε(x),ω, andΓ, solve the complex scattering problem Eq. (17) to obtain E(x,ω + iΓ).
3. Take the real part of L to get L.
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 [54], only one more linear system with the same operatorM (ε,ω) 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: 1. Solve the complex scattering problem to obtain A(x,ω + iΓ).
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.

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.

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][56][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 [58]. In contrast to direct solvers, iterative solvers (e.g., GMRES or BiCGStab) work quite well if one has a good preconditioner [36]. 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 [51]. We implemented the FDFD solver with the parallel sparse-matrix library PETSc [59][60][61] and the parallel sparse-direct solver PaStiX [62].

Optimization scheme
We use a free-software implementation [63] of standard gradient-based optimization algorithms, and we find that low-storage BFGS [64] and conservative convex separable approximations (CCSA) [65] 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 (∼ Q 3 , as shown in [51]) in the direction perpendicular to the ridge, and it is well known that optimization along narrow ridges is problematic [66]. 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 specifiedQ (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 forQ = 10, then use that result as a starting point for optimization at Q = 100, and so on.

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 , (n = √ 12.4). 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.

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 [14], 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 × 10 8 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 × 10 10 and mode volume V = 0.075(λ /n) 2 . Clearly, the optimization itself discovers a "radially periodic" structure, reminiscent of a Bragg onion [69] . We also run the optimization with vacuum as initial guess and obtain similar structure (Fig. 5) with Q = 1.30 × 10 9 and V = 0.075(λ /n) 2 .
In these two optimizations, we gradually increase the specifiedQ (we decrease the bandwidth (a). A square region for degrees of freedom.
(b). A thin strip region for degrees of freedom. 1/Q) from 10 to 10 5 . The optimization atQ = 10 actually gives a high Q cavity (almost the same radiation Q) with the resonance at about 1.003ω. The optimizations at higherQ 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 [14], because Q → ∞ for a finite modal volume V as the number of layers in a Bragg onion increases [69]. 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 [70]) 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 [73] 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.

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 (a). PhC cavity initial guess for 2d TM optimization. It has quality factor Q=1.41×10 8 and mode volume V = 0.097(λ /n) 2 .
(b). 2d TM optimized structure with Q=1.01×10 10 and V = 0.075(λ /n) 2 obtained from PhC initial guess. 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×10 8 and mode volume V = 0.092(λ /n) 2 . [Again, theQ = 10 gives an equally high-Q cavity with resonant frequency at 1.0007ω. The optimizations at higherQ = 10 2 to 10 5 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 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 [74]). 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 botĥ e 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 min j 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 [74] 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 C mv symmetry group [74] 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 [75] of introducing a "dummy" variable t and solving maxt subject to the nonlinear constraints t ≤ LDOS j at each angle j. (The resulting nonlinear-programming problem was solved via the CCSA algorithm [65].) 8.3. 2D optimization for different frequencies: max min n LDOS(ω 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][77][78][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 singlefrequency optimization in the previous sections, which is not surprising considering that Bragg mirrors tend to have gaps at integer frequency multiples [14] 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).

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 [14] 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][83][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 Q rad as a function of the requested bandwidthQ, and the result is shown in Fig. 9(a). For any given d, asQ is increased then Q rad is increased as well (as explained in Sec. 4), until it saturates at a maximum determined by two factors. First, the maximum Q rad is limited by d: a longer d gives more degrees of freedom and increases the maximum possible (a). Optimized structure for randomly polarized dipole. polarizations. Left: microcavities which maximize the minimum LDOS at two frequencies ω 1 andω 2 = 2ω 1 , e.g. for intra-cavity second-harmonic generation applications [79,81].
Confinement of such integer-multiple frequencies is physically enabled by the fact the concentric Bragg-onion "1d" bandgaps tend to occur at integer-multiple frequencies [14]. A more challenging case is a two-frequency cavity forω 2 = 1.5ω 1 , resulting in the more complicated structures shown at right. (All structures were optimized from vacuum initial guesses.) Q rad , as expected. Second, around Q rad ≈ 10 7 the improvement becomes limited by the numerical precision, which prevents the optimization from making progress even though higher Q rad 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 ′ ).]  (a). Q rad vsQ for 2d thin strip with different length (different DOF).
(b). An optimized structure for 2d thin strip from vacuum initial guess. Fig. 9. 2D TE optimization for thin strips with fixed width (geometry sketched in Fig. 3b). Fig. (a): Q rad vs.Q for 2d thin stirps with same width (λ ) but different length (d = λ , 2λ , 3λ , 5λ ). AsQ is increased in the optimization, higher Q rad are obtained until Q rad is limited by the degrees of freedom. As the degrees of freedom increase, Q rad first gets bigger, but becomes saturated at some level around 10 7 due to numerical precision. Fig. (b): An optimized 2d thin-strip structure with width λ and length d = 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.

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" [14] 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 A comparison with other large-or small-scale optimization work, such as 2.5d optimization [72], L 3 -type cavity [20] and H 0 -type cavity [1] 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 [1]. (Note that L 3 -type [20] and H 0 -type [1] cavities are designed from small-scale optimization and can be fabricated; while both 2.5d optimization [72] and our results are purely theoretical and computational investigation from a large-scale optimization perspective.) Table 1. Comparison of Q and V for structures from our result and the literature.

Optimization
Quality Factor Q Mode Volume V (λ /n) 3 2.5d optimization [72] 8000 0.32 L 3 -type cavity optimization [20]  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 [3] 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(λ /n Si ) 3 .
It is known that air-slot cavities [85-88] can have extremely small volumes. For example, Nomura [86] reported an air-slot cavity with Q = 4.8 × 10 6 and V = 0.015(λ /n air ) 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 (λ /n air ) 3 versus (λ /n Si ) 3 . In our units, their modal volume is 0.65(λ /n Si ) 3 . Second, these two types of cavities are for different applications: airslot 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 [85], after the introduction of an air-slot, the mode volume decreases by a factor of (n Si /n air ) 2 without changing Q. In our case, this factor is about 12.4, and the new mode volume is 4.8 × 10 −3 (λ /n Si ) 3 = 1.1 × 10 −4 (λ /n air ) 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 airslot 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.

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) [3], 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][6][7][8][92][93][94][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) [3]. 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].

Conclusion
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 cavityoptimization 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.