Maximum electromagnetic local density of states via material structuring

1 Introduction

The electromagnetic local density of states (LDOS), a measure of local response—the electric field from a dipolar current source—at a given point in space [1], plays a central role in many nanophotonic phenomena and applications, including spontaneous [2] and stimulated [3] emission, quantum information [4], surface enhanced Raman scattering [5, 6], photovoltaics [7], radiative heat transfer [8], nonlinear frequency conversion [9], scintillation [10], to name a few. Enhancing the LDOS by nanostructuring bulk media, a persistent theme in photonics design is often achieved through the creation of optical resonances [11]: a cavity or resonator supporting a single mode of quality factor Q and mode volume V enhances the LDOS in proportion to the Purcell factor Q/V. Many classic design schemes rely on maximizing Q (ring resonators [12]), minimizing V (plasmonic nanocavities [13, 14], bow-tie antennae [15]), or a combination thereof (photonic crystal cavities [16], metal–dielectric hybrid structures [17]). Going beyond the single mode picture, multi-mode and dispersion-based effects such as exceptional points [18, 19] and slow light devices [20] have also been explored. Recently, complex devices obtained by application of large-scale structural optimization appear to combine several confinement mechanisms [21, 22], providing further improvements on device performance.

Given its relevance in optics and the ever increasing structural freedom brought on by advances in nanofabrication [23], interest in assessing limits on achievable LDOS enhancement has grown over the past few decades [24]. Spectral sum rules [25–27] pin the frequency-integrated LDOS of any system but are of limited utility in practical settings involving finite and often narrow source bandwidths. Specialized results pertaining to achievable mode quality factors [28] or mode volumes in cavity settings [29] have proven useful in many applications but are not sufficiently general to account for multi-mode effects. Passivity requirements based on achievable material response were recently used to bind both single-frequency [30] and finite-bandwidth [31] LDOS. However, passivity alone cannot fully capture the wave nature of light as constrained by Maxwell’s equations, e.g., the necessity of phase matching to achieve resonance, leading to loose limits when contributions other than those coming from evanescent fields become relevant.

In this article, we generalize and lift several limitations imposed by these prior approaches to provide expressions and predictions for the largest spectral-integrated LDOS that may be achieved in a structured medium. The derived bounds incorporate wave and material constraints imposed by Maxwell’s equations over any desired length-scale and are geometry agnostic: besides specifying the material susceptibility χ, a frequency window of interest, and a bounding domain that the structured medium must reside within, no further assumptions are made on device topology. We consider sources both enclosed within and external to devices, obtaining bounds that generally come within an order of magnitude of complicated structures discovered through inverse methods [32]. Furthermore, we find that the mere requirement of energy conservation is sufficient to place tight limits on large devices; we exploit this fact to derive an analytical upper bound (13) along with integral expressions and asymptotic analysis for the particular scenario of a source above a semi-infinite, structured slab. The bounds show saturation to a finite value as χ → ∞ and varying power-law scalings with respect to source bandwidth Δω, with the maximum LDOS transitioning from a scaling ∝Δω⁻¹ to one ∝Δω^−1/4 as material absorption begins to limit the net enhancement that any one mode may contribute.

2 Problem formulation

Working in dimensionless units of ϵ₀ = μ₀ = 1, and considering only nonmagnetic materials, the partial LDOS at frequency ω and position x′ along the direction e ̂ can be shown to be, by Poynting’s theorem [33], directly proportional to the average power ρ(ω) emitted by a harmonic dipole source J ( r ) e − i ω t = e − i ω t δ ( x − x ′ ) e ̂ at the same location,

where the electric field generated by the J solves Maxwell’s equations,

Since real sources emit light over a finite bandwidth (e.g., the Planck spectrum of heated bodies or fluorescence/spontaneous emission linewidth of atoms [34]), we instead consider the more natural choice of a frequency average of ρ(ω) over a Lorentzian lineshape centered at ω₀ with bandwidth Δω_src = ω₀/2Q_src (corresponding source “quality factor” Q_src):

