All electromagnetic scattering bodies are matrix-valued oscillators

Scattering theory is the basis of all linear optical and photonic devices, whose spectral response underpins wide-ranging applications from sensing to energy conversion. Unlike the Shannon theory for communication channels, or the Fano theory for electric circuits, understanding the limits of spectral wave scattering remains a notoriously challenging open problem. We introduce a mathematical scattering representation that inherently embeds fundamental principles of causality and passivity into its elemental degrees of freedom. We use this representation to reveal strong constraints in the mathematical structure of scattered fields, and to develop a general theory of the maximum radiative heat transfer in the near field, resolving a long-standing open question. Our approach can be seamlessly applied to high-interest applications across nanophotonics, and appears extensible to general classical and quantum scattering theory.

Probing and harnessing the frequency dependence of electromagnetic interactions underlies atomic spectroscopy, molecular sensing, information and energy technologies, and more [1][2][3][4].Decades of research into resonant expansions and normal-and quasinormal-mode theories now enable complex scenarios to often be welldescribed by a small number of "physical oscillators" [5][6][7][8].These physical oscillators offer high-accuracy descriptive modeling, but they provide little prescriptive guidance: what lineshapes are physically possible, and what are the ultimate limits of corralling broadband radiation?
To address these foundational questions of wave physics, we introduce a new decomposition of the electromagnetic scattering process into an infinite set of fictitious mathematical "oscillators" that offer a general framework for prescriptive guidance and fundamental limits.Central to our approach is an arguably underappreciated scattering operator, the "T matrix" [9], which relates incident fields to the polarization fields they induce.The complete scattering response of any scattering body is determined by its T matrix, but what form can this matrix take?We show that causality and passivity impose strong constraints on its possible form.The culmination of this insight is a powerful and general representation: every scattering body's T matrix must be a superposition of lossless Drude-Lorentz oscillators with matrix-valued (spatially nonlocal) coefficients.Moreover, the oscillator coefficients are highly constrained, with three properties absent in other approaches: (1) they are Hermitian, even in lossy and open systems, (2) they have finite sum rules in many scenarios, and (3) they are positive-semidefinite in passive systems.As a whole, we arrive at a striking conclusion: the controllable degrees of freedom of any scattering body, for any application and spectral range, are only the (highly constrained) matrixvalued oscillator strengths.Such limited DOFs must im-ply strong constraints on scattering response; more generally, this new representation offers a general method for identifying fundamental limits to spectral control.
To demonstrate the power of this approach, we use our oscillator decomposition of T matrices to resolve a long-standing question in energy transport: what is the maximum rate at which two bodies can radiatively exchange heat in the near field?Going back many decades, it has been understood that radiative heat exchange in the near field can be substantially larger than its far-field counterpart [10][11][12], due to the enormous number of accessible evanescent channels in addition to propagating ones, yet the maximum extent of this enhancement-with ramifications for applications such as thermophotovoltaics [13,14], photonic refrigeration [15], and heat-assisted magnetic recording [16]-has been far less clear.Previous theoretical bounds [17][18][19][20][21] have suggested strong material-electron-density dependencies, unbounded response for low-loss materials, and ordersof-magnitude gaps from known designs (> 750X).Using our T-matrix representation, we offer a new simple, general theoretical limit that is within a small factor (5X) of the state-of-the-art.Encapsulated in our approach is a sum rule that explains why planar structures are better than sharp-tip patterns for large-area, broadband nearfield enhancements, and why unconventional plasmonic materials should offer the largest enhancements.

I. THE REFRACTIVE-INDEX PARADOX
At a microscopic level, accurate computations of the linear optical response (e.g., refractive index) of a material is sufficiently difficult to be at the frontier of modern electronic-structure methods [22].Yet modeling uncertainty does not imply unbounded possibility: there are well-known constraints to allowable refractive index [23][24][25].This understanding comes not from densityfunctional-theoretic constraints, but instead fundamental principles such as causality and passivity.At the macroscopic level of wave scattering, however, the situation is reversed.With fast-solver techniques pioneered over three decades, it is possible to compute the electromagnetic fields of structures that are thousands of wavelengths in size at machine precision [26].Yet our understanding of the extreme responses of those same structures is highly limited.Why has our understanding of limits at the macroscopic scale, where computational methods are more successful, not reached our understanding of limits at the microscopic scale?It is because the combined principles of causality and passivity have not been integrated into scattering-matrix representations.This shortcoming can be remedied in the scattering T matrix.

II. THE SCATTERING T MATRIX
In linear, time-invariant electrodynamics, there must be a linear operator that relates electromagnetic fields incident upon a scatterer to the polarization fields they induce; the standard name for this response function is the "T operator."For simplicity of notation and exposition, we assume any standard numerical discretization of sufficiently high accuracy such that the T operator becomes a T matrix; if we collate the incident fields E inc (x) into a vector e inc and the polarization fields P(x) into a vector p, then the frequency-domain relation (e −iωt sign convention) defining the corresponding T matrix is, p(ω) = T(ω)e inc (ω), (1) which is the discrete analog of the convolution equation P(x, ω) = ´T(x, x , ω)E inc (x , ω) dx .The T matrix can be derived from first principles via integral operators, as discussed in Ref. [9] and the SM.The key property of the T matrix, for our purposes, is that it is a causal response function: the polarization field at any x cannot be excited by an incident field at x until after the incident field has reached x .This is similar to the causality condition for material susceptibilities [27], and it has an analogous consequence.Using Fourier-based arguments that parallel those for material susceptibilities, we find a Kramers-Kronig-like equation that must be satisfied by any T matrix (detailed proof in SM): Equation ( 2) is a Kramers-Kronig (or Hilbert transform) relation for arbitrary scattering processes, for any underlying material.An important feature of Eq. ( 2) is that it is matrix-valued; the "Re" and "Im" operators denote the Hermitian and anti-Hermitian parts of their arguments, respectively.Equation (2) implies that once the anti-Hermitian part of the T matrix is known for all frequencies, its Hermitian part has been determined as well, and vice versa.
Next we identify sum rules for the T matrix.There are two special frequencies at which the left-hand side of Eq. ( 2) simplifies: zero frequency (electrostatics) and infinite frequency (where materials become transparent).Which is more useful depends on the frequency range of a given application; we include here the high-frequency sum rule, and leave the low-frequency sum rule for the SM.At infinite frequency, the electrons of a material can be regarded as free, and material susceptibilities must scale as χ(ω) → −ω 2 p /ω 2 , where ω 2 p is proportional to the total electron density of the material [27].In this limit, the first Born approximation is asymptotically exact, and the polarization field is given by P χE inc −(ω 2 p /ω 2 )E inc (in units where the free-space permittivity is 1), implying that the T matrix asymptotically approaches −(ω 2 p /ω 2 )I, where I is the identity matrix.Inserting this limit into the KK relation of Eq. ( 2) yields the sum rule, This sum rule constrains the total contributions from Im T(ω) over all frequencies, much like the f sum rule for material-susceptibility oscillator strengths [28][29][30].
The final piece of the puzzle is passivity.In passive scatterers, the polarization fields do no net work.Mathematically, the work done by the incident fields on the polarization currents J is 1  2 Re ´E * inc •J = ω 2 Im ´E * inc •P; positivity of this expression implies, at positive frequencies, the positive-semidefinite condition, Im T(ω) ≥ 0.

