Fundamental limits to radiative heat transfer: theory

TL;DR

This paper establishes fundamental, shape-independent limits on near-field radiative heat transfer, considering material and geometric constraints, and reveals that maximizing scattering rather than absorption leads to optimal heat transfer, impacting thermophotovoltaics and nanoscale cooling.

Contribution

It derives universal bounds on radiative heat transfer that incorporate realistic constraints and challenges previous assumptions about approaching Landauer limits.

Findings

Fundamental limits are shape-independent and include dissipation and size constraints.

Maximizing heat transfer correlates with maximizing scattering, not absorption.

Existing Landauer-based limits are overly optimistic compared to these new bounds.

Abstract

Near-field radiative heat transfer between bodies at the nanoscale can surpass blackbody limits on thermal radiation by orders of magnitude due to contributions from evanescent electromagnetic fields, which carry no energy to the far-field. Thus far, principles guiding explorations of larger heat transfer beyond planar structures have assumed utility in surface nanostructuring, which can enhance the density of states, and further assumed that such design paradigms can approach Landauer limits, in analogy to conduction. We derive fundamental shape-independent limits to radiative heat transfer, applicable in near- through far-field regimes, that incorporate material and geometric constraints such as intrinsic dissipation and finite object sizes, and show that these preclude reaching the Landauer limits in all but a few restrictive scenarios. Additionally, we show that the interplay of material response and electromagnetic scattering among proximate bodies means that bodies which maximize radiative heat transfer actually maximize scattering rather than absorption. Finally, we compare our new bounds to existing Landauer limits, as well as limits involving bodies maximizing far-field absorption, and show that these lead to overly optimistic predictions. Our results have ramifications for the ultimate performance of thermophotovoltaics and nanoscale cooling, as well as related incandescent and luminescent devices.

Equations

$$P={\int }_0^{\infty }[\Pi (\omega ,T_{\mathrm{B}})-\Pi (\omega ,T_{\mathrm{A}})]\Phi (\omega )\text{leavevmode}\mathrm{d}\omega ,$$

$$\begin{gathered} \Phi =\frac{2}{\pi }\mathrm{Tr}[{\mathbb{T}}_{\mathrm{B}}^{\star }({\mathbb{I}}_{\mathrm{B}}-{\mathbb{G}}_{\mathrm{BA}}^{\mathrm{vac}\star }{\mathbb{T}}_{\mathrm{A}}^{\star }{\mathbb{G}}_{\mathrm{AB}}^{\mathrm{vac}\star }{\mathbb{T}}_{\mathrm{B}}^{\star }{)}^{-1}\mathrm{Im}({\mathbb{V}}_{\mathrm{B}}^{-1\star })\times \\ {\mathbb{T}}_{\mathrm{B}}({\mathbb{I}}_{\mathrm{B}}-{\mathbb{G}}_{\mathrm{BA}}^{\mathrm{vac}}{\mathbb{T}}_{\mathrm{A}}{\mathbb{G}}_{\mathrm{AB}}^{\mathrm{vac}}{\mathbb{T}}_{\mathrm{B}}{)}^{-1}\times \\ {\mathbb{G}}_{\mathrm{BA}}^{\mathrm{vac}}{\mathbb{T}}_{\mathrm{A}}\mathrm{Im}({\mathbb{V}}_{\mathrm{A}}^{-1\star }){\mathbb{T}}_{\mathrm{A}}^{\star }{\mathbb{G}}_{\mathrm{AB}}^{\mathrm{vac}\star }], \end{gathered}$$

$$\Phi =\frac{2}{\pi }\mathrm{Tr}[{\mathbb{Y}}_{\mathrm{B}}^{\star }\mathrm{Im}({\mathbb{V}}_{\mathrm{B}}^{-1\star }){\mathbb{Y}}_{\mathrm{B}}{\mathbb{G}}_{\mathrm{BA}}^{\mathrm{vac}}{\mathbb{T}}_{\mathrm{A}}\mathrm{Im}({\mathbb{V}}_{\mathrm{A}}^{-1\star }){\mathbb{T}}_{\mathrm{A}}^{\star }{\mathbb{G}}_{\mathrm{AB}}^{\mathrm{vac}\star }],$$

$$\Phi =\frac{2}{\pi {\zeta }_{\mathrm{A}}{\zeta }_{\mathrm{B}}}{\Vert {\mathbb{Y}}_{\mathrm{B}}{\mathbb{G}}_{\mathrm{BA}}^{\mathrm{vac}}{\mathbb{T}}_{\mathrm{A}}\Vert }_{\mathrm{F}}^2$$

$$\Phi =\frac{2}{\pi }{\Vert \mathbb{Q}\Vert }_{\mathrm{F}}^2$$

$${\mathbb{Y}}_{\mathrm{B}}$$

$${\mathbb{G}}_{\mathrm{BA}}^{\mathrm{vac}}$$

$${\mathbb{T}}_{\mathrm{A}}$$

$$\Phi =\frac{2}{\pi {\zeta }_{\mathrm{A}}{\zeta }_{\mathrm{B}}}\sum _i(y_i{g}_i{\tau }_i{)}^2$$

$$\Phi =\frac{2{\zeta }_{\mathrm{B}}}{\pi {\zeta }_{\mathrm{A}}}\sum _i\frac{(g_i{\tau }_i{)}^2}{{\left[1+\frac{{\zeta }_{\mathrm{B}}}{{\zeta }_{\mathrm{A}}}(g_i{\tau }_i{)}^2\right]}^2}$$

$$\Phi_{\mathrm{opt}}=\sum_{i}\Big{[}\frac{1}{2\pi}\Theta(\zeta_{\mathrm{A}}\zeta_{\mathrm{B}}g_{i}^{2}-1)\\ +\frac{2}{\pi}\frac{\zeta_{\mathrm{A}}\zeta_{\mathrm{B}}g_{i}^{2}}{(1+\zeta_{\mathrm{A}}\zeta_{\mathrm{B}}g_{i}^{2})^{2}}\Theta(1-\zeta_{\mathrm{A}}\zeta_{\mathrm{B}}g_{i}^{2})\Big{]}$$

$${\Phi }_{\mathrm{L}}=\sum _i\frac{1}{2\pi },$$

$${\Phi }_{\mathrm{sc}}=\sum _i\frac{2}{\pi }\frac{{\zeta }_{\mathrm{A}}{\zeta }_{\mathrm{B}}{g}_i^2}{(1+{\zeta }_{\mathrm{A}}{\zeta }_{\mathrm{B}}{g}_i^2{)}^2}$$

$${\Phi }_{\mathrm{Born}}=\sum _i\frac{2}{\pi }{\zeta }_{\mathrm{A}}{\zeta }_{\mathrm{B}}{g}_i^2$$

$${\Phi }_{\mathrm{sc}}\le {\Phi }_{\mathrm{opt}}\le {\Phi }_{\mathrm{Born}},{\Phi }_{\mathrm{L}}.$$

$$\mathbb{A}$$

$$\mathbb{B}$$

$$\mathrm{Tr}[\mathbb{AB}]=\sum _{i=1}^N\sum _{j=1}^N{\sigma }_i{p}_{ij}{\tau }_j{q}_{ji}$$

$$|\mathbf{E}⟩$$

$$|\mathbf{P}⟩$$

$$|\mathbf{E}⟩$$

$${\mathbb{G}}^{\mathrm{vac}}$$

$${\mathbb{T}}^{-1}=\left[\begin{array}{cc}{\mathbb{T}}_{\mathrm{A}}^{-1}& -{\mathbb{G}}_{\mathrm{AB}}^{\mathrm{vac}}\\ -{\mathbb{G}}_{\mathrm{BA}}^{\mathrm{vac}}& {\mathbb{T}}_{\mathrm{B}}^{-1}\end{array}\right]$$

$$\begin{array}{rl}|\mathbf{E}⟩& ={\mathbb{G}}^{\mathrm{vac}}\mathbb{T}{\mathbb{V}}^{-1}|{\mathbf{P}}^{(0)}⟩\\ |\mathbf{P}⟩& =\mathbb{T}{\mathbb{V}}^{-1}|{\mathbf{P}}^{(0)}⟩\end{array}$$

$$\begin{array}{rl}{\mathbb{I}}_{\mathrm{A}}& =\left[\begin{array}{cc}{\mathbb{I}}_{\mathrm{A}}& 0\\ 0& 0\end{array}\right],{\mathbb{I}}_{\mathrm{B}}=\left[\begin{array}{cc}0& 0\\ 0& {\mathbb{I}}_{\mathrm{B}}\end{array}\right],\end{array}$$

$$\Phi =\frac{1}{2}\mathrm{Re}(⟨⟨{\mathbb{I}}_{\mathrm{B}}\mathbf{E},{\mathbb{I}}_{\mathrm{B}}\mathbf{J}⟩⟩)$$

$$\Phi =-\frac{\omega }{2}⟨⟨{\mathbf{P}}_{\mathrm{A}}^{(0)},{\mathbb{I}}_{\mathrm{A}}{\mathbb{V}}^{-1†}{\mathbb{T}}^{†}\text{leavevmode}\mathrm{Asym}({\mathbb{I}}_{\mathrm{B}}{\mathbb{G}}^{\mathrm{vac}})\mathbb{T}{\mathbb{V}}^{-1}{\mathbb{I}}_{\mathrm{A}}{\mathbf{P}}_{\mathrm{A}}^{(0)}⟩⟩$$

$$⟨|{\mathbf{P}}_{\mathrm{A}}^{(0)}⟩⟨{\mathbf{P}}_{\mathrm{A}}^{(0)}|⟩=\frac{4}{\pi \omega }\mathrm{Asym}({\mathbb{V}}_{\mathrm{A}})$$

$$\Phi =-\frac{2}{\pi }\mathrm{Tr}(\mathrm{Asym}({\mathbb{V}}_{\mathrm{A}}^{-1†}){\mathbb{I}}_{\mathrm{A}}{\mathbb{T}}^{†}\text{leavevmode}\mathrm{Asym}({\mathbb{I}}_{\mathrm{B}}{\mathbb{G}}^{\mathrm{vac}})\mathbb{T}{\mathbb{I}}_{\mathrm{A}})$$

$$\mathrm{Asym}({\mathbb{I}}_{\mathrm{B}}{\mathbb{G}}^{\mathrm{vac}})$$

$$\mathbb{T}{\mathbb{I}}_{\mathrm{A}}$$

$${\mathbb{I}}_{\mathrm{A}}{\mathbb{T}}^{†}$$

$$\Phi =\frac{2}{\pi }\mathrm{Tr}[\mathrm{Asym}({\mathbb{V}}_{\mathrm{A}}){\mathbb{G}}_{\mathrm{BA}}^{†}\mathrm{Asym}({\mathbb{V}}_{\mathrm{B}}){\mathbb{G}}_{\mathrm{BA}}]$$

$$\begin{gathered} \Phi (\omega )=\frac{2}{\pi }\mathrm{Tr}[\mathrm{Asym}({\mathbb{V}}_{\mathrm{A}}){\mathbb{V}}_{\mathrm{A}}^{-1†}{\mathbb{T}}_{\mathrm{A}}^{†}({\mathbb{S}}_{\mathrm{B}}^{\mathrm{A}}{\mathbb{G}}_{\mathrm{BA}}^{\mathrm{vac}}{)}^{†}{\mathbb{T}}_{\mathrm{B}}^{†}\times \\ {\mathbb{V}}_{\mathrm{B}}^{-1†}\mathrm{Asym}({\mathbb{V}}_{\mathrm{B}}){\mathbb{V}}_{\mathrm{B}}^{-1}{\mathbb{T}}_{\mathrm{B}}{\mathbb{S}}_{\mathrm{B}}^{\mathrm{A}}{\mathbb{G}}_{\mathrm{BA}}^{\mathrm{vac}}{\mathbb{T}}_{\mathrm{A}}{\mathbb{V}}_{\mathrm{A}}^{-1}] \end{gathered}$$

Full text