As pointed out in [32], such a modified figure of merit not only captures the spectral lineshape of many practical sources [34], but also offers a computational and conceptual advantage: instead of evaluating the frequency integral in (3a) directly, one can instead carry out a complex-ω contour integral over the upper half plane and exploit causality to convert, using the residue theorem, the spectral average ⟨ρ⟩ to the single complex-frequency ρ ( ω ̃ ) and an electrostatic (zero-frequency) contribution:

Here α = 1 2 R e p 0 ⋅ E 0 , where p₀ is a unit amplitude electrostatic dipole and E₀ the field it generates; this electrostatic contribution leads to an all-frequency integrated sum rule in the wide bandwidth limit [27, 31]. We note that (4) was previously derived in [31]; in this paper we will be focusing on the practically important case of moderate or narrow source bandwidth with Q_src ≫ 1 and thus will neglect the electrostatic term in (4) going forward.

Thus, the structural design problem for maximizing bandwidth-averaged LDOS can be formulated as

given the constraints

where V is a pre-specified design region that the structure resides within and χ ( ω ̃ ) the Fourier transform of the bulk material susceptibility evaluated at complex frequency ω ̃ . In later expressions we may suppress the ω ̃ dependence of time-harmonic variables for notation simplicity.

Due to the high dimensionality of the structural degrees of freedom ϵ(r) and the non-convex dependence of the field E(r) on ϵ(r), it is generally not possible to solve for the global optimum of (5c) [35]. However, it is possible to exploit an alternative parametrization of the problem in which the polarization density P(r) = (ϵ(r) − 1) E(r) rather than ϵ(r) serve as optimizations degrees of freedom, to obtain bounds on ρ ( ω ̃ ) applicable to any possible structure, given only the choices of design region V and material susceptibility χ [24]. In this context, it becomes convenient to decompose the total field into the vacuum field emitted by the source and scattered field emitted by the polarization currents:

with G ( r , r ′ ) denoting the vacuum dyadic Green’s function, the solution to

Note that our definition of the Green’s function has an extra global factor of ω ̃ 2 compared to the standard definition. Similarly, the net power extracted from the dipole source can be decomposed into constant (structure-independent) vacuum and scattered-field contributions:

where in the second line of (8c) we made use of the fact that J* = J (the global phase of dipole source is irrelevant) and the reciprocity relation G ( r , r ′ ) = G T ( r ′ , r ) . They key to formulating a shape-independent bound on ρ is to forgo the need for structural information by relaxing the requirement that P satisfy Maxwell’s equations everywhere, imposing instead a smaller subset of wave constraints [24]. Such constraints can be indeed be derived from Maxwell’s equations through a complex frequency generalization of the time-harmonic Poynting’s theorem (see [24] and SI), giving

where V _k ⊆ V is any spatial region within the design domain V and we have defined the composite operator U . Notably, (9) can be interpreted as a statement of the conservation of energy over each spatial region, requiring only specification of the allotted design footprint V and available susceptibility χ. Note also that geometric information is contained only implicitly in P, with material properties specified by the known complex scalar χ. The resulting optimization problem for the LDOS can be written as

such that ∀V _k ⊆ V

where I V k is the projection operator into the region V _k . Above, we employed Dirac bra-ket notation for brevity, with a denoting the vector field a(r) over the design domain V, ⟨a|b⟩ = ∫ _V a*(r) ⋅ b(r) dr is the conjugated inner product, and operator action A a = ∫ V A ( r , r ′ ) ⋅ a ( r ′ ) d r ′ . Equation (10a) is a quadratically constrained linear program for P. While direct optimization is still not guaranteed to find a global optimum due to the non-convexity of certain constraints in (10b), a bound on the problem can be computed efficiently via the Lagrange dual function [24, 36, 37], hereby referred to as the dual bound. Compared to prior bounds on ρ_sca based solely on passivity [31], the constraints in (10b) contain all relevant wave physics as encoded in the EM Green’s function G ; as will be seen in later sections this leads to qualitatively tighter limits.