III. OSCILLATOR REPRESENTATION
We now have three key ingredients: a Kramers-Kronig relation, a sum rule, and a positivity constraint.Together, we can use these to form a new spectral representation of any T matrix.Again, for simplicity of exposition, we work in a discrete basis, discussing the straightforward generalization to continuous frequencies in the SM.Consider a discrete set of frequencies ω i from 0 to ∞ at which we equate ω Im T(ω), with appropriate prefactors, to a set of frequency-independent Hermitian matrices T i : By passivity, all T i are positive semidefinite, T i ≥ 0, and by the sum rule, constrained to a finite sum: i T i ≤ I. Finally, the KK relation of Eq. ( 2) dictates that by encoding the anti-Hermitian part of T(ω) into matrix variables T i , the Hermitian part of T(ω) is completely determined; there are no more degrees of freedom.Moreover, the spectral lineshape of the Hermitian part around every frequency ω i is that of a lossless Drude-Lorentz oscillator.A few mathematical steps can make this connection more rigorous ([SM]), such that the three ingredients are encoded in a single, general representation: along with the conditions that T i ≥ 0 for all i and i T i = I.Equation ( 5) is the foundational result of our paper: the T matrix of every linear scattering body must be decomposable into a (fictitious) set of lossless Drude-Lorentz oscillators, with matrix-valued coefficients that are positive-semidefinite and constrained in total strength.This representation is contrasted with standard oscillator decompositions (coupled-mode theory, quasinormal modes) in the SM.In particular, we stress that this is not an eigendecomposition, and ω i are not eigenfrequencies of the operator.Given an electron density, the only degrees of freedom in the scattering process are the matrices T i .Moreover, the constraints on these matrices (T i ≥ 0, i T i ≤ I) are convex sets, and the T matrix is linear (and thereby convex) in these matrices.Hence, not only does Eq.( 5) appear to offer a strong constraint on possible spectral response, but its constituents are ideally suited for mathematical optimization and consequent fundamental limits.
As a first demonstration of the surprising structure implied by this representation, we consider broadband scattering from an elliptical dielectric cylinder.The scattered electric field (total field minus incident field) at various points within the scatterer, computed by full-wave simulations [SM], is shown in Fig. 1(b), but is hard to interpret due to its seemingly random undulations.The tremendous advances in quasinormal-mode (QNM) techniques suggest that one could accurately reproduce these fields with a relatively small number of QNMs [8], but that modeling capability does not imply an understanding of the extreme limits of what is possible.How many resonances can be excited?With what amplitudes, phases, and overlaps with power-carrying channels?
By contrast, consider the diagonal elements of the real (Hermitian) and imaginary (anti-Hermitian) parts of the T matrix, as depicted in Fig. 1(c).They show the characteristic features of oscillators!The lineshapes of the T-matrix elements mimic exactly the Drude-Lorentz-like behavior of electronic transitions, but they arise not from real material oscillators (the susceptibility is constant, Fig. 1(f)), but from complex wave-scattering behavior itself.Included in Fig. 1(d) are the matrix-valued coef-ficients at three frequencies, with the eigenvalues of all frequencies plotted in Fig. 1(e).They are all positive, as guaranteed by passivity.By uncovering this structure, we can now see that scattering processes must be decomposable into a collection of positive-semidefinite matrices with finite total strength.Such decompositions are ideally suited for identifying fundamental limits to spectral control.

IV. ULTIMATE LIMITS TO NFRHT
Now we apply our formulation to the question of maximal NFRHT.It has long been known that bringing two bodies closer than a characteristic thermal wavelength can lead to radiative transfer at rates far beyond the blackbody [10][11][12], via near-field evanescent tunneling.Yet, over all possible geometrical configurations, what is the maximum possible rate of energy transfer?NFRHT, as depicted in Fig. 2(a), faces prohibitive computational challenges-spatially and temporally incoherent, broadband thermal sources, exciting rapidly decaying near fields over large macroscopic areas-which has limited previous design efforts primarily to high-symmetry structures such as planar bodies [31][32][33].Numerous approaches have identified particular constraints with corresponding theoretical bounds [17][18][19][20][21], but as we show in Fig. 2, there are orders-of-magnitude differences between the best structures and the best bounds [20,21].We label the bounds by their distinguishing attributes: in Ref. [20] ("analyticity bound"), complex-analyticity played the central role for a finite bound, while in Ref. [21] ("channel bound"), a decomposition into power-carrying channels was the starting point.Last year, it was discovered that a set of unconventional plasmonic materials offer significant (10X) improvements over the previous best planar structures [34], but otherwise the field has been at an impasse: there has been no way to identify neither better structures nor a better theory of the upper limits.
The T matrix formulation resolves this impasse.The heat transfer coefficient (HTC) between two bodies is the net flux rate (per area and per degree K) of electromagnetic energy passing between bodies at temperatures T and T + ∆T , as measured by the integral of (1/2) Re (E × H * • n) through a separating plane with normal vector n.
The incoherent sources in body i with temperature T i and susceptibility χ i (ω), by the fluctuation-dissipation theorem [31], are given by J j (x, ω)J * k (x , ω) = (4ε 0 ω/π)Θ(ω, T i ) Im [χ i (ω)] δ jk δ(x − x ) at frequency ω, where Θ(ω, T i ) is the Planck spectrum.There are a variety of mathematical transformations that we make to this problem to make it more amenable to optimization, detailed in the SM, such as using reciprocity to move the sources out of the hotter body and onto the dividing surface, exploiting spatial symmetries of the bounding domains (two halfspaces, allowing for any patterning within), as well as a near-field generalization of the "op-tical theorem" [35].The key novelty, however, is our use of Eq. ( 5): once we have transformed the problem to an appropriate function of the two-body T matrix, we insert the representation theorem of T as a sum of (unknown) positive-semidefinite matrices with Drude-Lorentz lineshapes.NFRHT at moderate or high temperatures is dominated by low-frequency response (relative to the plasma frequency of common materials), which is proportional to an electrostatic constant α instead of the electron-density encoded in ω p .Use of the low-frequency sum rules simply requires replacement of ω 2 p with ω 2 i α in the representation theorem; moreover, the electrostatic constant is bounded above by the numerical value of 2 for high-contrast dielectric or metallic bodies enclosed in two halfspaces (SM).Once we insert the T-matrix representation into the NFRHT expression, the resulting optimization problem over the infinite set of matrix oscillator coefficients has an analytical upper bound.Straightforward algebraic manipulations (SM) lead to an ultimate limit to near-field radiative HTC given by where d is the minimum separation between the bodies, T the temperature of the cooler body, and β = 0.11(αk 2 B / ) = 3.8 × 10 5 Wnm 2 /m 2 /K 2 , a numerical constant.This limit cannot be surpassed by any geometric patterning, nor can exotic optical properties of any material alter its value.
Fig. 2(b) compares our theoretical limit with the current state-of-the-art, as well as the best known bounds.Whereas the gap between the optimal planar structures and the best previous bounds was at least 750X (but diverging to ∞ for some materials), the expression of Eq. ( 6) is only 5X larger the best design.The new bound has no material dependence, which resolves the problematic trend that if one orders the materials by their planar performance, as in Fig. 2(b), the previous bounds tended to predict worse maximal performance in the same direction.The resolution of this discrepancy is our use of the low-frequency sum rule, which encodes a constraint on the local density of states seen by thermal emitters that depends only on their gap separation, independent of material.Although it seemed possible that nano-structuring may lead to enhanced NFRHT through field-concentration (lightning-rod) effects, our sum rule explains why this is not the case: sharp tips can enhance the fields very close to a sharp tip, but not at the source location itself.The local density of states is proportional to the latter, and hence is not enhanced by lightning-rod effects.SM Sect.X illustrates this point, showing enormous field enhancements at points closer to a sharp tip than any source, but, surprisingly, smaller fields at the source itself, no matter how close the source is to the tip.The T-matrix approach predicts an optimal NFRHT frequency of ω max = 2.57 kBT , determined by the overlap of the Planck function with the Drude-Lorentz lineshape.The predictions are matched almost Heat-transfer coefficients of planar bodies comprising increasingly high-performance planar-body geometries (filled circles), in comparison with the best previous theoretical bounds [20,21] (open squares and triangles).The previous bounds diverged for some materials, while showing enormous gaps (> 750X) for others, and their trendline seems to decrease left-to-right, whereas planar-body performance increases.In black is the new theoretical bound offered by our T-matrix representation, very close to the best possible planar bodies.Panels (c,d) confirm the predictions of our spectral representation, showing that the state-of-the-art performance could be achieved precisely at the optimal frequencies identified by our approach.
exactly by computationally optimized planar Drude metals or 2D heterostructures, as shown in Fig. 2(c,d).For 300 K temperature, the spectra shown in Fig. 2(c) peak at almost exactly the optimal oscillator frequency, and the match persists across all relevant temperatures, as shown in Fig. 2(d).
The closeness of the arbitrary-structure bound of Eq. ( 6) to the best planar structures arises despite quite different mathematical routes to these results.The translational symmetry of planar bodies implies conserved wavevectors and thus a set of evanescent plane-wave channels that are independent, with Landauer-like transmissivities [33].Such an approach cannot describe patterned structures, and Eq. ( 6) culminates after using (generalized) reciprocity to move the sources from the hot body to the dividing surface, the sum rule to encapsulate the maximum densities of states seen by those sources, and the T-matrix representation to constrain the possible scattering lineshapes.The striking similarity of the two results suggests that even when confronted by spatial and temporal incoherence, rapidly decaying fields, and large areas, the oscillator representation compactly captures the key physics of maximal response in the near field.

