1. Introduction

Recently, we obtained theoretical upper bounds [1] to the (electric) local density of states (LDOS) $\rho (\mathbf {x},\omega )$, a key figure of merit for light–matter interactions (e.g. spontaneous emission) proportional to the power emitted by a dipole current at a position $\mathbf {x}$ and frequency $\omega$ [2–10]. For a resonant cavity with quality factor $Q$ (a dimensionless lifetime), LDOS is proportional to the “Purcell factor” $Q/V$ where $V$ is a modal volume [6,11], so LDOS is a measure of light localization in space and time. Our LDOS bounds $\sim |\chi |^{2}/\operatorname {Im}\chi / d^{3}$ (reviewed in Sec. 2) depend on the material used (described by the $\omega$-dependent susceptibility $\chi = \varepsilon - 1$) and the minimum separation $d$ between the emitter and the material, but are otherwise independent of shape, and hence give an upper limit to the localization attainable by any possible resonant cavity for $(\chi ,d,\omega )$. However, it is an open question to what extent these bounds are tight, i.e. is there any particular cavity design that comes close to the bounds? Initial investigations of a few simple resonant structures were often orders of magnitude from the upper bounds (except at the surface-plasmon wavelength for a given metal) [1,12,13]. In this paper, we perform computational optimization of 3D metallic cavities at many wavelengths and find that the bounds are much more nearly attainable than was previously known.

In particular, we perform many-parameter shape optimization of LDOS for cavities formed by voids in silver (since it has the largest $|\chi |^{2}/\operatorname {Im}\chi$, and therefore the largest bounds) at wavelengths $\lambda$ from 400–900 nm and a $d=50$ nm emitter–metal separation, depicted in Fig. 1. As described in Sec. 4, we obtain single-polarization LDOS values within a factor of $\approx 4$ of the theoretical upper bounds, and total (all-polarization) LDOS within a factor of $\approx 10$ of the bounds. Of course, real cavities would have a finite thickness of metal, but our goal is to attain the maximum possible LDOS—we find that a finite-thickness coating has slightly worse performance, but $> 95$% of the LDOS of the infinite metal is attained by $\approx 100$ nm thickness at $\lambda =500$ nm, and more generally we can theoretically bound [1] the improvement attainable with any additional structure of air voids outside of our cavity. Although our focus is on fundamental upper limits rather than manufacturable cavities, we find that simple shapes (optimized cylinders) are within $\approx 20$% of the LDOS of optimized many-parameter irregular shapes, analogous to results we obtained previously for optimized scattering and absorption [1]. (Our optimized cavities are deeply subwavelength along their shortest axes, very different from the non-plasmonic resonant regime where the diameter is much larger than the skin depth so that the walls simply act as mirrors.) Moreover, we find that optimizing for a single emitter polarization (the “polarized” LDOS) does nearly as well (within $\approx 10$%) as optimizing the total LDOS (power summed over all emitter polarizations), reminiscent of earlier results in 2D dielectric cavities where LDOS optimization arbitrarily picked one polarization to enhance [14]. We perform the shape optimization using an efficient boundary-element method [15,16] (which has unknowns only on the metal surface) coupled with adjoint sensitivity analysis [17,18] and an optimization algorithm robust to discretization errors [19] as described in Sec. 3. Although it is possible that even tighter LDOS bounds could be obtained in future results by incorporating additional physical constraints [20–23], we believe that our results show that the existing bounds are already closely related to attainable performance and provide useful guidance for optical cavity design.

 

Fig. 1. Schematic cavity-optimization problem: the shape of an air cavity in a metallic (silver) background is optimized to maximize the LDOS for emitters (dipoles) at the center $o$, constrained for a minimum separation $d$ (the metal lies outside of a sphere of radius $d$).

Download Full Size | PDF

In this work, we optimize the total electric LDOS, which is the sum of the absorbed and radiated powers from electric dipoles (Sec. 2) for comparison with the upper bounds, and in fact the power is entirely absorbed for an air-void cavity as in Fig. 1. However, the same theoretical procedure yielded a bound on the purely radiative power that was simply $1/4$ of the total-LDOS bound [1], and in fact the two results are closely related. In general, the addition of low-loss input/output channels can be analyzed via coupled-mode theory as small perturbations to existing cavity designs [24]. As reviewed in Appendix E, given an resonant absorptive cavity, one could modify it to radiate at most $\approx 1/4$ of the original LDOS by slightly perturbing it to add a radiative-escape channel (e.g. a small hole or thinning in the cavity wall) tuned to match the absorption-loss rate [25,26]. (The same channel could also be used to introduce input energy, e.g. for pumping an emitter.) Appendix D shows an example of this: by thinning the cavity walls, we achieve radiated power $\approx 1/4$ of the absorbed power in the purely absorbing cavity.

Many previous authors have computationally optimized the LDOS of cavities (or equivalent quantities such as the Purcell factor $Q/V$), including many-parameter shape or “topology” optimization [27–32], but in most cases these works did not compare to the recent upper bounds. In many cases, these works studied lossless dielectric materials where the bound diverges (though a finite LDOS is obtained for a finite volume [14,33] and/or a finite bandwidth [13,14]). Designs specifically for LDOS of metallic resonators that compared to the bounds initially yielded results far below the bounds except for the special case of a planar surface at the surface-plasmon frequency of the material [1,34], but recent topology optimization in two dimensions came within a factor of 10 of the 2D bound [33]. Semi-analytical calculations have also been published for resonant modes in spherical metallic voids [35], but did not calculate LDOS. Therefore, the opportunity remains for optimized metallic LDOS designs in three dimensions that approach the theoretical upper bounds. To come as close as possible to the bounds, we focus initially on the idealized case of an air void surrounded by metal filling the rest of space, so that there are no radiation losses; later in this paper, we consider the small corrections that arise due to finite metal thickness.

2. Local density of states (LDOS)

The (electric) LDOS is equivalent to the total power expended by three orthogonal dipole currents [2]:

From energy-conservation considerations, previous work found an upper bound for LDOS enhancement inside a cavity compared to vacuum electric LDOS ($\rho _0=\omega ^{2}/2\pi ^{2}c^{3}$ [7]), given any a material susceptibility $\chi$ and an emitter–material separation $d$ at a frequency $\omega = ck$, to be [1]:

It is important to emphasize that the derivation of Eq. (2) gives a rigorous upper bound to the LDOS, but does not say what structure (if any) achieves the bound. By actually solving Maxwell’s equations for various geometries, we can investigate how closely the bound can be approached (how “tight” the bound is). It is possible that incorporating additional constraints may lead to tighter bounds in the future [20–23], but our results below already show that Eq. (2) is achievable within an order of magnitude.

3. Cavity-optimization methods

To numerically compute the LDOS inside a metal cavity, we employed a free-software implementation [16] of the boundary element method (BEM) [15]. A BEM formulation only involves unknown tangential fields on the metal surface, leading to modest-size computations for 3D metallic voids (Fig. 1). The complex dielectric constant of silver was interpolated from tabulated data [37]. In addition, we implemented an adjoint method [17,18] to rapidly obtain the gradient of the LDOS with respect to the shape parameters described below. As reviewed in Appendix B, the gradient of LDOS with respect to all shape parameters simultaneously is obtained by the adjoint method using only two BEM simulations—the original problem and an adjoint problem (the same Maxwell problem with artificial “adjoint” sources). Since the adjoint problem is the same Maxwell/BEM operator, we need only form and factorize the BEM matrix a single time, and the computational cost to solve both the forward and adjoint problems is essentially equivalent to a single simulation. Validation against a semi-analytical solution for spheres [35] is discussed in Appendix C.

In order to parameterize an arbitrary cavity shape numerically, we use a level-set description [38,39], combined with a free-software surface-mesh generator CGAL [40,41]. In particular, we describe the radius of the cavity around the source point by a function $R(\theta ,\phi )$ (in spherical coordinates), which is expanded below in either spherical harmonics or other polynomials, and equivalently pass a level-set function $\Phi = r - R(\theta ,\phi )$ to CGAL (such that $\Phi = 0$ defines the surface). We considered various parameterizations of the shape function $R$. The simplest geometries considered were ellipsoids, cylinders, or rectangular boxes, described by two or three parameters. For many-parameter optimization with a minimum radius (separation) $d$, the function $R$ is expressed as an expansion in some basis functions $S_n(\theta ,\phi )$ as:

The level set is discretized for BEM (by the CGAL software) into a triangular surface mesh. We used $5\times$ greater resolution for surface points closer to the dipole source (radius $\lesssim 70$ nm), since the singularity of the fields at the source point leads to rapid variations nearby, for around $5000$ triangles overall. As we deformed the shape during optimization, we first deform the triangles smoothly as long as all angles remained between $30^{\circ }$ and $120^{\circ }$, after which point we triggered a re-meshing step. Unfortunately, re-meshing causes slight discontinuities in the objective function and its derivatives which tend to confuse optimization algorithms expecting completely smooth functions [42]. We tried various optimization algorithms designed to be robust to such “numerical noise” [19,43], and found that the Adam stochastic-optimization algorithm [19] seems to work best for our problem.

4. Cavity-optimization results

4.1 Total LDOS

We performed numerical shape optimization of the LDOS for cavities formed by voids in silver [37] at wavelengths $\lambda$ from 400–900 nm, for both simple geometries (cylinder, ellipsoid, and rectangular box) and complex many-parameter shape (spherical harmonics). To obtain a finite optimum LDOS, one must choose a lower bound on the emitter–metal separation distance $d$ [1]. We chose $d=50$ nm so that $kd < 1$ for all optimized wavelengths; this allows us to use Eq. (2) as the LDOS upper bound, neglecting additional far-field effects [1] (see also Appendix A), while a much smaller $d$ was inconvenient to model (due to extremely small feature sizes and even nonlocal effects [44] at such scales). Note that $d$ also sets a lengthscale for the region of the cavity with maximum LDOS: as long as one shifts the emitter location by $\ll d$ ($\lesssim 10$ nm), the optimized LDOS will be of the similar magnitude, but other (unoptimized) locations in the cavity will typically have a drastically different LDOS.

The results of the optimized LDOS as a function of the wavelength are displayed in Fig. 2(a). Note that each wavelength corresponds to a different structure optimized for that particular wavelength. If we fix the structure as the one optimized for $\lambda = 500$ nm, the resulting LDOS spectrum is shown in Fig. 2(b), exhibiting a peak at the optimized wavelength. For simple shapes (cylinder and ellipsoid), we swept the parameters through a large range and found the global optimum among several local optima. We also optimized rectangular boxes, but their performance was nearly identical to that of the cylinders (but slightly worse), so they are not shown. For the 16-parameter (spherical-harmonic) level-set optimization, we performed local optimization for $\sim 10$ random starting points and plotted the best result along with a few other typical local optima.

 

Fig. 2. (a) Total LDOS optima as a function of the wavelength $\lambda$ for a minimum separation $d=50$ nm, along with the upper bound (black line). A separately optimized structure is used for each wavelength, either optimized cylinders (orange line) and ellipsoids (green line) or general shape optimization via the optimized spherical-harmonic (SH) surfaces (blue dots) of Eq. (3). Several SH local optima are shown for each $\lambda$, whereas for cylinders and ellipsoids only the global optima are shown. (b) LDOS spectra of the spherical-harmonic (blue) and cylinder (orange) structures optimized for $\lambda =500$ nm, the the shapes (not to scale) inset (see supplementary Visualization 1 and Visualization 2 for 3D views). Also shown is the total-LDOS spectrum of a polynomial shape (dashed blue line) optimized for the polarized LDOS in Sec. 4.2, showing that optimizing for a single dipole orientation (polarized LDOS) is nearly equivalent in performance to optimizing for all orientations (total LDOS).

Download Full Size | PDF

We found that the optimized LDOS comes within a factor of 10 of the upper bound in the short-wavelength regions ($\lambda <550$ nm), and the optimized cylinders are surprisingly good (within $\approx 20$% of the many-parameter optima). The optimized cavity geometries at $\lambda =500$ nm are shown in the inset of Fig. 2(b). We can see that the optimized many-parameter shape has a three-fold rotational symmetry around one axis; consequently, it has equal polarized LDOS in two directions but $\sim 100\times$ smaller polarized LDOS in the third direction. The spherical-harmonic basis is unitarily invariant under rotations, so this means that the optimization of the total LDOS (an isotropic figure of merit) exhibits a spontaneous symmetry breaking: it chooses two directions to improve at the expense of the third. For the optimized cylinder and ellipsoid, the polarized LDOS is only large for one polarization (along the cylinder axis, the “short” axis). A similar spontaneous symmetry breaking was observed for optimization of LDOS in two dimensions [14], as well as in saturating the upper bounds for scattering and absorption [1,45] (where it was related to quasi-static sum rules constraining polarizability resonances [45]).

As discussed below, we found that we could approach the polarized LDOS bound at $\lambda = 500$ nm within a factor of $\sim 3$ for a single dipole orientation; the fact that the total LDOS optimization is worse compared to its bounds ($\approx 3\times {}$ polarized bound) reiterates the conclusion that it is probably not generally possible to maximize the polarized LDOS for all three directions simultaneously.

One possible shape that will have same polarized LDOS in all directions is a sphere. As a matter of fact, we found that at each wavelength $\lesssim 600$ nm, there exists a resonant sphere [35,46] such that the LDOS at the center of the spheres comes within $\approx 20\%$ of the corresponding-$d$ bound (see Appendix C). However, this resonant radius $d$ is relatively large ($\lambda /4 < d < \lambda /2$, in order to create a resonance at $\lambda$) leading to small LDOS and bound ($\rho / \rho _0 \sim {}$10–100), so saturating such a large-$d$ bound may have limited utility. Spheres at much smaller $d$ do not exhibit these “void resonances” and have much worse LDOS than the asymmetrical shapes in Fig. 2 for $d=50$ nm.

4.2 Polarized LDOS

Since the spontaneous symmetry-breaking in the previous section suggests that optimization favors maximizing LDOS in a single direction, we now consider optimizing the polarized LDOS. That is, we maximize the power expended by a dipole current with a single orientation (similar to previous work on cavity optimization in dielectric media [14,32,33]). As above, we performed few-parameter optimization of ellipsoids, cylinders, and rectangular boxes. For many-parameter optimization, we initially used spherical harmonics but observed that optimizing polarized LDOS naturally leads to structures that are rotationally symmetric around the dipole axis. To exploit this fact, we switched to simple polynomials in $\theta$ as described in Sec. 3. Specifically, we first performed a rough scan of degree-2 polynomials to obtain a starting point, then we performed a degree-5 optimization using the adjoint method (degree-10 gave similar results at greater expense). (Gradually increasing the number of degrees of freedom is “successive refinement,” a heuristic that has also been used in other work to avoid poor local minima [47,48].) The results are shown in Fig. 3. We only plotted the cylinder results (orange dashed line), because the ellipsoid and box results were worse.

 

Fig. 3. (a) Polarized LDOS optima (dashed lines) as a function of the wavelength $\lambda$ at a minimum separation $d=50$ nm, along with the upper bound (black line) and the shape-dependent bounds (black dots). Dashed lines are peak performance of separately optimized structures for each $\lambda$, either cylinders (orange) or the optimized polynomials (blue) of Eq. (4). Solid blue lines are spectra of the polarized LDOS for optimized polynomial designs at selected wavelengths $\lambda = 400, 500, 600, 700, 800$ nm, respectively. (b) Optimized polynomial (top) and cylinder (bottom) structures at the wavelength $\lambda = 400, 500, 600, 700$ nm (to scale; see supplementary Visualization 3, Visualization 4, Visualization 5, Visualization 6, Visualization 7, Visualization 8, Visualization 9, and Visualization 10 for 3D views).

Download Full Size | PDF

We obtain an optimized LDOS within a factor of about 4 of the polarized-LDOS bound in the short wavelength regions ($\lambda <550$ nm). At a wavelength of 400 nm, the optimized LDOS is only 2.5 times smaller than the bound, which greatly improves upon previous results that often came only within $10^{2}$–$10^{3}$ of the bound [1,13,34]. One interesting fact is that the optimized polarized LDOS is actually only slightly smaller ($\approx 10\%$) than the optimized total LDOS (blue dashed line in Fig. 2(b)), which is consistent with the spontaneous symmetry breaking we commented on above: optimizing total LDOS spontaneously chooses one or two directions to optimize at the expense of all others, and hence is often equivalent to optimizing polarized LDOS. (If an isotropic LDOS is required by an application, one approach is to maximize the minimum of three polarized LDOSes [14].)

The upper bounds can help us to answer another important question: how much additional improvement could be obtained by introducing additional void structures outside of our cavity? (For example, by giving the cavity walls a finite thickness.) An upper bound to this improvement is provided by computing a shape-dependent limit: we use the same bounding procedure, but evaluate the limit assuming the material lies outside our optimized shape rather than outside of a bounding sphere. This analysis, which is carried out in Appendix A, shows that our optimized polarized LDOS is nearly reaching this shape-dependent limit as shown by the black dots in Fig. 3(a). Therefore, little further improvement is possible using additional structures outside of the cavity, which justifies optimizing over simple voids in order to probe the bounds.

We also explicitly studied the effect of a finite thickness for the metallic walls. To study the cavity thickness effect, we implemented simulations of cylindrical shells using the optimized cylinder at $\lambda =500$ nm taken from Fig. 2(b). We found that the LDOS increases monotonically with the shell thickness (Fig. 5 in Appendix D). A shell thickness about 100 nm yields a polarized LDOS within 5% of the infinite-thickness result, which is not surprising considering that the skin depth [37,49] of silver is $< 33$ nm for $\lambda > 400$ nm.

5. Concluding remarks

In this work, we obtain LDOS values within a factor of $\approx 10$ of the total LDOS bound and a factor of $\approx 4$ of the polarized LDOS bound in a many-parameter metal-cavity optimization, showing that these upper bounds are much more nearly attainable than was previously known [1]. Unlike previous work on scattering/absorption by small particles [45], our optimized cavities do not appear to be in the quasi-static regime, since their largest axes are $\sim \lambda /2$, even while their smallest axes are deeply subwavelength and exhibit strong plasmonic effects. It is possible that further improvements could be obtained by a more extensive search of local optima, or by expanding the search to other classes of cavities beyond “star-shaped” structures that can be described by a $R(\theta ,\phi )$ level set, e.g. via full 3D topology optimization. We would also like to systematically explore the attainability of the finite-bandwidth bounds from Ref. [13], which are useful for lossless dielectrics (where the single-frequency LDOS bound diverges). Conversely, it is possible that incorporating additional constraints, such as considering a more complete form of the optical theorem, may lower the LDOS bounds [20–23]. There has also been recent interest in the magnetic LDOS, corresponding to magnetic-dipole radiation [50], and we expect that qualitatively similar results (albeit with different optimal shapes) would be obtained for the magnetic LDOS (or some combination of magnetic and electric, although trying to optimize both simultaneously would likely encounter difficulties similar to those for optimizing multiple polarizations). The magnetic LDOS would merely require one to replace our electric-dipole source with a magnetic-dipole source, along with a similar switch of the Green’s function that appears in the derivation of the bounds [1].

On a more practical level, a possible next step is to maximize LDOS (or similar figures of merit) for 3D geometries more amenable to fabrication, whereas our goal in the present paper was to probe the fundamental LDOS limits without concern for fabrication. Fortunately, our results show that relatively simple (constant cross-section) shapes such as cylinders can perform nearly as well as the irregular shapes produced by many-parameter shape optimization, and are relatively insensitive to small details (e.g. curved or flat walls). This is a hopeful sign for adapting such cavities to nano-manufacturing by lithography or other techniques. And, although infinite-thickness cavities completely absorb the emitted power, our computation of the radiated power in finite-thickness shells (Appendix D) agrees with the theoretical bound’s prediction that the optimal radiated power is $\approx 1/4$ of the total [1].

Appendix A: LDOS limit in metallic cavity

In this appendix, we briefly review the evaluation of the LDOS upper bounds described in Ref. [1]. In particular, the total (electric) LDOS limit can be evaluated as an integral over the entire scattering volume $V$ (the region containing the material $\chi$):

(5)$$\frac{\rho}{\rho_0} \leq 1+\frac{k^{3}}{4\pi}\frac{|\chi(\omega)|^{2}}{\operatorname{Im}\chi(\omega)}\int_V\left[\frac{3}{(kr)^{6}}+\frac{1}{(kr)^{4}}+\frac{1}{(kr)^{2}}\right]\mathrm{d}^{3}r \;,$$

where $\rho _0=\omega ^{2}/2\pi ^{2}c^{3}$ is the free-space electric LDOS [7] and $k=\omega /c$ is the wavenumber.

Ostensibly, this limit is dependent on the exact scattering geometry $V$ (leading to a shape-dependent limit). However, Eq. (5) is also an upper bound on any scatterer contained within $V$ [1]. In this paper, we are interested in a minimal separation $d$ as depicted in Fig. 1, so we take $V$ to be a spherical shell with inner radius $d$ and shell thickness $L$, with $L\rightarrow \infty$ for arbitrary thickness. The integral of Eq. (5) can then be evaluated as

(6)$$\frac{\rho}{\rho_0} \leq 1+\frac{|\chi(\omega)|^{2}}{\operatorname{Im}\chi(\omega)}\left[\frac{1}{(kd)^{3}}+\frac{1}{kd}+\mathrm{O}(kL)\right]\;,$$

where $\mathrm {O}(kL)$ is a “Big-O” asymptotic bound [51]. As discussed in Ref. [1], the $\mathrm {O}(kL)$ divergence as $L\rightarrow \infty$, which arises from far-field scattering, is unphysical and overly optimistic. The contribution of this term should be limited by the largest interaction distance over which polarization currents contribute to the LDOS, thus is generally small compared to the first two terms on the right-hand side of Eq. (6) and can be neglected at small separation distance $d$. Therefore, the total LDOS limit for a metallic cavity with minimum separation distance $d$ is

(7)$$\frac{\rho}{\rho_0} \leq 1+\frac{|\chi(\omega)|^{2}}{\operatorname{Im}\chi(\omega)}\left[\frac{1}{(kd)^{3}}+\frac{1}{kd}\right]\;.$$

Note that this limit, where $V$ is the exterior of a sphere, is about 8 times larger than the limit discussed in Ref. [1] where $V$ was a planar half-space. In practice, this factor-of-8 improvement may be difficult to realize since optimized cavities will typically have only a small surface area at the minimum separation $d$ (except for the resonant spheres discussed in Appendix C below).

For the polarized LDOS limit, the integral in Eq. (5) (squared Frobenius norm of the homogeneous Green’s function [1]) is replaced with the norm of the dipole polarization vector multiplied by the Green’s function norm [12]:

(8)$$\frac{\rho_\mathrm{p}}{\rho_0} \leq \frac{1}{3}+\frac{k^{3}}{8\pi}\frac{|\chi(\omega)|^{2}}{\operatorname{Im}\chi(\omega)}\int_V\left[a(r)+b(r)|\hat{\mathbf{p}}\cdot\hat{\mathbf{r}}|^{2}\right]\mathrm{d}^{3}r \;,$$

where $\hat {\mathbf {p}}$ is the unit vector in the polarization direction, and $a(r)$ and $b(r)$ are:

(9)$$a(r) = \frac{1}{(kr)^{6}}-\frac{1}{(kr)^{4}}+\frac{1}{(kr)^{2}} $$

(10)$$b(r) = \frac{3}{(kr)^{6}}+\frac{5}{(kr)^{4}}-\frac{1}{(kr)^{2}} \;. $$

Similar to the total LDOS limit analysis, we can use a spherical bounding surface of radius $d$ to derive a general upper bound (also excluding the diverging $\mathrm {O}(kL$) term):

(11)$$\frac{\rho_\mathrm{p}}{\rho_0} \leq \frac{1}{3}+\frac{|\chi(\omega)|^{2}}{3\operatorname{Im}\chi(\omega)}\left[\frac{1}{(kd)^{3}}+\frac{1}{kd}\right]\;,$$

which is exactly $1/3$ of the total LDOS limit. That is, the total LDOS bound is equivalent to assuming that the polarized LDOS bound can be attained for all three polarizations simultaneously, which our results show to be unlikely.

To compute the shape-dependent polarized-LDOS bound (for a given optimized shape) in Sec. 4.2, we performed numerical integration of Eq. (8) over spherical angles (with the $r$ integral performed analytically), but excluding the $1/(kr)^{2}$ radiative term that yields the $\mathrm {O}(kL)$ divergence.

Appendix B: LDOS gradient from adjoint method

The gradient of the LDOS with respect to many shape parameters can be computed by solving Maxwell’s equations a single additional time (for “adjoint” fields) via the adjoint method, a key algorithm for large-scale photonics optimization [17,18]. The specific case of a boundary perturbation is reviewed in Ref. [52], which shows that the variation of an objective function $F$ in response to small shape deformations $\delta R$ (the surface displacement in the normal direction) over the surface $\partial \Omega$ is

(12)$$\delta F = 2\operatorname{Re}\iint_{\partial\Omega} \delta R(\mathbf{x}^{\prime}) \left[(1 -\varepsilon) \mathbf{E}_{\parallel}(\mathbf{x}^{\prime})\cdot\mathbf{E}^{\mathrm A}_{\parallel}(\mathbf{x}^{\prime}) + \left(\frac{1}{\varepsilon}-1\right)\mathbf{D}_{\perp}(\mathbf{x}^{\prime})\cdot\mathbf{D}^{\mathrm A}_{\perp}(\mathbf{x}^{\prime})\right]\mathrm{d}S,$$

where $\varepsilon$ is the electric permittivity of the metal, $\mathbf {E}_{\parallel }$ is the surface-parallel electric field, $\mathbf {D}_{\perp }$ is the surface-parallel displacement field, and the superscript “$\mathrm {A}$” denotes the adjoint field excited by an adjoint current source $\mathbf {J} = \partial F/\partial \mathbf {E}$.

In the case of LDOS, a further simplification arises. The objective function $F=\rho$ at position $\mathbf {x}_0$ can be expressed as [2]

(13)$$\rho=\operatorname{Im}\left[\frac{\epsilon_0}{\pi\omega}\sum_{j=1}^{3}\hat{\mathbf{s}}_j\cdot \mathbf{E}_j(\mathbf{x}_0)\right],$$

where $\mathbf {E}_j$ denotes the field excited by a unit dipole source at $\mathbf {x}_0$ polarized in the $\hat {\mathbf {s}}_j$ direction, and the sum over $j$ accounts for all three possible dipole orientations. We can see from Eq. (13) that the LDOS is proportional to the electric field, leading to an adjoint field that is also proportional to the original problem for each orientation $j$,

(14)$$\mathbf{E}^{\mathrm A}_{j}(\mathbf{x}) = \operatorname{Im}\left[\frac{\epsilon_0}{\pi\omega}\mathbf{E}_{j}\right]\;.$$

Inserting Eq. (14) into Eq. (12) gives us the LDOS gradient (first-order variation) with respect to any shape deformation:

(15)$$\delta \rho = \frac{\epsilon_0}{\pi\omega}\operatorname{Im}\sum_j\iint_{\partial\Omega} \delta R(\mathbf{x}^{\prime}) \left[(1 -\varepsilon) \mathbf{E}_{j\parallel}(\mathbf{x}^{\prime})^{2}+ \left(\frac{1}{\varepsilon}-1\right)\mathbf{D}_{j\perp}(\mathbf{x}^{\prime})^{2}\right]\mathrm{d}S\;.$$

Appendix C: Resonant sphere

For a void sphere cavity, the resonant electromagnetic surface modes can be analytically obtained by solving the equation [35,46] (after correcting a typographical error in Ref. [35]):

(16)$$\epsilon_\mathrm{m} (\omega) H_\ell(k_\mathrm{m}a)\left[k_\mathrm{d}aJ_\ell(k_\mathrm{d}a)\right]^{\prime} = \epsilon_\mathrm{d} J_\ell(k_\mathrm{d}a)\left[k_\mathrm{m}aH_\ell(k_\mathrm{m}a)\right]^{\prime} \;,$$

where $a$ corresponds to the void radius, $\ell$ is the (integer) index denoting the angular momentum, $k_{\mathrm {m}}=\sqrt {\epsilon _{\mathrm {m}}}k$ and $k_{\mathrm {d}} = \sqrt {\epsilon _{\mathrm {d}}}k$ are wave vectors in metal and air/vacuum respectively, $J_\ell$ and $H_\ell$ are spherical Bessel and Hankel functions of the first kind, and the prime denotes differentiation with respect to $k_{\mathrm {d}} a$ or $k_{\mathrm {m}} a$. Since the excitation source in our LDOS problem is a dipole at the center of the sphere, only an $\ell =1$ mode can be excited. Therefore, for each wavelength $2\pi /k$, there is a minimal resonant sphere: a minimal radius $d=a_{\mathrm {res}}$ satisfying Eq. (16) for $\ell = 1$.