In principle, the maximum LDOS achievable in a structured medium is known to grow indefinitively in the single-frequency limit of Δω_src → 0. For instance, in lossless media, the whispering gallery modes of ring resonators and defect modes of photonic crystals exhibit exponential scaling in Q with increasing system size [16, 38], while their mode volumes either increase polynomially or remain constant, respectively, leading to exponential growth in the Purcell factor. Geometries with sharp tips can also give rise to field singularities and hence vanishing mode volumes [39]. However, in practice, finite source bandwidths, device footprint, and material losses limit the utility of diverging lifetimes, while fabrication imperfections and atomic-scale effects preclude realization of arbitrarily small features. As will be seen below, the imposition of finite Δω_src and a minimum vacuum gap d between source and medium regularizes such divergences, paving the way for investigations of LDOS growth characteristics in realistic settings. We also note that for finite Δω_src the vacuum LDOS contribution ρ_vac diverges due to contributions from high frequencies. This is an artifact of the Lorentzian lineshape chosen and has little bearing on practical applications; the results presented will thus focus on the structural contribution ρ_sca which does remain finite.

To make explicit the scale invariance of Maxwell’s equations, all lengths are given in units of the center wavelength λ₀ = 2πc/ω₀ and ρ_sca normalized by the single-frequency vacuum dipole radiation ρ₀ ≡ ρ_vac(ω₀).

3 Numerical results

The dual bound can be obtained numerically for arbitrary domains via a suitable representation of the Green’s function and Maxwell operator. We do so for two standard settings: an external–dipole configuration where the dipole is adjacent to the structure, relevant for instance to gratings and solid-state defect couplers [40, 41]; and an enclosed–dipole configuration in which the dipole is surrounded by the structure, relevant to photonic crystal cavities [16], bullseye gratings [42], and bowtie antennas [15]. For computational convenience, all calculations are performed in 2d for both the out-of-plane (scalar) TM and in-plane (vectorial) TE electric-field polarizations. The finite-difference frequency-domain method is used to represent all relevant fields and operators (SI). Progressively tighter bounds were obtained by gradually increasing the number of constraint sub-regions V _k down to the computational pixel level. Inverse designs were also obtained to compare against the bounds, following the topology optimization (TO) approach detailed in [32].

Results pertaining to maximum LDOS for both external and enclosed configurations are shown in Figure 1 as a function of the design footprint. Since the vacuum LDOS diverges when integrated against a Lorentzian lineshape [25, 31, 32], this and subsequent results pertain only to the structural contribution ρ_sca, defining the enhancement factor in relation to the single-frequency vacuum emission ρ₀ ≡ ρ_vac(ω₀). In both settings the medium is assumed to be lossless, Im χ = 0. For a fixed Q_src, the bounds are seen to grow exponentially before saturating with increasing domain size L. Conversely, for fixed L the bounds grow linearly before saturating with increasing Q_src. Both observations are consistent with a resonant enhancement mechanism: since the finite vacuum gap d precludes field divergences in the vicinity of the dipole, thereby imposing an upper bound on its coupling to any mode, one would expect improvements in the Purcell factor to be mainly driven by growth in the modal lifetimes. Supposing a system with response dominated by a single mode of quality factor Q_mode = ω₀/2Δω_mode, (3a) yields a frequency-integrated LDOS that scales as

Thus, this modal picture predicts linear growth ∝Q_src in (3a) so long as the system supports a sufficiently long-lived mode, Q_mode ≫ Q_src, eventually saturating to a value ∝Q_mode whenever Q_mode ≪ Q_src. In the absence of material dissipation, the highest achievable Q_mode is constrained solely by radiative losses and, as confirmed by our bound calculations, known to scale at least exponentially with L in ring-resonator and photonic-crystal cavities. Comparing our bounds (solid lines) against inverse designs (squares), one observes remarkable alignment, with bounds and designs often coming within an order of magnitude of each other. The biggest performance gaps occur at high Q_src in the external configuration and can at least be partially attributed to TO getting trapped in local optima that, regardless of exhaustive initial conditions, underperform simple ring resonator geometries. Even starting with ring resonators as initial seeds, TO was only able to make modest improvements. In contrast, Bragg–onion type structures discovered by TO in the cavity setting are seen to tightly approach corresponding bounds. Comparing the high degree of precision needed to specify the dimensions of the ring resonator against the relative robustness of photonic bandgap confinement [43], we hypothesize that attaining high-performing designs for the external configuration at narrow bandwidths requires resonances based on sensitive interference cancellation. This leads to an ill-behaved optimization landscape with many subpar local optima to get stuck in, with ring resonators as an example of a class of high-performing designs not easily discoverable via brute-force optimization.