V. A GENERAL FRAMEWORK
In this article we have introduced a new viewpoint on electromagnetic scattering.The example of Fig. 1 reveals a surprising structure that has previously been hidden underneath the complexity of scattering dynamics.Our application of this framework to NFRHT offers clear guidance for the fundamental limits of radiative heat transfer and the physical mechanisms underlying them.The generality of our T matrix representation offers tantalizing prospects for wide-ranging applications across nanophotonics.Metasurfaces [36,37], for example, offer a new, compact form factor for optics.A central question is the extent to which metasurfaces can control incoming waves [38][39][40], across varying frequency and angular bandwidths, for applications from lenses to virtual and augmented reality.Similarly, techniques for imaging through opaque media have flourished with modern spatial light modulators [41], with a key open question being ultimate limits to spectral control.In photovoltaics as well as photodetection, the quest for ever-thinner devices must ultimately contend with fundamental limits.And similarly across almost every application of nanophotonics.There has long been a need to quantify the ultimate limits to spectral control; our approach offers a theory to do so.
Our approach also dovetails seamlessly with a recent flurry of activity in understanding the limits controlling spatial degrees of freedom in nanophotonic systems [42][43][44][45].Transforming the typical Maxwell differential equations into a set of local conservation laws in space, for real and reactive power flows, leads to a mathematical form of the design problem that is amenable to systematic approaches to computational bounds.For a single frequency (or a small number of them [46]), such conservation laws have shown powerful capabilities for identifying fundamental limits to spatial control.In these approaches, the degrees of freedom of the system are typically encoded not in the electric and magnetic fields, but rather in the electric and magnetic polarization currents that they induce.Those polarization fields are exactly those that are determined by the T matrix, which means that our spectral expansion of the T matrix should be seamlessly compatible with the spatial conservation laws proposed in Ref. [42,43].Together, the two approaches may enable a complete understanding of the spatio-spectral limits of electromagnetic systems.
More broadly, the insight at the foundation of our framework, about the mathematical properties of scattering T matrices, can be directly applied to any classical wave equation.These techniques should be readily extensible to linear scattering problems in acoustics, elasticity, fluid dynamics, and beyond.The mathematical structure of the wave equation is similar in each case, and the resulting T matrices should therefore have similar representations.An interesting twist may arise in acoustic scattering theory, where materials with higher-thanvacuum speeds of sound lead to "non-causal" scattering processes [47] that have prevented the development of classical sum rules, and would appear to prohibit a corresponding T matrix representation.Yet the T matrix itself may offer a new route to complex-analytic response functions in exactly such scenarios.The reason higher sound speeds lead to "non-causal" response is that the scattered field appears at a location within the scatterer earlier than the incident wave itself.Hence, locally, the process appears non-causal.Yet the nonlocal nature of the T matrix may be precisely what is needed to resolve this paradox.A T matrix isolates the response at any point x to the contributions from the wave incident at each point x in the scatterer; each of which, individually, must be causal.Hence, not only should the T matrix be extensible to such scenarios; it may further resolve impediments that had previously stymied even simple sum rules in these fields.
Finally, we speculate that the approach described here may even be extensible to quantum scattering.In the frequency domain, the key difference between quantum and classical scattering is the analytic structure of the governing equations.In classical wave equations, second derivatives in space are proportional to second derivatives in time, which lead to poles in the lower half of the complex-frequency plane and analyticity in the upper half.In quantum scattering, second derivatives in space are proportional to first derivatives in time, which leads to bounds states for negative real energies and branch cuts on the positive real axis.Our standard semicircular contours likely need to be replaced by "keyhole" contours [27], with the open question of whether there are meaningful sum rules that can be derived (perhaps dependent on bound-state properties, as in Levinson's theorem [48,49] for spherically symmetry potentials).If such sum rules could be derived, it is likely that an infinite-oscillator description could be used to identify fundamental limits for quantum scattering as well.

I. PHYSICAL OSCILLATORS VS. MATHEMATICAL OSCILLATORS
In the main text we contrasted our T-matrix "mathematical-oscillator" representation with the well-known "physical-oscillator" decompositions in use today.Here we detail the similarities and differences in these approaches.At the highest level, physical-oscillator approaches are meant to be efficient for simulation and modeling: for a given structure, can one identify a small number of parameters (e.g.normal-or quasinormal-mode coefficients, etc.) that accurately model the complete response of the system?Typically these models are highly nonlinear in the unknown parameters, but there are standard computational methods for finding them.Yet for problems of design, these representations are difficult or impossible to work with: there is not a single given structure, anymore, but instead a large class of structures.It is not known a priori how many resonances or modes may contribute; to be safe, very large numbers must be used.So one is left with large, highly nonlinear models to optimize over, without any beneficial mathematical structure.The "mathematical-oscillator" theory developed in the main text is well-suited to this scenario.
The use of lossless (and hence narrow-linewidth) oscillators naturally also leads to a large numbers of parameters (matrix-valued oscillator coefficients), but these parameters have ideal properties from an optimization perspective: they are positive-definite, constrained in sum, and linear in the only degrees of freedom.In fact, this decomposition is quite illsuited for modeling: one needs to invert a large, dense matrix at every frequency to get the corresponding coefficients, which is computationally prohibitive except for small structures.
But for design, one never needs to do any matrix inversion, and the mathematical structure of the decomposition is far superior for optimization.Below, we provide the mathematical expressions supporting these qualitative assertions.
In electromagnetic scattering simulations, there are two primary classes of "physicaloscillator" approaches: coupled-mode theory (CMT) [1][2][3][4][5], and quasinormal-mode (QNM) theories [6][7][8][9][10][11][12].It can be shown that the former can be derived from the latter in the limit of isolated, high-Q resonances with negligible non-resonant scattering contributions [13].We will describe both of these constructions.We start with coupled-mode theory.In coupledmode theory, there is a basis of resonant modes described by a matrix Ω, whose diagonal terms are complex-valued resonance frequencies of the modes, and whose off-diagonal terms describe coupling rates between each pair of modes.Each resonance has "'overlap coeffients" with each outgoing-wave channel, the matrix containing these elements is often referred to as K (or D, which is identical to K in reciprocal systems).Finally, there is typically a nontrivial background scattering matrix S bg that contains non-resonant contributions to the scattering process.In full, in any general (reciprocal) coupled-mode theory, the scattering matrix is given by the expression [1,13] subject to reciprocity and unitarity conditions given by One can immediately see that a CMT model constructed from these three equations will be impossible to optimize over.The degrees of freedom are the matrices Ω, S bg , and K, none of which are Hermitian (let alone positive definite).Moreover, one cannot even presuppose any finite, constrained size of the matrices, as there is no sum rule constraining any norm of the entries.
One route towards using CMT models for understanding limits is to remove much of the complexity from Eqs. ( 1)-( 3).If one assumes background scattering cannot occur (although note that even in Mie resonators it plays a quite important role [13]), then S bg = I and K * = −K, such that one can rewrite the scattering matrix relation as a form of the S-matrix that also arises in nuclear scattering theory [14,15].Next, one can assume that none of the modes are coupled, such that Ω is a diagonal matrix.(This condition will typically conflict with the requirement that K † K = 2 Im Ω, but we ignore that for simplicity.)Finally, special quantities such as absorption (given by I − S † S [16]) can be written as: after repeated use of From Eq. ( 5), one can integrate over all frequencies to simplify the interior matrix product involving frequencies (which is a diagonal matrix with Lorentzians along the diagonal) to a constant.Finally, one is left with an expression involving only the loss rates of each resonator and the number of resonances [17][18][19].But what are these values?How large can they be?One is always left with more free unconstrained parameters.And the extreme limits of these models, where one wants to operate for fundamental limits, are precisely where the assumptions mentioned above (high-Q resonances, uncoupled resonances, frequencyindependent K matrix, isolated resonances, no background processes, etc.) break down.
To move beyond the assumptions of CMT, expansions via quasinormal modes (QNMs) have become more popular in recent years.There are various expansion techniques, many of which can be shown to be equivalent [7].We will assume a Maxwell equation of the form where we assume a sufficiently high-resolution discretization of Maxwell's equations, for matrix M and diagonal matrix ε, unknown electric-field vector e, and free-current source vector j.The boundary conditions or PMLs are assumed to be encoded in the matrix M , as well as the curl-curl operator.The pair of matrices M and ε form generalized eigenproblem pairs according to M U = εU Λ, where U are the eigenfields and Λ the squared eigenfrequencies.We can assume reciprocity, in which case U −1 = U T .Inserting this eigendecomposition into our Maxwell equation and solving for the electric field yields One can interpret Eq. ( 7) intuitively: U T ε −1 j is a decomposition of normalized free currents into modal fields, (Λ − ω 2 ) −1 is the resonant enhancement associated with real frequencies close to the resonant frequencies, and the final U on the left converts from the modal basis back to the original (real-space) basis.Superficially, Eq. ( 7) actually looks quite similar to the coupled-mode scattering-matrix equation of Eq. ( 1): a resonant amplification term inversely proportional to the differences between the complex-valued resonant frequencies and the real excitation frequencies, and frequency-independent matrices surrounding the resonant-amplification term.However, to connect to the "scattering channels" that bring energy into or out of such systems, one would need to pre-and post-multiply these matrices with Green's-function matrices that are highly frequency-dependent.Then, again, one is left with a complex set of degrees of freedom: the number of resonances (the size of Λ and number of columns of U ), the locations of the resonant poles in the complex plane (the values of Λ), and the resonant field patterns (the values of the columns of U ).There are almost no constraints on these degrees of freedom, except that the field patterns must be orthogonal in the unconjugated inner product, corresponding to U T U = I.There is no meaningful way to convert this representation to upper bounds or fundamental limits.Again, however, this representation is quite useful for modeling, with a number of exemplary successes over the past decade [6][7][8][9][10][11][12]26].
To summarize: Eqs.(1,7) are the key "physical-oscillator" descriptions of classical scattering processes.At a glance, they share similarities with each other and with the T-matrix representation of the main text: at a coarse level, each has a resonant-enhancement term and one or more matrices that can be described as a "coupling" matrix.With more granularity, however, there are crucial mathematical differences between the two expressions of Eqs.(1,7) with the T matrix expression.In the CMT and QNM approaches, all of the degrees of freedom (the resonant pole locations and the coupling matrices) are complex-valued quantities without any Hermiticity or positive-definiteness qualities.Moreover, the number of resonances can never be constrained for the arbitrarily patterned nanophotonic systems of interest.By contrast, in the T matrix expansion, all of the "resonant poles" are approaching the real axis, there is an infinite set (one need not limit the number of "resonances"), the degrees of freedom (the scattering-oscillator strengths) are positive semidefinite, Hermitian matrices, and their sum is constrained, thanks to sum rules.Because their poles are not related to the normal-or quasinormal-mode eigenfrequencies, T matrix expansions are computationally expensive for a given structure.But in optimizations over all possible geometries, their mathematical structure is unique, and pays significant dividends.

