Fundamental limits to near-field optical response, over any bandwidth

We develop an analytical framework to derive upper bounds to light--matter interactions in the optical near field, where applications ranging from spontaneous-emission amplification to greater-than-blackbody heat transfer show transformative potential. Our framework connects the classic complex-analytic properties of causal fields with newly developed energy-conservation principles, resulting in a new class of power--bandwidth limits. These limits demonstrate the possibility of orders-of-magnitude enhancement in near-field optical response with the right combination of material and geometry. At specific frequency and bandwidth combinations, the bounds can be closely approached by canonical plasmonic geometries, with the opportunity for new designs to emerge away from those frequency ranges. Embedded in the bounds is a material"figure of merit,"which determines the maximum response of any material (metal/dielectric, bulk/2D, etc.), for any frequency and bandwidth. Our bounds on local density of states (LDOS) represent maximal spontaneous-emission enhancements, our bounds on cross density of states (CDOS) limit electromagnetic-field correlations, and our bounds on radiative heat transfer (RHT) represent the first such analytical rule, revealing fundamental limits relative to the classical Stefan--Boltzmann law.

Conceptually, the bounds arise because LDOS and related near-field quantities are given by the real (or imaginary) parts of causal linear-response functions. We use the complex-analytic properties of such functions to transform bandwidth-averaged response to that at a single, complexvalued frequency, where we develop generalized energyconservation constraints, ultimately leading to bounds over arbitrary bandwidths. A distinctive feature of our arbitrarybandwidth approach is that it predicts a simple material figure of merit (FOM) that determines the maximum possible response of any material (metal, dielectric, 2D, 3D, etc.), for any frequency and bandwidth. In the case of RHT, this FOM provides insight into which materials can facilitate optimal heat transfer for any temperature. There is significant ongoing debate about whether a plasmonic or an all-dielectric approach is better and in which scenarios 2D materials might be better than conventional bulk materials. Unlike all previous bounds and sum rules [24][25][26][27][28][29][30][31][32][33][34][35][36][37][38][39][40][41][42][43], the material figure of merit we derive here enables general quantitative answers to these questions. In a frequencybandwidth phase space, we map out which materials are optimal and where the critical thresholds, from dielectric to plasmonic and bulk to 2D, occur. The techniques developed herein for LDOS, CDOS, and radiative heat transfer should be extensible to other near-field effects ranging from engineered Lamb shifts [44,45] and Förster resonance energy transfer [46,47] to nonlinear (Raman) or fluctuation-induced (Casimir) phenomena.
Near-field electromagnetism, in which localized sources interact with scatterers separated by less than an optical wavelength, offers transformative potential for wideranging applications. Quantum emitters that only weakly couple to the radiation continuum can be dramatically amplified by near-field engineering: Optical antennas offer prospects for imaging single molecules [13,14,48,49] or for designing nanoscale light-emitting diodes that are faster than lasers [50]. Nonlinear emitters such as Raman-active molecules [51,52] experience even more dramatic enhancements: Surface-enhanced Raman scattering (SERS) [12][13][14], for example, scales with the square of the spontaneous-emission enhancement rate. Thermal emission can be accessed and controlled for the productive transfer of heat energy, at rates many orders of magnitude beyond classical blackbodies [18][19][20][21][22][23][24]. The emission can be stimulated by the vacuum itself: The field of Casimir physics is probing a vast expanse of materials and structures to explore how vacuum forces can be controlled and manipulated [7][8][9][53][54][55][56][57][58].
For such a broad scope of applications, there is a fundamental theoretical question that remains unanswered: For a given bandwidth of interest, what is the maximum near-field response that is possible? Sum rules enable at least a partial answer. Relying on the analytic properties embedded within particular response functions-such as susceptibilities and cross sections-sum rules relate integrated response over all frequencies to behavior at a single frequency, and they have been derived in a variety of classical and quantum frameworks [25][26][27][28][29][30][31][32][33]. In the near field, there is a well-known sum rule for spontaneous emission [34] that suggests spontaneous-emission enhancements average out to zero over all frequencies. Yet, this sum rule neglects the longitudinal component of the Green's function, thus neglecting the near-to-far-field coupling that is crucial for spontaneous-emission engineering (and hence recovering a far-field refractive-index sum rule). Very recently, a near-field sum rule for electric LDOS was derived [35], which represents a specialized version of a sum rule that we derive in Eq. (4). However, sum rules make no predictions as to how a finite available bandwidth affects maximal response. The difference between infinite and finite bandwidth is stark for "dielectric" (Re ϵ > 0) materials because infinite-bandwidth sum rules include technologically irrelevant high-frequency contributions that dominate the integrated response. For example, when applied to Silicon, the sum rule for plane-wave scattering is dominated by contributions at energies on the order of 100 eV [59], thus providing little insight into maximal response over typical bandwidths of interest (infrared, visible, etc.).
At the other extreme, single-frequency limits to power extinction and other physical observables have been discovered in both the near field [24,36,37,43] and far field [38][39][40][41][42]60] based on energy-conservation principles, but they necessarily fail to account for the effects of nonzero bandwidth. (As an example, they predict infinite maximal response for any lossless material. Such a prediction is in fact correct-it is possible to make LDOS arbitrarily large [61]. But the bounds developed herein show that the average response over any nonzero bandwidth has a finite upper bound). Thus, previous approaches do not provide any meaningful metric for lossless dielectrics at optical frequencies, either for quantitative comparisons amongst each other or to plasmonic and other metallic systems.
The key idea of our work is that two seemingly independent ideas-causality for sum rules and energyconservation principles for single-frequency bounds-can be unified into a single framework that yield bounds for any bandwidth, as illustrated in Fig. 1. In this framework, single-frequency bounds and all-frequency sum rules emerge as asymptotic limits of a more general arbitrarybandwidth approach. Our bounds over arbitrary bandwidths, which we term "power-bandwidth limits," arise by connecting the properties that enable sum rules to those that enable energy-conservation principles. Sum rules for power quantities (such as optical cross section) require one to be able to compute the quantity by taking the imaginary (or real) part of some amplitude-for extinction, the optical theorem [33,[62][63][64] guarantees such a form in terms of the scattering amplitude. The amplitude is a causal linearresponse function and thus analytic in the upper-half of the complex-frequency plane [65]. With suitable boundary conditions, a Hilbert transform (i.e., a Kramers-Kroniglike transform) then enables a sum rule, relating integrated response over all frequencies to that of a single frequency. Conversely, energy-conservation bounds-recognized primarily within the past decade [24,[36][37][38][39][40][41][42][43]60]-exploit the power-quantity-by-amplitude form in a different way. In writing a power quantity as the imaginary part of an amplitude, the amplitude itself is linear in the electromagnetic fields and/or currents (holding the incident field fixed). By contrast, many power quantities-absorption, scattering, radiation, etc.-are explicitly quadratic functionals of the fields and/or currents. If it can be shown that the linear quantity must be larger than the quadratic one, then an upper bound can be derived. Remarkably, it is the same optical theorem that provides such a constraint. Thus, we see that sum rules and single-frequency bounds both start with particular response functions that can be written as the imaginary part of an amplitude, but they diverge in their approaches thereafter. (For response functions that do not admit such expressions, the power-bandwidth limits established above can be applied by transformation to the real or imaginary part of a causal linear-response function. ) We connect the sum-rule and single-frequency approaches through the use of a "window function," an averaging function over a prescribed bandwidth. In general, such a function will have one (or multiple) poles in the upper-half plane (UHP), and thus the typical contourintegral analysis of a given power quantity requires the computation of residues not at a single real frequency (as in sum rules) but at a discrete set of complex frequencies. At this juncture, we identify the energy-conservation constraints at those complex frequencies, deriving bounds on how large they can be. This multistep procedure (fleshed out in detail below) thus provides perhaps the first approach to arbitrary-bandwidth bounds. For maximal clarity, we start with local density of states-the prototypical optical response-in Sec. II. We first derive near-field sum rules for LDOS (Sec. II A), showing that near-field response integrated over all frequencies must equal a new electrostatic constant, α LDOS . Then, we use geometric perturbation theory to prove a monotonicity theorem that enables us to bound the electrostatic constant itself in terms of only the material permittivity and near-field separation distance (Sec. II B). We introduce the window function in Sec. II C, combining it with the complex-frequency energyconservation idea to develop arbitrary-bandwidth bounds and show how closely they can be approached for specific choices of frequency and bandwidth by canonical structures. Having established the theoretical bound framework, we derive general bounds for cross density of states (Sec. III) and near-field radiative heat transfer (Sec. IV). Emerging within all these bounds is a common material figure of merit, and in Sec. V, we discuss the physical intuition of the FOM and compare a wide variety of materials across frequency and bandwidth. Finally, in Sec. VI, we discuss extensions of our formulation to near-field phenomena such as Lamb shifts, Raman scattering, Casimir forces, and more.

II. LOCAL DENSITY OF STATES
In this section, we develop a theoretical framework for upper bounds to near-field optical-response functions. The prototypical near-field interaction is the alterationand potentially dramatic enhancement-of spontaneous emission from a two-level dipolar transition in a quantum emitter by an inhomogeneous environment. The power radiated by such an emitter, and thus the spontaneousemission rate enhancement, is proportional to the LDOS [5,61,[66][67][68][69][70]. It has long been recognized that changing the electromagnetic environment of an emitter alters its spontaneous-emission rate [71][72][73]; applications where such enhancements play an important role include singlemolecule imaging [13,14], micro-LED design [74], and photovoltaics [75]. Mathematically, the spontaneousemission rate is determined by the imaginary part of the total Green's function [70,76,77]. To avoid unwieldy expressions and derivations, we use six-vector notation for fields and currents, treating the electric and magnetic fields, and electric and magnetic dipolar transitions, on equal footing. (We take the background to be vacuum throughout this paper and work in dimensionless units such that ε 0 ¼ μ 0 ¼ 1, with generalization to nonvacuum background in the Supplemental Material (SM) [78].) We denote the fields ψ, the polarization currents ϕ, and dipolar sources ξ: Then, the spontaneous-emission rate of randomly oriented electric (p) and magnetic (m) dipoles at a point x 0 is FIG. 1. Near-field optics exhibits phenomena ranging from spontaneous-emission enhancements through LDOS engineering, fieldcorrelation phenomena as measured by CDOS, RHT, Casimir effects, Smith-Purcell radiation, and Raman scattering. Our theoretical framework, connecting causality principles with energy-conservation constraints, yields bounds over any arbitrary bandwidth. In the limit of zero bandwidth, we obtain recently discovered single-frequency bounds [36,37]; in the infinite-bandwidth limit, we arrive at a sum rule for integrated all-frequency response.
where Γ denotes the 6 × 6 dyadic Green's function, the "s" subscripts denote scattered-field contributions (thus, Γ s is the total Green's function minus the free-space Green's function), and we use the convention [68] that LDOS is the sum of electric-dipole and magnetic-dipole contributions. It is important to subtract the free-space rate and consider only the scattered-field contribution, to ensure sufficiently fast decay at high frequency and thus convergence of integrals over frequency. The random dipole orientation (for dipoles of unit amplitude, i.e., ξ † ξ ¼ 1) is encoded in the summation over directions j ¼ fx; y; zg and, ultimately, the trace of the Green's function. In Eq. (2), we denote a term sðωÞ (suppressing the implicit position dependence), which we identify as a near-field scattering amplitude, as measured at the location of the emitter. It is this term that enables the sum rule and the power-bandwidth limits. Maxwell's equations do not prevent us from taking the frequency to be complex. In the complex-ω plane, we can define the complex extension of the near-field scattering amplitude as where we have made the typical assumption that the dipole amplitudes are real valued, such that ξ † ¼ ξ T . (For complex dipole amplitudes, a few additional steps in the derivations below are needed, but the results remain unchanged.) Since ξ is constant (analytic everywhere), and the scattered field ψ s is a causal linear-response function [65], the amplitude sðωÞ is analytic in the upper half of the complex-ω plane. This is analogous to the classical result that quantities such as refractive index and far-field scattering amplitudes are analytic in the upper half plane (UHP) [65]. On the real line, sðωÞ has a pole at the origin due to the singularity in the 1=πω prefactor and the fact that ξ T ψ s (usually) has a nonzero value in the zero-frequency limit.

A. Sum rules
Thus, a sum rule can be derived for ρðx 0 ; ωÞ through contour integration of the scattering amplitude sðωÞ in the UHP. If we enclose the UHP in a typical contour that is semicircular going to infinity and follows the real line with a "bump" at the origin [see Fig. 2(a)], then analyticity ensures that the total integral is zero. For local, linear susceptibilities, sðωÞ falls off sufficiently rapidly as ω → ∞ (because the free-space contribution was subtracted out), such that the real-line integral is well defined and the contour at infinity does not contribute (shown explicitly in the SM [78]). Thus, there are two contributions to the integral: the (principal value of the) integral over the real line and the residue of the simple pole at zero. Because negative-real-frequency fields are conjugates of their positive counterparts [79], ρð−ωÞ ¼ ρðωÞ, such that the real-line integral can be reduced to only positive frequencies, after which a few algebraic steps FIG. 2. (a) Contour of integration in the complex-ω plane used to obtain the sum rule for LDOS, which has a simple pole at the origin and is analytic everywhere on the upper half plane. For high-symmetry geometries such as a half-space or planar sheet, the sum-rule constant is known analytically. (b) Electric and (c) magnetic LDOS for Ag, Al, Au (half-spaces), and graphene (planar sheet). Although they all have resonant peaks at different frequencies with varying amplitudes and widths, their integrated response converges to the same constant, 1=16πd 3 , as shown in the inset. The emitter-scatterer distance d is set to 10 nm, and the Fermi level for graphene is set to 0.6 eV.
(SM [78]) yield a general expression for the value of ρðx 0 ; ωÞ integrated over all frequencies: where α LDOS ¼ 1 2 Re½TrΓ s ðx 0 ; x 0 Þj ω¼0 . The electrostatic constant α LDOS measures the scattered field at the position of the dipole source and is shown in the SM [78] to always be positive. Equation (4) marks our first result: Over all frequencies, integrated LDOS must equal an electrostatic constant. [An electric-only specialization of Eq. (4) was very recently discovered [35], albeit without the analytical bounds to follow.] For materials with an electrical conductivity (e.g., metals), α LDOS is independent of the value of the conductivity, for any geometry. More generally, Eq. (4) applies to any material (whose susceptibility decays in the limit ω → ∞), including the wide array of newly emerging 2D materials.
Some care is required with Eq. (4) in singular situations. At high-symmetry points near high-symmetry scatterers, e.g., at the center of a hemispherical bowl, LDOS may appear to not decay at frequencies going to infinity because the optical rays reflect off the perfect spherical interface and constructively interfere at the origin. However, any random deviation from a hemisphere, no matter how small, would destroy such interference at high enough frequencies, restoring the natural, sufficiently rapid decay. Hence, the correct approach to regularizing such singularities is to compute the sum rule for a hemisphere with imperfections and then take the limit as the imperfections go to zero, such that Eq. (4) still applies. (Such a procedure is a geometrical analog of the limiting absorption principle [80], defining "lossless" materials in the limit as loss goes to zero from above.) For any scatterer, the constant α LDOS can be found from an electrostatic calculation. In high-symmetry geometries, the calculation may be analytically tractable. Consider a half-space of permittivity ε and permeability μ. (In this paper, we consider materials that have scalar material susceptibilities, with tensor-valued generalizations in the SM [78].) The value of α LDOS at a separation d can be computed via the image charge within the half-space, leading to an expression (derived in the SM [78]) for α LDOS , and thus of the frequency-integrated LDOS, of where εð0Þ and μð0Þ are the zero-frequency (electrostatic or magnetostatic) permittivities and permeabilities, respectively. For metals and any material with a conductivity, the permittivity and/or permeability diverges in the zerofrequency limit, such that the corresponding term in square brackets in Eq. (5) simplifies to 1. This result is also the case for any 2D conductive sheet, which at zero frequency represents the same perfect-conductor boundary condition as a conductive 3D half-space. In the SM [78], we also derive the simple α LDOS expression for conductive spheres. An interesting feature of the LDOS sum rule is that it depends only on zero-frequency behavior, where electric and magnetic fields decouple. For the nonmagnetic materials that are ubiquitous at optical frequencies, this implies very different behavior for electric-dipole sources (i.e., electric-dipole transitions) as compared to magnetic-dipole sources (magnetic-dipole transitions). To illustrate the difference, one can separately define electric LDOS ρ E as the component arising from the electric sources only, and magnetic LDOS ρ H as arising from the magnetic sources only: where G EE and G HH are the electric and magnetic dyadic Green's functions, respectively. The two terms in the square brackets in Eq. (5) correspond to these individual LDOS constituents. For a nonmagnetic half-space (which we take as conductive just for simplicity), the electric and magnetic LDOS sum rules are The magnetic LDOS must average out to zero because, in the zero-frequency limit, a magnetic dipole does not interact with a nonmagnetic medium. Figures 2(b) and 2(c) illustrate the generality of the electric and magnetic LDOS sum rules for bulk metals (Ag, Al, Au) [81] as well as 2D materials such as graphene (material model from Refs. [82,83]). In our figures, we normalize electric and magnetic LDOS by the free-space electric (or magnetic) LDOS, ρ 0 ¼ ω 2 =ð2π 2 c 3 Þ which is half that from Ref. [68] as they consider total free-space LDOS. Each of these materials supports surface plasmonpolaritons [84], which are excited by near-field sources. These materials exhibit very different resonant frequencies and line widths, as seen in Figs. 2(b) and 2(c). In terms of electric LDOS in Fig. 2(b), graphene exhibits a very large and narrow-band response at infrared frequencies, whereas metals exhibit varying levels of maximum response, with inversely proportional bandwidths, at visible or ultraviolet energies. In contrast to the large order-of-magnitude enhancements for electric LDOS, the magnetic LDOS in Fig. 2(c) shows only limited response-note the scale of the y axis in Fig. 2(c) relative to Fig. 2(b). The modest, fluctuating magnetic LDOS arises because of the small electric field generated by the magnetic source or, equivalently, because the magnetic source cannot induce any magnetization in nonmagnetic media. Across the wide variations of response seen for both electric and magnetic LDOS, for systems with different materials and dimensionality, the all-frequency response must converge to the sum rules of Eqs. (4)- (6), as shown in the insets of Figs. 2(b) and 2(c).

B. All-frequency bounds
Equation (4) is an equality for any geometry. For structures without a high degree of symmetry, the electrostatic constant α LDOS would typically require an electrostatic computation. In this section, we use perturbation theory to derive a "monotonicity theorem," showing that if one material body (with static permittivity and permeability greater than 1) encloses a second body of the same material, the electrostatic constant α LDOS must be larger for the first than for the second. With this result, we can bound the allfrequency integrated response for any geometry in terms of the analytically known α LDOS for high-symmetry enclosures, yielding general analytical bounds.
Any outward deformation can be broken down into a sequence of outward perturbations. Thus, if one can show that any outward perturbation of a geometry must increase some response function, a "mononicity theorem" has been proved, guaranteeing that such a function increases for any outward deformation. Such theorems are known for electrostatic polarizability under plane-wave excitations [31,85,86]. To understand how α LDOS changes under geometrical perturbations, we use a variational-calculus approach applicable to any frequency and then isolate the electrostatic behavior. Within the variational-calculus approach, quantities incorporating the displacement fields D and B as well as the scalar permittivity and permeability will be necessary, for which we define the six-vector field Ψ and tensor ν: where I is the 3 × 3 identity matrix (and as discussed above, generalizations to anisotropic materials are included in the SM [78]). Consider a scattering FOM such as LDOS. Under geometrical perturbations, the variation in the FOM can be written as an overlap integral between two fields: (1) a "direct" field, which is the response of the unperturbed geometry to the original source (e.g., a nearby dipolar emitter), and (2) an "adjoint" field, which is the response of the same unperturbed geometry to a source whose phase, amplitude, and position depend on the precise FOM [87,88]. For any FOM, if we write the displacement of the boundary in the normal direction as Δh n ðxÞ, the variation in the FOM can generally be written as an overlap integral of the direct and adjoint fields over the geometrical boundary [88]: where (similarly for μ), Δν ¼ ν 1 − ν 0 , and ν 1 and ν 0 represent the material properties of the scatterer and its surroundings, respectively, while the "adj" superscript denotes adjoint field solutions. Implied in the above integral over the boundary is multiplication with an area element dA along the boundary, where ψ k and ψ ⊥ denote the field components parallel and perpendicular to the locally flat boundary, respectively. We have explicitly used the electrostatic constant α LDOS for the figure of merit since it is that constant for which monotonicity will apply. For any figure of merit f, the adjoint fields are a solution with source currents given by the functional derivatives ∂f=∂ψ T . From Eq. (3) and the discussion following Eq. (4), we know that α LDOS is given by which means that for any given dipole orientation j, the adjoint source field is given by ∂f=∂ψ T j ¼ 1 2 ξ j ðx 0 Þ. This shows a unique circumstance: The dipolar sources for the adjoint field are exactly half of ξ j ðx 0 Þ, which were the original LDOS dipolar excitations, such that ψ ðadjÞ ¼ 1 2 ψ and Ψ ðadjÞ ¼ 1 2 Ψ. Moreover, at zero frequency, without any material or radiative losses, the fields can be chosen to be real valued. Thus, for materials with positive static permittivities and permeabilities that are greater than those of their surroundings at zero frequency (ν 0 , ν 1 , and Δν positivedefinite), Eq. (8) can be written as the integral of a positive quantity, Equation (10) ensures that Δα LDOS is positive for any outward deformation (Δh n > 0 everywhere on the boundary). Regardless of the size and shape of a given scatterer By connecting this monotonicity theorem to the sum rule in Eq. (4), one can see that the integrated LDOS near a scatterer must increase with the size or shape of that scatterer. Note that although our derivation does not strictly apply to the limiting case of a conductive material with Δν → ∞ (because ε → ∞ or μ → ∞) in the zero-frequency limit, it does apply for any arbitrarily large but finite material susceptibility, and in the SM [78], we provide an alternative proof that confirms the validity of Eq. (11) for conductive materials.
The use of Eq. (8) assumes a smooth perturbation of the boundary. In the case of a region with "kinks," or sharp corners, such a boundary represents the limit of smooth deformations. Since the fields are finite and the discontinuity is a region of zero measure, it would not contribute to first order [87], and monotonicity would hold. (It is not clear whether monotonicity would hold for fractal surfaces.) Although not covered explicitly by our derivation, monotonicity would also hold for a gradient-index (at zero frequency) medium, if the increase in index is non-negative everywhere (and positive over some region) across the medium.
A simple application of the monotonicity theorem is to enclose the scatterer within some larger half-space, which is possible as long as there is a separating plane between the emitter and the scatterer. The electrostatic constant of any half-space is given in Eq. (5), involving only the separation distance and the zero-frequency material parameters. Since the material factors ðε − 1Þ=ðε þ 1Þ and ðμ − 1Þ=ðμ þ 1Þ are bounded above by 1 (for static material constants larger than 1), we can replace them with 1 in the upper bounds. By the monotonicity theorem, the all-frequency integrated LDOS for any structure enclosed by a half-space (which is separated from the emitter by a minimum distance d) must obey the general bounds for any material: and, for nonmagnetic materials: If the scatterer in question can be more tightly enclosed by another shape, such as a sphere or two half-spaces, one can replace the rhs of Eqs. (12) and (13) with the respective electrostatic constants to obtain a tighter bound on α LDOS . Figure 3 shows the electric LDOS for an emitter at the center of an aluminum double cone (similar to a bowtie antenna [89]), computed with an open-source software implementation [90,91] of the boundary element method (BEM) [92]. Figure 3(a) shows that for an emitter-antenna separation of d ¼ 10 nm (a cone-cone separation of 20 nm), thousandfold enhancements in electric LDOS are possible, at resonant frequencies determined by the geometry and the opening angle θ. Enlarging the opening angle represents a way to increase the size of the scatterer, and Fig. 3(b) demonstrates a monotonic increase in integrated electric LDOS, as expected from the monotonicity theorem. In conjunction, the sum rule, Eq. (4), and the monotonicity theorem, Eq. (11), suggest a critical takeaway: Isolated sharp tips do not occupy enough of the near field to maximize electric LDOS; instead, large-area structured surfaces offer significantly greater potential.
Our identification of the causal linear-response function sðωÞ defined in Eq. (2) as the underpinning of near-field sum rules ultimately yields a relation between all-frequency response and the single pole at the origin. Such relations form the crux of all sum rules [65], where the pole is almost always chosen at the origin or in the limit ω → ∞ (on the real line) because the response at those two particular frequencies can often be simplified: ω ¼ 0 is the regime of electrostatics, while electromagnetism in the ω → ∞ limit is perturbative, as material susceptibilities converge to zero. From a theoretical viewpoint, of course, a pole can be introduced anywhere on the real line (not just ω → 0; ∞), or even anywhere in the UHP. Typically, however, one cannot make general statements about the response at arbitrary frequencies. In the next section, we show how to employ recently developed energy-conservation techniques to derive general bounds at such frequencies, moving beyond the single-frequency or all-frequency dichotomy to a framework that works for any bandwidth.

C. Power-bandwidth limits
In this section, we introduce two ideas that transform the sum-rule approach of Sec. II A to an approach that bounds the response over any bandwidth: (1) We use the notion of a window function to connect the average response over some bandwidth to discrete frequencies in the upper half plane (at the window function's poles), and (2) we show how energy-conservation and passivity constraints can be applied at those complex frequencies, yielding analytical bounds on the bandwidth-averaged response.
Bandwidth plays a key role in any electromagnetic scattering problem, whether arising from the intrinsic line width of a source or as a primary technological constraint. For example, enhancements in LDOS over a broad spectrum could enable simultaneous imaging of multiple molecular species. Also, broadband LDOS enhancements provide a key criterion in designing optimal photovoltaic structures capable of maximal absorption enhancements [93].
There are many ways in which one might want to average the response over bandwidth (equal weighting, Lorentzian weighting, etc.), and one can accommodate almost any by prescribing a window function H ω 0 ;Δω ðωÞ that serves as a weighting function-it is concentrated around a center frequency ω 0 , defined by a frequency width Δω, and normalized ( R ∞ −∞ H ω 0 ;Δω ðωÞdω ¼ 1). Then, the average LDOS over that range of frequencies, which we denote hρi, can be defined by For the remainder of this section, we choose a Lorentzian function for H ω 0 ;Δω : where Δω is the half width at half maximum. We use a Lorentzian for simplicity: H ω 0 ;Δω ðωÞ extended into the UHP has only a single pole at ω ¼ ω 0 þ iΔω. Other window functions can be used with the simple addition of extra (or higher-order [61]) poles. For example, one can approximate a rect function (where H ω 0 ;Δω is constant over the bandwidth of interest and zero elsewhere) by the generalization H ω 0 ;Δω ðωÞ ¼ ½cðΔωÞ 2n−1 =½ðω − ω 0 Þ 2n þ ðΔωÞ 2n , which has n poles in the UHP. The bounds below at a single pole would then become a linear combination of the bounds at the new poles, with slightly modified numerical factors but the same physical ramifications for material and structural design.
Inserting the Lorentzian of Eq. (15) into the average LDOS [Eq. (14)] and writing the LDOS in terms of the near-field scattering amplitude, ρ ¼ Im sðωÞ, the product ρðωÞHðωÞ in the averaging integrand is given by Im½sðωÞHðωÞ. Outside of the lower-half plane, the function sðωÞHðωÞ has two simple poles: one at the origin, which was responsible for the sum rule from Sec. II A, and another at ω 0 þ iΔω, from the Lorentzian, as shown in Fig. 4(a). One can integrate over the contour in Fig. 4(a) and use Cauchy's residue theorem to equate the allfrequency integral of Eq. (14) to the evaluation of two complex-frequency quantities: Contour of integration in the complex-ω plane used to obtain the average LDOS, which contains singularities at the origin ("LDOS pole") that is intrinsic to LDOS and at a complex frequency ("Lorentzian pole") determined by the parameters chosen for the Lorentzian window function H ω 0 ;Δω . Apart from these two poles, the product ρH ω 0 ;Δω is analytic everywhere on the upper half plane.
(b) Average electric and (c) magnetic LDOS near a Ag half-space centered at different frequencies around its peak (3.65 eV), compared to their respective bounds. Taking the bandwidth to zero gives the single-frequency limit found in earlier works [36,37]. In the opposite limit of infinitely large bandwidth that includes the entire LDOS spectrum, our bounds reproduce the electric and magnetic sum rules in Sec. II A. The emitter-scatterer distance d is set to 10 nm.
where sðω 0 þ iΔωÞ is the near-field scattering amplitude evaluated at the single complex frequency ω 0 þ iΔω, and α LDOS is the electrostatic constant appearing in the sum rule defined in Eq. (4). In the limit of zero bandwidth, H ω 0 ;Δω ð0Þ equals zero, and Eq. (16) comprises only the first term, which represents single-frequency LDOS: ρ ¼ Im sðωÞ; conversely, as the bandwidth goes to infinity, Im sðω 0 þ iΔωÞ decays rapidly (as we show below) and the second term converges to the sum rule [Eq. (4)]. Between these extremes, the two terms comprising hρi combine to represent a bandwidth-averaged response. At the complex frequency ω 0 þ iΔω, the positive imaginary part of the wave number k means that the incident field emanating from the dipolar source is exponentially decaying, as can be understood from the outwardgoing wave e ikr =r → e iω 0 r=c e −Δωr=c . The decaying source is mathematically equivalent to a scattering problem in which the frequency is real valued but material loss is increased in both the scatterer and the background [94]. Through either viewpoint, one can see that large broadbandwidth response is inherently more difficult to achieve than large single-frequency response due to an inherent bandwidth-induced dissipation.
By expressing the weighted integral of LDOS in terms of residues evaluated at single complex frequencies, Eq. (16) is now conceptually similar to single-frequency response functions at real frequencies for which we have developed an energy-conservation or passivity-based approach to upper bounds [36,37]. The key idea as applied here is that Eq. (16) is the imaginary part of a function that is linear in the scattering amplitude sðω 0 þ iΔωÞ, while representing the total (bandwidth-averaged) power lost by the dipolar source, to either far-field radiation or near-field absorption. By contrast, absorption itself is dissipation in the medium, computed with the imaginary part of the susceptibility and the field intensity jEj 2 , a function that is quadratic in the fields. Since absorption must be smaller than total LDOS (absorption þ radiation), this implies that the quadratic functional must be smaller than the linear functional, a convex optimization constraint that necessarily bounds how large the fields and induced currents can be. We defer to Ref. [36] for a more detailed discussion of such single-frequency optimization, and we emphasize below the new developments in the case of a complex frequency.
To bound the first term, sðω 0 þ iΔωÞ, in Eq. (16) (the second term is the known electrostatic constant), it is helpful to rewrite the near-field scattering amplitude not in terms of the field at the source location but instead in terms of the fields within the scatterer, at the complex frequency ω ¼ ω 0 þ iΔω. The amplitude can be written most succinctly in terms of the material susceptibility χðωÞ, which is the difference between the scatterer permittivity or permeability and that of the background: (In the SM [78], we treat the most general scenario in which the susceptibility is a 6 × 6 tensor that can be magnetic, anisotropic, nonreciprocal, and spatially inhomogeneous.) If we consider the LDOS for a single dipole orientation j and momentarily drop the j notation for simplicity, the scattering amplitude is, per Eq.
The scattered field at the dipole location, ψ s ðx 0 Þ, is given by the convolution of the free-space Green's function Γ with the polarization currents ϕ ¼ χψ in the scatterer; then, reciprocity [95] can connect the free-space Green's function to the incident field from the dipole itself (a procedure we followed at real frequencies in Ref. [36]). After defining a modified incident field,ψ inc ¼ ðE inc −H inc Þ T , the near-field scattering amplitude can be written as sðωÞ ¼ ð1=πωÞ R Vψ T inc χψ. Finding a convex constraint that encodes energy conservation-requiring absorbed power (properly normalized) to be smaller than total power expended-requires some care at complex frequencies. One cannot analytically continue the absorbed and scattered powers into the UHP, as they are not analytic everywhere (originating from their quadratic field dependence). To recover the notion of an energy-conservation constraint, we start with passivity, which states that everywhere in the UHP, the product ImðωχÞ must be positive-definite [96]: where ImðωχÞ ¼ ½ωχ − ðωχÞ † =2i, and Eq. (18) extends into the complex plane the notion that passivity implies positive(-definite) imaginary susceptibilities. From passivity, we define two positive functionals: an integral of the (positive) quantity Im½ψ † ðωχÞψ within the scatterer volume, and an integral of the (positive) quantity Im½ψ † s ðωχÞψ s outside the scatterer volume. Through repeated application of the divergence theorem and the complex-frequency Maxwell equations, we can define two new functionals, φ A ðωÞ and φ E ðωÞ. At real frequencies, the functionals equal absorption and extinction, but we define them here in the UHP: In the SM [78], we show that these functionals indeed satisfy absorption and extinction-like constraints everywhere in the UHP: φ A ðωÞ < φ E ðωÞ for Im ω > 0. This constraint is precisely the type of convex constraint needed for a bound, providing a mechanism to derive one at complex frequencies. Given the expression above for sðωÞ, and the constraint φ A < φ E , we formulate the upper bound as the solution of a (convex) optimization problem at ω ¼ ω 0 þ iΔω: maximize ψðx;ωÞ Equation (20) has a unique, globally optimal solution. The optimal field distribution ψðωÞ and scattering amplitude sðωÞ can be found through variational calculus, following a similar procedure to that developed in Ref. [36] and detailed in the SM [78]. A crucial term that emerges is a material-dependent material figure of merit, fðωÞ. For bulk (non-2D), nonmagnetic materials with scalar electric susceptibilities χðωÞ, the material FOM is given by The optimal field is proportional to this material function, as well as to the conjugate of the incident field. Now, reintroducing the average over dipole orientation j, the optimal field yields a frequency-averaged LDOS bound (SM [78]): where the complex frequency ω ¼ ω 0 þ iΔω encodes the center frequency and bandwidth of interest. Equation (22) shows that a near-field LDOS, averaged in a half-width-at-half-max bandwidth Δω around a center frequency ω 0 , is fundamentally limited by the field of a dipole in free space and by the frequency-dependent material composition of the scatterer(s). The volume integral of the incident field can be further simplified by enclosing the scatterer within some bounding shape of high symmetry over which the integral can be calculated analytically. A typical example is that of an emitter above a structured (or randomly textured [97,98]) surface, in which case the scatterer can be enclosed in a half-space that is a separation distance d from the emitter. Near-field interactions (jωjd=c ≪ 1) are dominated by the rapidly decaying evanescent fields emanating from the sources, which implies that the overall shape of the scatterer, aside from its dimensionality, has little effect on the volume integral in Eq. (22). Enclosing any structure by a half-space and keeping only the dominant near-field term in the integral in Eq. (22), we obtain a simple, shape-independent, analytical expression. (All remaining terms, which are nondivergent and typically small, are included in the SM [78].) The bound scaling is very different when the incident field is generated by a magnetic dipole rather than an electric one, and thus we can separately derive for total LDOS ρ, electric LDOS ρ E , and magnetic LDOS ρ H , the general bandwidth-averaged bounds (SM [78]): where we have defined a complex-valued wave number k ¼ ω=c, the function fðωÞ is the material FOM from Eq. (21), and for the second term on the first line, we have inserted the half-space bounds for α LDOS from Eq. (13).
Equations. (22) and (23) are foundational results of our paper. No geometrical engineering of resonances or coupling can overcome their limits. For any structure and bandwidth, the bound of Eq. (23) depends only on the frequency range of interest, the material properties at those frequencies, and the emitter-scatterer separation. Figure 4 compares the LDOS near a silver half-space to the bounds of Eq. (23). Center frequencies ranging from ω 0 ¼ 3.65 eV=ℏ to ω 0 ¼ 4.2 eV=ℏ (with corresponding wavelengths from 340 nm to 295 nm) are considered near the surface-plasmon frequency of silver. One can show analytically that in the zero-bandwidth, near-field (kd ≪ 1) limit, the electric LDOS above a nonmagnetic, surfaceplasmon-resonant interface should approach the bound of Eq. (23), while the magnetic LDOS above the same interface should approach its respective bound within a factor of 2. (This can be shown starting from asymptotic expressions in Ref. [68].) Such close approaches in the zero-bandwidth limit are visible in both Figs. 4(b) and 4(c). In the large-bandwidth limit, the bounds converge to the sum rules of Sec. II A, ensuring that they are "tight" (i.e., there is no smaller upper bound) in that regime as well. To simplify the ultrahigh-bandwidth computations (Δω ≳ 100 eV=ℏ), we use a five-pole Drude-Lorentz multi-oscillator model for silver that closely approximates tabulated susceptibility data [81] (a comparison is included in the SM [78]). There is an interesting peak in the bound at moderate bandwidths that arises due to the large, lossy permittivity of silver at about 5 eV; standard quasistatic theory [99] would predict that a half-space is nonoptimal for such a bandwidth but that perhaps another structure is optimal. A key question prompted by these bounds is whether nonplanar, designed structures-or perhaps even randomly corrugated structures-can approach the bounds at frequencies away from the surface-plasmon resonance. In Fig. 3(c), it was observed that double cones approach within almost of factor of 2 of their upper bounds [using Eq. (22) for their specific geometry] over a large range of bandwidths, further suggesting that such structural design should be possible.
The convergence of the electric-magnetic LDOS bounds to their respective sum rules, in the infinitebandwidth limit, can be verified directly from Eq. (23). The first term in the ρ E bound goes to zero as Δω → ∞ due to the e −2dΔω=c factor, in which case one can rewrite the bound as where the last term is precisely the average electric LDOS in the large-bandwidth limit [since HðΔω → ∞Þ ¼ 1=πΔω]. Similarly, the ρ H bound only contains the e −2dΔω=c factor for nonmagnetic materials, and hence the bound tends to hρ H i ≤ 0, in agreement with the sum rule. By construction, the bounds agree with their respective sum rules for large bandwidths.
The power-bandwidth limit of Eq. (22) applies equally well to 2D materials characterized by a spatial conductivity σðωÞ, with the substitution ωχðωÞ → iδ S ðxÞσðωÞ, where δ S ðxÞ is a delta function on the surface of the (not necessarily planar [83]) 2D material. In doing so, the delta function transforms the volume integral in Eq. (22) to a surface integral over the incident field. We can enclose the 2D scatterer in a high-symmetry enclosure: For a 2D plane enclosure, and keeping only the highest-order terms in the near-field limit (jωjd=c ≪ 1), we find that the LDOS above the 2D material is bounded above by (SM [78]) where for 2D materials the material FOM is (In SI units, there would be an additional factor of the freespace impedance Z 0 multiplying jσj 2 =Re σ.) There are two distinct features that emerge for 2D materials: The material FOM is jσðωÞj 2 =Re σðωÞ, and the electric-and magnetic-LDOS bounds have terms that scale as 1=d 4 and 1=d 2 , instead of 1=d 3 and 1=d in the bulk-material bounds. The different distance scaling is a natural consequence of integrating the norm of the Green's function over an area instead of a volume. Yet, there is an interesting contrast embedded within the bounds for hρi and hρ E i: Their first term, dominant over narrow bandwidths, scales as 1=d 4 , whereas the second term, dominant over wider bandwidths and corresponding to the sum rule, scales as 1=d 3 . The faster scaling with 1=d of the first term suggests a scenario in which the average LDOS over a narrow bandwidth may be larger than the sum rule would seem to allow. One possibility is that the bound is "loose" and that the 1=d 4 scaling is artificial, but in multiple previous studies [37,100] of singlefrequency behavior, 1=d 4 scaling has been observed in the LDOS near 2D materials. Another possibility is that such large response is only possible over a narrow bandwidth, though the connection of broadband response to singlecomplex-frequency response would seem to suggest that large response is likely not restricted to single frequencies.
Finally, perhaps the most likely possibility is that such a bound is achievable and simply requires negative (scattered) LDOS at frequencies outside the range of interest. There is no requirement that scattered LDOS be positive at all frequencies since a scatterer can suppress all modes and reduce the total LDOS to nearly zero. The idea of exploiting such suppression to achieve anomalously large response over some desired bandwidth is intriguing. Our derivations provide general insight into optimal structures that would reach the bounds. First, as noted above, it is critical to have as much material as possible in the near-field region of the source-that material enables the polarization currents that ultimately drive the large response. Sharp tips, though potentially exhibiting strong resonances, are not ideal. Second, the convexity-based optimization provides not only the maximal bandwidthaveraged response but also the optimal fields ψ that would generate such response. In all cases, those optimal fields, throughout the volume of the scatterer, are proportional to the incident fields. Hence, one would want to generate, perhaps via computational design [61,[101][102][103], resonances with the same phase and amplitude profile as the fields emanating from the dipolar sources. Finally, through a volume-integral-equation (VIE) framework [36], one can reinterpret our bounds as upper limits that occur when the incident field couples only to a single VIE mode at the optimal resonance location. Thus, if one can completely avoid exciting other modes, then the upper bound is in fact guaranteed to be achieved.

III. CROSS DENSITY OF STATES
The previous section developed the theoretical framework for power-bandwidth limits in the context of LDOS. We translate that framework to other near-field optical response functions, starting with CDOS [104], which measures field correlations between two points of a structured environment. In addition to fundamental interest as a correlation function, CDOS is also the critical term in the frequency integrand for resonant near-field dipolar energy transfer, as in Förster energy transfer [46,47,105], as well as for quantum entanglement and super-radiative coupling between qubits [106][107][108][109][110]. Whereas LDOS is given by the Green's function for identical source and measurement points, CDOS is given by the Green's function between different source and measurement points, where i and j are the measurement and source polarizations, respectively, and we again use the "s" subscript to denote the scattered-field contribution (subtracting off the known free-space contribution). By analogy to Eq. (4), there is a sum-rule constant for CDOS, which we denote α CDOS (defined as 1 2 Re½Γ s;ij ðx; x 0 Þj ω¼0 ); by the similarity to LDOS, one can follow a similar procedure to derive power-bandwidth limits. The key new feature is that instead of the single separation distance between the emitter and scatterer controlling the bound, now there are two relevant separation distances: the distance between the emitter and the scatterer, denoted d 1 , and the distance between the scatterer and the measurement point, d 2 .
To bound Eq. (26) averaged over any bandwidth, we first rewrite the Green's function in terms of the polarization currents of the scatterer(s). This leads to an overlap integral between the polarization currents induced by the field incident from the source position with a (parity-reversed) secondary field incident from the measurement position. Ideally, the polarization response is maximally aligned to both fields (SM [78]), in which case the maximal response is proportional to the square root of the energy of each "incident" field, ψ inc;1 and ψ inc;2 . Applying the bandaveraged bound approach modified as described above, we arrive at the following bound on average CDOS for bulk, nonmagnetic materials (SM [78]): where the complex frequency ω ¼ ω 0 þ iΔω encodes the center frequency and bandwidth of interest. In the nearfield regime, (jωjd=c ≪ 1), Eq. (27) further simplifies (SM [78], neglecting α CDOS for small to moderate bandwidths): where we separate the electric-and magnetic-source contributions to the CDOS. Now, the electric bounds depend on the source-scatterer and measurement-scatterer separation distances to the three-halves power, instead of the cubic dependence for LDOS when the source and measurement points are identical. Note that the material dependence of the bound is encoded in the same material figure of merit, fðωÞ, as defined in Eq. (21), suggesting the universal role it may play in determining the maximal broadband response of any material.

IV. RADIATIVE HEAT TRANSFER
Near-field radiative heat transfer (NFRHT) can be substantially larger than far-field radiative heat transfer, via photon-tunneling evanescent-wave energy transfer, and has generated much interest for applications such as thermophotovoltaics [11,[111][112][113][114]. In RHT, there are two bodies at temperatures T 1 and T 2 , with minimal separation distance d. The net radiative heat flux between the two bodies [115] is given by H 1→2 ðωÞ ¼ ΦðωÞ½Θðω; T 1 Þ − Θðω; T 2 Þ, where Θðω; TÞ denotes the mean energy of the harmonic oscillator (without the zero-point energy ℏω=2) and Φ is a temperature-independent flux rate from incoherent sources in body 1 radiating to body 2.
Since Θ is positive for all frequencies, we can bound the difference ΘðT 1 Þ − ΘðT 2 Þ by its maximum value, which is simply ΘðT 1 Þ (taking T 1 > T 2 ), i.e., (Φ is non-negative at all frequencies for passive media): where T ¼ T 1 and the equality holds if body 2 is at absolute zero.
In order to apply contour-integration techniques (described in Sec. II C) to bound Eq. (30), one would need to extend the mean energy Θ to negative frequencies, and Θ has to decay fast enough such that the integral in Eq. (30) is finite for all real frequencies. While it is bounded and decays for large positive frequencies, it diverges when extended to negative frequencies. To avoid such asymmetry, one could add the vacuum energy ℏω=2 and work with Θ v ¼ 1 2 ℏω cothðℏω=2k B TÞ [115] (k B denotes the Boltzmann constant), which is symmetric about the origin. However, Θ v diverges linearly for large frequencies and would enforce dramatic restrictions on the flux rate Φ to fall rapidly at high frequencies. Even if convergence were not an issue, Θ v contains infinitely many singularities along the imaginary axis. This would reduce a contour integral evaluation of Eq. (30) (albeit with Θ v instead of Θ) to a sum of infinitely many residues (at the "Matsubara frequencies" [116]), which is cumbersome to handle. We can avoid all of these issues (asymmetry, nonconvergence, and the Matsubara sum) in a single stroke by virtue of the following fact: For any temperature T, the mean-energy spectrum Θðω; TÞ is simultaneously bounded above and closely approximated ( [78]). [Note that unlike the spectrum of a blackbody, whose peak wavelength is nonzero and scales inversely with temperature, the mean-energy spectrum ΘðωÞ peaks at zero frequency.] Since the integral in Eq. (31) is the Lorentzian-averaged flux rate, it might be tempting to close the contour in the UHP to relate the integral to a sum of residues, in the spirit of Sec. II C. However, unlike LDOS, the flux rate Φ is not directly given by the real or imaginary part of a scattering amplitude. But we can transform the problem by generalized reciprocity [117], recasting the flux rate from a surface integral of fields generated by volume sources to a volume integral of fields generated by sources along a surface, revealing a surprising similarity to LDOS and ultimately leading to NFRHT bounds.
The heat flux between the two bodies is given by the power flow through a separating surface S: Þ is a real symmetric matrix. The fields can be expressed [118] as convolutions of the system Green's function Γðx; x 0 Þ with the thermal sources ϕðx 0 Þ: ψðxÞ ¼ R V Γðx; x 0 Þϕðx 0 Þdx 0 . The incoherence of the thermal sources can be incorporated via the fluctuationdissipation theorem in a standard substitution [115]. The key step is that we then use generalized reciprocity [117] (cf. SM [78]) to interchange the source positions in body 1 with the measurement position on the surface S separating the bodies. This transforms the problem to sources over a surface in free space (or a homogeneous background) radiating into one of the bodies, and the net flux rate is an energylike quadratic form evaluated in the body, for a specific linear combination of electric and magnetic dipoles (cf. SM [78]). In the ideal scenario, the dipoles are only emitted into body 1; in any case, the radiation into body 1 is bounded above by the total radiation. In the near field, the total radiation is essentially exactly equal to the scattered-field LDOS, such that we can directly bound (via Cauchy-Schwarz arguments) the flux rate at any frequency by the product of the electric and magnetic LDOS (SM [78]): To bound the bandwidthaveraged flux rate hΦi, the product of the electric and magnetic LDOS prevents direct identification of a complexanalytic quantity, but we can again use Cauchy-Schwarz for a bound in terms of the individually bandwidth-averaged ρ E and ρ H (ideally, they would exhibit the same line shape, in which case the Cauchy-Schwarz substitution would be an equality), such that we can write (SM [78]) In the case of near-field RHT, the center frequency and bandwidth are fixed by the temperature as discussed above. Using Eq. (33), we can bound NFRHT in Eq. (31): Equation (34) is the culmination of our transformations: The product of flux and oscillator energy in the integrand of Eq. (30), which is not easily extensible to negative frequencies nor analytic over the upper-half plane, is bounded above by (and, for optimal designs, equal to) an integral over a bounding surface between the bodies of the geometric mean of the bandwidth-averaged electric and magnetic LDOS. Then, for bulk, nonmagnetic materials, we can directly insert the LDOS bounds, Eq. (23) (with the second term slightly modified to allow for a two half-spaces enclosure), into Eq. (34) to obtain a near-field bound on net radiative power transfer (with suitable generalizations for magnetic and/or 2D materials). The material figure of merit of Eq. (21) again plays a key role; in the case of near-field RHT, the temperature determines the bandwidth. Moreover, because the center frequency is zero, the material figure of merit is evaluated at the purely imaginary frequency i ffiffi ffi 2 p k B T=ℏ; since susceptibilities are real and positive for imaginary frequencies in the UHP [79], the material FOM (now written as a function of temperature) can be greatly simplified: If we denote r 1 and r 2 as the distances from a surface point to bodies 1 and 2, respectively, and f 1 ðTÞ and f 2 ðTÞ as the corresponding material figures of merit, then the bound for an arbitrary separating surface S is where the constant terms proportional to 4.21 arise from the LDOS constant α LDOS for two half-spaces (SM [78]). For a planar bounding surface halfway between the two bodies at d=2 for a minimal separation d, the integral in Eq. (37) can be done analytically, simplifying the bound to maximum heat transfer per unit surface area A: where we have dropped constant terms (which are small relative to f 1 and f 2 for all practical materials and temperatures of interest). Equation (38) represents the first general bound to near-field radiative heat transfer. We can more easily interpret Eq. (38) by recasting the expression as a product of dimensionless enhancement factors with the far-field blackbody limit, H BB ¼ σT 4 , where σ is the Stefan-Boltzmann constant. Using the thermal de Broglie wavelength, λ T ¼ π 2=3 ℏc=k B T, algebraic manipulations yield an equivalent alternative, where β ¼ 120=π 16=3 ≈ 0.268 is of order 1. Equation (39) succinctly identifies two maximum possible enhancements beyond the blackbody limit. First, a distance-dependent enhancement ðλ T =dÞ 2 emerges, which accounts for the increased amplitudes of evanescent waves at shorter separations. This enhancement factor is intermediate between that appearing in the bounds for electric and magnetic LDOS (1=d 3 and 1=d, respectively), as RHT is equivalent to a combination of electric and magnetic dipolar radiation in free space. The thermal wavelength λ T (which is about 10 μm at room temperature) sets the threshold for the nearfield regime, highly sensible since the blackbody radiation limit holds only when the length scales involved are greater than λ T . The second enhancement factor is a materialdependent factor f 1 ðTÞ þ f 2 ðTÞ, which accounts for material-based resonant enhancements. Per the material FOM of Eq. (35), the larger the susceptibility is at the complex frequency set by the temperature of the emitter, the larger the possible response is. The bounds of Eqs. (38), and (39) cannot be overcome by any metamaterial, metasurface, or other design approaches.

V. OPTIMAL MATERIALS
Embedded throughout the LDOS, CDOS, and near-field RHT bounds is a material metric fðωÞ, defined in Eqs. (21) and (25), that indicates the intrinsic capability of any material to exhibit large optical response over a frequency bandwidth Δω around a center frequency ω 0 . This material metric enables comparison of any material-dielectric and metal, 2D and 3D, lossless and lossy-many of whose capabilities cannot be understood through single-frequency bounds or sum rules.
Sum rules, such as Eq. (5), typically have little-to-no dependence on material parameters, suggesting that different materials only alter resonant bandwidths, without impacting total optical response. Yet, this is misleading on two fronts: (1) It only applies over infinite bandwidth; over any finite bandwidth, material properties play an important role in maximal response, and (2) sum rules require susceptibilities that satisfy Kramers-Kronig relations, diminishing to zero at high frequencies. The decay-tozero requirement, though physically reasonable, means that even "dielectric" media (semiconductors, insulators, etc.) have a plasmalike response at large enough frequencies.
Such response contributes to sum rules, often in a large way due to the negative susceptibility. This contribution obscures the behavior of, for example, a transparent dielectric at optical frequencies, by accounting for transitions that occur at UV and x-ray frequencies. Thus, sum rules miss finitebandwidth effects and dramatically overestimate dielectricmaterial interactions. At the other end of the continuum, single-frequency bounds [24,36,37,42,43] apply at any given frequency, but they use material loss as the intrinsic system limitation and thereby diverge for materials with vanishingly small imaginary susceptibilities (such as dielectrics). The material FOM embedded in the power-bandwidth limits does not have any of these limitations: It can account for finite bandwidths, it does not require susceptibilities that asymptotically approach zero at large frequencies, and it provides a finite bound for lossless materials for any nonzero bandwidth. Hence, fðωÞ is a simple expression that enables comparison among the multitude of possible optical materials.
To gain intuition about the material FOM, we consider the small-bandwidth limit in which Δω ≪ ω 0 . We delineate two types of (bulk, 3D) materials: lossy materials, with a nonzero Im χðωÞ in the small-bandwidth limit, and lossless materials, with Im χðωÞ ≈ 0 [and Im χðωÞ ≪ Δω=ω 0 even in the small-bandwidth limit], as would characterize many transparent materials at optical frequencies. In the small-bandwidth limit, the material FOM is approximately given by where we have retained the full material FOM for 2D materials since it is already simple. For high-index lossless materials, the expression would simplify even further: fðωÞ ≈ ω 0 Δω χðωÞ lossless; high-index: For small to moderate bandwidths, a natural dichotomy emerges: Lossy materials are inherently restricted by material loss in Im χðωÞ, whereas lossless materials are inherently limited by the relative bandwidth Δω=ω 0 . Intuitively, in simple single-mode interactions, one could interpret the figures of merit as dictating that lossy materials have maximum responses proportional to jχðωÞj, over a bandwidth proportional to jχðωÞj=Im χðωÞ, whereas lossless materials have maximum responses proportional to χðωÞ, over bandwidths proportional to ω 0 =Δω. The intuitive interpretation about lossy-material bandwidths is supported by previous results in quasistatic plasmonic frameworks [119,120]. Of course, we make no assumption of singlemode or quasistatic behavior, and our scattering framework is valid for any number of resonances as well as more complex phenomena such as Fano interactions [121] and exceptional points [122,123]. And perhaps more importantly, it enables consideration of lossless and lossy media on equal footing. As discussed in the Introduction, the maximum response of lossless media has been impossible to  accurately capture with either the sum-rule or the singlefrequency-bound approaches known today. In the complexfrequency approach, bandwidth naturally adds a form of "loss" to the system, yielding finite bounds that vary smoothly with bandwidth. Figure 5 compares the material FOM for a large variety of materials at optical frequencies. To evaluate the material susceptibilities and conductivities at complex frequencies, we use analytic models (e.g., Lorentz-Drude oscillators) that can be continued into the complex plane, and we ensure that they are accurate over the range of bandwidths considered. On the left side of the figure, we model the material FOM for canonical material types: (a) a Drude metal, χðωÞ ¼ −ω 2 p =ðω 2 þ iγωÞ, for plasma frequency ω p and loss rate γ, (d) a lossless, constant-susceptibility [χðωÞ ¼ 9] material, and (g) a Drude 2D material, with conductivity σðωÞ ¼ iω p =ðω þ iγÞ. One can see that these three material types show very different characteristic dependencies of their FOM on frequency and bandwidth. The Drude-metal FOM is nearly independent of small-tomoderate bandwidths, as expected from Eq. (40)-for metals, intrinsic loss is the limiting factor. The FOM of a Drude metal increases with the center wavelength (of the frequency band of interest) λ 0 since the increasing wavelength increases the magnitude of the susceptibility. By contrast, a constant-permittivity "dielectric" has nearly opposite dependencies. The figure of merit is independent of center wavelength and highly dependent on the bandwidth. Because the bandwidth is the source of loss, there is a trade-off between the average response and bandwidth. Finally, 2D Drude conductors are somewhere in between.
Loss originates from both the material parameter γ and the bandwidth, with increasing FOM towards the lower-righthand corner of Fig. 5(g): small bandwidth and large wavelength (for a large conductivity). These simplified metal, dielectric, and 2D conductor profiles capture well the key dependencies of the FOM for real materials: The plots in Figs. 5(b), 5(c), and 5(h) follow the same trends as those in Figs. 5(a), 5(d), and 5(g): Metal [81] FOM increases with wavelength, whereas dielectrics (Si [124,125] and SiC [126,127]) and polaritonic materials (SiO2 [128][129][130] and TiO2 [131,132]) that support surface phonon-polaritons at mid-IR frequencies [84] do not depend appreciably on wavelength. Conversely, the plots in Figs. 5(e), 5(f), 5(i) show the effects of increasing bandwidth, with metal material FOMs nearly unchanged but those of the dielectrics and polaritonic materials decreasing nearly linearly. The material FOM of 2D conductors increases with both wavelength and smaller bandwidths. We consider the 2D conductivities of graphene for various Fermi levels [82], magnetic biasing [133], and AA-type bilayer stacking (BLG) [134], hBN [135], and metals Ag, Al, and Au, with conductivities set by a combination [136] of bulk properties and interlayer atomic spacing.
An intriguing prediction that emerges from the LDOS and CDOS power-bandwidth limits is that the 2D-material bounds increase more rapidly for smaller separations (∼1=d 4 ) than for bulk materials (∼1=d 3 ), suggesting that 2D materials should overtake bulk materials as optimal, with the precise transition depending on the bound prefactors and, crucially, the relative 2D and bulk material figures of merit. In Fig. 6, we consider Drude models for both a 2D conductivity [σ ¼ iω p =ðω þ iγÞ] and a bulkmaterial susceptibility [χ ¼ −ω 2 p =ðω 2 þ iγωÞ] and plot isocontours for the material FOMs of each in panel (a). In Fig. 6(b), we trace out the region of frequency and bandwidth for which the bulk 3D material has a larger maximal response and the region for which the 2D material offers larger maximal response. This line will be different for every 2D/bulk-material pair, and it is determined by Eqs. (23) and (24) and their CDOS analogs.
The above discussion about the material FOM as a function of bandwidth is particularly relevant for LDOS and CDOS, where different source configurations lead to different center frequencies and bandwidths. For near-field RHT, as shown in Sec. IV, the center frequency is fixed (at zero), the bandwidth is uniquely determined by the temperature, and the material FOM takes a particularly simple form, fðTÞ ¼ χði ffiffi ffi 2 p k B T=ℏÞ. At practical temperatures of interest, the material FOM is determined by the susceptibility at an infrared frequency, except along the imaginary frequency axis instead of the real line. Dielectrics have near-constant susceptibilities over this range; surprisingly, polar dielectrics, which have resonant transitions in the infrared, corresponding to surface-phonon-polariton modes, have no resonant peaks along the imaginary axis, and exhibit dielectriclike permittivities and material FOMs. For plasmonic materials, we can use a Drude model to describe the material FOM since k B T=ℏ is typically an infrared frequency. For a Drude model with plasma frequency ω p and loss rate γ, the material FOM is Materials with large plasma frequency and small γ are ideal for maximum heat transfer. Figure 7 shows the material FOM for common plasmonic materials, all exhibiting a decrease with temperature. Although this may appear counterintuitive, it does not imply that maximal NFRHT decreases with temperature, as Eq. (39) has a separate T 4 multiplicative factor. Instead, such behavior reflects the fact that increasing temperature has the same effect as increasing material loss due to the larger bandwidth over which NFRHT occurs. We can understand this through two alternative vantage points: The larger bandwidth increases the imaginary part of the complex frequency at which the equivalent scattering problem occurs, and moving higher into the UHP increases loss; alternatively, in Sec. II A, we showed that LDOS quantities are subject to sum rules, and thus increasing bandwidth necessarily reduces the average response. For nondispersive dielectrics, the material FOM does not depend on temperature, and it takes on small values relative to metals (which usually have large susceptibilities at small frequencies), as shown in the inset of Fig. 7.
We can derive a particularly simple form of the NFRHT bound, Eq. (39), for plasmonic materials with f given by Eq. (42). For many such materials, at practical temperatures of interest, the material loss rate γ for many materials is much larger than γ ≫ k B T=ℏ. In the near field, the separation distance d is much smaller than λ T , such that the exponential factor in Eq. (39) is approximately 1. Assuming that both bodies consist of the same material, we can plug Eq. (42) into Eq. (39) to arrive at the following bound: Thus, we see that there are three enhancements relative to a blackbody: the near-field distance enhancement, an enhancement from the ratio of the plasma frequency to the loss rate, and an enhancement from the ratio of the plasma frequency to an effective thermal frequency. The material FOM extends in a natural way to anisotropic, magnetic, and even spatially inhomogeneous media, as shown in the SM [78]. Nonlocality, wherein the polarization FIG. 7. Comparison of material FOM fðTÞ in the context of NFRHT for a variety of conventional metals. The left inset shows the orders-of-magnitude difference between silver and SiC, a polar dielectric that supports surface-phonon polaritons at infrared frequencies yet has a small material FOM due to its non-Drude-like permittivity. The right inset shows the mean energy spectrum Θ (normalized to its maximum value) at different temperatures. The decaying FOM reflects increasing loss from bandwidth, which broadens with temperature. For a practical range of temperatures as shown here, metals are superior under this metric compared to dielectrics. Generally, fðTÞ favors materials with large plasma frequency ω p and small material loss rate γ. field at a position x depends on the electromagnetic field at another x 0 , can also be incorporated for certain hydrodynamic models [137][138][139]. An intriguing open question is whether density functional theory models can be bounded in a similar fashion. Such bounds could motivate and clarify the search for new "quantum materials."

VI. EXTENSIONS AND SUMMARY
We have established a framework for identifying upper bounds to near-field optical response over any frequency bandwidth of interest, with an emergent material FOM that enables quantitative comparisons of any material. We derived bounds for three optical-response functions: the LDOS, a measure of the spontaneous-emission enhancement for any electric or magnetic dipole (or atomic dipolar transition), CDOS, a field-correlation function, and radiative heat transfer, a measure of energy transfer from thermal fluctuations. The property of these response functions that is critical to our framework is the fact that they can be related to the imaginary part of a function that is analytic in the upper half of the complex-frequency plane. Here, we explore how our complex-analytic framework can be extended to other optical response functions.
There are a few near-field quantities that map closely to LDOS. First, atomic Lamb shifts due to inhomogeneous electromagnetic environments [44,45] are given by frequency integrals of Im Γ ij ðx; x; ωÞ; ð44Þ multiplied by frequency-dependent prefactors that include the atomic frequencies and position matrix elements. Hence, the sum rules and power-bandwidth limits derived here can be directly extended to the emitter-environment coupling rate in the Lamb shift. Second, Raman scattering [140] is a process in which a pump wave interacts with a molecular transition and subsequent emission that is potentially enhanced by the electromagnetic environment. It appears possible to bound this interaction above by the product of the LDOS at the separate pump and emission frequencies, in which case the framework herein can be applied for sum rules and bounds. We discuss the derivation and bounds to this process in a separate publication [141]. A more complex case is that of free-electron radiation (e.g., Smith-Purcell, Cherenkov, etc.), in which a free electron at high speed (of order c) interacts with a structured medium to generate electromagnetic radiation. The incident electromagnetic field of a free electron is proportional to a modified Bessel function. One difficulty that arises in considering sum rules and power-bandwidth limits in this case is that the modified Bessel function has a logarithmic frequency dependence at the origin, rendering it difficult to apply standard contour-integral techniques as we have done here. In two dimensions, a constant-velocity free electron emits (evanescent) plane waves, and sum rules and bounds appear to emerge in a straightforward way. The three-dimensional case may be more difficult, however.
Another complication emerges when the dipolar sources are embedded within the scatterers of interest. This is typical of the Casimir force, which is a momentum-transport quantity arising from vacuum-induced fluctuations. It appears possible to derive sum rules for such a quantity by exploiting the same generalized reciprocity [117] that we used for near-field radiative heat transfer. We will consider sum rules, powerbandwidth limits, and interesting physical consequences for Casimir physics in an upcoming publication [142].
Finally, we consider extensions of this framework to cases when the incident field is not generated by a localized dipolar or free-electron source but instead by a plane wave. At first glance, it would appear that the conventional optical theorem [33,[62][63][64] provides a simple analytic quantity to serve as the basis for the contour-integral and energyconservation approaches developed here. Yet, the bounds derived by such an approach yield a term that grows exponentially with bandwidth [the opposite of the exponential decay seen in Eq. (22)]. Such a dependence is not physical-known sum rules [25] would contradict it-and is instead an indication that the energy constraints developed herein, based on the positive-definite quantities in Eq. (19), may not be optimal for plane-wave sources. Modified constraints may be required to develop powerbandwidth limits for plane-wave scattering.
The bounds derived here, and those suggested above, suggest a tremendous design opportunity in near-field nanophotonics. For various canonical structures, there are frequency ranges at which they come close to reaching the bounds, but there are also wide frequency ranges at which there is a sizable gap. New designs, and new approaches such as large-scale computational "inverse design" [61,[101][102][103], offer the prospect for overcoming the gap and revealing the physical principles underlying optimal operation.

I. SUM RULES AND POWER-BANDWIDTH LIMITS FOR NON-VACUUM BACKGROUND
If the emitter and scatterer are embedded in a background material with a scalar, non-dispersive permittivity ε b and permeability µ b , LDOS is modified as follows: where ν b is a diagonal 6 × 6 matrix: Then, to update the sum rules and power-bandwidth limits throughout the paper, ν b (or ε b and µ b individually) is directly inserted into Eqs. (3,5,9,17,21) and the expression for α LDOS after Eq. (4) of the main text.

II. HIGH-FREQUENCY SCALING BEHAVIOR OF LDOS
This section analyzes the high-frequency asymptotic behavior of the frequency-resolved LDOS, to show that its allfrequency integral converges to a finite value. In order to do so, we need to rewrite LDOS in terms of the polarization currents induced in the scatterer V by a point dipole source in its vicinity. By encoding the material properties of the scatterer in the induced polarization P instead of the scattered field E s , we can better understand the high-frequency scaling behavior of LDOS and hence the validity of its sum rule as will be demonstrated here. Starting from Eq. (2) of the main text, LDOS can be expressed in terms of the free-space Green's function G EE , omitting the subscript in this section for clarity (we only consider electric LDOS in the presence of a nonmagnetic material for simplicity, but the high-frequency asymptotics are unchanged even if we include magnetic LDOS and/or consider magnetic materials): where G(x, x 0 )p j (x 0 ) is the electric field incident on x from an electric dipole oriented along p j at x 0 . In the derivation, we have also taken the transpose of the integrand (which does not affect ρ, a scalar) and then used the relation G T (x 0 , x) = G(x, x 0 ), which follows from reciprocity [1]. Note that Eq. (S3) represents the electric component of the expression in Eq. (20) in the main text, and the steps leading to Eq. (S3) show explicitly how the expression in Eq. (20) is derived as well. In the high-frequency limit, the polarization field must decay towards zero (the bound charges cannot respond to such high frequencies), and on physical grounds [2] the decay must occur in proportion to 1/ω 2 . Conventionally, the decay constant is chosen to be a "plasma frequency" ω p that is physically meaningful for metals but applies to dielectrics as well. Because the scatterer becomes transparent at high frequencies, the Born approximation applies and the polarization field will be directly proportional to the incident field: Inserting Eq. (S4) into Eq. (S3) results in a term of the form Im´V E 2 inc , which can be solved using the free-space Green's function G by noticing that the conventional definition for LDOS averages over all dipole orientations: 1/ω FIG. S1: High-frequency behavior of electric LDOS for a Drude metal with loss rate γ = 0.1ω p and distance d = 0.1c/ω from halfspace. For frequencies much larger than the plasma frequency ω p , it falls off as 1/ω (in an oscillatory manner).
The incident field is given by Green's function multiplied by the dipole polarization (with unit magnitude). The Green's function is [3]: where a = kr, k = ω/c. As a side remark, s(ω) (the complex extension of LDOS) decays exponentially in the upper half complex-ω plane due to the term e ikr , leaving us to only worry about real frequencies to confirm the validity of LDOS sum rule. In the limit as ω → ∞, we can keep only the highest order term in ω for G 2 ij , resulting in the following asymptotic expansion (dropping uninteresting numerical prefactors -we only care about the scaling behavior of ω): An interesting remark about the sin(2ωr/c) term is that it arises from taking the square of individual components of G instead of its Frobenius norm [4], which would get rid of the oscillatory term e ikr and hence sin(2ωr/c) in Eq. (S7). This is closely related to the fact that we take our dipole source to be real (p = p), as it allows us to rewrite LDOS as a volume-integral expression in terms of E inc in Eq. (S3). For various computations, it appears that´V sin(2ωr/c) tends to scale as sin(ω)/ω 2 , rendering ρ E (ω → ∞) ∼ sin(ω)/ω. While the high-frequency scaling behavior of LDOS depends on the scatterer geometry, we can still investigate under what conditions ρ(ω → ∞) falls off "sufficiently rapidly" for´∞ 0 ρ(ω) dω to be finite. Recalling that´∞ 0 sin(x)/x p dx is convergent for p > 0, our requirement for a valid sum rule is as follows (assuming an oscillatory term): What Eq. (S8) says is that´V sin(2ωr/c) should fall off as sin(ω)/ω or faster for LDOS to decay at large enough frequencies. As a prototypical example, Fig. S1 shows that electric LDOS falls off as 1/ω for a halfspace. In fact, our sum rule for LDOS appears to hold for all geometries (subject to the limiting procedure for singular geometries discussed in the main text).

III. SUM RULE DERIVATION
In this section, we lay out in detail how the sum rule in Eq (4) can be derived starting from Eq (2) and Eq (3) of the main text. Eq (3) defines s(ω) as a complex-valued function that is analytic in the upper half of the complex-ω plane. Our goal is to construct a closed contour C that includes the real line and is analytic inside C, so that we can apply Cauchy's integral theorem to conclude that the integral over C vanishes. However, s(ω) has a simple pole at ω = 0, which forces us to avoid the origin and instead integrate over the infinitesimally small upper semicircle ("bump" at the origin as shown in Fig. 2(a) of the main text). Since the scattered fields decay exponentially at positive imaginary frequencies, the integral over the infinite semicircle vanishes. This leaves us with the integral over the real line (except at the singularity at ω = 0) and that over the "bump." The latter can be straightforwardly evaluated in polar coordinates, resulting in the following (principal value of the) integral of s(ω) over all real frequencies (assuming real-valued dipole amplitudes as in the main text): The zero-frequency limit follows from the fact that our "bump" can be chosen arbitrarily small in radius. Recalling from Eq. (2) that ρ(ω) = Im s(ω) on the real axis and that ρ(ω) = ρ(−ω), we obtain the sum rule for LDOS after taking the imaginary part on both sides of Eq. (S9): Equation (S10) is our sum rule for LDOS, identical to Eq. (4) in the main text since j ξ T j (x 0 )ψ s,j (x 0 ) = Tr Γ s (x 0 , x 0 , ω) for summation over j = {x, y, z}.
In the main text, we have assumed a scalar permittivity ε as it reduces α LDOS = 1 2 Re j [ξ T j (x 0 )ψ s,j (x 0 )]| ω=0 to simple analytical form for certain geometries. However, we did not make such an assumption in deriving Eq. (S10), and so the definition of α LDOS is valid for tensor permittivities as well.

IV. POSITIVITY OF αLDOS
While not immediately obvious, the sum-rule constant α LDOS is a positive quantity for arbitrary geometries and materials. For simplicity, we only consider electric fields-the result presented here carries over to magnetic LDOS after appropriate replacements (ε → µ, E → H). To prove the positivity of α LDOS , it is useful to express it in terms of the polarization currents induced, and hence the field, in the scatterer V (following the procedure in Eq. (S3)): In what follows, the fields are all evaluated in the electrostatic regime. Expressing E as a sum of the incident and scattered components E inc and E s respectively, where V + denotes the entire volume outside the scatterer (where the susceptibility is zero): In going from the first to the second line, we have used the following relation, which can be derived by expressing the integral over V + ,´V + E T s εE s , in terms of the fields inside the scatterer V through repeated use of the divergence theorem:ˆV Rewriting Eq. (S13) in a more suggestive form, > 0 (S17) since χ and ε are assumed to be hermitian and positive-definite (they also commute, since χε = (ε − ε 0 I) ε = ε (ε − ε 0 I) = εχ and similarly for their inverses and squares). In going from the first to the second line, we used the fact that χ − χ 2 ε −1 = χ I − (ε − ε 0 I) ε −1 = ε 0 χε −1 . Equation (S17) thus shows that α LDOS > 0. An intuitive way of seeing that α LDOS is positive is to consider an infinitesimally small volume of a scatterer. Since we know that α LDOS is by definition zero in free space, the presence of an arbitrarily small scatterer should also have α LDOS arbitrarily close to zero and positive (and hence positive for any nonzero volume). Otherwise, the presence of a scatterer with negative α LDOS would imply from our monotonicity theorem that further shrinking its volume (until it becomes negligible in spatial extent) will cause α LDOS to become more and more negative-in contradiction with the fact that α LDOS = 0 in the absence of a scatterer.

V. DERIVATION OF αLDOS FOR CANONICAL GEOMETRIES
If the scattered fields are known at zero frequency, we can use the expression from Eq. (S10) to directly evaluate α LDOS : where the dipole source is at x 0 and ψ s,j , here evaluated at zero frequency, denotes the scattered field due to ξ j , a unit dipole polarized alongĵ. In what follows, we use the image charge method to exactly determine the electro/magnetostatic scattered field for canonical geometries. Except for a halfspace, these geometries do not admit closed-form expressions for α LDOS for finite ε, and so for non-halfspace geometries we assume conductive materials (that are perfect conductors at zero frequency). For simplicity, we focus on the electric LDOS in this section. Duality simplifies the corresponding computations for magnetic LDOS.

A. Halfspace
Consider a plane interface separating the two semi-infinite half spaces-vacuum with permittivity ε 0 = 1 at z > 0 and a dielectric with (electrostatic) permittivity ε at z < 0. In the electrostatic regime, a dipole p j a distance d above the interface will create an image dipole with some orientation at a distance d below the interface. Since the LDOS due to randomly oriented dipoles equals that for dipoles each along one of the x, y, z-axis, where p j is the image dipole due to p j andr the unit vector pointing from p j to p j (r =ẑ here). The effect of the dielectric medium with permittivity ε at z < 0 is to introduce the prefactor ε−1 ε+1 (which simplifies to 1 for perfect conductors, as expected) in the image charges [5] in order to satisfy the requirement of continuous tangential fields at z = 0. From p x = − ε−1 ε+1 p x , p y = − ε−1 ε+1 p y , p z = ε−1 ε+1 p z , we obtain the final result in Eq. (S19). Adding the corresponding term for magnetic LDOS, we arrive at α LDOS for a halfspace in Eq. (5) of the main text.

B. Sphere
A dipole p j is placed at a distance h from the center of a perfectly conducting sphere of radius R such that it is outside of the sphere (h > R). For simplicity, assume that the dipole lies along the z-axis with the origin coinciding with the center of the sphere. For vanishing tangential fields on the spherical surface, the image dipole p j will be at a distance R 2 h from the center of sphere (also along z-axis) such that p Denoting the distance between the source and image dipoles as (S20) Since h = R + d where d is the distance from the dipole to the nearest surface of the sphere, we arrive at the following result: While we assumed the dipole to be on the z-axis, Eq. (S21) holds for any dipole a distance d away from the surface (as long as we sum over random orientations) from spherical symmetry.

C. Two halfspaces
The derivation of α LDOS for two perfectly conducting halfspaces closely mimics the halfspace case, but now there are two boundary conditions to safisfy-the tangential fields must vanish at both interfaces. Taking the two halfspaces to be at z < 0 and z > d (and vacuum with permittivity ε 0 in between the two halfspaces), consider a dipole p j placed in between the two halfspaces at z = d 2 . This will create image dipoles on each halfspace a distance d from the dipole source. However, the image dipole in one halfspace will cause a non-vanishing tangential field at the interface of the other halfspace, so that these image charges leads to another set of image charges (now a distance 2d from the dipole source on both sides). Iterating this procedure, we obtain the following expressions for the parallel (z-axis) and perpendicular (x, y-axis) components: where the alternating/constant sign for the perpendicular fields follow from the fact that the image dipole parallel/orthogonal to the interfaces is opposite/same in sign to the dipole that created it (which itself is an image dipole until one arrives at the original dipole source). As expected on physical grounds, the infinite summations

VI. GENERALIZATION OF MONOTONICITY THEOREM FOR ANISOTROPIC MATERIALS
In deriving the monotonicity theorem in Sec. IB of the main text, we have taken the permittivity and permeability to be isotropic. In this section, we generalize Eq. (8) of the main text for the case of anisotropic media, following the prescription given in [7]. For simplicity, we consider electric fields only-extension to the magnetic case is straightforward after appropriate replacements (ε → µ, E → H, D → B). Working in a local coordinate system for each point on the geometrical boundary dividing regions of permittivity ε 1 (scatterer) from ε 0 (background), we define F = (D 1 , E 2 , E 3 ) such that F 1 = D ⊥ and F 2,3 = E . We can then write ∆α LDOS under outward deformation of the scatterer:ˆ∆ where τ is a 3 × 3 matrix given by: For scalar, isotropic permittivities ε 1 and ε 0 , Eq. (S24) reduces to Eq. (8) of the main text. Thus, a monotonicity theorem holds in the tensor-material case if the matrix τ (ε 1 ) − τ (ε 0 ) is positive-definite.

VII. MONOTONICITY THEOREM FOR CONDUCTIVE MATERIALS
The monotonicity theorem derived in the main text applies for any finite ε, and the fact that it works for arbitrarily large ε suggests that such a theorem may apply in the perfect-conductor (at zero frequency case). In this section, we provide an alternative derivation of the monotonicity theorem for the perfect-conductor case. From Sec. IB of the main text, we wish to show that ∆α LDOS = [ 1 2 ξ(x 0 ) · ∆ψ s (x 0 )]| ω=0 > 0. Working at zero frequency throughout this section, imagine a small change in the geometry of the scatterer such that the LDOS dipole source induces a polarization current P ind in a volume δV = δx n dA small compared to V , where δx n indicates the size of deformation normal to the scatterer boundary. Under this condition, P ind = σn = ε 0 E wheren denotes the normal to the conducting surface and σ the surface charge density. This follows from the fact that for electric conductors, the electric field is normal to the conducting surface such that E = σ ε0n . Likewise, M ind = µ 0 H for magnetic conductors. Based on these relations, we arrive at the following expression (absorbing ε 0 and µ 0 in the definition of G, which represents the Green's function in the presence of scatterer V ): Using Eq. (S26) and taking transpose on both sides (to exploit reciprocity [1]), we can derive δF as follows: Equation (S27) proves our monotonicity theorem for perfect conductors. Notice that the derivation presented here closely mirrors that of Eq. (S3) in that we take the transpose of the integrand, use reciprocity (at zero frequency, the "off-diagonal" 3 × 3 matrices G EH and G HE vanish and we do not have to worry about sign issues associated with them in taking the transpose of Γ), and the relation ψ(x) = Γ(x, x 0 )ξ(x 0 ) (but now Γ is not the free-space Green's function, but one in the presence of scatterer V before shape deformation).

VIII. CONVEX CONSTRAINTS IN THE UHP
In this section, motivate our definition of ϕ A and ϕ E in Eq. (19) of the main text, the complex extensions of the non-radiative and total power respectively. To this end, we start from the complex-frequency Maxwell's equations: where Θ is Hermitian and complex-symmetric. If we separate the fields into incident and scattered components, ψ inc and ψ s respectively, with the incident field as the solution of Θψ inc + iων = J (ν 0 is the background medium), then the scattered field is the solution to where χ = ν − ν 0 . At complex frequencies, passivity requires and, consequently, Im (ωε) > 0 as well. We can use this simple fact to show that two quantities (identical to the absorbed and scattered powers at real frequencies) are positive-definite. We partition the space into two regions, the scatterer occupying the volume V , and the external region V out . We assume the outer region is finite, which can be accomplished using perfectly matched layers [8,9]. Then by the passivity condition, two quantities that must be greater than zero are: where ϕ A and ϕ s are complex extensions of absorbed (non-radiative) and scattered power in that they reduce to the corresponding real-valued powers at real frequencies.
Through two applications of the divergence theorem we can rewrite ϕ s in terms of quantities over the interior volume V , where the overline denotes complex conjugation andn a unit normal pointing outwards from the volume V. Note that in the first line we used the fact that χ = 0 in the exterior volume. This is the key to the derivation: the scattered power is the power radiated to the exterior by currents in the interior volume, and this separation of sources and power flow ensures positivity. If we define ϕ E = ϕ A + ϕ s (motivated by extinction), then we have: Finally, with some notational simplification (using double prime to denote the imaginary part of a matrix), we have We can see that when the frequency is real and the background is lossless (ν 0 = 0), the last three terms are zero and we are left with the first term, which is identical to the extinction at real frequencies. Again, this motivates the interpretation of ϕ E as representing a complex extension of total power. Combining Eq. (S31) and Eq. (S34) ,We can write the constraint ϕ A ≤ ϕ E as or, equivalently, In the main text, ϕ A and ϕ E are defined as the left-and right-hand sides of Eq. (S36). While this differs slightly from the definition given in Eq. (S31) (where ϕ E = ϕ A + ϕ s ) , such difference does not matter -complex extension of power need not be uniquely defined as long as the different extensions agree at real frequencies. In fact, the two ways of defining ϕ A and ϕ E give the same constraint in Eq. (S36).

IX. FORMULATION OF POWER-BANDWIDTH BOUNDS
Returning to the optimization problem of Im s(ω 0 +i∆ω) appearing in average LDOS (see Eq. (16) of the main text), we wish to formulate it in terms of a quantity linear in (total) field, which is then amenable to convex optimization. To this end, we follow the procedure in Eq. (S3) and define a modified incident field ψ inc = E inc −H inc T (the minus sign in front of H inc appears from using reciprocity to the magnetic Green's function). After using the following relation between polarization currents φ and fields ψ, φ = χψ, we obtain the following expression (where ω = ω 0 + i∆ω): In light of Eq. (S37), maximizing Im s(ω 0 + i∆ω) subject to energy-conservation constraints in the UHP reduces to the following optimization problem (evaluated at a complex frequency): where the volume integrals are implicit in the vector-vector multiplications and f = 1 πω 2 ψ inc . For electric LDOS, ψ inc only contains electric fields, which greatly simplifies the problem. But, we keep the term as it is for sake of generality, simplifying only after the optimal solution is derived. Before solving the optimization problem, it is helpful to rewrite it in terms of generic variables: and b = b † (guaranteeing that ψ † bψ is real). A straightforward Lagrange multiplier approach yields the optimal current ψ in terms of a Lagrange multiplier λ: Enforcing the constraint in Eq. (6), per the KKT conditions [10], requires the Lagrange multiplier to satisfy We defer the calculations of the integrals appearing in Eq. (S49) to the next section, where we will see how we can use various approximations to simplify their expressions for canonical geometries. We can also extend Eq. (S49) to the case of magnetic materials as well as magnetic LDOS by modifying the free-space Green's function to accommodate for magnetic fields and/or dipoles.

X. POWER-BANDWIDTH BOUNDS
In our power-bandwidth LDOS bound (Eq. (S49)), we end up with an integral that has 5 terms: where a = 2 Im ω/c is proportional to the bandwidth over which averaging occurs. The 1/r n contribution to the integrand in the first three terms decays sufficiently quickly that the exponential term is not required for convergence.
The term e −ar can simply be bounded above by e −ad and brought out of the integral; for (∆ω)d/c 1, such a simplification is also nearly tight. Thus, we can use bounds for the first term: and the second term:ˆe −ar where E 1 (x) = − Ei(−x) and Ei(−x) is the exponential integral.
For the halfspace, for example, we can put all of this together to bound Eq. (S49) above by For a spherical shell, we get the bound For a 2D plane: Finally, for a 2D spherical surface: For simplicity in the above equations, we have dropped the last term in Eq. (S49) coming from α LDOS , which is negligible for bandwidths ∆ω that are smaller than typical frequencies |ω| considered.

XI. GREEN'S-FUNCTION-RELEVANT INTEGRALS: HALFSPACES AND SPHERES
We consider integrals of various-order 1/r n and e −ar /r n terms, that arise repeatedly in our various bound formulations.

A. Planar integrals
Here we consider integrals from a point to either a 2D planar surface or a 3D halfspace. For the halfspace, the integral of some function f (x) can be written Equation (S60) can be applied to a 2D planar surface simply by removing the integral over z on the right-hand side, and thus the 2D-plane integrals below emerge as intermediate steps of the halfspace solution.

1.´1 r 6 dV
The integral of 1/r 6 iŝ 1 The integral of 1/r 5 iŝ Note that for ad small (ad 1), the halfspace bound simplifies to πe −ad /d, which is exactly what we get when we pull e −ad out of the integrand and multiply it by the integral of 1/r 4 . Hence for the typical case of ad small, pulling the exponential out of the integrand hardly loosens the bound.

B. Spherical integrals
Here we consider integrals from a point to either a 2D spherical surface or a 3D spherical shell of infinite thickess. Then, the integral of some function f (x) can be written  Average electric LDOS at various center frequencies, which show close agreement at center frequencies of 3.65eV and 4.2eV. At frequencies between these two, the difference in average LDOS reflects the fact that the DL model does not match Palik data very well at these frequencies. The emitter-scatterer distance d is set to 10nm, and the scatterer here is Ag halfspace. In Fig. 4 of the main text, we have used the Drude-Lorentz (DL) multi-oscillator model to describe the average LDOS and bounds for Ag, since it behaves nicely for complex frequencies in the UHP. In principle, however, we can analytically continue Palik data [11] for permittivity of Ag to the UHP (but using DL model is simpler as it explicitly provides an analytic expression for permittivity). Figure S2(a) shows how the DL model fit compares to LDOS for bulk Ag (halfspace). In Fig. S2(b), we compute the average LDOS for both the DL model and Palik data. Notice that for the latter, since we do not have data for complex frequencies, we have numerically integrated LDOS (weighted by the Lorentzian in Eq. (15) of main text) over real frequencies. Although the two results are extremely close to each other at center frequencies of 3.65eV and 4.2eV, with greater deviation at 3.8eV and 4.0eV due to the error incurred in approximating the Palik data by a DL model.

XIII. LORENTZIAN APPROXIMATION TO THE MEAN PLANCK ENERGY Θ(ω)
In bounding near-field RHT, we have chosen to replace the mean Planck energy Θ(ω) by a Lorentzian, which is more amenable to contour-integration techniques (involved in bounding RHT). As long as the Lorentzian strictly is larger than Θ for all frequencies, the Lorentzian approximation is indeed a bound for the actual RHT. However, a more ideal bound is one that minimizes their difference, and in Fig. S3 we compare how closely our Lorentzian approximates the mean energy. As Fig. S3(a) illustrates, there is close agreement between the two curves (normalized such that they both have unit magnitude at the origin). In Fig. S3(b), the integrated values are not fully converged, but one can show analytically that the area for the (peak-normalized-to-unity) Lorentzian and the mean energy is √ 2π/2 ≈ 2.22 and π 2 /6 ≈ 1.64, respectively. Thus, the Lorentzian has about 1.35 times more area than for the mean energy Θ as quoted in the main text. Although their high-frequency asymptotic behavior differs qualitatively (polynomial vs. exponential decay), the Lorentzian approximation provides a tight bound to Θ and thus comes in handy for RHT bounds.

XIV. POWER-BANDWIDTH BOUNDS FOR NEAR-FIELD RHT
Given object 1 at temperature T 1 and object 2 at T 2 , we can write the radiative heat flux between the two bodies as H 1→2 (ω) = Φ(ω) [Θ(ω, T 1 ) − Θ(ω, T 2 )], where Θ(ω, T ) denotes the Planck spectrum and Φ is a temperature independent flux rate from incoherent sources in body 1 radiating to body 2. In order to derive upper limits to the heat flux, we need to bound the structure-dependent flux rate Φ(ω) (since Θ(ω, T ) is a known, broad spectrum).
Starting from the heat flux between two bodies expressed in terms of the electromagnetic Poynting vector, we can relate the fields to the fluctuation currents in volume 1, V 1 , in terms of the Green's functions (in index notation): The fluctuation-dissipation theorem [12] gives the ensemble-average of the currents (taking ε 0 = µ 0 = 1 for simplicity): where FD indicates that we are averaging over the fluctuation currents (not to be confused with averaging over frequency bandwidth). Since the Planck spectrum is known, we will simply work with the energy flux Φ(ω) in what follows.
Inserting the average fluctuations into the power, taking the transpose, and using generalized reciprocity [13] (Γ(x s , x) = Γ T (x, x s ) where denotes flipping the sign of 3 × 3 off-diagonal matrices G EH and G HE ) plus Cholesky decomposition on the Hermitian, positive-definite matrix Im χ (expressed as LL † ) leads to where {i, j} indicates the allowed dipole orientations by Λ ji , which is nonzero only for certain combinations of i, j and whose structure is reflected in the coefficients c ij . Note that if the fields were from the same dipole (remember Γ kj (x, x s ) represents the field due to a φ j up to a sign), then the integral over V 1 in Eq. (S78) would correspond to power absorbed in V 1 due to that dipole (note also the Im χ 1 prefactor). Through the Cauchy-Schwarz inequality, we can bound Eq. (S78) exactly in terms of such dissipation: where ψ ξi(xs) (x) denotes the field at x due to a dipolar source ξ i at x s . We can relate this to power absorption due to each individual dipole by inserting an ω/2 prefactor: where the terms in parentheses are now exactly the powers absorbed in V 1 by each dipole. Thus we can write Φ(ω) ≤ 2 πω 2 {i,j}ˆS |c ij | P V1,ξi(xs) (ω)P V1,ξj (xs) (ω).
The power radiated by a given dipole into V 1 is less than or equal to the total power radiated by that dipole. If we denote the polarization-resolved LDOS as ρ i = (6/πω) Im φ † i ψ, with the factor of 6 introduced because we are no longer averaging over dipole orientation, then Here we will make our first approximation: that in the near-field cases of interest, the LDOS ρ i is much larger than its free-space value, i.e. ρ i ρ 0 , in which case we can approximate ρ i as the scattered-field LDOS. Henceforth ρ i will denote the scattered-field LDOS. The inequality becomes Φ(ω) ≤ 1 3 {i,j}ˆS |c ij | ρ ξi(xs) (ω)ρ ξj (xs) (ω).
Now, we want to integrate over frequency and multiply by a bandwidth weighting function H ω0,∆ω (ω). We will hold off on the weight function for a moment, and again use Cauchy-Schwarz to simplify the integrals: Written out in full form, we can bound average CDOS for nonmagnetic materials in the presence of an electric dipole (where r 1 and r 2 respectively denote the distance from the scatter to the source and measurement positions): B n e −ar1 |k| n r n 1 ˆV 6 n=4 B n e −ar2 |k| n r n 2 + 2H ω0,∆ω (0)α CDOS .
For a halfspace enclosure, we can exactly calculate all the integrals appearing in Eq. (S90) as was done in Sec. XI to arrive at the following bound (dropping the last term in Eq. (S90) coming from α CDOS , which is negligible for bandwidths ∆ω small compared to typical frequencies |ω|): where d 1 and d 2 denote the minimum value of r 1 and r 2 respectively. We can also extend Eqs. (S90,S91) for magnetic materials as well as magnetic LDOS by modifying the free-space Green's function to accommodate for magnetic fields and/or dipoles. It can also be extended for 2D materials by modifying the volume to a surface integral.