Figure 1:

Bounds and inverse designs with TM polarization, χ = 4, Q_src ∈ {10², 10⁴, 10⁶}, and finite design domain size L for (a) dipole on the side of a square design region and (b) dipole in the center of a square design region. Solid lines are bounds that are converged with respect to increasing the number of constraint subregions V _k . Dotted lines are bounds with just the global constraint V _k = V. Squares are topology optimized structures from random initializations. The hollow circle is a ring resonator with inner and outer radii approximately 1.85 and 2.05, respectively, with the exact resonant radius requiring accuracy up to 8 significant figures to achieve the plotted enhancement (SI); the filled circle is a topology optimized structure starting from that ring resonator. LDOS enhancement as a function of Q_src for fixed L = 1.6 is shown beneath the main figures, along with designs corresponding to filled shapes.

Finally, the observed dependence of the dual bound on the number and size of subregion constraints (10b) reveals two important features: first, the necessity of imposing many subwavelength constraint regions (SI) for the bounds to exhibit saturation as Q_src → ∞; second, the fact that a single, global energy-conservation constraint appears sufficient to produce tight bounds when either device sizes or source bandwidths become sufficiently large. This suggests that subregion constraints are necessary to adequately “resolve” wave effects, including phase-matching restrictions, which limit mode confinement in finite systems; conversely, when system sizes no longer limit the highest achievable Q_mode in relation to Q_src, the latter and not the former becomes a bottleneck for further enhancements.

4 Large-size semi-analytics

We now exploit the remarkable success of global energy conservation in describing large system behavior to obtain a semi-analytical bound on the maximum LDOS above a semi-infinite structure, i.e., the L → ∞ limit of the external–dipole configuration depicted in Figure 1a. As shown in SI, the solution of the dual problem under both resistive and reactive energy conservation constraints (10b) can be mapped to a family of solutions involving a unitary phasor parametrized by the single constraint angle θ, leading to the modified problem:

such that for P ∈ V,

where Asym ( A ) = ( A − A † ) / 2 i is the Hermitian anti-symmetric component of the operator A . Owing to its simplicity, the solution to this dual bound for any given θ can be written explicitly as

which notably, depends only on the analytically known quantities E_vac, χ, and U . The tightest dual bound consistent with (13) can thus be obtained numerically by carrying out an additional optimization over θ, restricted so that Asym ( e i θ U ) is positive definite (SI). Further analytical insight can be gleaned by expanding the operators and fields in a spectral basis conforming to the symmetry of the design domain: for a half-space enclosure, the natural choice is a Fourier basis e i k ‖ x ‖ parametrized by the wavevector k_‖ parallel to the half-space surface. Carrying out this expansion for the incident field E vac yields,

where k ⊥ = ω ̃ 2 − k ‖ 2 and Im(k_⊥) ≥ 0. For each conserved k_‖, the inverse operator image Asym ( e i θ U ) k ‖ − 1 ⋅ e i k ⊥ x ⊥ can be similarly expanded and evaluated explicitly as the sum of two complex sinusoids, leading to

where the closed-form complex coefficients R_1,2(k_‖; θ) and decay constants r_1,2(k_‖; θ) are given in the SI.

The main challenge in evaluating (13) comes from the need to compute the inverse image of Asym ( e i θ U ) = I m e i θ / χ * + Asym ( e − i θ G ) − 1 , which contains constraints on wave propagation via its dependence on the vacuum Green’s function G . In that regard, it is useful to consider an expansion of (13) in orders of Asym ( e − i θ G ) (after using the Cauchy–Schwartz inequality to relax the term involving E vac * (SI)):