The concept of a blackbody, derived from electromagnetic reciprocity (or detailed balance), has provided a benchmark of the largest emission rates that can be achieved by a heated object: through nanoscale texturing, gray objects can be designed in myriad ways to mimic the response of a blackbody at selective wavelengths Jin et al. (2019); Costantini et al. (2015); Rodriguez et al. (2011), with implications for a variety of technologies, including high-efficiency solar cells, selective emitters, and thermal sensors Asano et al. (2016). Over the past few decades, much effort has gone toward understanding analogous limits to enhancements of near-field radiative heat transfer (RHT) Pendry (1999); Biehs et al. (2010); Ben-Abdallah and Joulain (2010); Miller et al. (2015), supported by a rich and growing number of experimental Shen et al. (2009); St-Gelais et al. (2014); Cui et al. (2017); Kloppstech et al. (2017) and theoretical Luo et al. (2004); Otey et al. (2014); Polimeridis et al. (2015); Rodriguez et al. (2013); Krüger et al. (2012) investigations, and motivated by potential applications to thermophotovoltaics Lenert et al. (2014); Karalis and Joannopoulos (2016), nanoscale cooling Guha et al. (2012), and thermal microscopy Boudreau et al. (1997); Jones and Raschke (2012). A key principle underlying further near-field RHT enhancements is the use of materials supporting bound (plasmon and phonon) polaritons in the infrared, where the Planck distribution peaks at typical temperatures probed in experiments. This leads to strong subwavelength responses tied to corresponding enhancements in the density of states Volokitin and Persson (2001); Domingues et al. (2005); Volokitin and Persson (2007); Song et al. (2015); consequently, the amplified near-field RHT spectrum exhibits a narrow lineshape, justifying focus on selective wavelengths. However, while the properties of such polaritons, particularly their resonance frequencies, associated densities of states, and scattering characteristics can be modified through nanoscale texturing, only recently have computational methods Reid et al. (2013); Otey et al. (2014); Rodriguez et al. (2013); Polimeridis et al. (2015) arisen to model RHT between bodies of arbitrary shapes beyond those with high symmetry Krüger et al. (2012); Pérez-Madrid et al. (2008); Messina et al. (2013). Furthermore, the challenge of gaining simultaneous control over the scattering properties of large numbers of contributing surface waves has generally precluded general upper bounds on RHT.

RHT between two bodies A and B in vacuum is given as

[TABLE]

in terms of their local temperatures $T_{\mathrm{A}}$ and $T_{\mathrm{B}}$, entering the Planck function $\Pi(\omega,T)=\hbar\omega/[\mathrm{exp}(\hbar\omega/(k_{\mathrm{B}}T))-1]$ (and it has been assumed, without loss of generality, that $T_{\mathrm{B}}>T_{\mathrm{A}}$ so $P>0$), and the spectral function $\Phi(\omega)$, which can be enhanced by changing material and geometric properties through the creation of resonances and changes in the electromagnetic density of states. In particular, nanostructuring metallic surfaces or polar dielectrics makes it possible to shift resonant frequencies from the visible or far infrared into the infrared, such that the peak of the spectrum $\Phi$ may coincide with the peak of the Planck distribution near room temperature. It remains an open question, however, to what extent the peak value of $\Phi$ itself may be enhanced through appropriate geometric and material choices, as well as what such optimal structures should be.

Previous attempts at deriving bounds on RHT have primarily focused on extended media Pendry (1999); Bimonte (2009); Biehs et al. (2010); Ben-Abdallah and Joulain (2010), showing that at least for translationally invariant structures, $\Phi$ can be expressed as the trace of a “transmission” matrix whose singular values (corresponding to evanescent Fourier modes) each contribute a finite flux, bounded above by a Landauer limit in analogy with conduction Datta (1995); Klöckner et al. (2016). Aside from being restricted to planar geometries, these bounds turn out to be either pessimistic Ben-Abdallah and Joulain (2010), ignoring the large densities of states that can arise in nanostructured and low-loss materials, or too optimistic Pendry (1999); Biehs et al. (2010), ignoring any constraints imposed by Maxwell’s equations and assuming instead that all such Fourier modes, up to an unrealistic cutoff on the order of the atomic scale, can saturate the flux Pendry (1999). From a design perspective, Landauer limits present a hurdle as they rely on ad-hoc estimates of the number and relative contribution of radiative modes/channels, which depend on specific material and geometric features. More recent works have derived complementary material limits on electromagnetic absorption in subwavelength regimes Miller et al. (2016), showing that absorbed power in a medium of susceptibility $\chi$ increases in proportion to an “inverse resistivity” figure of merit, $|\chi|^{2}/\operatorname{Im}\chi$, in principle diverging with increasing indices of refraction and decreasing dissipation. Saturation of these bounds for a subwavelength absorber in the quasistatic regime can generally be achieved through the strong polarization currents arising in resonant media supporting surface plasmon or phonon polaritons. These arguments have been extended to near-field RHT Miller et al. (2015) by exploiting energy conservation and reciprocity, finding the upper bound of $\Phi$ at a polariton resonance to scale quadratically with $|\chi|^{2}/\operatorname{Im}\chi$, corresponding to enhanced absorption and emission in both objects. While near-field RHT between dipolar objects can attain these bounds in a dilute limit, such a universal scaling has yet to be observed in large-area structures. This naïvely suggests room for improvement in $\Phi$ through nanostructuring via enhancements in the density of states or equivalently, via saturation of modal contributions, yet trial-and-error explorations and optimization procedures Jin et al. (2017); Fernández-Hurtado et al. (2017) have failed to produce nanostructured geometries that bridge this gap, leading to the alternative possibility that existing bounds are too loose.

In this paper, we derive new algebraic bounds on RHT, valid in the near-, mid-, and far-field regimes, through analysis of the singular value decompositions of relevant response quantities. In contrast to prior limits, our bounds incorporate constraints imposed by material losses and multiple scattering, and are therefore tighter; moreover, they are formulated to be independent of object shapes while simultaneously accounting for finite size effects. In particular, every channel of energy transmission is shown to be generally prohibited from saturating its Landauer limit, in contrast to predictions based on modal decompositions Pendry (1999); Bimonte (2009); Biehs et al. (2010); Ben-Abdallah and Joulain (2010) that neglect material properties and are most applicable in the ray optics regime. Furthermore, the growth of RHT with decreasing material dissipation is shown to be strongly limited by radiative losses, in contrast to predictions based on energy-conservation limits to material response Miller et al. (2015) that neglect finite-size scattering effects and are thus most applicable in the quasistatic regime. Finally, we discuss the difficulties in optimizing each individual energy transmission channel and show how more restricted bounds may be obtained by approximating each body to maximally absorb all incident fields in isolation. In upcoming papers closely related to this one, we apply these bounds to various scenarios of interest, providing predictions of the maximum RHT achievable in compact and extended geometries Venkataram et al. (2019), and deriving related bounds on far-field thermal emission from single bodies in isolation Molesky et al. (2019).

I Heat transfer definitions

For two bodies A and B in vacuum [Fig. 1], the spectral function $\Phi$ appearing in (1) represents the average power absorbed in B due to fluctuating current sources in A , depicted in and is reciprocal (invariant under interchange of A and B). Using formal operator notation, this average absorbed power can be written in terms of the susceptibilities $\mathbb{V}_{p}$, the vacuum Maxwell Green’s function $\mathbb{G}^{\mathrm{vac}}_{pq}$, and scattering T-operators $\mathbb{T}_{p}$, for $p,q\in\{\mathrm{A},\mathrm{B}\}$. Each susceptibility is written as $\mathbb{V}_{p}=\chi_{p}\mathbb{I}_{p}$, where each $\chi_{p}$ is assumed to be homogeneous, local, and isotropic. The vacuum Maxwell Green’s function $\mathbb{G}^{\mathrm{vac}}$ solves $[(c/\omega)^{2}\nabla\times(\nabla\times)-\mathbb{I}]\mathbb{G}^{\mathrm{vac}}=\mathbb{I}$ in all space, and its blocks are denoted as $\mathbb{G}^{\mathrm{vac}}_{pq}$ for sources in body $q$ propagating fields to body $p$. Finally, the T-operators $\mathbb{T}_{p}=(\mathbb{V}_{p}^{-1}-\mathbb{G}^{\mathrm{vac}}_{pp})^{-1}$ represent the total induced polarization moment in body $p$ due to a localized field of unit magnitude incident upon it. All of these quantities are reciprocal, so they are equal to their unconjugated transposes in position space: $\mathbb{V}_{p}=\mathbb{V}_{p}^{\top}$, $\mathbb{T}_{p}=\mathbb{T}_{p}^{\top}$, and $\mathbb{G}^{\mathrm{vac}}_{pq}=(\mathbb{G}^{\mathrm{vac}}_{qp})^{\top}$. This means that Hermitian conjugation is equivalent to complex conjugation: $\mathbb{V}_{p}^{\dagger}=\mathbb{V}_{p}^{\star}$, $\mathbb{T}_{p}^{\dagger}=\mathbb{T}_{p}^{\star}$, and $(\mathbb{G}^{\mathrm{vac}}_{pq})^{\dagger}=\mathbb{G}_{qp}^{\mathrm{vac}\star}$. Additionally, all of these quantities depend on frequency $\omega$, though this will be suppressed for the sake of notational brevity.

Given these definitions and relations (see appendices for more details), the RHT spectrum can be written as Krüger et al. (2012)

[TABLE]

where $\operatorname{Im}(\mathbb{A})=(\mathbb{A}-\mathbb{A}^{\star})/(2\mathrm{i})$ and $\operatorname{Asym}(\mathbb{A})=(\mathbb{A}-\mathbb{A}^{\dagger})/(2\mathrm{i})$ for any operator $\mathbb{A}$; if $\mathbb{A}$ is reciprocal, then $\operatorname{Asym}(\mathbb{A})=\operatorname{Im}(\mathbb{A})$. This expression is manifestly reciprocal in A and B, and treats the T-operators of A and B on an equal footing, linked only by the Green’s function $\mathbb{G}^{\mathrm{vac}}_{\mathrm{BA}}$ propagating fields in vacuum from one body to the other. However, it is possible to write this spectrum more suggestively in terms of operator combinations that hide this reciprocity in order to more strongly link this expression to absorbed and emitted powers. In particular, $\Phi$ may be rewritten as

[TABLE]

in terms of the reciprocal operator $\mathbb{Y}_{\mathrm{B}}=\mathbb{T}_{\mathrm{B}}\mathbb{S}^{\mathrm{A}}_{\mathrm{B}}$, which is in turn written in terms of the scattering operator $\mathbb{S}^{\mathrm{A}}_{\mathrm{B}}=(\mathbb{I}_{\mathrm{B}}-\mathbb{G}^{\mathrm{vac}}_{\mathrm{BA}}\mathbb{T}_{\mathrm{A}}\mathbb{G}^{\mathrm{vac}}_{\mathrm{AB}}\mathbb{T}_{\mathrm{B}})^{-1}$: essentially, $\mathbb{Y}_{\mathrm{B}}$ is a new “dressed T-operator” describing absorption and scattering in B in the presence of A, just as the bare T-operators $\mathbb{T}_{p}$ describe absorption and scattering from each body $p\in\{\mathrm{A},\mathrm{B}\}$ in isolation.

Having assumed the susceptibility in each body $p\in\{\mathrm{A},\mathrm{B}\}$ to be homogeneous, uniform, and isotropic, makes it possible to identify $\operatorname{Im}(\mathbb{V}_{p}^{-1\star})=\frac{\operatorname{Im}(\chi_{p})}{|\chi_{p}|^{2}}\mathbb{I}_{p}$. For convenience, we denote $\zeta_{p}=\frac{|\chi_{p}|^{2}}{\operatorname{Im}(\chi_{p})}$ as the material response factor. Using this, we write the RHT spectrum as

[TABLE]

where $\left\lVert\mathbb{A}\right\rVert_{\mathrm{F}}^{2}=\operatorname{Tr}(\mathbb{A}^{\dagger}\mathbb{A})$ denotes the Frobenius norm for any operator $\mathbb{A}$.

As we show in the appendices, the RHT spectrum may alternatively be written as

[TABLE]

