Approaching the upper limits of the local density of states via optimized metallic cavities

By computational optimization of air-void cavities in metallic substrates, we show that the local density of states (LDOS) can reach within a factor of $\approx 10$ of recent theoretical upper limits, and within a factor $\approx 4$ for the single-polarization LDOS, demonstrating that the theoretical limits are nearly attainable. Optimizing the total LDOS results in a spontaneous symmetry breaking where it is preferable to couple to a specific polarization. Moreover, simple shapes such as optimized cylinders attain nearly the performance of complicated many-parameter optima, suggesting that only one or two key parameters matter in order to approach the theoretical LDOS bounds for metallic resonators.


Introduction
Recently, we obtained theoretical upper bounds [1] to the (electric) local density of states (LDOS) ρ(x, ω), 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 x and frequency ω [2][3][4][5][6][7][8][9][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 ∼ | χ| 2 /Im χ/d 3 (reviewed in Sec. 2) depend on the material used (described by the ω-dependent susceptibility χ = ε − 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 ( χ, d, ω). 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 | χ| 2 /Im χ, and therefore the largest bounds) at wavelengths λ 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 ≈ 4 of the theoretical upper bounds, and total (all-polarization) LDOS within a factor of ≈ 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 ≈ 100 nm thickness at λ = 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 ≈ 20% 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).

O d
of the LDOS of optimized many-parameter irregular shapes, analogous to results we obtained previously for optimized scattering and absorption [1]. (Even a sphere is nearly optimal, but only for a specific separation d equal to a resonant-sphere radius for a given λ.) Moreover, we find that optimizing for a single emitter polarization (the "polarized" LDOS) does nearly as well (within ≈ 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][21][22][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.
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 [24][25][26][27][28][29], 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,30] 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,31], but recent topology optimization in two dimensions came within a factor of 10 of the 2d bound [30]. Semi-analytical calculations have also been published for resonant modes in spherical metallic voids [32], 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.