II. INTEGRAL-OPERATOR DEFINITION OF THE T MATRIX
In the main text, we used linearity as a sufficient condition to argue that the polarization fields P(x) induced by an incident field E inc (x) must be related through a linear operator that one can call "T:" or, in our vector notation, In this section, we show how to construct (or define) the T matrix from known integralequation operators.
For any scattering problem, one can solve for all field degrees of freedom within only the body of the scatterer, from which all other fields are computable.One starts by writing the total field as the sum of the incident and scattered fields, where E inc (x) is known (for a given application).The scattered field is given by the convolution of the free-space (or background) Green's function G 0 with the polarization fields: Within the scatterer, the total field E(x) equals the polarization field P(x) divided by the susceptibility χ(x) (which is nonzero everywhere by definition of "within the scatterer"), so that we can rewrite the decomposition of the electric field in Eq. (10) in terms of the polarization fields or, in vector notation, which in electromagnetism is known as the volume integral equation, or Lippmann-Schwinger equation [27].We can simply invert the matrices in square brackets, to find the definition of T: Hence the T matrix is the inverse of well-known integral-equation matrices, which can be computed via standard methods [27].

III. T MATRIX KRAMERS-KRONIG RELATION: RECIPROCAL CASE
In this section, we fill in the details for the derivation of our "Kramers-Kronig" T matrix equation, Eq. ( 2) of the main text.Our derivation mirrors the derivation of conventional KK relations for susceptibilities, cf.[28,29].
We start with the time-domain equation defining the T matrix as the linear matrix connecting incident fields to induced polarization fields: This definition enables us to encode causality in the T: at any point x in the scatterer, the polarization field at that point cannot be generated from the incident field at any other point x until after the incident field has reached that point x .We can define the first time at which the incident field arrives at any point in the scatterer as the time origin, in which case causality implies that where the right-hand side should be interpreted as the zero matrix, i.e., every entry is zero.
The frequency-domain T matrix is the Fourier transform of the time-domain T matrix: where the integrand starts at zero in the second line because of the causality condition of Eq. ( 17).We assume the integral in Eq. ( 18) is well-behaved for all real frequencies.Then, if one inserts a complex-valued frequency into the expression of Eq. ( 18), it is guaranteed to converge (under suitable conditions [28]), which leads to the usual property of complexanalyticity in the upper half of the complex-frequency plane, now for each entry of the T matrix.
Once we have complex-analyticity, we can consider an integral of the form where C is the typical semicircular contour in the complex plane with a small deviation around the frequency ω.The integral in Eq. ( 19) has a simple pole at ω. Assuming sufficient decay at infinite frequencies in the upper-half plane (where physical materials become transparent and the T matrix goes to zero), the semicircular arc in the upper-half plane does not contribute to the integral, and one is left only with the principal-valued integral along the real line at all points except ω, and the residue term at ω given by: The total integral in Eq. ( 19) equals zero, so separating the principal-value and residue terms and equating their negatives, we have (assuming principal-value integrals where needed): This is a matrix-valued equation that applies entrywise on each side.We can "take the imaginary part of each side," really meaning that on each side we take the difference between the current term and its conjugate transpose, and divide by 2i.This "imaginary part" of Eq. ( 21) then reads (moving the terms on the left to the right and vice versa), where, as in the main text, "Re" and "Im" denote the Hermitian and anti-Hermitian parts of their matrix arguments.
Finally, we can use symmetry of the T operator around the origin.For an electric field or a polarization field, the usual symmetry relation is E(−ω) = E * (ω) and P(−ω) = P * (ω), which arises via the real-valued nature of the time-domain fields and the Fourier-transform definition of the frequency-domain fields.By exactly the same reasoning for the T matrix, we have: Because of reciprocity, T ij = T ji , we can equivalently write T(−ω) = T † (ω), which implies a symmetry in the matrix imaginary part of T: Hence the integral over the real line can be simplified to positive frequencies: Now we have our "KK relation:" IV. T MATRIX KK RELATION AND REPRESENTATION THEOREM: NONRE-