Using Eq. (16), we can directly compute this minimal resonance radius $a_{\mathrm {res}}$ and then evaluate the corresponding LDOS. Furthermore, since a semi-analytical solution for the resonant modes is known for a void sphere [35], the LDOS can also be evaluated analytically via Eq. (1). We find that our BEM results are within $2\%$ of the analytical values using a surface mesh of $\approx 5000$ triangles, validating our numerical solver. The resulting LDOS values are shown in Fig. 4, which shows that the LDOS of the resonant sphere is very close to the theoretical limit (within $\approx 10$%) for the corresponding minimal separation $d = a_{\mathrm {res}}$. These strong results verify the limit at least at the resonance combinations of $d$ and $\lambda$.

 

Fig. 4. LDOS of a resonant air sphere in silver as a function of the wavelength $\lambda$ (blue line), where for each $\lambda$ we choose the smallest radius $a_{\mathrm {res}}$ for which we couple to a resonant mode at $\lambda$. The black line is the corresponding upper bound from Eq. (7), setting the minimum separation distance $d = a_{\mathrm {res}}$. The LDOS slightly exceeds the bound at small wavelengths where the radius becomes so large that one would need to include the $\mathrm {O}(kL)$ term that we dropped in Eq. (6).

Download Full Size | PDF

Notice that the resonant-sphere LDOS seems to actually slightly exceed the theoretical limit obtained with Eq. (7) in short wavelength range. There is no contradiction however: this is simply the effect of the $\mathrm {O}(kL)$ we dropped in Eq. (6). In particular, the resonant radius here is relatively large compared to the wavelength. For example, at $\lambda =700$ nm we get $a_{\mathrm {res}}=273$ nm, for which $kd=2.45$ and $1/(kd)^{3}+1/(kd)=0.47$. As a matter of fact, $1/(kd)^{3}+1/(kd)<1$ for all resonant spheres, thus the $\mathrm {O}(kL)$ term in Eq. (6) will have a non-negligible influence on the bound, causing the actual bound to be slightly higher than Eq. (7).

Appendix D: Cavity thickness effect

Here, we study the effect of a finite thickness of the metallic walls, replacing the infinite metallic regions of Fig. 1. To do this, we took the optimized cylinder at $\lambda =500$ nm from Fig. 2(b) and modified the silver walls to have finite thickness with the same inner surface. The LDOS as a function of the shell thickness is shown in Fig. 5. We observe that the LDOS increases monotonically with the shell thickness, and that a shell thickness of about 100 nm (about 3.7 times the skin depth) yields a polarized LDOS within 5% of the infinite-thickness result.

For a finite-thickness shell, some of the expended power (total LDOS) is absorbed and some “leaks” through the finite thickness to radiate away, and it is interesting to consider the radiative LDOS $\rho _{\mathrm {rad}}$ defined as the latter radiated power for the same dipole source (green line in Fig. 5). As a function of shell thickness, $\rho _{\mathrm {rad}}$ exhibits a peak: too thin and the resonance is too weak to enhance LDOS, but too thick and no power escapes to radiate (all power is absorbed). The theoretical limit for $\rho _{\mathrm {rad}}/\rho _0 - 1$ is $1/4$ of the limit for the total (absorbed+radiated) LDOS [1]. Correspondingly, we plot $4\rho _{\mathrm {rad}}/\rho _0$ in Fig. 5 (dashed green line) and see that, at the optimum $\rho _{\mathrm {rad}}$, the radiative LDOS is approximately $1/4$ of the total ($4\rho _{\mathrm {rad}} \approx \rho$), agreeing with the prediction of polarization-maximization in Ref. [1].

 

Fig. 5. LDOS of optimized cylindrical cavity a function of the wall thickness (blue line) at a wavelength $\lambda =500$ nm and a minimum separation $d=50$ nm, compared to the infinite-thickness LDOS (orange line). Also shown are the radiative LDOS $\rho _{\mathrm {rad}}/\rho _0$ (solid green line: the radiated/non-absorbed power), as well as $4\times \rho _{\mathrm {rad}}/\rho _0$ (dashed green line) because the theoretical bounds predict that the maximum $\rho _{\mathrm {rad}}$ is $1/4$ of the total LDOS [1].

Download Full Size | PDF

Appendix E: Coupled-mode theory for radiative LDOS

Given a purely absorbing resonant cavity such as the ones optimized in this paper, the addition of a radiation-loss channel (e.g. by coupling the cavity to a waveguide, permitting radiation through a small hole in the cavity walls, or simply thinning a portion of the wall as in Appendix D) can be analyzed quantitatively using the technique of temporal coupled-mode theory (TCMT) as long as the lifetime of the cavity remains long enough to be treated as an isolated resonant mode [24,53]. TCMT yields a result identical to that of Ref. [1]: the optimal radiated power is $1/4$ of the original absorbed power. We describe that straightforward analysis in this appendix, because it helps to connect the results in this paper with applications to radiative cavities.

In particular, consider a purely absorptive cavity with a quality factor [24,53] $Q_a \gg 1$, so that the Purcell enhancement of the LDOS for a dipole source is proportional to $Q_a/V$ where $V$ is a corresponding modal volume [6,24]. Now, suppose that one perturbs the cavity to add a radiative loss channel with a corresponding quality factor $Q_r$. As long as the radiation loss is low ($Q_r \gg 1$, as is necessary to retain a well-defined resonance), it can be treated as a small perturbation and the effect on $Q_a$ and $V$ can be neglected as higher-order in $1/Q_r$ [24]. With the addition of this channel, the total cavity $Q$ becomes $Q = Q_a Q_r / (Q_a + Q_r)$, corresponding to a total nondimensionalized loss rate of $1/Q = 1/Q_a + 1/Q_r$ [24], so the Purcell enhancement factor is reduced to $Q/V$. However, only a fraction $Q/Q_r$ of the dipole power goes into radiation (vs. absorption), so the resonant enhancement of the radiated power is proportional to

(17)$$\frac{Q^{2}}{VQ_r} = \frac{Q_a^{2} Q_r}{(Q_a + Q_r)^{2} V}.$$

Eq. (17) is maximized when $Q_r = Q_a$, i.e. when the absorptive and radiative loss rates are matched, similar to results obtained previously [25,26]. For $Q_r = Q_a$, Eq. (17) becomes $Q_a/4V$, or exactly $1/4$ of the $Q_a/V$ Purcell enhancement in the purely radiative case.

This result is consistent with the fact that our radiative-LDOS bound is $1/4$ of the total LDOS bound [1], but is more far-reaching. It prescribes how any high-$Q$ absorbing cavity can be converted to a radiating cavity with about $1/4$ the radiated power, similar to the results we obtained numerically in Appendix D.

Funding

Air Force Office of Scientific Research (FA9550-17-1-0093); Defense Advanced Research Projects Agency (HR0011-20-90016); Army Research Office (W911NF-13-D-0001).

Disclosures

The authors declare no conflicts of interest.

References

1. O. D. Miller, A. G. Polimeridis, M. T. H. Reid, C. W. Hsu, B. G. DeLacy, J. D. Joannopoulos, M. Soljačić, and S. G. Johnson, “Fundamental limits to optical response in absorptive systems,” Opt. Express 24(4), 3329–3364 (2016). [CrossRef]  

2. K. Joulain, R. Carminati, J.-P. Mulet, and J.-J. Greffet, “Definition and measurement of the local density of electromagnetic states close to an interface,” Phys. Rev. B 68(24), 245405 (2003). [CrossRef]  

3. L. Novotny and B. Hecht, Principles of Nano-Optics (Cambridge University, 2012), 2nd ed.

4. O. J. F. Martin and N. B. Piller, “Electromagnetic scattering in polarizable backgrounds,” Phys. Rev. E 58(3), 3909–3915 (1998). [CrossRef]  

5. 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(5), 057601 (2004). [CrossRef]  

6. A. Oskooi and S. G. Johnson, Electromagnetic wave source conditions, Advances in FDTD Computational Electrodynamics: Photonics and Nanotechnology, A. Taflove, A. Oskooi, and S. G. Johnson, eds. (Artech, Boston, 2013), chap. 4.

7. K. Joulain, J.-P. Mulet, F. Marquier, R. Carminati, and J.-J. Greffet, “Surface electromagnetic waves thermally excited: Radiative heat transfer, coherence properties and casimir forces revisited in the near field,” Surf. Sci. Rep. 57(3-4), 59–112 (2005). [CrossRef]  

8. R. R. Chance, A. Prock, and R. Silbey, Molecular Fluorescence and Energy Transfer Near Interfaces (John Wiley & Sons, Ltd, 2007), pp. 1–65.

9. F. Wijnands, J. B. Pendry, F. J. Garcia-Vidal, P. M. Bell, P. J. Roberts, and L. M. Moreno, “Green’s functions for maxwell’s equations: application to spontaneous emission,” Opt. Quantum Electron. 29(2), 199–216 (1997). [CrossRef]  

10. Y. Xu, R. K. Lee, and A. Yariv, “Quantum analysis and the classical analysis of spontaneous emission in a microcavity,” Phys. Rev. A 61(3), 033807 (2000). [CrossRef]  

11. M. Agio and D. M. Cano, “The purcell factor of nanoresonators,” Nat. Photonics 7(9), 674–675 (2013). [CrossRef]  

12. J. Michon, M. Benzaouia, W. Yao, O. D. Miller, and S. G. Johnson, “Limits to surface-enhanced raman scattering near arbitrary-shape scatterers,” Opt. Express 27(24), 35189–35202 (2019). [CrossRef]  

13. H. Shim, L. Fan, S. G. Johnson, and O. D. Miller, “Fundamental limits to near-field optical response over any bandwidth,” Phys. Rev. X 9(1), 011043 (2019). [CrossRef]  

14. X. Liang and S. G. Johnson, “Formulation for scalable optimization of microcavities via the frequency-averaged local density of states,” Opt. Express 21(25), 30812–30841 (2013). [CrossRef]  

15. P. K. Banerjee and R. Butterfield, Boundary element methods in engineering science (McGraw-Hill Book Co., 1981).

16. M. T. H. Reid and S. G. Johnson, “Efficient computation of power, force, and torque in bem scattering calculations,” IEEE Trans. Antennas Propag. 63(8), 3588–3598 (2015). [CrossRef]  

17. S. Molesky, Z. Lin, A. Y. Piggott, W. Jin, J. Vucković, and A. W. Rodriguez, “Inverse design in nanophotonics,” Nat. Photonics 12(11), 659–670 (2018). [CrossRef]  

18. G. Strang, Computational Science and Engineering (Wellesley, 2007).

19. D. P. Kingma and J. Ba, “Adam: A method for stochastic optimization,” (2014). Https://arxiv.org/abs/1412.6980.

20. S. Molesky, P. S. Venkataram, W. Jin, and A. W. Rodriguez, “Fundamental limits to radiative heat transfer: Theory,” Phys. Rev. B 101(3), 035408 (2020). [CrossRef]  

21. M. Gustafsson, K. Schab, L. Jelinek, and M. Capek, “Upper bounds on absorption and scattering,” New Journal of Physics (2020).

22. S. Molesky, P. Chao, and A. W. Rodriguez, “$\mathbb {T}$ operator bounds on electromagnetic power transfer: Application to far-field cross sections,” (2020). Https://arxiv.org/abs/2001.11531.

23. Z. Kuang, L. Zhang, and O. D. Miller, “Maximal single-frequency electromagnetic response,” (2020). Https://arxiv.org/abs/2002.00521.

24. J. D. Joannopoulos, S. G. Johnson, J. N. Winn, and R. D. Meade, Photonic Crystals: Molding the Flow of Light (Princeton University Press, 2008), 2nd ed.

25. Z. Ruan and S. Fan, “Superscattering of light from subwavelength nanostructures,” Phys. Rev. Lett. 105(1), 013901 (2010). [CrossRef]  