in terms of the transmission operator $\mathbb{Q}=\operatorname{Im}(\mathbb{V}_{\mathrm{B}})^{1/2}\mathbb{G}_{\mathrm{BA}}\operatorname{Im}(\mathbb{V}_{\mathrm{A}})^{1/2}$, which in turn depends on the total Green’s function $\mathbb{G}_{\mathrm{BA}}=\mathbb{V}_{\mathrm{B}}^{-1}\mathbb{T}_{\mathrm{B}}\mathbb{S}^{\mathrm{A}}_{\mathrm{B}}\mathbb{G}^{\mathrm{vac}}_{\mathrm{BA}}\mathbb{T}_{\mathrm{A}}\mathbb{V}_{\mathrm{A}}^{-1}$ connecting dipole sources in body A to total fields in body B and accounting for multiple scattering to all orders within and between both bodies. This obeys a Landauer limit, as the singular values of $\mathbb{Q}^{\dagger}\mathbb{Q}$ do not exceed $1/4$, so including the prefactor $2/\pi$, the contribution of each mode/channel in the trace expression to $\Phi$ does not exceed $\frac{1}{2\pi}$. We will shortly explain the conditions under which the Landauer bounds for each of these contributions may be saturated.

II Singular value bounds

Each of the operators $\mathbb{Y}_{\mathrm{B}}$, $\mathbb{G}^{\mathrm{vac}}_{\mathrm{BA}}$, and $\mathbb{T}_{\mathrm{A}}$ may be decomposed into their respective singular values and vectors. In general, the singular vectors of each of these operators will not bear any relation to one another. However, we prove in the appendices that an upper bound to $\Phi$ can be found by ensuring that the right singular vectors of $\mathbb{Y}_{\mathrm{B}}$ and the left singular vectors of $\mathbb{T}_{\mathrm{A}}$ are respectively the same as the left and right singular vectors of $\mathbb{G}^{\mathrm{vac}}_{\mathrm{BA}}$. Intuitively, these conditions can be interpreted as representing perfect coupling and propagation of fields due to sources in A which are perfectly absorbed in body B (while accounting for multiple scattering from the presence of body A). Formally, we write the statements

[TABLE]

where as a reminder, the singular values $\tau_{i}$, $g_{i}$, and $y_{i}$ are nonnegative real numbers, and the vectors $|\mathbf{a}_{i}\rangle$ and $|\mathbf{b}_{i}\rangle$ are orthonormal among themselves so $\langle\mathbf{a}_{i},\mathbf{a}_{j}\rangle=\delta_{ij}$ and $\langle\mathbf{b}_{i},\mathbf{b}_{j}\rangle=\delta_{ij}$. Our explicit conjugation of $|\mathbf{b}_{i}^{\star}\rangle$ in the definition of $\mathbb{Y}_{\mathrm{B}}$ and of $\langle\mathbf{a}_{i}^{\star}|$ in the definition of $\mathbb{T}_{\mathrm{A}}$ ensure that they still respect reciprocity. This allows us to write a bound on the RHT spectrum as

[TABLE]

which is maximized when $y_{i}$ and $\tau_{i}$ attain their maximum values as functions of $g_{i}$; the singular values $g_{i}$ depend only on the geometry and topology of the system as they enter $\mathbb{G}^{\mathrm{vac}}_{\mathrm{BA}}$, not the material properties, so we take those as fixed.

As we show in the appendices, nonnegative far-field scattering from the combined system of bodies A and B together implies that the operator $\operatorname{Im}(\mathbb{Y}_{\mathrm{B}})-\mathbb{Y}_{\mathrm{B}}\left(\frac{1}{\zeta_{\mathrm{B}}}\mathbb{I}_{\mathrm{B}}+\frac{1}{\zeta_{\mathrm{A}}}\mathbb{G}^{\mathrm{vac}}_{\mathrm{BA}}\mathbb{T}_{\mathrm{A}}\mathbb{T}_{\mathrm{A}}^{\star}\mathbb{G}_{\mathrm{AB}}^{\mathrm{vac}\star}\right)\mathbb{Y}_{\mathrm{B}}^{\star}$ is not only Hermitian but also positive-semidefinite. This implies that if $\mathbb{Y}_{\mathrm{B}}$ had a nontrivial Hermitian part, any bounds on the singular values $y_{i}$ arising from this total nonnegative far-field scattering constraint must be more restrictive than if $\mathbb{Y}_{\mathrm{B}}$ were purely anti-Hermitian, meaning $\mathbb{Y}_{\mathrm{B}}=\mathrm{i}\operatorname{Im}(\mathbb{Y}_{\mathrm{B}})$. Thus, we take $\mathbb{Y}_{\mathrm{B}}$ to be anti-Hermitian, which also means that $\langle\mathbf{b}_{j},\mathbf{b}_{j}^{\star}\rangle=\mathrm{i}$ for every channel $j$. Additionally, we note that $\mathbb{G}^{\mathrm{vac}}_{\mathrm{BA}}\mathbb{T}_{\mathrm{A}}\mathbb{T}_{\mathrm{A}}^{\star}\mathbb{G}_{\mathrm{AB}}^{\mathrm{vac}\star}=\sum_{i}g_{i}^{2}\tau_{i}^{2}|\mathbf{b}^{(i)}\rangle\langle\mathbf{b}^{(i)}|$ follows from our singular value decompositions, so it remains to be proved that the operator $\sum_{j}\frac{y_{j}}{2\mathrm{i}}\left(|\mathbf{b}_{j}^{\star}\rangle\langle\mathbf{b}_{j}|-|\mathbf{b}_{j}\rangle\langle\mathbf{b}_{j}^{\star}|\right)-\sum_{j}y_{j}^{2}\left(\frac{1}{\zeta_{\mathrm{B}}}+\frac{1}{\zeta_{\mathrm{A}}}(g_{j}\tau_{j})^{2}\right)|\mathbf{b}_{j}^{\star}\rangle\langle\mathbf{b}_{j}^{\star}|$ is positive-semidefinite. To do this, we take the inner product of this operator with $|\mathbf{b}_{i}^{\star}\rangle$ on the right and $\langle\mathbf{b}_{i}^{\star}|$ on the left for a single fixed index $i$ (not to be confused with the imaginary unit), which yields the condition $y_{i}-y_{i}^{2}\left(\frac{1}{\zeta_{\mathrm{B}}}+\frac{1}{\zeta_{\mathrm{A}}}(g_{i}\tau_{i})^{2}\right)\geq 0$. This inequality is saturated to give the largest possible $y_{i}=\frac{\zeta_{\mathrm{B}}}{1+\frac{\zeta_{\mathrm{B}}}{\zeta_{\mathrm{A}}}(g_{i}\tau_{i})^{2}}$ in terms of $g_{i}\tau_{i}$, which in turn gives the largest possible flux,

[TABLE]

in terms of $g_{i}\tau_{i}$.

The above expression is maximized for each channel $i$, and for the sum in turn, if $\frac{\zeta_{\mathrm{B}}}{\zeta_{\mathrm{A}}}(g_{i}\tau_{i})^{2}=1$. This is solved to yield $\tau_{i}=\frac{1}{\sqrt{\frac{\zeta_{\mathrm{B}}}{\zeta_{\mathrm{A}}}}g_{i}}$, and plugging this into the expression for $y_{i}$ yields $y_{i}=\frac{\zeta_{\mathrm{B}}}{2}$, while plugging this into the expression for $\Phi$ yields a contribution $\frac{2}{\pi\zeta_{\mathrm{A}}\zeta_{\mathrm{B}}}(y_{i}g_{i}\tau_{i})^{2}=\frac{2}{\pi}\times\frac{1}{4}=\frac{1}{2\pi}$. We interpret this to mean that to obtain optimal heat transfer, the T-operator of body A in isolation must be engineered in a way that depends on the presence of body B, due to both the presence of the material enhancement factor $\zeta_{\mathrm{B}}$ and the dependence on the singular values $g_{i}$ of $\mathbb{G}^{\mathrm{vac}}_{\mathrm{BA}}$ propagating electromagnetic fields in vacuum from A to B, though no higher-order scattering effects come into play. In turn, the expression $y_{i}=\frac{\zeta_{\mathrm{B}}}{2}$ means that the effective T-operator of body B dressed by scattering from body A must actually exhibit maximal scattering, and not absorption, in the presence of body A Miller et al. (2016), though it is more difficult to extract information about the implications for the T-operator of body B in isolation; we clarify that maximal scattering includes both far-field scattering from body B in the presence of body A, as well as absorption from body A in the presence of body B. If these two conditions can be met simultaneously for the given channel $i$, which is effectively a rate-matching condition relating the absorption and scattering rates of each body in the presence of the other, then the contribution $\frac{1}{2\pi}$ is exactly the Landauer upper bound on energy transmission.

However, this analysis has so far ignored constraints on $\tau_{i}$. In particular, nonnegativity of far-field scattering from body A in isolation requires that $\operatorname{Im}(\mathbb{T}_{\mathrm{A}})-\frac{1}{\zeta_{\mathrm{A}}}\mathbb{T}_{\mathrm{A}}\mathbb{T}_{\mathrm{A}}^{\star}$ be a positive-semidefinite Hermitian operator Miller et al. (2016). Following similar steps as above, we may see that for the purposes of bounding $\tau_{i}$, the loosest constraints on $\tau_{i}$ arise if $\mathbb{T}_{\mathrm{A}}$ is purely anti-Hermitian, so taking the inner product of the above operator with respect to $|\mathbf{a}_{i}\rangle$ on the right and $\langle\mathbf{a}_{i}|$ on the left yields the condition $\tau_{i}\leq\zeta_{\mathrm{A}}$. The expression $\tau_{i}=\frac{1}{\sqrt{\frac{\zeta_{\mathrm{B}}}{\zeta_{\mathrm{A}}}}g_{i}}$ is consistent with this inequality only if the inequality $g_{i}\geq\frac{1}{\sqrt{\zeta_{\mathrm{A}}\zeta_{\mathrm{B}}}}$ also holds for the given index $i$. If not, and instead $g_{i}<\frac{1}{\sqrt{\zeta_{\mathrm{A}}\zeta_{\mathrm{B}}}}$, then $\tau_{i}=\zeta_{\mathrm{A}}$ must be used to maximize the contribution at index $i$ to $\Phi$, saturating the bound on the response of A in isolation. This yields $y_{i}=\frac{\zeta_{\mathrm{B}}}{1+\zeta_{\mathrm{A}}\zeta_{\mathrm{B}}g_{i}^{2}}\geq\frac{\zeta_{\mathrm{B}}}{2}$ and a contribution of $\frac{2}{\pi}\frac{\zeta_{\mathrm{A}}\zeta_{\mathrm{B}}g_{i}^{2}}{(1+\zeta_{\mathrm{A}}\zeta_{\mathrm{B}}g_{i}^{2})^{2}}\leq\frac{1}{2\pi}$ to $\Phi$. We interpret this to mean that if the singular value $g_{i}$ of $\mathbb{G}^{\mathrm{vac}}_{\mathrm{BA}}$ falls below a threshold involving the two material enhancement factors, then the optimal T-operator of body A in isolation corresponds to maximal absorption, the optimal effective T-operator of body B dressed by body A evinces the effects of multiple scattering with A, and the contribution to $\Phi$ similarly shows the effects of multiple scattering between the two bodies and is unable to saturate the Landauer bound for that channel.

To summarize, the bound on RHT may be written as

[TABLE]

where $\Theta$ is the Heaviside step function. As the material response (encoded in $\zeta_{\mathrm{A}}\zeta_{\mathrm{B}}$) increases, progressively more channels may saturate the Landauer limit per channel, so that the Landauer limit (summed over all channels) is reached asymptotically as $\zeta_{\mathrm{A}}\zeta_{\mathrm{B}}\to\infty$. However, the rate at which this divergence occurs depends on the general topology of the problem, as that determines how the singular values $g_{i}$ of $\mathbb{G}^{\mathrm{vac}}_{\mathrm{BA}}$ depend on the index $i$.

We emphasize that while the singular values $g_{i}$ of $\mathbb{G}^{\mathrm{vac}}_{\mathrm{BA}}$ are technically restricted to the domains of the objects to give the tightest bound on heat transfer, such a restriction is less than ideal given the explicit dependence on the shapes of the objects. However, as we prove in the appendices, the singular values $g_{i}$ of $\mathbb{G}^{\mathrm{vac}}_{\mathrm{BA}}$ are domain monotonic, meaning that they increase monotonically as the volumes of regions A and B increase; consequently, $\Phi_{\mathrm{opt}}$ is domain monotonic, as it is monotonically non-decreasing with respect to $g_{i}$ for each $i$. Separately from this, the regions containing only the material degrees of freedom of each body can be replaced by larger regions that fully enclose each body, as the T-operators of each body will commute with projections into the smaller subspaces corresponding to the actual material degrees of freedom. Thus, these bounds can be slightly loosened to be independent of body shapes, and can then be evaluated subject to constraints on topology and domain volumes as determined by the desired application [Fig. 1], e.g. ellipsoids with prescribed aspect ratios or films of prescribed thicknesses representative of compact or extended geometries, respectively. Essentially, the effective rank of $\mathbb{G}^{\mathrm{vac}}_{\mathrm{BA}}$, which determines the number of modes that could participate in RHT, is largely determined by the size and topology of the choice of bounding surface, which represents a general and fundamental geometric constraint on the bounds of RHT analogous to the general material constraints imposed by $\zeta_{p}$ for each body $p\in\{\mathrm{A},\mathrm{B}\}$.

III Comparison to alternative bounds

The bound for the RHT spectrum $\Phi$ in (8) may be compared to a number of other bounds. Strictly speaking, this bound is not necessarily the tightest general bound that could be formulated. In particular, using the relation $\mathbb{T}_{\mathrm{A}}^{-1}=\mathbb{V}_{\mathrm{A}}^{-1}-\mathbb{G}^{\mathrm{vac}}_{\mathrm{AA}}$ allows for writing (3) in terms of $\mathbb{T}_{\mathrm{A}}$ and $\operatorname{Im}(\mathbb{G}^{\mathrm{vac}}_{\mathrm{AA}})$ without reference to $\mathbb{V}_{\mathrm{A}}$. Such a procedure, in analogy with bounds on thermal emission which we detail in an upcoming paper Molesky et al. (2019), would more explicitly capture far-field radiative losses from bodies of finite size, which becomes more relevant at large separations where such losses may compete with RHT itself. However, as we show in the appendices, we find the resulting bound to be intractable, requiring self-consistent solution of systems of nonlinear equations to find the optimal singular values of $\mathbb{T}_{\mathrm{A}}$. Therefore, we do not further consider such a bound, and henceforth refer only to (8).

With respect to prior work, the most obvious point of comparison is the Landauer bound Pendry (1999), namely

[TABLE]

which simply depends on the number of modes participating in RHT, without any reference to separation, geometric constraints, multiple scattering, or even material constraints, let alone their interplay; consequently, in contrast to our bounds, there is no metric to evaluate how many participating modes can actually saturate the limit $\frac{1}{2\pi}$. Even modal analyses that technically do not necessarily assume saturation of the Landauer limits for every mode Pendry (1999); Bimonte (2009); Biehs et al. (2010); Ben-Abdallah and Joulain (2010) tend to neglect material effects, so the purely geometric arguments are valid only in the ray-optical regime where blackbody limits are reproduced.

An alternative bound can be derived by modifying the form of (8) such that the contributions from saturating the singular value bounds for $\tau_{i}$ are used for all indices $i$, not only those where $\zeta_{\mathrm{A}}\zeta_{\mathrm{B}}g_{i}^{2}<1$. We may write this as

[TABLE]

and term this the “scalar approximation” with the subscript “sc” because (as shown in the appendices) it can equivalently be derived by plugging into (2) the T-operators $\mathbb{T}_{p}=\mathrm{i}\tau_{p}\mathbb{I}_{p}$ for each body $p\in\{\mathrm{A},\mathrm{B}\}$, corresponding to uniform singular values $\tau_{p}$ for each $\mathbb{T}_{p}$; in the appendices, we further derive that this bound on RHT is maximal when each body exactly satisfies the condition for perfect isolated absorption Miller et al. (2016) for every incident field, meaning $\tau_{p}=\zeta_{p}$. As each contribution to $\Phi_{\mathrm{sc}}$ depends on the singular values $g_{i}$ in a nonmonotonic way, domain monotonicity is not obvious, though arguments based on the embedding of the T-operators in larger spaces allow for consideration of general bounding surfaces as geometric constraints, making this a useful bound restricted to the class of T-operators with uniform singular values. However, as we show in an upcoming paper, $\Phi_{\mathrm{sc}}$ overall can be shown to be domain monotonic in the near-field.

Another such bound neglects multiple scattering effects between the two bodies, and can be derived in several ways. As we show in the appendices, this bound follows from the scalar approximation $\Phi_{\mathrm{sc}}$ by neglecting the scattering operator $\mathbb{S}^{\mathrm{A}}_{\mathrm{B}}\to\mathbb{I}_{\mathrm{B}}$, or alternatively by ignoring the denominators $(1+\zeta_{\mathrm{A}}\zeta_{\mathrm{B}}g_{i}^{2})^{2}$ in $\Phi_{\mathrm{sc}}$ encoding the effects of multiple scattering within each channel. Regardless of the particular derivation, the end result may be written as

[TABLE]

which we term the “Born bound” as its neglect of multiple scattering effects between the two bodies is characteristic of a Born approximation to the full solution of Maxwell’s equations. The domain monotonicity of this bound trivially follows from that of $g_{i}$ for each channel $i$, and the validity of the embedding argument with respect to $\mathbb{T}_{p}$ still holds, meaning that $\Phi_{\mathrm{Born}}$ is also a useful bound for the RHT spectrum when considering bounding surfaces of arbitrary size and topology. We note that Miller et al Miller et al. (2015) equivalently derive this bound by limiting the singular values of $\mathbb{Y}_{\mathrm{B}}$ such that only far-field scattering of body B in the presence of body A needs to be nonnegative, meaning $\operatorname{Im}(\mathbb{Y}_{\mathrm{B}})-\mathbb{Y}_{\mathrm{B}}^{\star}\operatorname{Im}(\mathbb{V}_{\mathrm{B}}^{-1\star})\mathbb{Y}_{\mathrm{B}}$ should be positive-semidefinite; this leads to the bound $y_{i}\leq\zeta_{\mathrm{B}}$, so maximization of (7) subject to that constraint as well as $\tau_{i}\leq\zeta_{\mathrm{A}}$ simply requires saturation of both of these constraints. However, we have shown that positive-semidefiniteness of $\operatorname{Im}(\mathbb{Y}_{\mathrm{B}})-\mathbb{Y}_{\mathrm{B}}^{\star}\operatorname{Im}(\mathbb{V}_{\mathrm{B}}^{-1\star})\mathbb{Y}_{\mathrm{B}}$, corresponding to nonnegative scattering from body B in the presence of body A, is a looser constraint on $y_{i}$ than nonnegativity of scattering from the system as a whole, corresponding to positive-semidefiniteness of $\operatorname{Im}(\mathbb{Y}_{\mathrm{B}})-\mathbb{Y}_{\mathrm{B}}\left(\operatorname{Im}(\mathbb{V}_{\mathrm{B}}^{-1\star})+\mathbb{G}^{\mathrm{vac}}_{\mathrm{BA}}\mathbb{T}_{\mathrm{A}}\operatorname{Im}(\mathbb{V}_{\mathrm{A}}^{-1\star})\mathbb{T}_{\mathrm{A}}^{\star}\mathbb{G}_{\mathrm{AB}}^{\mathrm{vac}\star}\right)\mathbb{Y}_{\mathrm{B}}^{\star}$, and that only the latter produces a bound on $y_{i}$ that explicitly accounts for multiple scattering and absorption losses through body A. The Born bound not only neglects multiple scattering but also ignores far-field radiative losses, which means that while it is technically a looser bound, its usefulness is limited to deeply subwavelength systems exchanging heat in the near-field.

We are now in a position to relate $\Phi_{\mathrm{opt}}$, $\Phi_{\mathrm{sc}}$, $\Phi_{\mathrm{Born}}$, and $\Phi_{\mathrm{L}}$ to each other via inequalities. First, the discussion of (8) immediately makes clear that because some of the contributions to $\Phi_{\mathrm{opt}}$ equal corresponding contributions to $\Phi_{\mathrm{L}}$ while all of the others equal corresponding contributions to $\Phi_{\mathrm{sc}}$, then it must be the case that $\Phi_{\mathrm{sc}}\leq\Phi_{\mathrm{opt}}\leq\Phi_{\mathrm{L}}$. Additionally, because the contributions to $\Phi_{\mathrm{opt}}$ from each channel $i$ never exceed $\frac{1}{2\pi}$, and in particular because the saturation of the Landauer bound for each channel occurs for $\zeta_{\mathrm{A}}\zeta_{\mathrm{B}}g_{i}^{2}\geq 1$, then it must be the case that $\Phi_{\mathrm{sc}}\leq\Phi_{\mathrm{opt}}\leq\Phi_{\mathrm{Born}}$, because in the expression for $\Phi_{\mathrm{Born}}$, the contributions $\zeta_{\mathrm{A}}\zeta_{\mathrm{B}}g_{i}^{2}$ will always trivially exceed $1/4$ when $\zeta_{\mathrm{A}}\zeta_{\mathrm{B}}g_{i}^{2}\geq 1$, and will always exceed $\frac{\zeta_{\mathrm{A}}\zeta_{\mathrm{B}}g_{i}^{2}}{(1+\zeta_{\mathrm{A}}\zeta_{\mathrm{B}}g_{i}^{2})^{2}}$ when $\zeta_{\mathrm{A}}\zeta_{\mathrm{B}}g_{i}^{2}<1$. In general, it is not possible to write an inequality relation between $\Phi_{\mathrm{Born}}$ and $\Phi_{\mathrm{L}}$ in all situations, because $\Phi_{\mathrm{Born}}$ may have some contributions $\frac{2}{\pi}\zeta_{\mathrm{A}}\zeta_{\mathrm{B}}g_{i}^{2}$ which fall above or below $\frac{1}{2\pi}$, and the geometry determining $g_{i}$ would have to be known in order to know how many fall above or below. Thus, we may write the overall inequalities as

[TABLE]

IV Concluding remarks

We have determined bounds for the RHT spectrum $\Phi$ based purely on algebraic arguments. In particular, we show that there is a tension between optimizing transmission channels and material/geometric constraints placed on each object in isolation as well as in in the presence of the other; as a result, some but not all channels can saturate previously derived Landauer bounds, while others are restricted by the aforementioned constraints. By virtue of domain monotonicity, these bounds can be applied in a shape-independent manner, so while they can be evaluated analytically in highly symmetric bounding surfaces, they can just as easily be evaluated numerically in more complicated domains depending on specific design constraints [Fig. 1]. Similarly, the dependence on the material response factor $\zeta=|\chi|^{2}/\operatorname{Im}\chi$ does not make explicit reference to a particular frequency or material model. In comparison, the Landauer bounds yield overly optimistic predictions, while choosing a scalar response for each object corresponding to maximal absorption of every incident field in isolation yields overly pessimistic predictions. Additionally, we find that previous work by Miller et al Miller et al. (2015) also yields overly optimistic predictions compared to our current bounds, because those derivations neglect the geometric effects of multiple scattering between the two bodies and consequently overestimate the optimal response of one body in the presence of the other. We point out that while our bounds are always tighter than Landauer and Born limits for any given bounding domain, they say nothing about which domains may yield the tightest per-volume limits given material constraints, or whether they may in fact be attained by physically realizable structures.In summary, while Born and Landauer limits are technically upper bounds on RHT, their neglect of multiple scattering and material losses, respectively, render them loose compared to the bounds presented here.

We emphasize that in contrast to Born limits Miller et al. (2015), our bounds are valid from the near-field all the way through the far-field for bodies of arbitrary size, as we never needed to exploit nonretarded or quasistatic approximations. We also point out that in the near-field, $\mathbb{G}^{\mathrm{vac}}_{\mathrm{BA}}$ is real-valued in position space, so arguments involving the anti-Hermiticity of $\mathbb{G}^{\mathrm{vac}}_{\mathrm{BA}}\mathbb{T}_{\mathrm{A}}\mathbb{G}^{\mathrm{vac}}_{\mathrm{AB}}$ imply the anti-Hermiticity of $\mathbb{T}_{\mathrm{A}}$, as $\operatorname{Im}(\mathbb{G}^{\mathrm{vac}}_{\mathrm{BA}}\mathbb{T}_{\mathrm{A}}\mathbb{G}^{\mathrm{vac}}_{\mathrm{AB}})\to\mathbb{G}^{\mathrm{vac}}_{\mathrm{BA}}\operatorname{Im}(\mathbb{T}_{\mathrm{A}})\mathbb{G}^{\mathrm{vac}}_{\mathrm{AB}}$. In this context, Born bounds only yield sensible results in the near-field, as they tend to diverge beyond the near-field regime for objects of infinite extent in at least one spatial dimension faster than expected.

In a complementary paper Venkataram et al. (2019), we analyze these bounds in the near-field in specific geometries of interest, particularly dipolar and extended (infinite area) bodies. We take particular note of the scalar approximation because, though representing a restricted version of our general bounds and resulting in comparably pessimistic predictions, it may prove challenging if not impossible to optimize every individual channel to saturate its Landauer limit when the material constraints are not saturated. While this mere observation does not preclude the possibility of larger RHT (as allowed by our general bounds), from a design perspective it might prove more feasible to consider bodies capable of maximally absorbing every incident field in isolation, as could occur near a polaritonic resonance. Moreover, as we show in this follow-up work, the simplicity of the scalar approximation in the near field lends itself to various mathematical manipulations which lead to the following suggestive results. First, the scalar bound is a local stationary point with respect to perturbations of the T-operators of each body, away from the perfect isolated absorption condition. Second, $\Phi_{\mathrm{sc}}$ can be shown to be domain-monotonic, meaning that for fixed $\zeta_{p}$ for each body $p\in\{\mathrm{A},\mathrm{B}\}$, the near-field RHT spectrum will always increase with increasing object volumes. Together, these findings suggest that the role of nanostructuring in enhancing near-field RHT will be limited, and this has important implications for the theoretical and experimental design of devices for cooling, heat dissipation, and energy generation.

Acknowledgments

The authors would like to thank Riccardo Messina and Pengning Chao for helpful discussions. This work was supported by the National Science Foundation under Grants No. DMR-1454836, DMR 1420541, DGE 1148900, the Cornell Center for Materials Research MRSEC (award no. DMR-1719875), and the Defense Advanced Research Projects Agency (DARPA) under agreement HR00111820046. The views, opinions and/or findings expressed are those of the authors and should not be interpreted as representing the official views or policies of the Department of Defense or the U.S. Government.

Appendix A Notation

We briefly discuss the notation used through the main text and the appendices. A vector field $\mathbf{v}(\mathbf{x})$ will be denoted as $|\mathbf{v}\rangle$. The conjugated inner product is $\langle\mathbf{u},\mathbf{v}\rangle=\int\leavevmode\nobreak\ \mathrm{d}^{3}x\leavevmode\nobreak\ \mathbf{u}^{\star}(\mathbf{x})\cdot\mathbf{v}(\mathbf{x})$. An operator $\mathbb{A}(\mathbf{x},\mathbf{x}^{\prime})$ will be denoted as $\mathbb{A}$, with $\int\leavevmode\nobreak\ \mathrm{d}^{3}x^{\prime}\leavevmode\nobreak\ \mathbb{A}(\mathbf{x},\mathbf{x}^{\prime})\cdot\mathbf{v}(\mathbf{x}^{\prime})$ denoted as $\mathbb{A}|\mathbf{v}\rangle$. The Hermitian conjugate $\mathbb{A}^{\dagger}$ is defined such that $\langle\mathbf{u},\mathbb{A}^{\dagger}\mathbf{v}\rangle=\langle\mathbb{A}\mathbf{u},\mathbf{v}\rangle$. The anti-Hermitian part of a square operator (whose domain and range are the same size) is defined as the operator $\operatorname{Asym}(\mathbb{A})=(\mathbb{A}-\mathbb{A}^{\dagger})/(2\mathrm{i})$. Finally, the trace of an operator is $\operatorname{Tr}(\mathbb{A})=\int\leavevmode\nobreak\ \mathrm{d}^{3}x\leavevmode\nobreak\ \operatorname{Tr}(\mathbb{A}(\mathbf{x},\mathbf{x}))$. Through this paper, unless stated explicitly otherwise, all quantities implicitly depend on $\omega$, and such dependence will be notationally suppressed for brevity.

Appendix B Proof of trace maximization via shared singular vectors

In this section, we prove the following lemma: if operators $\mathbb{A}_{n}$ for $n\in\{1,2,\ldots,N\}$ have fixed singular values labeled $\sigma^{(n)}_{i}$, then the singular vectors that maximize $\operatorname{Tr}[\mathbb{A}_{1}\mathbb{A}_{2}\ldots\mathbb{A}_{N}]$ are common between operators multiplied together. That is, the singular value decomposition of $\mathbb{A}_{n}$ should follow $\mathbb{A}_{n}=\sum_{i}\sigma^{(n)}_{i}|\mathbf{a}^{(n)}_{i}\rangle\langle\mathbf{a}^{(n+1)}_{i}|$ for $n\in\{1,2,\ldots,N-1\}$, with $\mathbb{A}_{N}=\sum_{i}\sigma^{(N)}_{i}|\mathbf{a}^{(N)}_{i}\rangle\langle\mathbf{a}^{(1)}_{i}|$, where the vectors $|\mathbf{a}^{(n)}_{i}\rangle$ are orthonormal for each $n$ such that $\langle\mathbf{a}^{(n)}_{i},\mathbf{a}^{(n)}_{j}\rangle=\delta_{ij}$. This lemma will hold even if each $\mathbb{A}_{n}$ is not square, as long as $\mathbb{A}_{n}\mathbb{A}_{n+1}$ forms a valid nontrivial operator product, as these can be embedded in larger spaces padded with more vanishing singular values. Thus, we restrict our consideration to square operators. Moreover, associativity means $\mathbb{A}_{n}\mathbb{A}_{n+1}\mathbb{A}_{n+2}=(\mathbb{A}_{n}\mathbb{A}_{n+1})\mathbb{A}_{n+2}$, and the trace of a product of operators is invariant under cyclic permutations, so we ultimately only consider maximizing the trace of a product of two operators, as maximization of the trace of products of more than two operators follows inductively from this.

To maximize $\operatorname{Tr}[\mathbb{A}\mathbb{B}]$, we start by writing

[TABLE]

where $N$ is the size of the space; this may be larger than the rank of either $\mathbb{A}$ or $\mathbb{B}$, but the point is moot because the singular values are fixed, whether they vanish or not, and it has already been assumed that $\mathbb{A}$ and $\mathbb{B}$ are square. We also assume that the singular values are ordered such that $\sigma_{i}\geq\sigma_{i+1}$ for all $i\in\{1,2,\ldots,N-1\}$ and $\tau_{j}\geq\tau_{j+1}$ for all $j\in\{1,2,\ldots,N-1\}$. This allows for writing

[TABLE]

in terms of $p_{ij}=\langle\mathbf{v}_{i},\mathbf{w}_{j}\rangle$ and $q_{ji}=\langle\mathbf{y}_{j},\mathbf{u}_{i}\rangle$. As the singular vectors are orthonormal, then $p_{ij}$ and $q_{ji}$ are the elements of unitary matrices, satisfying $\sum_{j=1}^{N}|p_{ij}|^{2}=\sum_{i=1}^{N}|p_{ij}|^{2}=\sum_{j=1}^{N}|q_{ji}|^{2}=\sum_{i=1}^{N}|q_{ji}|^{2}=1$. We consider only real nonnegative $\operatorname{Tr}[\mathbb{A}\mathbb{B}]$, so that is maximized when $\sigma_{i}p_{ij}\tau_{j}q_{ji}$ are all nonnegative; this means the singular vectors can be chosen without loss of generality such that $p_{ij}$ and $q_{ji}$ are real and nonnegative, implying $p_{ij}$ and $q_{ji}$ are the elements of real-valued orthogonal matrices.

We use induction to prove that maximizing the trace requires that $\{|\mathbf{u}_{i}\rangle\}$ be the duals of $\{\langle\mathbf{y}_{j}|\}$, and that $\{\langle\mathbf{v}_{i}|\}$ by the duals of $\{|\mathbf{w}_{j}\rangle\}$. The case $N=1$ is trivial, as all quantities are scalars. For $N=2$, we use orthogonality to note that $p_{1,2}=p_{2,1}$, $q_{1,2}=q_{2,1}$, $p_{1,1}=p_{2,2}=\sqrt{1-p_{1,2}^{2}}$, and $q_{1,1}=q_{2,2}=\sqrt{1-q_{1,2}^{2}}$. As a result, we may write $\operatorname{Tr}[\mathbb{A}\mathbb{B}]=\sqrt{(1-p_{1,2}^{2})(1-q_{1,2}^{2})}(\sigma_{1}\tau_{1}+\sigma_{2}\tau_{2})+p_{1,2}q_{1,2}(\sigma_{1}\tau_{2}+\sigma_{2}\tau_{1})$. As the first term in parentheses is larger than the second term in parentheses by the nonnegative value $(\sigma_{1}-\sigma_{2})(\tau_{1}-\tau_{2})$ given the ordering of singular values, having the left singular vectors of one operator not be duals of the right singular vectors of the other and vice versa could only increase the trace if $\sqrt{(1-p_{1,2}^{2})(1-q_{1,2}^{2})}+p_{1,2}q_{1,2}>1$, but this leads to the impossible condition $0>(p_{1,2}-q_{1,2})^{2}$, so we can only have $p_{1,2}=q_{1,2}$ for the trace to be maximized, implying the duality result must hold.

The inductive step assumes an arbitrary $N-1$ and moves from there to proving the statement for $N$. Without loss of generality, we consider first the contribution of the largest singular value $\tau_{1}$ of $\mathbb{B}$, namely $\tau_{1}|\mathbf{w}_{1}\rangle\langle\mathbf{y}_{1}|$, interacting with $\mathbb{A}=\sum_{i}\sigma_{i}|\mathbf{u}_{i}\rangle\langle\mathbf{v}_{i}|$ in the trace. This yields the contribution $\sum_{i}(p_{i,1}-q_{1,i})^{2}=2\left(1-\sum_{i}p_{i,1}q_{1,i}\right)\geq 0$ using the fact that $\sum_{i}p_{i,1}^{2}=\sum_{i}q_{1,i}^{2}=1$, so this in turn gives the condition $\sum_{i}p_{i,1}q_{1,i}\leq 1$. The trace can be seen to be maximal when the above condition is saturated, so $\sum_{i}p_{i,1}q_{1,i}=1$, which implies $p_{i,1}=q_{1,i}$ for every $i$. As this also holds when the roles of $\mathbb{A}$ and $\mathbb{B}$ are interchanged, and as this can be progressively carried out for each successively smaller singular value given orthonormality of the singular vectors, then the duality condition must hold, completing the proof.

Appendix C Derivation of radiative heat transfer formulas

In this section, we derive the formula for the radiative heat transfer spectrum between two bodies, valid even without assumptions about retardation, homogeneity, locality, or isotropy. The formula depends on individual T-operators and the vacuum Green’s function, and follows a previous derivation Venkataram et al. (2018) which considered energy transfer by fluctuating volume currents. Through further derivation, we also equate this formula to another formula involving the susceptibilities and the full Maxwell Green’s function, and use that to recast the heat transfer spectrum in a Landauer form, whence we prove that the singular values of the Landauer transmission operator for RHT do not exceed $1/4$. Finally, we prove that the formula for thermal emission of a single body in isolation can be derived from the formula for RHT between two bodies in vacuum, by taking the second body to fully enclose the first and to be perfectly absorbing, thus taking on the role of a perfectly absorbing medium (vacuum).

C.1 T-operator formula

Our derivation of the heat transfer spectrum from the fluctuation–dissipation theorem for dipole sources in each body follows prior work Venkataram et al. (2018), which we reproduce here for clarity. Consider two bodies A and B in vacuum with general susceptibilities $\mathbb{V}_{p}$ for $p\in\{\mathrm{A},\mathrm{B}\}$ which may be inhomogeneous, nonlocal, or anisotropic. Maxwell’s equations may be written in integral form as

[TABLE]

for the fields $|\mathbf{E}\rangle$ and total polarizations $|\mathbf{P}\rangle$ in terms of the polarization sources $|\mathbf{P}^{(0)}\rangle$, after defining

[TABLE]

in block form for the material degrees of freedom constituting each object. By defining the total T-operator via

[TABLE]

where $\mathbb{T}_{p}^{-1}=\mathbb{V}_{p}^{-1}-\mathbb{G}^{\mathrm{vac}}_{pp}$, then Maxwell’s equations can be formally solved to yield

[TABLE]

obtained by applying formulas for the block matrix inverse to compute $\mathbb{T}$. We also define the projection operators,

[TABLE]

such that (abusing notation) $\mathbb{I}_{p}$ is the projection onto the material degrees of freedom of body $p$.

We consider the energy flow from fluctuating dipole sources only in body A into material degrees of freedom in body B, noting that reciprocity would yield the same heat transfer if the roles of bodies A and B were interchanged. This means $|\mathbf{P}^{(0)}\rangle=\begin{bmatrix}|\mathbf{P}^{(0)}_{\mathrm{A}}\rangle\\ 0\end{bmatrix}$ defines the fluctuating sources in body A. The heat transfer spectrum is the ensemble-averaged work, denoted by $\langle\cdots\rangle$, done by the field,

[TABLE]

where $|\mathbf{J}\rangle=-\mathrm{i}\omega|\mathbf{P}\rangle$. Using the Hermiticity and idempotence of $\mathbb{I}_{p}$ yields $\Phi=-\frac{\omega}{4\mathrm{i}}(\langle\langle\mathbb{I}_{\mathrm{B}}\mathbf{P},\mathbf{E}\rangle\rangle-\langle\langle\mathbf{E},\mathbb{I}_{\mathrm{B}}\mathbf{P}\rangle\rangle)$, and using the results of (21) gives

[TABLE]

in terms of the fluctuating sources $|\mathbf{P}^{(0)}_{\mathrm{A}}\rangle$. As these fluctuations are thermal in nature, their correlations are given by the fluctuation–dissipation theorem

[TABLE]

(suppressing the Planck function $\Pi$ as it has already been factored to be separate from $\Phi$), yielding

[TABLE]

as the dressed radiative heat transfer spectrum.

To prove equivalence of this expression for $\Phi$ to that involving only $\mathbb{G}^{\mathrm{vac}}$ and $\mathbb{T}_{p}$, it is useful to explicitly invoke reciprocity: $\mathbb{V}_{p}^{\top}=\mathbb{V}_{p}$, $\mathbb{T}_{p}^{\top}=\mathbb{T}_{p}$, and $(\mathbb{G}^{\mathrm{vac}}_{pq})^{\top}=\mathbb{G}^{\mathrm{vac}}_{qp}$ for $p,q\in\{\mathrm{A},\mathrm{B}\}$, so this also means $\mathbb{V}_{p}^{\dagger}=\mathbb{V}_{p}^{\star}$, $\mathbb{T}_{p}^{\dagger}=\mathbb{T}_{p}^{\star}$, and $(\mathbb{G}^{\mathrm{vac}}_{pq})^{\dagger}=\mathbb{G}^{\mathrm{vac}\star}_{qp}$; this means $\operatorname{Asym}(\mathbb{A})=\operatorname{Im}(\mathbb{A})$ for $\mathbb{A}\in\{\mathbb{V}_{p},\mathbb{T}_{p},\mathbb{G}^{\mathrm{vac}}_{pp}\}$. This allows for writing the operators

[TABLE]

in block matrix form, where the projection onto A allows for truncation to the appropriate block column or row for notational convenience; note that $\mathbb{I}_{\mathrm{A}}\mathbb{T}^{\dagger}$ should actually be a row vector, but has been written as a column for ease of reading. Multiplying these matrices together, it can be noted that $\mathbb{G}^{\mathrm{vac}}_{\mathrm{BA}}(\mathbb{T}_{\mathrm{A}}^{-1}-\mathbb{G}^{\mathrm{vac}}_{\mathrm{AB}}\mathbb{T}_{\mathrm{B}}\mathbb{G}^{\mathrm{vac}}_{\mathrm{BA}})^{-1}=\mathbb{G}^{\mathrm{vac}}_{\mathrm{BA}}(\mathbb{I}_{\mathrm{A}}-\mathbb{T}_{\mathrm{A}}\mathbb{G}^{\mathrm{vac}}_{\mathrm{AB}}\mathbb{T}_{\mathrm{B}}\mathbb{G}^{\mathrm{vac}}_{\mathrm{BA}})^{-1}\mathbb{T}_{\mathrm{A}}=(\mathbb{I}_{\mathrm{B}}-\mathbb{G}^{\mathrm{vac}}_{\mathrm{BA}}\mathbb{T}_{\mathrm{A}}\mathbb{G}^{\mathrm{vac}}_{\mathrm{AB}}\mathbb{T}_{\mathrm{B}})^{-1}\mathbb{G}^{\mathrm{vac}}_{\mathrm{BA}}\mathbb{T}_{\mathrm{A}}=\mathbb{S}^{\mathrm{A}}_{\mathrm{B}}\mathbb{G}^{\mathrm{vac}}_{\mathrm{BA}}\mathbb{T}_{\mathrm{A}}$ using the definition of the scattering operator $\mathbb{S}^{\mathrm{A}}_{\mathrm{B}}=(\mathbb{I}_{\mathrm{B}}-\mathbb{G}^{\mathrm{vac}}_{\mathrm{BA}}\mathbb{T}_{\mathrm{A}}\mathbb{G}^{\mathrm{vac}}_{\mathrm{AB}}\mathbb{T}_{\mathrm{B}})^{-1}$. Additionally, using the definition $\mathbb{T}_{p}^{-1}=\mathbb{V}_{p}^{-1}-\mathbb{G}^{\mathrm{vac}}_{pp}$, it is easy to prove that $\operatorname{Im}(\mathbb{T}_{\mathrm{B}})-\mathbb{T}_{\mathrm{B}}^{\star}\operatorname{Im}(\mathbb{G}^{\mathrm{vac}}_{\mathrm{BB}})\mathbb{T}_{\mathrm{B}}=\mathbb{T}_{\mathrm{B}}^{\star}\operatorname{Im}(\mathbb{V}_{\mathrm{B}}^{-1\star})\mathbb{T}_{\mathrm{B}}$, and likewise $\operatorname{Im}(\mathbb{T}_{\mathrm{A}})-\mathbb{T}_{\mathrm{A}}\operatorname{Im}(\mathbb{G}^{\mathrm{vac}}_{\mathrm{AA}})\mathbb{T}_{\mathrm{A}}^{\star}=\mathbb{T}_{\mathrm{A}}\operatorname{Im}(\mathbb{V}_{\mathrm{A}}^{-1\star})\mathbb{T}_{\mathrm{A}}^{\star}$. Thus, the result is (2) in the main text, as expected.

C.2 Derivation of Green’s function heat transfer formula

Our derivation of the bounds in the main text relies on the relationship between the heat transfer spectrum $\Phi$ written in terms of the vacuum Green’s function and the T-operators of individual objects, to the heat transfer formula Jin et al. (2019)

[TABLE]

where $\mathbb{G}=(\mathbb{G}^{\mathrm{vac}-1}-\mathbb{V}_{\mathrm{A}}-\mathbb{V}_{\mathrm{B}})^{-1}$ is the full Maxwell Green’s function in the presence of both bodies, with the block $\mathbb{G}_{\mathrm{BA}}$ representing the fields in body B due to dipole sources in body A. We start with the T-operator by rewriting

[TABLE]

after using $\operatorname{Asym}(\mathbb{V}_{\mathrm{B}}^{-1\dagger})=\mathbb{V}_{\mathrm{B}}^{-1\dagger}\operatorname{Asym}(\mathbb{V}_{\mathrm{B}})\mathbb{V}_{\mathrm{B}}^{-1}$ and $\operatorname{Asym}(\mathbb{V}_{\mathrm{A}}^{-1\dagger})=\mathbb{V}_{\mathrm{A}}^{-1}\operatorname{Asym}(\mathbb{V}_{\mathrm{A}})\mathbb{V}_{\mathrm{A}}^{-1\dagger}$ along with invariance of the trace under cyclic permutations of operator products. From this, it can be seen that the two expressions for $\Phi(\omega)$ are guaranteed to be the same if the operator $\mathbb{G}_{\mathrm{BA}}$ is the same as $\mathbb{V}_{\mathrm{B}}^{-1}\mathbb{T}_{\mathrm{B}}\mathbb{S}^{\mathrm{A}}_{\mathrm{B}}\mathbb{G}^{\mathrm{vac}}_{\mathrm{BA}}\mathbb{T}_{\mathrm{A}}\mathbb{V}_{\mathrm{A}}^{-1}$. We use the fact that $\mathbb{V}_{\sigma}^{-1}=\mathbb{T}_{\sigma}^{-1}+\mathbb{G}^{\mathrm{vac}}_{\sigma\sigma}$ to say that

[TABLE]

must hold. To prove that this is equal to $\mathbb{G}_{\mathrm{BA}}$, we use the definition

[TABLE]

in conjunction with definitions of $\mathbb{G}^{\mathrm{vac}}$ and $\mathbb{T}$ as $2\times 2$ block matrices in (20) to write

[TABLE]

for this system. Performing this matrix multiplication, recognizing that $\mathbb{G}^{\mathrm{vac}}_{\mathrm{BA}}\mathbb{T}_{\mathrm{A}}\mathbb{G}^{\mathrm{vac}}_{\mathrm{AB}}(\mathbb{T}_{\mathrm{B}}^{-1}-\mathbb{G}^{\mathrm{vac}}_{\mathrm{BA}}\mathbb{T}_{\mathrm{A}}\mathbb{G}^{\mathrm{vac}}_{\mathrm{AB}})^{-1}=\mathbb{S}^{\mathrm{A}}_{\mathrm{B}}-\mathbb{I}_{\mathrm{B}}$, using the fact definition of $\mathbb{S}^{\mathrm{A}}_{\mathrm{B}}$, and collecting and canceling terms leads to the proof of the equality $\mathbb{G}_{\mathrm{BA}}=\mathbb{V}_{\mathrm{B}}^{-1}\mathbb{T}_{\mathrm{B}}\mathbb{S}^{\mathrm{A}}_{\mathrm{B}}\mathbb{G}^{\mathrm{vac}}_{\mathrm{BA}}\mathbb{T}_{\mathrm{A}}\mathbb{V}_{\mathrm{A}}^{-1}$.

C.3 Landauer bounds on heat transfer singular values

We now prove that radiative heat transfer between arbitrarily shaped bodies can also be expressed as the trace of a transmission matrix whose singular values can be bounded above, similar to previously derived bounds in planar media. This relation intuitively connects the finite value of the RHT bounds and approximate low rank of $\mathbb{G}_{\mathrm{BA}}$, and can be proved as follows. For this, we use the cyclic property of the trace to define

[TABLE]

where $\mathbb{Q}=\operatorname{Im}(\mathbb{V}_{\mathrm{B}})^{1/2}\mathbb{G}_{\mathrm{BA}}\operatorname{Im}(\mathbb{V}_{\mathrm{A}})^{1/2}$ is the heat transmission operator.

The definition $\mathbb{G}^{-1}=\mathbb{G}^{\mathrm{vac}-1}-(\mathbb{V}_{\mathrm{A}}+\mathbb{V}_{\mathrm{B}})$ along with the fact that the vacuum Maxwell operator $\mathbb{G}^{\mathrm{vac}-1}$ is real-valued in position space allows for writing

[TABLE]

which relates dissipation in polarization currents and electromagnetic fields in equilibrium. Additionally, the fact that $\operatorname{Asym}(\mathbb{V}_{p})$ is a Hermitian positive-definite operator for each body $p\in\{\mathrm{A},\mathrm{B}\}$ means it has a unique square root $\operatorname{Asym}(\mathbb{V}_{\mathrm{p}})^{1/2}$. Rearranging the above equation, multiplying both sides by $2\leavevmode\nobreak\ \operatorname{Asym}(\mathbb{V}_{\mathrm{A}})^{1/2}$, and adding $\mathbb{I}_{\mathrm{A}}$ to both sides gives

[TABLE]

where we recognize the equality $\mathbb{Q}^{\dagger}\mathbb{Q}=\operatorname{Asym}(\mathbb{V}_{\mathrm{A}})^{1/2}\mathbb{G}^{\dagger}\operatorname{Asym}(\mathbb{V}_{\mathrm{B}})\mathbb{G}\operatorname{Asym}(\mathbb{V}_{\mathrm{A}})^{1/2}$. Following this substitution, this may be factored as

[TABLE]

where $\mathbb{G}$ has been replaced by its blocks $\mathbb{G}_{\mathrm{AA}}$ and $\mathbb{G}_{\mathrm{AB}}$ due to multiplications on each each side by $\operatorname{Asym}(\mathbb{V}_{p})^{1/2}$ for $p\in\{\mathrm{A},\mathrm{B}\}$ (and likewise for $\mathbb{G}^{\dagger}$). This expression is the sum of two Hermitian positive-semidefinite operators equal to the identity; though this has been done for body A, reciprocity of heat transfer yields a similar expression in terms of the operators for body B. Consequently, the singular values of the operator $\mathbb{Q}^{\dagger}\mathbb{Q}$ entering the trace expression for $\Phi(\omega)$ must all be less than or equal to $1/4$. We emphasize that this derivation is valid for compact or extended structures of arbitrary geometry, without any need to expand heat transfer in terms of incoming and outgoing plane waves specific to translationally symmetric systems Pendry (1999).

C.4 Single-body thermal radiation from two-body radiative heat

transfer

In this section, we prove that the formula for thermal emission of a single body in isolation in vacuum can be derived by starting from the formula for heat transfer between two bodies in vacuum under the following conditions. We take body A to be the thermal emitter in question, while body B is taken to fully surround body A as a shell of inner radius $r_{\mathrm{B}}$ and outer radius $R_{\mathrm{B}}$ with susceptibility $\mathbb{V}_{\mathrm{B}}=\chi_{\mathrm{B}}\mathbb{I}_{\mathrm{B}}$, and take the simultaneous limits $\operatorname{Im}(\chi_{\mathrm{B}})\to 0$ and $\omega r_{\mathrm{B}}/c\to\infty$ constrained by $\omega(R_{\mathrm{B}}-r_{\mathrm{B}})/c\to\infty$ and $\omega(R_{\mathrm{B}}-r_{\mathrm{B}})\operatorname{Im}(\chi_{\mathrm{B}})/c\to 1$, in which case body B takes on the role of a perfectly absorbing medium.

To start, we note that $\mathbb{T}_{\mathrm{B}}^{-1}=\mathbb{V}_{\mathrm{B}}^{-1}-\mathbb{G}^{\mathrm{vac}}_{\mathrm{BB}}$ in conjunction with reciprocity allows for writing (2) as

[TABLE]