CIPROCAL CASE
In this section we identify the KK relation, and general representation theorem, for a scattering T matrix in the case when reciprocity does not apply.
First, we want a "KK relation."In this case, we can just use Eq. ( 22): This, already, is nearly an oscillator representation, including both positive and negative oscillator frequencies.We want the expression in the numerator to be positive semidefinite, but in this case, that is only true at positive frequencies.However, ω Im T(ω ) is positive semidefinite across all frequencies, positive and negative.So we can multiply by ω to find: This will form the basis for our oscillator representation, and can be regarded as a mathematical-oscillator representation in nonreciprocal (as well as reciprocal) systems.
Next we need a passivity condition and sum rules.For passive systems, the polarization fields do no net work, which implies that at all (positive and negative) frequencies.
Lastly, we need sum rules.The zero-frequency sum rule can be derived directly from Eq. ( 28); if we denote Re T(ω = 0) = α, then A high-frequency sum rule is harder to identify.The intuition is as follows: we can decompose any T matrix into a "reciprocal" component, that is complex-symmetric, and a "nonreciprocal" (or "anti-reciprocal" component), which is skew-symmetric, at every frequency.The reciprocal component satisfies the matrix analogs of the usual symmetry rules around the imaginary frequency axis, while the nonreciprocal component satisfies the opposite.Over all negative and positive frequencies, the reciprocal component will comprise the entire sum rule, while the nonreciprocal component will contribute zero.Let's work through this mathematically.
We start with the symmetry relation about the origin, Eq. ( 23).The key point about Eq. ( 23) is that it is entrywise, i.e., each entry satisfies its own symmetry relation.This prohibits direct symmetry relations for the Hermitian and anti-Hermitian parts of T in the nonreciprocal case.(In the reciprocal case, we could use the transpose operator to turn the conjugation into a conjugate transpose operation.) The key, then, is to break the T matrix into its complex-symmetric and skew-symmetry parts, which we can refer to as its "reciprocal" and "nonreciprocal" (or really, "antireciprocal") parts: where This is useful because now each component obeys matrix symmetry relations: which, in matrix form, can be written: If we look for a sum rule for the quantity then we can insert the symmetric and skew-symmetric parts of T into this expression, and use the symmetry relations of Eq. ( 35) to find The first integrand, ω Im X(ω), is symmetric around the origin, while the second is antisymmetry and thus integrates to zero.Hence, Then we only need to find a sum rule for the symmetric (reciprocal) part of T(ω) to get a sum rule over all frequencies for ω Im T(ω).
To find a sum rule for X(ω), it is useful to find a separate KK relation for this component.
We can do this from the KK relation for T(ω), Eq. ( 21).We can write this down for two entries of the T matrix: the ij entry and the ji entry: We can add these two equations together to cancel the skew-symmetric parts, leaving only This is a KK relation for X(ω); moreover, X(ω) obeys the anticipated matrix symmetry relations around the imaginary axis, so we can write Next, we use the high-frequency asymptote of T(ω), T(ω) → −ω 2 p /ω 2 I as ω → ∞ to write Inserting this into the KK relation for X, Eq. ( 44), we get, without any assumption of reciprocity: We know this integral determines the integral of ω Im T(ω), and indeed can insert it into Eq.( 38) to find: Now, over negative and positive frequencies, we have our three needed items: a positive semidefinite quantity, a Kramers-Kronig relation in terms of that positive-semidefinite quantity, and a sum rule (both zero-and high-frequency).We can compile here in one place: At first blush, it might seem from the sum rule, which now has πω 2 p instead of πω 2 p /2 on its right-hand side, like you should be able to do a factor of two (or more) better in the nonreciprocal case than in the reciprocal case.But the symmetry condition on T, that each entry equal the conjugate of its negative-frequency counterpart, imposes a strict condition on ω Im T(ω): Hence, we must always put a concomitant oscillator at negative frequencies when an oscillator is placed at positive frequencies.Hence, for any optimization problem over all possible oscillator strengths, one should also impose this symmetry condition.
But for just about any objective, there is an even simpler optimization trick we can play to enforce this symmetry constraint.Instead of enforcing the constraint in the oscillator strength itself, we can enforce it in the objective.The fields E, H, P, etc., all satisfy the usual symmetry condition s(−ω) = s * (ω), and any objective made of these functions appears to be symmetric about the origin.We can use Poynting flux as an example: Its negative-frequency counterpart is Because of this, instead of optimizing a positive-frequency objective subject to the constraint of symmetry in the oscillator strength, we can drop the constraint of symmetry in the oscillator strength, and instead maximize any linear combination of the negative-and positive-frequency objectives: Typically it may be sufficient to simply choose α = 1/2.Hence there are no additional degrees of freedom in the negative frequencies.
One more version of a representation theorem in the nonreciprocal case can be specified, and this case most neatly allows comparison with the representation theorem in the reciprocal case.We showed above that there is a KK relation for X(ω), the complex-symmetric part of T, Eq. ( 44).Similarly there is a KK relation for Y(ω); it follows the same derivation, but now has the opposite symmetry around the origin, so that the final relation becomes Notice that the term in the numerator of the fraction in the integrand has ω now, and not ω .We can put the KK relations for X and Y together to get a positive-frequency-only representation of T(ω): where we have intentionally written the numerator of the integrand in terms of the component parts of ω Im T(ω ) = ω Im X(ω ) + ω Im Y(ω ).Let us define two new matrices U and V such that which inherit the sum rules and symmetries of X and Y, such that U is symmetric about the frequency origin while V is anti-symmetric about the frequency origin.Moreover V is real-symmetric (in space) while V is anti-symmetric.We do not necessarily know that whether each individual matrix is positive semidefinite, but we know that their sum is, at all positive and negative frequencies.Hence, for a positive frequency ω, The second condition, using the symmetry relations about the origin, can be converted to Together, these two conditions imply that U is positive semidefinite, U(ω) ≥ 0. So we now have: For nonreciprocal problems, then, we get to choose two matrices at every frequency: one realsymmetry and one skew-symmetric.The real-symmetric matrix must be positive definite and has an all-positive-frequency sum rule.The skew-symmetry part has no sum rule but must satisfy two positivity conditions, in tandem with the symmetric part, at all frequencies.
The skew-symmetric matrix is the only new degree of freedom in nonreciprocal systems.

V. LOW-FREQUENCY SUM RULE FOR THE T MATRIX
In the main text, we discussed the simple method for deriving a high-frequency sum rule for ω Im T(ω), by taking the pole frequency to infinity and using the transparency of materials at high frequencies to simplify the analysis.The other special frequency at which the analysis is simplified is zero, where the scattering problem becomes electrostatic.In this section, we derive the low-frequency sum rule, and the corresponding representation theorem.
As a reminder, the KK relation for To derive a low-frequency sum rule, we simply evaluate this expression at ω = 0: The matrix T(ω = 0) at zero frequency is a generalized polarizability matrix (since it relates the induced dipole density to the incident field), which we can denote by α.This matrix is real-symmetric per the discussion in the previous section: decomposing T(ω) into its symmetric (X) and anti-symmetric (Y) parts, the latter has zero Hermitian part at zero frequency.Now our sum rule reads which further shows that the generalized polarizability tensor is positive-semidefinite (since the left-hand side is).
As discussed in Sec.VIII, this generalized polarizability tensor obeys "domain monotonicity," meaning that the difference between the polarizability tensor of a domain and any second domain that is enclosed by the first must be positive-definite.Because of this, we can typically choose a high-symmetry enclosure (e.g.half-spaces) for which the polarizability is a scalar α multiplied by the identity matrix, and all possible geometric domains inside of them will have sums that satisfy: From the sum rule of Eq. ( 65), we can form a representation of T similar to the highfrequency representation of the main text.If we start by choosing Im T(ω)/ω to be a summation of delta functions, then we find a natural representation of ω Im T(ω) as where we used the fact that ω 2 δ(ω − ω i ) = ω 2 i δ(ω − ω i ).Now we can see that this oscillator representation is exactly the same as in the high-frequency case, except with the replacement ω 2 p → αω 2 i .Hence the representation theorem will be the same, but with the coefficients Importantly, note that the Fresnel coefficients are independent of k x and k y , and therefore after the inverse Fourier transform, we have E = tE inc , and α = χt.Similarly, the arguments expressing the fields with Fourier basis apply when we consider two parallel half-spaces separated by d 0 , but the transmission coefficient needs to be substituted by that of two interfaces Using p ω 2 at ω → 0, one can obtain t 2hs (ω = 0) = 2 ε 2 and α 2hs = χt 2hs = 2. Therefore the electrostatic T matrix for the bounding volume of two half-spaces is T(ω = 0) = α 2hs I where α 2hs = 2.

VI. DRUDE-LORENTZ REPRESENTATION OF THE T MATRIX
Here we discuss some of the mathematical aspects of our Drude-Lorentz oscillator representation in the main text.We started with three identities, From these three identities, we represented the imaginary parts of the scattering matrix via delta functions at a discrete set of points along the real line, from 0 to arbitrarily high frequencies: Then we claimed that this leads to a single unified representation, where the T i satisfy T i ≥ 0 and i T i ≤ I. First, we can simply validate that this expression is indeed consistent with the three identities above.The positivity and sum-rule conditions are clearly met.We next want to ensure that the real and imaginary parts are correct.In the discrete representation of Eq. ( 78), the resulting real part is given by Re Clearly this is consistent with Eq. (79) in the limit as γ → 0. (We will see another validation in a moment.)Is the imaginary part of Eq. ( 79) consistent with the delta function representation?One way to see this would be to separate the real and imaginary parts of Eq. ( 79), see that the linewidth of the imaginary part is going to zero but that its integral is fixed at 1.An alternative route is to use partial fractions to rewrite where "PV" denotes principal value, and to reach the second line we used the well-known Sokhotski-Plemelj identity [28] lim We can insert the expression of Eq. (82) into the expression for T(ω), Eq. (79), to find: where to simplify the last term we used the identities δ(ω i − ω)/ω i = δ(ω i − ω)/ω and exactly as required.Hence we have shown that the oscillator representation of Eq. ( 79) encodes all of the identities we derived.
For simplicity of exposition and intuition we used the basis of delta functions.One can expect arbitrarily high accuracy using this approach when working with objectives that satisfy two conditions: first, integration over bandwidths at least as large as the (infinitesimal) oscillator strengths γ, and second, an oscillator spacing well below any other characteristic x 5 = (0.93, 0.31)a, using the center of the ellipse as the origin.The incident field is a plane wave propagating along the y direction with the electric field polarized perpendicular to the plane.

VIII. DOMAIN MONOTONICITY: WHY SHARP TIPS DO NOT ENHANCE LDOS OR NFRHT
In this section, we explain why sharp tips, or any other geometric patterning, do not enhance the NFRHT bound.As discussed in the main text, for any scattering body, the only degrees of freedom are the scattering oscillator strengths T i .The oscillator strengths are positive semidefinite, and their key constraint is the sum rule of Eq. ( 64), which as discussed in Sec.V is given by T(ω = 0), the electrostatic T matrix, which we also call the "generalized polarizibility matrix" α.The question, then, is what structure leads to the largest generalized polarizability / electrostatic T matrix?In this section, we derive a "domain montonicity theorem" for the electrostatic T matrix, which shows that enclosing any domain with another will always "increase" (in a matrix-valued sense made precise below) the electrostatic T matrix.Then, among all structures with a minimum separation d, the largest electrostatic T matrix will be that of the domain of two planar half-spaces.
We prove the mathematics of this theorem below, and then discuss the physics: sharp tips enhance fields infinitely close to the tips themselves, but not back at the location of the source itself.The latter is the key quantity for local densities of states and related quantities such as NFRHT.
To prove domain monotonicity, we need to prove that quantities of the form x † Tx increase, for all x = 0, when the domain of the T matrix increases.We can interpret the multiplication of T with x as the polarization field induced by an "incident field" x, and then multiplication on the left by x takes the overlap of that incident field with the polarization that it induces.
Hence we will label our arbitrary vectors as e inc instead of x, for clarity in the mathematical relations to follow, though we impose no constraints on the "incident field" and indeed allow it to be an arbitrary vector.In computing the response to such a vector, however, we can use a few important physical consequences of electromagnetism.We will assume a scalar (and therefore reciprocal) electric susceptibility; the generalization to the arbitrary tensor case follows the same process as below in tandem with the modifications described in Sec.
VI of the SM of Ref. [41].In electrostatics the fields (and T matrix) can be chosen to be real-valued, so that we can consider the objective as x T Tx, without any conjugation.
When an incident field e inc interacts with a scatterer, for example through its T scattering matrix, there is a scattered field produced.The volume-equivalence principle allows us to write the scattered field in two forms.The first is as the convolution of the induced polarization field p with the background (vacuum) Green's function operator G 0 , i.e.G 0 p.
The second is a convolution of the full Green's function operator, G, with an overlap of the susceptibility χ with the incident e inc : e scat = Gχe inc .The total field is the incident field plus the scattered field, We are interested in the quantity F = e T inc Te inc = e T inc p, and how it changes when the domain changes.We will consider only continuous, increasing changes in susceptibility: ∆χ(x) ≥ 0 everywhere.Hence a variation in F can be written for any increases in the domain size or shape; since this is true for any vector e inc , then variations in the electrostatic T matrix must themselves be monotonic.This means that given a scatterer Ω 1 of any size and shape whose static T matrix is T (1) (ω = 0), any other scatterer Ω 2 whose volume encloses that of Ω 1 must have a T (2) (ω = 0) greater than T (1) (ω = 0), i.e.: T (2) (ω = 0) > T (1) (ω = 0), when the scatterer domain Ω 2 entirely encloses the scatter domain Ω 1 .Recalling from Eq. (64) that the electrostatic T matrix determines the constant in the low-frequency sum rule for the corresponding structure, the domain monotonicity theorem dictates that the NFRHT bound obtains a similar monotonicity: the bound only increases for structures with larger bounding volumes.For any given gap sizes of NFRHT, two planar half-spaces provide the largest T(ω = 0), which we use as the sum rule constant to derive the upper bound for NFRHT.Therefore, our bound applies to any structure, regardless of their size and shape, as they must satisfy the gap separation and therefore be contained by the two half-spaces.
Nanopatterned sharp tips, despite potentially providing local field enhancement near the tip, must have smaller upper bounds due to smaller bounding volumes, according to the monotonicity theorem.
To further illustrate why nano-patterning such as sharp-tip structure is inferior, we design a numerical experiment.The key objective function in the above analytical proof is the quantity e T inc p, i.e. the overlap of an incident field with the polarization field it induces.In the case of NFRHT, the incident field arises from point sources along a separating plane between the two bodies.An incident field emanating from a dipolar point source p 0 can be written e inc = G 0 p 0 , in which case the overlap can be written e T inc p = p T 0 G T 0 p = p T 0 G 0 p = p T 0 e scat , where e scat is the scattered field, emanating from the polarization fields, back at the location of the source on the separating plane.Moreover, the overlap measures the zerofrequency response proportional to the component of E scat (x 0 ) parallel to the polarization of the dipole source p(x 0 ).We should pinpoint the static E scat (x 0 ) as its magnitude is proportional to the electrostatic T matrix.
We simulated the static E scat (x 0 ) for a selection of representative structures to verify the monotonicity theorem.The simulations below uses the finite element method (FEM) of a commercial software, COMSOL Multiphysics.Point dipoles with unit amplitude are placed at the center of gaps of size 2d between identical wedges or half-spaces.The maximum size of triangular meshes are 0.01d in the central region of the air gap and around the tip, and 0.1d for the rest of region away from the dipole.E scat (x 0 ) is computed for vertical "⊥" (perpendicular to the separating plane) and horizontal " " (parallel to the separating plane) polarized p(x 0 ) for 8 different structures with inner angles from 0.125π to π, essentially from a pair of sharp and pointed wedges to two planar half-spaces.Due to the intrinsic difficulty in measuring E scat (x 0 ) which is at the point of dipole source location, it is instead obtained by interpolation of E scat (x) at locations near the dipole source.The numerical data points for β = π (two half-spaces) agree well with analytical prediction through method of images which gives E ⊥ scat (x 0 ) = π 24 0 d 2 for vertical polarization and E scat (x 0 ) = π 48 0 d 2 for horizontal polarization.
The results are shown in the figure Fig. 1(c) as dark red lines with vertical markers representing vertical polarization and square markers representing horizontal polarization.
Additionally, scattered field at a point close to the tip (0.1d away from the tip along the line connecting two tips) are shown as grey lines.As wedge angle increases, the vertically polarized scattered field decreases near the tip and increases at source location, the horizontally polarized scattered field increases both near the tip and at source location.In other words, although there is field enhancement near the tip with small wedge structures for the vertical polarization, which is the reasoning behind nano-patterning to improve near-field interaction and increase NFRHT efficiency, the quantity that is much more relevant, the scattered field right at the source point E scat (x 0 ), is decreased with nano-patterning.The largest E scat (x 0 ) is achieved with the largest wedge angle β = π, namely two unpatterned half-spaces.As the increase of wedge angle essentially increases scatterer domain size, this numerical result provides quantitative confirmation that static E scat (x 0 ) is increasing as domain size monotonically increases.

IX. DERIVATION OF THE NFRHT BOUND
To investigate radiative heat transfer from object 1 (bottom) to object 2 (top), we first break down the problem to power integrations at every frequency.The power flowing in the positive z direction across the middle separating plane (perpendicular to z) between the two objects is: where the superscripts denote the polarizations of the current sources in the bottom object, whose amplitudes are dictated by the fluctuation-dissipation theorem: where the susceptibility of object 1 is χ Then the field correlations are found by rewriting the fields in terms of the combining Green's function and thermal source correlations Our bound will not distinguish between the x and y directions (which are symmetric in the bounding domain, even though they of course are not for many allowable patterns), in which case the upper bounds on either of the two terms in power integration in Eq. ( 99) are identical: Hence the maximum flux S z (ω, T ) equals the maximum of the function By Lorentz reciprocity, we can equate the fields at r s produced by sources at r v with fields at r v produced by sources at r s .This will be helpful to our bound formalism which would essentially be a volume integration of field operators.In light of the correlations for currents sources inside the volume, Eq. (100), we can define the correlations for current sources on the middle flux plane as The amplitude ω is chosen so that E Jx inc (r v ) * E My inc (r v ) is independent of frequency, which will be important later.Now we can rewrite Eq. ( 101) as and Eq.(103) as Suppressing all position arguments of the relevant fields in a vector notation, where O is a diagonal matrix whose only nonzero diagonal entries that equal 1 correspond to the bottom body (where the thermal sources of interest come from) and the rest of the diagonal entries that equal 0 correspond to the top body.Both T matrix and E inc vectors are defined across both the top and bottom bodies.This O matrix implies that the outer spatial integral is over the bottom body, instead of both bodies.The matrix E = Re e My inc e Jx inc † is a rank-2 matrix that can be decomposed into one positive eigenvalue term and one negative eigenvalue term: One can now see that our choice of source amplitudes in Eq. ( 106) leads to frequencyindependent incident-field eigenvalues.
To bound the expression of Eq. ( 112), we will relax it in a few ways.(Interestingly, intensive numerical optimizations using manifold-optimization techniques [42,43] directly on Eq. ( 112) lead to the same upper limits that we derive below, suggesting that these "relaxations" are minimal and do not loosen the analysis given the information, such as sum rules, that we use.)First, the E matrix defined by the two renormalized incident fields has one positive and one negative eigenvalue, per Eq. ( 113).Physically, we can interpret the negative sign of the second eigenvalue via the power expression of Eq. ( 112) containing E, as the difference in powers absorbed for the two renormalized incident fields.This is of course bounded above by the absorption of only the first incident field, dropping the subtracted term, leaving only the contribution of the single positive eigenvalue of E. Thus we have: My inc e Jx inc † = T † OTE ≤ T † OT λ 1 q 1 q † 1 . (117) Next, we note that O indicates absorption only in the lower body; of course this quantity is bounded above by the total absorption in both bodies.Mathematically, the O matrix is bounded above by the identity matrix, I. Finally, the absorption in both bodies is less than the net extinction of the two bodies (their far-field scattered powers are positive, and essentially zero in the near-field case, so that this relaxation is negligible), giving: Tr T † OT λ 1 q 1 q † 1 = λ 1 Tr q † 1 T † OTq 1 (118) which is quadratic in T. We can use "optical theorem" constraint to bound this quadratic absorption-like quantity with a linear extinction-like quantity.The idea is that absorption must be smaller than extinction: For our problem, scattering power is negligible for the entire system, as the near-field heat exchange dominates over the scattering to outside the system.In other words, P abs ≈ P ext .
Thus we can write without introducing much relaxation.
We can now rewrite Eq. (108) as Surprisingly, the various transformations to this point have removed all explicit dependencies on material susceptibility χ 1,2 (ω), with the only implicit dependence embedded in ImT.We will now focus on the upper bound for HTC, and the upper bound for RHT can be found by taking similar steps.To switch from the RHT to HTC bound computation, we just need to take the temperature derivative of the last expression to get where, for simplicity of notation, we can write x = ω k B T so that ∂Θ(ω, T ) ∂T = k B x 2 e x (e x − 1) 2 . (125) The oscillator representation of the T matrix, per Eq. ( 67), is: where T i are the positive-definite oscillator amplitudes and f i (ω) = Now we will integrate with respect to frequency, and it will then be obvious that a single T i oscillator at the optimal frequency maximizes the total HTC.For each T i oscillator term, the frequency integral is where x i = ω i k B T .The maximum Ω i occurs for x i = 2.57, which corresponds to an optimal oscillator frequency of which is exactly the near-field Wien frequency that we found from HTC optimization for planar, unpatterned geometries [44]!The maximum Ω i at this frequency is Apart from the frequency integral, everything else is in the trace, where T i is constrained by T i > 0 and i T i = I.The answer to achieving maximum ´∞ 0 dω S HTC (ω, T ) is obviously allocating all possible T i to the optimal oscillator frequency.It is easy to find the maximal value of the trace term, as Max Tr T i q 1 q † 1 = 1 (134) when T i = I or q 1 q † 1 .Then the bound will be the product of the frequency integral term and the trace term, both of this optimal oscillator: HTC ≤ 0.21 where β = 3.8 × 10 5 Wnm 2 /m 2 /K 2 .For T = 300 K and d = 10 nm, HTC ≤ 1.1 × 10 6 W/m 2 /K, which is 5X the optimal planar performance.Hence this theoretical framework offers a close prediction to the best known designs, it predicts the optimal resonance frequency where the oscillator-strength should be concentrated, and it explain why previous material-dependent predictions were incorrect.

FIG. 1 .
FIG.1.Broadband scattering of incident plane waves from an elliptical cylinder.A schematic depiction is provided in (a).Whereas the scattered fields (at points 1-5 of part a) exhibit seemingly random variations, as depicted in (b) for a single incident angle, the T matrix has far more spectral structure.Part (c) shows the real and imaginary parts of the diagonal elements of the T matrix, which look reminiscent of the structured oscillations of a material susceptibility.Yet this cylinder has a constant susceptibility χ = 4, as shown in (f).The undulations in the T matrix elements are not material oscillators, but a new type of matrix-valued, nonlocal scattering oscillator.Parts (d) and (e) show the matrix-valued scattering-oscillator coefficients Ti, and their eigenvalues, respectively, showcasing the extent to which passivity and causality restrict the possible form of the scattering T matrix.

FIG. 2 .
FIG. 2. (a)NFRHT between two closely separated bodies.(b) Heat-transfer coefficients of planar bodies comprising increasingly high-performance planar-body geometries (filled circles), in comparison with the best previous theoretical bounds[20, 21]  (open squares and triangles).The previous bounds diverged for some materials, while showing enormous gaps (> 750X) for others, and their trendline seems to decrease left-to-right, whereas planar-body performance increases.In black is the new theoretical bound offered by our T-matrix representation, very close to the best possible planar bodies.Panels (c,d) confirm the predictions of our spectral representation, showing that the state-of-the-art performance could be achieved precisely at the optimal frequencies identified by our approach.
δF = e T inc δp.(92) The change in the polarization field p = χe across the scatterers has two contributions: a change in the permittity multiplied by the field, and the permittivity multiplied by the change in field: δp = (δχ) e + χ(δe).(93) The change in field equals the convolution of the Green's function with the newly induced current, which is (δχ)e.Hence we have δp = (δχ) e + χG(δχ)e, (94) and, δF = e T inc [(δχ) e + χG(δχ)e] = e T inc (I + χG) δχe = e T δχe ≥ 0, (95) where to proceed from the second to the third line we used the previous scattering-matrix equality, Eq. (91), in tandem with the assumption of scalar (reciprocal) materials.The fourth line follows from the third because the fields are real and the change in susceptibility is positive, so e T δχe = ˆV δχ(x)|E(x)| 2 ≥ 0. (96) Hence we have shown that e † inc δTe inc ≥ 0 (97) FIG. 1. Domain monotonicity of static E scat (x 0 ).(a) Structures patterned with sharp cones v.s.(b) planar structures without any patterning.(c) Scattered fields from a dipole placed in the middle of a pair of conducting wedges of different inner angles from 0.125π to π. Square markers signify horizontal dipole orientation, and vertical markers signify vertical dipole polarization.Despite there is monotonic increase of scattered field near the tip of the wedges with decreasing inner angles for the vertical polarization according to the top grey line, the two red lines, directly proportional to static E scat (x 0 ), show monotonic increase of scattered field back at the source location with increasing inner angles.Largest static E scat (x 0 ) is when β = π, namely the structure is two planar half-spaces.

− 1 .
and Θ(ω, T ) is the Planck distribution, Θ(ω, T ) = The subscripts on the position vectors r s and r v indicate whether their domain is the surface (separating flux plane) or volume (of the emitter/absorber).

αω 2 i ω 2 i
−ω 2 −iωγ are the oscillator lineshapes.The imaginary part of T matrix reads