26. R. E. Hamam, A. Karalis, J. D. Joannopoulos, and M. Soljačić, “Coupled-mode theory for general free-space resonant scattering of waves,” Phys. Rev. A 75(5), 053801 (2007). [CrossRef]  

27. J. Vučković, M. Pelton, A. Scherer, and Y. Yamamoto, “Optimization of three-dimensional micropost microcavities for cavity quantum electrodynamics,” Phys. Rev. A 66(2), 023808 (2002). [CrossRef]  

28. J. Jensen and O. Sigmund, “Topology optimization for nano-photonics,” Laser Photonics Rev. 5(2), 308–321 (2011). [CrossRef]  

29. A. Mazaheri, H. R. Fallah, and J. Zarbakhsh, “Application of ldos and multipole expansion technique in optimization of photonic crystal designs,” Opt. Quant. Electron. 45(1), 67–77 (2013). [CrossRef]  

30. E. J. R. Vesseur, F. J. G. de Abajo, and A. Polman, “Broadband purcell enhancement in plasmonic ring cavities,” Phys. Rev. B 82(16), 165419 (2010). [CrossRef]  

31. T. W. Saucer and V. Sih, “Optimizing nanophotonic cavity designs with the gravitational search algorithm,” Opt. Express 21(18), 20831–20836 (2013). [CrossRef]  

32. F. Wang, R. E. Christiansen, Y. Yu, J. Mørk, and O. Sigmund, “Maximizing the quality factor to mode volume ratio for ultra-small photonic crystal cavities,” Appl. Phys. Lett. 113(24), 241101 (2018). [CrossRef]  

33. R. E. Christiansen, J. Michon, M. Benzaouia, O. Sigmund, and S. G. Johnson, “Inverse design of nanoparticles for enhanced raman scattering,” Opt. Express 28(4), 4444–4462 (2020). [CrossRef]  

34. O. D. Miller, O. Ilic, T. Christensen, M. T. H. Reid, H. A. Atwater, J. D. Joannopoulos, M. Soljačić, and S. G. Johnson, “Limits to the optical response of graphene and two-dimensional materials,” Nano Lett. 17(9), 5408–5415 (2017). [CrossRef]  

35. T. A. Kelf, Y. Sugawara, R. M. Cole, J. J. Baumberg, M. E. Abdelsalam, S. Cintra, S. Mahajan, A. E. Russell, and P. N. Bartlett, “Localized and delocalized plasmons in metallic nanovoids,” Phys. Rev. B 74(24), 245415 (2006). [CrossRef]  

36. F. Wijnands, J. B. Pendry, F. J. Garcia-Vidal, P. M. Bell, P. J. Roberts, and L. M. Moreno, “Green’s functions for Maxwell’s equations: Application to spontaneous emission,” Opt. Quantum Electron. 29(2), 199–216 (1997). [CrossRef]  

37. E. D. Palik, Handbook of Optical Constants of Solids (Elsevier Science, 1998).

38. S. Osher and R. Fedkiw, Level Set Methods and Dynamic Implicit Surfaces (SpringerVerlag, New York, 2003).

39. N. P. van Dijk, K. Maute, M. Langelaar, and F. van Keulen, “Level-set methods for structural topology optimization: a review,” Struct. Multidisc. Optim. 48(3), 437–472 (2013). [CrossRef]  

40. The CGAL Project, CGAL User and Reference Manual (CGAL Editorial Board, 2020), 5.0.2 ed.

41. J.-D. Boissonnat and S. Oudot, “Provably good sampling and meshing of surfaces,” Graph. Model. 67(5), 405–451 (2005). [CrossRef]  

42. D. A. Tortorelli and P. Michaleris, “Design sensitivity analysis: Overview and review,” Inverse Probl. Eng. 1(1), 71–105 (1994). [CrossRef]  

43. M. Heinkenschloss and L. N. Vicente, “Analysis of inexact trust-region sqp algorithms,” SIAM J. Optim. 12(2), 283–302 (2002). [CrossRef]  

44. S. Raza, S. I. Bozhevolnyi, M. Wubs, and N. A. Mortensen, “Nonlocal optical response in metallic nanostructures,” J. Phys.: Condens. Matter 27(18), 183204 (2015). [CrossRef]  

45. O. D. Miller, C. W. Hsu, M. T. H. Reid, W. Qiu, B. G. DeLacy, J. D. Joannopoulos, M. Soljačić, and S. G. Johnson, “Fundamental limits to extinction by metallic nanoparticles,” Phys. Rev. Lett. 112(12), 123903 (2014). [CrossRef]  

46. A. D. Boardman, Electromagnetic Surface Modes (Wiley-Interscience, 1982).

47. V. Sima, Algorithms for Linear-Quadratic Optimization (CRC Press, 1996).

48. F. Palacios-Gomez, L. Lasdon, and M. Engquist, “Nonlinear optimization by successive linear programming,” Manag. Sci. 28(10), 1106–1120 (1982). [CrossRef]  

49. J. D. Jackson, Classical Electrodynamics (Wiley, 1998), 3rd ed.

50. D. G. Baranov, R. S. Savelev, S. V. Li, A. E. Krasnok, and A. Alù, “Modifying magnetic dipole spontaneous emission with nanophotonic structures,” Laser Photonics Rev. 11(3), 1600268 (2017). [CrossRef]  

51. T. H. Cormen, C. E. Leiserson, R. L. Rivest, and C. Stein, Introduction to Algorithms (MIT, Cambridge, 2009), 3rd ed.

52. O. D. Miller, “Photonic design: From fundamental solar cell physics to computational inverse design,” Ph.D. thesis, Berkeley (2012).

53. W. Suh, Z. Wang, and S. Fan, “Temporal coupled-mode theory and the presence of non-orthogonal modes in lossless multimode cavities,” IEEE J. Quantum Electron. 40(10), 1511–1518 (2004). [CrossRef]