The zeroth-order approximation given in (16) is especially simple to evaluate as it only depends on χ, and will hereafter be referred to as a material bound. A special yet sub-optimal value of θ in (16) recovers prior LDOS bounds based on passivity (SI) [31]. Notably, this zeroth order approximation becomes loose whenever I m e i θ / χ * − 1 Asym ( e − i θ G ) is large as compared to I ; this is the case for large |χ| or large design domains (e.g., the half-space) where Asym G dominates. Intuitively, terms containing the vacuum Green’s function capture physical limitations imposed by multiple scattering, screening, and the finite speed of light, which limit achievable LDOS.

Figure 2 shows upper bounds on the LDOS for both TM and TE sources, obtained by evaluating (15). In the near-field limit of d → 0, the TM bounds are found to asymptote to a constant while the TE bounds scale ∝1/d², in agreement with a detailed asymptotic analysis of the evanescent k_‖ → ∞ behavior of the integrand in (15) (SI), and matching prior results based on passivity [31]. The constant asymptote observed for TM sources as they approach the device reflects a known artifact of 2d scalar electromagnetism, which precludes non-integrable field singularities at sharp corners [44]. In contrast, the 1/d² scaling of TE sources confirms the known divergence in the energy concentration of a vector dipole in its near field, with the exponent of d tied to the number of spatial dimensions; note that the maximum LDOS for a 3d source grows ∝1/d³ as d → 0 (SI). While Figure 2 shows data only for lossless dielectrics, the asymptotics in the SI apply equally as well to metals with Re(χ) < 0, which exhibit the same scaling with d.

Figure 2:

Maximum LDOS enhancement near a half-space design region as a function of (a) separation distance and (b) real material susceptibility χ for lossless dielectrics. While (b) plots bounds for lossless dielectrics, the material independent limit as χ → ∞ is a bound for general lossy χ as well (SI).

In the opposite far-field limit of d → ∞ and for finite Q_src, the bounds ultimately tend toward zero, reflecting a kind of space–bandwidth constraint on the ability of structuring to affect far-field emission over a finite bandwidth. One can, however, observe plateaus occurring at intermediate regimes, ωd/c ≪ Q_src, wherein structuring can efficiently reflect traveling planewaves back onto the dipole position. Achievable enhancements in this regime can be shown to decay ∝ e − d / Q src , manifesting as plateaus in the various plots of Figure 2b. Only in the strictly single-frequency limit Q_src → ∞ do far-field LDOS contributions tend toward constants, ρ₀/2 and ρ₀ for TE and TM sources, respectively, independently of separation. The importance of properly capturing relevant wave interference effects in these far-field regimes becomes evident when analyzing corresponding predictions for material bounds [31], which exhibit drastically inflated ∝ Q src 2 scaling in this half-space setting (SI).

Asymptotic analysis of (15) supported by Figure 2b reveal that LDOS enhancements saturate to finite values as χ → ∞ for both dielectrics and metals (SI), in contrast to prior bounds which grow indefinitely with increasing material response [31]. This is somewhat surprising given that larger χ implies a larger (potentially infinite) density of states within the material itself; ultimately, multiple scattering and screening effects lead to restrictions on the possible field localization at the source location that cannot be overcome with clever structuring. Similar conclusions have been reached concerning the efficacy of hyperbolic metamaterials [45] at enhancing the LDOS, though we emphasize that our results are more general and make no assumption on the particulars of the device structure. It is also worth noting that the saturation characteristics of the full bounds as a function of d are distinct for the TM and TE polarizations. For TM, reducing d increases the relative advantage of stronger materials; ρ(χ → ∞) ∝ log(1/d) as d → 0 whereas ρ(χ < ∞) is finite as d → 0 (SI). In contrast, for TE the saturation behavior for large |χ| scales uniformly with d; decreasing separation does not increase the relative advantage of stronger materials.

With regard to bandwidth scaling, the broad structural freedom afforded by an infinite design domain coupled with the lack of material dissipation allows for creation of resonances with arbitrarily small radiative loss rates, Q_mode → ∞, leading to the same ∝Q_src dependence observed in Figure 1. The situation changes in dissipative media where mode lifetimes Q mode = ( 1 Q abs + 1 Q rad ) − 1 contain both radiative Q_rad and absorptive Q_abs contributions, with Q_abs ∝ 1/Im(χ) denoting the proportion of stored energy dissipated in the medium per unit cycle. Hence, under finite absorption, Q_abs sets a bound on the highest achievable mode lifetime, with (11) suggesting a saturating LDOS ∝Q_abs as Q_src → ∞.

In contrast, Figure 3 shows that the bounds exhibit a transition from the linear ∝Q_src scaling of a wide-bandwidth source coupled to a single resonance toward diverging response ∝ Q src 1 / 4 as Q_src → ∞, with the transition taking place as Q_src → Q_abs. Intuitively, in this regime wherein the lifetime of a single resonance is capped, one might expect the optimal design strategy to shift toward exploiting degeneracies. One class of geometry capable of achieving high modal degeneracies are waveguides, e.g., coupled cavities or photonic crystal gratings, which at least in 1d are known to generate van Hove singularities in the density of states [46]. Generally, a spectral singularity of the form |ω − ω₀|^−α, 0 < α < 1 will yield a Lorentzian-averaged LDOS that scales as

Figure 3:

Maximum LDOS enhancement near a halfspace design region as a function of the source bandwidth for dielectric Re(χ) = 5. Squares are the performance of TM inverse designs for χ = 5 + 1i, with points between A and B (inclusive) over a 10 × 10 design region and C, D over a 3 × 200 design region. Insets A and B show the structures found via TO for Q_src = 1 and Q_src = 100, respectively. Additional data pertaining to metallic χ is included in the SI.

The bounds for the 2d half-space configuration thus suggest the possibility of near-field absorbers supporting quartic “band-edge” dispersions capable of efficiently extracting evanescent fields, resulting in singularities of the form |ω − ω₀|^−1/4 (anomalous quartic ”slow” group-velocity dispersions have been studied, for instance, in the context of waveguide solitons [47]). Indeed, TO for χ = 5 + 1i discovers structures resembling adiabatic gratings (Figure 3 insets) with performances that approach the bounds and continue growing with increasing Q_src beyond what is expected of single mode enhancement given Q_abs on the order of 10. The remaining gap between inverse design and bound for narrow bandwidths is partially due to the finite size of the TO design domain, with finite domain bounds saturating as Q_src → ∞ (SI); further investigation is needed to clarify the tightness of the bound and the feasibility of engineering semi-infinite structures with the requisite quartic root spectral singularity.

5 Summary

We proposed an improved framework for evaluating upper bounds on the LDOS in structured media, and study in detail the effects of source bandwidth, material susceptibility, and device footprint. The bounds provide a strong top–down complement to bottom–up design in at least two ways: First, they are useful in quantifying the optimality of existing “best” structures while diagnosing potential areas of improvement: as seen in narrow bandwidth regimes, existing formulations of inverse design fail to converge upon structures with performance close to the global optimum, often outperformed by traditional ring resonators. Second, they allow studies of fundamental scaling characteristics without prior assumptions on device topology.

For finite bandwidths, the quick saturation of the bounds with design size indicate that a device footprint of a few wavelengths is generally enough to achieve near-optimal performance. The observed saturation with increasing susceptibility has direct implications for material choice: a weak material response can be mitigated given sufficient design size and judicious structuring; conversely there are diminishing returns associated with seeking large absolute values of the susceptibility, raising the importance of other concerns such as material loss and dispersion characteristics. The impact of material loss was also evaluated in the context of maximum LDOS above a semi-infinite structure, showing a transition from the intuitive linear ∝Q_src scaling expected of single-mode resonant enhancement to a less obvious ∝ Q src 1 / 4 dependence given absorptive materials, inviting further studies into the associated dispersion engineering mechanisms making such scaling possible.

While only results for dielectrics (Re(χ) > 0) where shown in the figures, the framework can handle metals (Re(χ) < 0) just as well. Metals may have an edge over dielectrics for electrically small design domains due to material resonances allowing for highly subwavelength confinement at the cost of high loss [48]. For large design domains such as the half-space, the bounds for metals have the same scaling with regards to source bandwidth and separation, and the existence of a material-independent limit blunts the advantage of the large χ of metals relative to their drawback of high loss.

Code for the computation of halfspace bounds makes use of the Python arbitrary precision arithmetic package mpmath [49] and is given at [50].