In the aforementioned size and susceptibility limits for body B, $\mathbb{T}_{\mathrm{B}}\to 0$ so $\mathbb{S}^{\mathrm{A}}_{\mathrm{B}}\to\mathbb{I}_{\mathrm{B}}$, and $\mathbb{T}_{\mathrm{B}}\mathbb{S}^{\mathrm{A}}_{\mathrm{B}}\operatorname{Asym}(\mathbb{V}_{\mathrm{B}}^{-1\dagger})\mathbb{S}^{\mathrm{A}\dagger}_{\mathrm{B}}\mathbb{T}_{\mathrm{B}}^{\dagger}\to\operatorname{Asym}(\mathbb{V}_{\mathrm{B}})$. Using reciprocity, this yields $\Phi=\frac{2}{\pi}\operatorname{Tr}[\mathbb{G}^{\mathrm{vac}}_{\mathrm{AB}}\operatorname{Im}(\mathbb{V}_{\mathrm{B}})\mathbb{G}_{\mathrm{BA}}^{\mathrm{vac}\star}(\operatorname{Im}(\mathbb{T}_{\mathrm{A}})-\mathbb{T}_{\mathrm{A}}^{\star}\operatorname{Im}(\mathbb{G}^{\mathrm{vac}}_{\mathrm{AA}})\mathbb{T}_{\mathrm{A}})$. For a system with a general susceptibility $\mathbb{V}$ and Maxwell Green’s function $\mathbb{G}=(\mathbb{G}^{\mathrm{vac}-1}-\mathbb{V})^{-1}$, the relations $\mathbb{G}\operatorname{Im}(\mathbb{V})\mathbb{G}^{\star}=\mathbb{G}^{\star}\operatorname{Im}(\mathbb{V})\mathbb{G}=\operatorname{Im}(\mathbb{G})$ will always hold. Considering B in isolation, the simultaneous constrained limits of infinite size and infinitesimal susceptibility mean that $\mathbb{G}^{\mathrm{vac}}_{\mathrm{AB}}\operatorname{Im}(\mathbb{V}_{\mathrm{B}})\mathbb{G}_{\mathrm{BA}}^{\mathrm{vac}\star}\to\operatorname{Im}(\mathbb{G}^{\mathrm{vac}}_{\mathrm{AA}})$. Finally, this yields the emission formula

[TABLE]

in agreement with the formula derived by Krüger et al Krüger et al. (2012), where the subscripts A have been dropped as there is only one material body under consideration given that body B has effectively vanished.

Appendix D Nonnegative far-field scattering

In this section, we derive conditions on operators describing the response of bodies $p\in\{\mathrm{A},\mathrm{B}\}$ for far-field scattering to be nonnegative. In general, given a susceptibility $\mathbb{V}$ and an associated T-operator $\mathbb{T}$, the far-field scattered power from a given incident field $|\mathbf{E}^{\mathrm{inc}}\rangle$ is $\frac{\omega}{2}\langle\mathbf{E}^{\mathrm{inc}},(\operatorname{Im}(\mathbb{T})-\mathbb{T}^{\star}\operatorname{Im}(\mathbb{V}^{-1\star})\mathbb{T})\mathbf{E}^{\mathrm{inc}}\rangle$, and for this to be nonnegative, the operator $\operatorname{Im}(\mathbb{T})-\mathbb{T}^{\star}\operatorname{Im}(\mathbb{V}^{-1\star})\mathbb{T}$ must be positive-semidefinite. This must hold true for each body in isolation, meaning $\operatorname{Im}(\mathbb{T}_{p})-\mathbb{T}_{p}^{\star}\operatorname{Im}(\mathbb{V}_{p}^{-1\star})\mathbb{T}_{p}$ must be positive-semidefinite for each $p\in\{\mathrm{A},\mathrm{B}\}$. However, more conditions can be derived when the two bodies are proximate to each other.

As a first step, we show that the operator $\mathbb{Y}_{\mathrm{B}}=(\mathbb{T}_{\mathrm{B}}^{-1}-\mathbb{G}^{\mathrm{vac}}_{\mathrm{BA}}\mathbb{T}_{\mathrm{A}}\mathbb{G}^{\mathrm{vac}}_{\mathrm{AB}})^{-1}=\mathbb{T}_{\mathrm{B}}\mathbb{S}^{\mathrm{A}}_{\mathrm{B}}$ is an effective T-operator for body B dressed by the proximity of body A. In particular, the definitions (19) and (20) can be plugged into (21) and rearranged in order to write

[TABLE]

and if only sources in A are relevant, then we may set $|\mathbf{P}^{(0)}_{\mathrm{B}}\rangle\to 0$ and define an effective incident field $|\mathbf{E}^{\mathrm{inc}}\rangle=\mathbb{G}^{\mathrm{vac}}_{\mathrm{BA}}\mathbb{T}_{\mathrm{A}}\mathbb{V}_{\mathrm{A}}^{-1}|\mathbf{P}^{(0)}_{\mathrm{A}}\rangle$ which depends on multiple scattering within A but not on any properties of B apart from projection onto its volumetric degrees of freedom. This means $|\mathbf{P}_{\mathrm{B}}\rangle=\mathbb{Y}_{\mathrm{B}}|\mathbf{E}^{\mathrm{inc}}\rangle$, which is interpreted to mean that the total induced polarization in B arises from the response of B dressed in the presence of A, namely $\mathbb{Y}_{\mathrm{B}}$, acting on the effective incident field $|\mathbf{E}^{\mathrm{inc}}\rangle$ accounting only for body A; this is analogous to $\mathbb{T}_{\mathrm{B}}$ which relates the total polarization induced in B to incident fields in vacuum.

Given this and the fact that $|\mathbf{E}_{\mathrm{B}}\rangle=\mathbb{V}_{\mathrm{B}}^{-1}|\mathbf{P}_{\mathrm{B}}\rangle$ after setting $|\mathbf{P}^{(0)}_{\mathrm{B}}\rangle\to 0$, the scattered power only from body B (in the presence of body A) may be written as the difference between extinction and absorption powers only from body B (in the presence of body A), namely $\frac{\omega}{2}\left(\operatorname{Im}(\langle\mathbf{E}^{\mathrm{inc}},\mathbf{P}_{\mathrm{B}}\rangle)-\langle\mathbf{E}_{\mathrm{B}},\mathbf{P}_{\mathrm{B}}\rangle\right)=\frac{\omega}{2}\langle\mathbf{E}^{\mathrm{inc}},(\operatorname{Im}(\mathbb{Y}_{\mathrm{B}})-\mathbb{Y}_{\mathrm{B}}^{\star}\operatorname{Im}(\mathbb{V}_{\mathrm{B}}^{-1\star})\mathbb{Y}_{\mathrm{B}})\mathbf{E}^{\mathrm{inc}}\rangle$. Nonnegativity of this quantity for any $|\mathbf{P}^{(0)}_{\mathrm{A}}\rangle$, or more generally any $|\mathbf{E}^{\mathrm{inc}}\rangle$, means that $\operatorname{Im}(\mathbb{Y}_{\mathrm{B}})-\mathbb{Y}_{\mathrm{B}}^{\star}\operatorname{Im}(\mathbb{V}_{\mathrm{B}}^{-1\star})\mathbb{Y}_{\mathrm{B}}$ must be positive-semidefinite.

However, nonnegativity of far-field scattering from the system in general means that upon evaluating the inverse of (20), the operator $\operatorname{Im}(\mathbb{T})-\mathbb{T}^{\star}\operatorname{Im}(\mathbb{V}^{-1\star})\mathbb{T}$ must be positive-semidefinite, which means in turn that each of its diagonal blocks must be positive-semidefinite. Manipulating operators allows for showing that the bottom-right block is $\operatorname{Im}(\mathbb{Y}_{\mathrm{B}})-\mathbb{Y}_{\mathrm{B}}\left(\operatorname{Im}(\mathbb{V}_{\mathrm{B}}^{-1\star})+\mathbb{G}^{\mathrm{vac}}_{\mathrm{BA}}\mathbb{T}_{\mathrm{A}}\operatorname{Im}(\mathbb{V}_{\mathrm{A}}^{-1\star})\mathbb{T}_{\mathrm{A}}^{\star}\mathbb{G}_{\mathrm{AB}}^{\mathrm{vac}\star}\right)\mathbb{Y}_{\mathrm{B}}^{\star}$, which involves another positive-semidefinite operator (as $\operatorname{Im}(\mathbb{V}_{\mathrm{A}}^{-1\star})$ is positive-semidefinite and is multiplied on its left and right by operators which are Hermitian adjoints of each other, leaving the positive-semidefiniteness unaffected) subtracted from the operator $\operatorname{Im}(\mathbb{Y}_{\mathrm{B}})-\mathbb{Y}_{\mathrm{B}}^{\star}\operatorname{Im}(\mathbb{V}_{\mathrm{B}}^{-1\star})\mathbb{Y}_{\mathrm{B}}$. Therefore, the positive-semidefiniteness of $\operatorname{Im}(\mathbb{Y}_{\mathrm{B}})-\mathbb{Y}_{\mathrm{B}}\left(\operatorname{Im}(\mathbb{V}_{\mathrm{B}}^{-1\star})+\mathbb{G}^{\mathrm{vac}}_{\mathrm{BA}}\mathbb{T}_{\mathrm{A}}\operatorname{Im}(\mathbb{V}_{\mathrm{A}}^{-1\star})\mathbb{T}_{\mathrm{A}}^{\star}\mathbb{G}_{\mathrm{AB}}^{\mathrm{vac}\star}\right)\mathbb{Y}_{\mathrm{B}}^{\star}$ yields stronger bounds on the singular values of $\mathbb{Y}_{\mathrm{B}}$ than the positive-semidefiniteness of $\operatorname{Im}(\mathbb{Y}_{\mathrm{B}})-\mathbb{Y}_{\mathrm{B}}^{\star}\operatorname{Im}(\mathbb{V}_{\mathrm{B}}^{-1\star})\mathbb{Y}_{\mathrm{B}}$ alone, because the former also subtracts the absorption in body A in the presence of body B, whereas the latter does not.

Appendix E Alternative derivations of scalar approximation and Born

bound

In this section, we provide alternative derivations of the scalar approximation $\Phi_{\mathrm{sc}}$ and the Born bound $\Phi_{\mathrm{Born}}$. We start by writing the RHT spectrum in (2) as

[TABLE]

in terms of the scattering operator $\mathbb{S}^{\mathrm{A}}_{\mathrm{B}}$, after using the fact that $\operatorname{Asym}(\mathbb{V}_{p}^{-1\dagger})=\frac{1}{\zeta_{p}}\mathbb{I}_{p}$ for each body $p\in\{\mathrm{A},\mathrm{B}\}$. The scalar approximation follows from choosing the T-operator for each body to be $\mathbb{T}_{p}=\mathrm{i}\tau_{p}\mathbb{I}_{p}$ for some choice of uniform singular values $\tau_{p}\in[0,\zeta_{p}]$. This gives $\Phi_{\mathrm{sc}}=\frac{2\tau_{\mathrm{A}}^{2}\tau_{\mathrm{B}}^{2}}{\pi\zeta_{\mathrm{A}}\zeta_{\mathrm{B}}}\sum_{i}\frac{g_{i}^{2}}{(1+\tau_{\mathrm{A}}\tau_{\mathrm{B}}g_{i}^{2})^{2}}$, and it can be seen that the contribution from each channel $i$ monotonically increases with increasing $\tau_{\mathrm{A}}\tau_{\mathrm{B}}$. Thus, the only limit comes from the constraints on $\tau_{p}$, which therefore must saturate at $\tau_{p}=\zeta_{p}$. This is equivalent to directly plugging in $\mathbb{T}_{p}=\mathrm{i}\zeta_{p}\mathbb{I}_{p}$, which immediately yields $\Phi_{\mathrm{sc}}=\frac{2}{\pi}\zeta_{\mathrm{A}}\zeta_{\mathrm{B}}\left\lVert\mathbb{S}^{\mathrm{A}}_{\mathrm{B}}\mathbb{G}^{\mathrm{vac}}_{\mathrm{BA}}\right\rVert_{\mathrm{F}}^{2}$ where $\mathbb{S}^{\mathrm{A}}_{\mathrm{B}}=(\mathbb{I}_{\mathrm{B}}+\zeta_{\mathrm{A}}\zeta_{\mathrm{B}}\mathbb{G}^{\mathrm{vac}}_{\mathrm{BA}}\mathbb{G}^{\mathrm{vac}}_{\mathrm{AB}})^{-1}$ in the scalar approximation. Using the singular value decomposition of $\mathbb{G}^{\mathrm{vac}}_{\mathrm{BA}}$, particularly the facts that $\langle\mathbf{a}^{(j)},\mathbf{a}^{(j)\star}\rangle=\langle\mathbf{b}^{(j)},\mathbf{b}^{(j)\star}\rangle\mathrm{i}$ for each channel $j$, leads to the expression in the main text in terms of the singular values $\rho_{j}$ of $\mathbb{G}^{\mathrm{vac}}_{\mathrm{BA}}$.

The Born bound can be derived from $\Phi_{\mathrm{sc}}=\frac{2}{\pi}\zeta_{\mathrm{A}}\zeta_{\mathrm{B}}\left\lVert\mathbb{S}^{\mathrm{A}}_{\mathrm{B}}\mathbb{G}^{\mathrm{vac}}_{\mathrm{BA}}\right\rVert_{\mathrm{F}}^{2}$ by taking the scattering operator $\mathbb{S}^{\mathrm{A}}_{\mathrm{B}}\to\mathbb{I}_{\mathrm{B}}$, which is exactly a Born approximation. This exactly gives $\Phi_{\mathrm{Born}}=\frac{2}{\pi}\zeta_{\mathrm{A}}\zeta_{\mathrm{B}}\left\lVert\mathbb{G}^{\mathrm{vac}}_{\mathrm{BA}}\right\rVert_{\mathrm{F}}^{2}$, and substituting the singular values of $\mathbb{G}^{\mathrm{vac}}_{\mathrm{BA}}$ in the Frobenius norm squared gives the result in the main text. Alternatively, the operator $\mathbb{S}^{\mathrm{A}}_{\mathrm{B}}$ contributes the denominators $(1+\zeta_{\mathrm{A}}\zeta_{\mathrm{B}}g_{i})^{2}$ to each channel $i$, so the replacement $\mathbb{S}^{\mathrm{A}}_{\mathrm{B}}\to\mathbb{I}_{\mathrm{B}}$ makes those denominators disappear. In any case, the result in the main text is recovered by neglecting multiple scattering between the two bodies.

Appendix F Proof of domain monotonicity of singular values of

$\mathbb{G}^{\mathrm{vac}}_{\mathrm{BA}}$

In this section, we prove that the singular values $g_{i}$ of $\mathbb{G}^{\mathrm{vac}}_{\mathrm{BA}}$ are domain monotonic. The singular values of $\mathbb{G}^{\mathrm{vac}}_{\mathrm{BA}}$ are the eigenvalues of $\mathbb{G}^{\mathrm{vac}}_{\mathrm{BA}}(\mathbb{G}^{\mathrm{vac}}_{\mathrm{BA}})^{\dagger}$. We consider the effects of a perturbative addition of volume only to body A; a perturbative effect on body B can be considered through reciprocity, and the proof will remain the same. This allows for writing the block row vector of operators

[TABLE]

where $\mathbb{G}^{\mathrm{vac}}_{\mathrm{BA}_{0}}$ is the operator propagating fields in vacuum from the unperturbed volume $\mathrm{A}_{0}$ to body B, and $\mathbb{G}^{\mathrm{vac}}_{\mathrm{B}\Delta\mathrm{A}}$ is the operator propagating fields in vacuum from the perturbative volume $\Delta\mathrm{A}$ to body B. Using reciprocity, we may then write

[TABLE]

for which the first term $\mathbb{G}^{\mathrm{vac}}_{\mathrm{BA}_{0}}\mathbb{G}^{\mathrm{vac}\star}_{\mathrm{A}_{0}\mathrm{B}}$ is the Hermitian positive-semidefinite unperturbed operator, the terms $\mathbb{G}^{\mathrm{vac}}_{\mathrm{BA}_{0}}\mathbb{G}^{\mathrm{vac}\star}_{\Delta\mathrm{AB}}$ and $\mathbb{G}^{\mathrm{vac}}_{\mathrm{B}\Delta\mathrm{A}}\mathbb{G}^{\mathrm{vac}\star}_{\mathrm{A}_{0}\mathrm{B}}$ vanish because the projections onto the volume $\mathrm{A}_{0}$ and $\Delta\mathrm{A}$ are orthogonal to each other, and $\mathbb{G}^{\mathrm{vac}}_{\mathrm{B}\Delta\mathrm{A}}\mathbb{G}^{\mathrm{vac}\star}_{\Delta\mathrm{AB}}$ is the Hermitian positive-semidefinite perturbation. From Rayleigh-Schrödinger perturbation theory, if $\rho_{i}$ is an unperturbed singular value of $\mathbb{G}^{\mathrm{vac}}_{\mathrm{BA}_{0}}$ with $|\mathbf{b}_{i}\rangle$ being the corresponding normalized right singular vector, then the perturbation to $\rho_{i}$ is $\langle\mathbf{b}_{i},\mathbb{G}^{\mathrm{vac}}_{\mathrm{B}\Delta\mathrm{A}}\mathbb{G}^{\mathrm{vac}\star}_{\Delta\mathrm{AB}}\mathbf{b}_{i}\rangle$, which is nonnegative by virtue of the positive-semidefiniteness of $\mathbb{G}^{\mathrm{vac}}_{\mathrm{B}\Delta\mathrm{A}}\mathbb{G}^{\mathrm{vac}\star}_{\Delta\mathrm{AB}}$. Therefore, any increase in the volume of a body will increase the singular values of $\mathbb{G}^{\mathrm{vac}}_{\mathrm{BA}}$.

Appendix G Alternative bounds incorporating far-field radiative losses through

$\operatorname{Im}(\mathbb{G}^{\mathrm{vac}}_{\mathrm{AA}})$

In this section, we derive an alternative bound to (8) that involves $\operatorname{Im}(\mathbb{G}^{\mathrm{vac}}_{\mathrm{AA}})$, thus capturing constraints on scattering losses purely from finite object sizes rather than through multiple scattering. Starting from $\mathbb{T}_{\mathrm{A}}^{-1}=\mathbb{V}_{\mathrm{A}}^{-1}-\mathbb{G}^{\mathrm{vac}}_{\mathrm{AA}}$, operator manipulations allow for writing $\mathbb{T}_{\mathrm{A}}^{-1}\operatorname{Im}(\mathbb{V}_{\mathrm{A}}^{-1\star})\mathbb{T}_{\mathrm{A}}^{\star}=\operatorname{Im}(\mathbb{T}_{\mathrm{A}})-\mathbb{T}_{\mathrm{A}}\operatorname{Im}(\mathbb{G}^{\mathrm{vac}}_{\mathrm{AA}})\mathbb{T}_{\mathrm{A}}^{\star}$. This allows for rewriting (3) as

[TABLE]

which now hides reciprocity, as no similar transformation has been made to eliminate terms giving rise to $\zeta_{\mathrm{B}}$. Using the singular value decomposition from (6) for $\mathbb{Y}_{\mathrm{B}}$ and $\mathbb{G}^{\mathrm{vac}}_{\mathrm{BA}}$ but leaving $\mathbb{T}_{\mathrm{A}}$ general, we find that the constraint on $y_{i}$ is saturated when $y_{i}=\left(\zeta_{\mathrm{B}}^{-1}+\zeta_{\mathrm{A}}^{-1}g_{i}^{2}\langle\mathbf{a}_{i},\mathbb{T}_{\mathrm{A}}\mathbb{T}_{\mathrm{A}}^{\star}\mathbf{a}_{i}\rangle\right)^{-1}$. Completing the square, one may write

[TABLE]

This strongly suggests that the optimal $\mathbb{T}_{\mathrm{A}}$ should be diagonalized in the same basis as $\operatorname{Im}(\mathbb{G}^{\mathrm{vac}}_{\mathrm{AA}})$, so if we write

[TABLE]

then one may also write $\mathbb{T}_{\mathrm{A}}=\mathrm{i}\sum_{i}\tau_{i}|\mathbf{q}_{i}\rangle\langle\mathbf{q}_{i}|$. This implies that the constraint on the singular values of $\mathbb{Y}_{\mathrm{B}}$ becomes

[TABLE]

so $y_{i}$ depends on $\tau_{j}$ for every channel $j$, not just $j=i$. Consequently, we arrive at the following bound:

[TABLE]

for which finding the optimal values of $\tau_{i}$ for each channel $i$ requires self-consistently solving a large set of nonlinear equations subject to the constraint $\tau_{i}\leq\zeta_{\mathrm{A}}$ for each $i$. While this expression should yield tighter bounds on RHT owing to the incorporation of constraints on scattering losses for both objects in isolation (in addition to multiple scattering), it appears to be analytically intractable and must therefore be evaluated numerically, which we leave to future work; that said, numerical solution of the optimal values of $\tau_{i}$ for evaluating this bound requires solving a large set of nonlinear polynomial equations, which is generally computationally easier than brute-force optimization of the RHT spectrum due to the non-polynomial dependence of $\Phi$ on the T-operators in general. We point out that in the nonretarded quasistatic limit, $\operatorname{Im}(\mathbb{G}^{\mathrm{vac}}_{\mathrm{AA}})\to 0$, so all of its singular values may be taken to vanish as well; this means that its singular vectors become arbitrary, allowing for choosing $|\mathbf{q}_{i}\rangle=|\mathbf{a}_{i}\rangle$. Doing so, the above expression simplifies to $\Phi=\frac{2}{\pi}\sum_{i}\frac{\zeta_{\mathrm{B}}g_{i}^{2}\tau_{i}}{(1+\zeta_{\mathrm{A}}^{-1}\zeta_{\mathrm{B}}g_{i}^{2}\tau_{i}^{2})^{2}}$, from which it can be seen that if the material bound $\tau_{i}\leq\zeta_{\mathrm{A}}$ is saturated, each contribution is identical to the corresponding contribution from (8); if the material bound is not saturated, then the optimal $\tau_{i}=\frac{1}{\sqrt{3\zeta_{\mathrm{A}}^{-1}\zeta_{\mathrm{B}}}g_{i}}$ leads to a larger contribution than what we find in (8), yielding an overall looser bound. This corroborates the notion that our present bounds, which do not start with an explicit expression in terms of $\operatorname{Im}(\mathbb{G}^{\mathrm{vac}}_{\mathrm{AA}})$, do not fully account for dissipation via far-field radiative losses, while consideration of the singular values $\rho_{i}$ will remedy this.