The local density of states (LDOS)
The (electric) LDOS is equivalent to the total power expended by three orthogonal dipole currents [2]: where 0 is the vacuum electric permittivity, E j denotes the field produced by a frequency-ω unit-dipole source at x 0 polarized in theŝ j direction, and the sum over j accounts for all three possible dipole orientations. We refer to the power Im[ 0 πωŝ j · E j (x 0 )] expended by only a single dipole current as the "polarized" LDOS.
From energy-conservation considerations, previous work found an upper bound for LDOS enhancement inside a cavity compared to vacuum electric LDOS (ρ 0 = ω 2 /2π 2 c 3 [7]), given any a material susceptibility χ and an emitter-material separation d at a frequency ω = ck, to be [1]: The details of this bound are reviewed in Appendix A. Two details in Eq. (2) require some comment. First, the bounding surface lying between the dipole source and the material for Eq. (2) is a sphere of radius d around the source (Fig. 1). If the bounding surface is a separating plane, there would be a factor of 1/8 multiplying the | χ| 2 /Im χ term as well as a small modification to the 1/kd term [1]. Second, the separation distance should be small compared to the wavelength, otherwise a third term O(k L) can have non-negligible contribution to the bound (also discussed in Appendix C). Finally, the polarized LDOS limit is 1/3 of the total limit in Eq. (2) for the same spherical bounding surface.
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][21][22][23], but our results below already show that Eq. (2) is achievable within an order of magnitude.

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 [33]. 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.
In order to parameterize an arbitrary cavity shape numerically, we use a level-set description [34,35], combined with a free-software surface-mesh generator CGAL [36,37]. In particular, we describe the radius of the cavity around the source point by a function R(θ, φ) (in spherical coordinates), which is expanded below in either spherical harmonics or other polynomials, and equivalently pass a level-set function Φ = r − R(θ, φ) to CGAL (such that Φ = 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 (θ, φ) as: ( For the basis functions S n , we used either spherical harmonics Y m (θ, φ) (for arbitrary asymmetrical "star-shaped" cavities) or simple polynomials in θ (to impose azimuthal symmetry round the z axis and a z = 0 mirror plane): The level set is discretized for BEM (by the CGAL software) into a triangular surface mesh. We used 5× greater resolution for surface points closer to the dipole source (radius 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 • and 120 • , 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 [38]. We tried various optimization algorithms designed to be robust to such "numerical noise" [19,39], and found that the Adam stochastic-optimization algorithm [19] seems to work best for our problem.

Total LDOS
We performed numerical shape optimization of the LDOS for cavities formed by voids in silver [33] at wavelengths λ from 400-900 nm, for both simple geometries (cylinder, ellipsoid, and rectangular box) and complex many-parameter shape (spherical harmonics). We chose an emitter-metal separation distance of 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 long-range effects [1] (see also Appendix A).
The results of the optimized LDOS as a function of the wavelength are displayed in Fig. 2a. Note that each wavelength corresponds to a different structure optimized for that particular wavelength. If we fix the structure as the one optimized for λ = 500 nm, the resulting LDOS spectrum is shown in Fig. 2b, 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 ∼ 10 random starting points and plotted the best result along with a few other typical local optima.
We found that the optimized LDOS comes within a factor of 10 of the upper bound in the short-wavelength regions (λ < 550 nm), and the optimized cylinders are surprisingly good (within ≈ 20% of the many-parameter optima). The optimized cavity geometries at λ = 500 nm are shown in the inset of Fig. 2b. 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 nearly-zero 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 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,40].
As discussed below, we found that we could approach the polarized LDOS bound at λ = 500 nm within a factor of ∼ 3 for a single dipole orientation; the fact that the total LDOS optimization is worse compared to its bounds (≈ 3 × 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 600 nm, there exists a resonant sphere [32,41] such that the LDOS at the center of the spheres comes within ≈ 20% of the corresponding-d bound (see Appendix C). However, this resonant radius d is relatively large (λ/4 < d < λ/2, in order to create a resonance at λ) leading to small LDOS and bound (ρ/ρ 0 ∼ 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.

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,29,30]). 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 θ 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 [42,43].) The results are shown in Fig. 3. We only plotted the cylinder results (orange dashed line), because the ellipsoid and box results were worse.
We obtain an optimized LDOS within a factor of about 4 of the polarized-LDOS bound in the short wavelength regions (λ < 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,31]. One interesting fact is that the optimized polarized LDOS is actually only slightly smaller (≈ 10%) than the optimized total LDOS (blue dashed line in Fig. 2b), 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. 3a. 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 λ = 500 nm taken from Fig. 2b. 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 [33,44] of silver is < 33 nm for λ > 400 nm.

Concluding Remarks
In this work, we obtain LDOS values within a factor of ≈ 10 of the total LDOS bound and a factor of ≈ 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. 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(θ, φ) level set, e.g. via full 3D topology optimization. Conversely, it is possible that incorporating additional constraints may lower the LDOS bounds [20][21][22][23]. 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 ≈ 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 χ): where ρ 0 = ω 2 /2π 2 c 3 is the free-space electric LDOS [7] and k = ω/c is the wavenumber. Ostensibly, this limit is dependent on the exact scattering geometry V (leading to a shapedependent 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 → ∞ for arbitrary thickness. The integral of Eq. (5) can then be evaluated as where O(k L) is a "Big-O" asymptotic bound [45]. As discussed in Ref. 1, the O(k L) divergence as L → ∞, 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 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]: wherep is the unit vector in the polarization direction, and a(r) and b(r) are: 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 O(k L) term): 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 O(k L) 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. 46, which shows that the variation of an objective function F in response to small shape deformations δR (the surface displacement in the normal direction) over the surface ∂Ω is where ε is the electric permittivity of the metal, E is the surface-parallel electric field, D ⊥ is the surface-parallel displacement field, and the superscript "A" denotes the adjoint field excited by an adjoint current source J = ∂F/∂E. In the case of LDOS, a further simplification arises. The objective function F = ρ at position x 0 can be expressed as [2] ρ = Im Fig. 4. LDOS of a resonant air sphere in silver as a function of the wavelength λ (blue line), where for each λ we choose the smallest radius a res for which we couple to a resonant mode at λ. The black line is the corresponding upper bound from Eq. (7), setting the minimum separation distance d = a res . The LDOS slightly exceeds the bound at small wavelengths where the radius becomes so large that one would need to include the O(k L) term that we dropped in Eq. (6).
where E j denotes the field excited by a unit dipole source at x 0 polarized in theŝ 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, Inserting Eq. (14) into Eq. (12) gives us the LDOS gradient (first-order variation) with respect to any shape deformation: