$\mathbb{T}$-operator bounds on angle-integrated absorption and thermal radiation for arbitrary objects

TL;DR

This paper establishes fundamental bounds on the angle-integrated absorption and thermal radiation of arbitrary objects using the scattering T-operator, capturing the transition from volume to area scaling in different regimes.

Contribution

It introduces a unified theoretical framework for bounds on absorption and radiation that incorporates material, geometric, and passivity constraints using the T-operator.

Findings

Bounds accurately reflect the transition from volume to area scaling.

Comparison with prior limits shows improved accuracy.

Topology optimization results approach the theoretical bounds.

Abstract

We derive fundamental per-channel bounds on angle-integrated absorption and thermal radiation for arbitrary bodies---for any given material susceptibility and bounding region---that simultaneously encode both the per-volume limit on polarization set by passivity and geometric constraints on radiative efficiencies set by finite object sizes through the scattering $\mathbb{T}$-operator. We then analyze these bounds in two practical settings, comparing against prior limits as well as near optimal structures discovered through topology optimization. Principally, we show that the bounds properly capture the physically observed transition from the volume scaling of absorptivity seen in deeply subwavelength objects (nanoparticle radius or thin film thickness) to the area scaling of absorptivity seen in ray optics (blackbody limits).

Tables

Table 1. Table 1: Comparison of 𝕋 𝕋 \mathbb{T} response values for exemplary structure. Comparison of the values of | ⟨ q i | 𝕋 | q i ⟩ | quantum-operator-product subscript q 𝑖 𝕋 subscript q 𝑖 \left|\langle\textbf{q}_{i}\left|\mathbb{T}\right|\textbf{q}_{i}\rangle\right| (where | q i ⟩ ket subscript q 𝑖 |\textbf{q}_{i}\rangle are the singular vectors of Im ⁡ [ 𝔾 vac ] Im superscript 𝔾 vac \operatorname{Im}\left[\mathbb{G}^{\text{vac}}\right] ) for the structure of the rightmost green square, χ = 20 + 4 ​ i 𝜒 20 4 𝑖 \chi=20+4i ( ζ = 104 𝜁 104 \zeta=104 ), of Fig. 2 to the magnitude of the ideal τ i subscript 𝜏 𝑖 \tau_{i} values set by Φ opt subscript Φ opt \Phi_{\text{opt}} . Although the two set of values agree to great extent, pronounced differences between the ideal response dictated by the bound and the actual 𝕋 𝕋 \mathbb{T} -operator of the optimized structure are seen for higher spherical harmonics ℓ ℓ \ell -numbers when the magnitude of ideal value of ⟨ q i | 𝕋 | q i ⟩ quantum-operator-product subscript q 𝑖 𝕋 subscript q 𝑖 \langle\textbf{q}_{i}\left|\mathbb{T}\right|\textbf{q}_{i}\rangle approaches ζ 𝜁 \zeta .

$\ell$-Type	Ideal Value	Average Ratio	Weight
$1-1$	$0.32$	$0.94$	$0.05$
$1-2$	$0.40$	$1.32$	$0.05$
$2-1$	$0.81$	$1.29$	$0.08$
$2-2$	$0.50$	$1.19$	$0.08$
$3-1$	$4.18$	$0.83$	$0.11$
$3-2$	$1.40$	$0.67$	$0.11$
$4-1$	$36.25$	$0.53$	$0.13$
$4-2$	$7.48$	$0.60$	$0.13$
$5-2$	$66.23$	$0.04$	$0.17$

Equations

$$P_{\text{sct}}\left(\omega \right)=\frac{k_o}{2Z}⟨{\text{E}}_{\text{inc}}|{\mathbb{T}}^{†}\mathrm{Im}\left[{\mathbb{G}}^{\mathrm{vac}}\right]\mathbb{T}|{\text{E}}_{\text{inc}}⟩$$

$$=\frac{k_o}{2Z}⟨{\text{E}}_{\text{inc}}|\mathrm{Im}\left[\mathbb{T}\right]-{\mathbb{T}}^{†}\mathrm{Im}\left[{\mathbb{V}}^{-1†}\right]\mathbb{T}|{\text{E}}_{\text{inc}}⟩,$$

$$H\left(\omega ,T\right)=\Pi \left(\omega ,T\right)\Phi \left(\omega \right)$$

$$\Phi =\frac{2}{\pi }\mathrm{Tr}\left[\mathrm{Im}\left[{\mathbb{G}}^{\mathrm{vac}}\right]\left(\mathrm{Im}\left[\mathbb{T}\right]-{\mathbb{T}}^{†}\mathrm{Im}\left[{\mathbb{G}}^{\mathrm{vac}}\right]\mathbb{T}\right)\right].$$

$$\Vert \mathbb{T}\Vert \le \zeta =\frac{{\left|\chi \left(\omega \right)\right|}^2}{\text{Im}\left[\chi \left(\omega \right)\right]},$$

$$\mathrm{Im}\left[{\mathbb{G}}^{\mathrm{vac}}\right]=\sum _i{\rho }_i\left|{\text{q}}_i⟩⟨{\text{q}}_i\right|,$$

$${\Phi }_{\text{opt}}=\frac{2}{\pi }\sum _i{\tau }_i{\rho }_i-{\left({\tau }_i{\rho }_i\right)}^2$$

$${\Phi }_{\text{opt}}$$

$${\Phi }_{\text{qs}}\left(\omega \right)=\sum _i\zeta {\rho }_i=\zeta \underset{V}{\int }d\text{r}\mathrm{Im}\left[{\mathbb{G}}^{\mathrm{vac}}\left(\text{r},\text{r}\right)\right]$$

$$\displaystyle\operatorname{Im}\left[\mathbb{G}^{\mathrm{vac}}\right]\left(\textbf{r}\right)=\frac{k_{o}^{3}}{4\pi r}\Bigg{[}\left(\text{sin}\left(r\right)+\frac{\text{cos}\left(r\right)}{r}-\frac{\text{sin}\left(r\right)}{r^{2}}\right)\overline{\text{I}}-$$

$$\displaystyle\left(\text{sin}\left(r\right)+\frac{3~{}\text{cos}\left(r\right)}{r}-\frac{3~{}\text{sin}\left(r\right)}{r^{2}}\right)\hat{\textbf{r}}\otimes\hat{\textbf{r}}\Bigg{]},$$

$${\rho }_{\ell }^{\left(1\right)}$$

$$\displaystyle\hskip 14.45377pt+\frac{\ell}{2\ell+1}\left(J_{\ell+\frac{3}{2}}^{2}\left(R\right)-J_{\ell+\frac{1}{2}}\left(R\right)J_{\ell+\frac{5}{2}}\left(R\right)\right)\Bigg{)},$$

$${\rho }_{\ell }^{\left(2\right)}$$

$$\text{ln}\left(\frac{\zeta }{2}\right)\ge \left(2\ell +1\right)\text{ln}\left(\frac{\ell }{e\pi R}\right),$$

$$\text{E}={\text{E}}_{\text{inc}}+{\mathbb{G}}^{\mathrm{vac}}\mathbb{V}\text{E},$$

$${\mathbb{T}}^{-1}={\mathbb{V}}^{-1}-{\mathbb{G}}^{\mathrm{vac}},\text{J}=\mathbb{T}{\mathbb{V}}^{-1}{\text{J}}_{\text{inc}},$$

$$\displaystyle\Phi=\frac{2}{\pi}\operatorname{\text{Tr}}\Big{\{}$$

$$\displaystyle\left(\left(\mathbb{G}_{ba}^{\mathrm{vac}}\mathbb{T}_{a}\right)\operatorname{Im}\left[\mathbb{V}_{a}^{-1*}\right]\left(\mathbb{G}_{ba}^{\mathrm{vac}}\mathbb{T}_{a}\right)^{\dagger}\right)\Big{\}}$$

$${\mathbb{M}}_a={\left({\mathbb{I}}_a-{\mathbb{T}}_a{\mathbb{G}}_{ab}^{\mathrm{vac}}{\mathbb{T}}_b{\mathbb{G}}_{ba}^{\mathrm{vac}}\right)}^{-1}$$

$${\mathbb{M}}_b={\left({\mathbb{I}}_b-{\mathbb{T}}_b{\mathbb{G}}_{ba}^{\mathrm{vac}}{\mathbb{T}}_a{\mathbb{G}}_{ab}^{\mathrm{vac}}\right)}^{-1}.$$

$${\mathbb{O}}_b={\mathbb{M}}_b{\mathbb{W}}_b\mathrm{Im}\left[{\mathbb{V}}_b\right]{\mathbb{W}}_b^{†}{\mathbb{M}}_b^{†}$$

$${\mathbb{O}}_a={\mathbb{G}}_{ba}^{\mathrm{vac}}\left(\mathrm{Im}\left[{\mathbb{T}}_a\right]-{\mathbb{T}}_a^{\ast }\mathrm{Im}\left[{\mathbb{G}}_a^{\mathrm{vac}}\right]{\mathbb{T}}_a\right){\mathbb{G}}_{ab}^{{}^{\mathrm{vac}}\ast },$$

$${\mathbb{O}}_b\to \mathrm{Asym}\left\{{\mathbb{V}}_b\right\}$$

$$\underset{V_b{V}_c}{\iint }{\mathbb{G}}^{\ast }\left({\text{r}}_a,{\text{r}}_b\right)\mathrm{Im}\left[\mathbb{V}\right]\left({\text{r}}_b,{\text{r}}_c\right)\mathbb{G}\left({\text{r}}_c,{\text{r}}_d\right)=\mathrm{Im}\left[\mathbb{G}\right]\left({\text{r}}_a,{\text{r}}_d\right),$$

$${\mathbb{G}}^{\ast }\mathrm{Im}\left[\mathbb{V}\right]\mathbb{G}=\mathrm{Im}\left[\mathbb{G}\right].$$

$$\underset{{\mathbb{V}}_b\to 0}{\lim }{\mathbb{G}}_{ab}^{\ast \mathrm{vac}}\mathrm{Im}\left[{\mathbb{V}}_b\right]{\mathbb{G}}_{ba}^{\mathrm{vac}}=\mathrm{Im}\left[{\mathbb{G}}_a^{\mathrm{vac}}\right].$$

$$\displaystyle\Phi=\frac{2}{\pi}\operatorname{\text{Tr}}\Bigg{\{}\mathbb{G}_{ab}^{\mathrm{vac}*}\operatorname{Im}\left[\mathbb{V}_{b}\right]\mathbb{G}_{ba}^{\mathrm{vac}}\left(\operatorname{Im}\left[\mathbb{T}_{a}\right]-\mathbb{T}_{a}^{*}~{}\operatorname{Im}\left[\mathbb{G}_{a}^{\mathrm{vac}}\right]\mathbb{T}_{a}\right)\Bigg{\}}$$

$$\displaystyle=\frac{2}{\pi}\Bigg{(}\operatorname{\text{Tr}}\Bigg{\{}\operatorname{Im}\left[\mathbb{G}_{a}^{\mathrm{vac}}\right]\left(\operatorname{Im}\left[\mathbb{T}_{a}\right]-\mathbb{T}_{a}^{*}~{}\operatorname{Im}\left[\mathbb{G}_{a}^{\mathrm{vac}}\right]\mathbb{T}_{a}\right)\Bigg{\}}-\lim_{\mathbb{V}_{a}\rightarrow 0}$$

$$\displaystyle\operatorname{Tr}\left\{\mathbb{G}_{a}^{\mathrm{vac}*}\operatorname{Im}\left[\mathbb{V}_{a}\right]\mathbb{G}_{a}^{\mathrm{vac}}\left(\operatorname{Im}\left[\mathbb{T}_{a}\right]-\mathbb{T}_{a}^{*}~{}\operatorname{Im}\left[\mathbb{G}_{a}^{\mathrm{vac}}\right]\mathbb{T}_{a}\right)\right\}\Bigg{)}$$

$$\displaystyle=\frac{2}{\pi}\operatorname{\text{Tr}}\Bigg{\{}\operatorname{Im}\left[\mathbb{G}_{a}^{\mathrm{vac}}\right]\left(\operatorname{Im}\left[\mathbb{T}_{a}\right]-\mathbb{T}_{a}^{*}~{}\operatorname{Im}\left[\mathbb{G}_{a}^{\mathrm{vac}}\right]\mathbb{T}_{a}\right)\Bigg{\}},$$

$$\mathrm{Im}\left[{\mathbb{G}}_{{}_U}^{\mathrm{vac}}\right]=\sum _i{u}_i\left|{\text{u}}_i⟩⟨{\text{u}}_i\right|$$

$$\mathrm{Im}\left[{\mathbb{G}}_{{}_D}^{\mathrm{vac}}\right]=\sum _i{d}_i\left|{\text{d}}_i⟩⟨{\text{d}}_i\right|$$

$$\underset{V}{\iint }d\text{r}d{\text{r}}^{\prime }{\text{f}}^{\ast }\left(\text{r}\right)\mathrm{Im}\left[{\mathbb{G}}^{\mathrm{vac}}\right]\left(\text{r},{\text{r}}^{\prime }\right)\text{f}\left({\text{r}}^{\prime }\right)$$

$$\begin{gathered} \mathrm{Im}\left[{\mathbb{G}}^{\mathrm{vac}}\right]\left({\text{r}}_a,{\text{r}}_a^{\prime }\right)=k_o^3\sum _{\ell =1}^{\infty }\sum _{m=-l}^{\ell }{\left(-1\right)}^m \\ \sum _{j=1}^2{\text{rS}}_{\ell ,m}^{\left(j\right)}\left(r,\theta ,\varphi \right)\otimes {\text{rS}}_{\ell ,-m}^{\left(j\right)}\left(r^{\prime },{\theta }^{\prime },{\varphi }^{\prime }\right) \end{gathered}$$

$${\text{rS}}_{\ell m}^{\left(1\right)}\left(r,\theta ,\varphi \right)=\sqrt{\frac{2\ell +1}{4\pi \ell \left(\ell +1\right)}}{j}_{\ell }\left(r\right){\text{V}}_{\ell m}^{\left(3\right)}\left(\theta ,\varphi \right)$$

Full text

Motivated by the increasing control of light offered by micro and nanoscale structuring Koenderink et al. (2015); Molesky et al. (2018), impetus to find bounds analogous to the blackbody limit for geometries that violate the assumptions of ray optics (nanoparticles Du and Yu (2017), thin films Tong et al. (2015), photonic crystals Zhu et al. (2015); Ilic et al. (2016), etc.) has steadily grown over the past few decades. It is now well established that the absorption (radiative thermal emission) cross-sections of a compact object can be much greater than its geometric area Bohren and Huffman (2008); Tribelsky (2011); Ruan and Fan (2011); Biehs and Ben-Abdallah (2016); Fernández-Hurtado et al. (2018); Thompson et al. (2018) (“super-Planckian” emission), and that deeply subwavelength films can achieve near unity absorptivity via surface texturing Kats and Capasso (2016); Dyachenko et al. (2016). Limits applicable to all length scales and materials could both provide insight into these representative phenomena and guide efforts in related application areas such as integrated and meta-optics Kruk and Kivshar (2017); Khorasaninejad and Capasso (2017); Hampson et al. (2018), photovoltaics Yablonovitch (1982); Atwater and Polman (2010); Zhang et al. (2016); Jariwala et al. (2017), and photon sources Thompson et al. (2013); Galfsky et al. (2015); Somaschi et al. (2016).

Development of bounds for arbitrary objects have primarily followed two overarching strategies: modal decompositions based on quasi-normal, singular value, Fourier and/or multipole expansions McLean (1996); Hamam et al. (2007); Kwon and Pozar (2009); Yu et al. (2010a, b); Ruan and Fan (2010); Hugonin et al. (2015); Jia et al. (2015); Yang et al. (2017), relating absorption cross-section to the number of excitable optical modes (channels); or material/region bounds, utilizing energy Callahan et al. (2012); Miller et al. (2016) and/or spectral sum rules Fuchs and Liu (1976); Miller (2000, 2007); Miller et al. (2014); Shim et al. (2019) to constrain achievable properties. Separately, each of these approaches present challenges for photonic design. Modal decompositions incorporate the specific size and shape characteristics of a body through expansion coefficients, and hence, inherently, require some enumeration and characterization of the participating modes to determine the range of values these coefficients can take Hamam et al. (2007); Ruan and Fan (2012); Alpeggiani et al. (2017). Although fundamental considerations (transparency, energy, size, etc.) can and have been used in this regard Miller (2000); Sohl et al. (2007); Yu et al. (2010a); Yang et al. (2017), such cut-offs have yet to tightly bound potential coefficient values for arbitrary compact geometries, particularly when applied to metallic nanoparticles and antennas Gustafsson et al. (2007); Kwon and Pozar (2009); Miller et al. (2016). Conversely, material bounds set by intrinsic dissipation naturally reproduce the volumetric scaling of absorptivity characteristic of deeply subwavelength objects (and are highly accurate for the special case of weak polarizability in this regime Miller et al. (2016)). However, because such approaches intrinsically suppose an optimally large response field existing at all points within an arbitrary object for any incident field, the same volumetric scaling persists for all length scales. Consequently, material bounds can rapidly become too loose beyond quasi-static settings, yielding unphysical divergences with both increasing object size and material response.

In this letter, we derive bounds on thermal radiation and absorption that combine these two approaches, linking the impact of material response with the influence of an object’s geometry through the scattering $\mathbb{T}$-operator. This leads to a per-channel limit on integrated absorption capturing both extraction and radiative (scattering) loss processes through the singular values of the imaginary part of the vacuum Green function. The result is applicable to objects of any size, exhibiting a smooth transition in absorptivity from the volume scaling achievable in the quasi-static (deeply-subwavelength) regime to the area scaling limit of macroscopic ray optics. Further, the bounds always asymptotically approach the ray optics limit (when all characteristic lengths are large) and diverge sub-logarithmically (rather than linearly) with material quality for objects of finite extent, significantly reducing cross-section limits for typical optical media even when all characteristic lengths are small. Throughout, we compare the present results to prior bounds as well as structures discovered using topology optimization, realizing a variety of examples (metallic and dielectric) that nearly achieve the predicted limits.

**Derivation—**From the relations of scattering theory, both the power scattered from an incident field ($\left|\textbf{E}_{\text{inc}}\right>$), and the thermal radiation emitted at temperature $T$, can be expressed in terms of the scattering $\mathbb{T}$-operator of an object and the vacuum Green function $\mathbb{G}^{\mathrm{vac}}$ Krüger et al. (2012) as

[TABLE]

and,

[TABLE]

Here, $\omega$ is the angular frequency, $k_{o}=2\pi/\lambda$ is the wavenumber, $Z$ is the impedance of free-space, $\Pi\left(\omega\right)=\hbar\omega/\left(\text{Exp}\left(\beta\hbar\omega\right)-1\right)$ with $\left(\beta=1/k_{\text{B}}T\right)$ is the Planck energy of a harmonic oscillator, $\text{Tr}\left[\ldots\right]$ denotes the trace, $\operatorname{Im}\left[\mathbb{T}\right]=\left(\mathbb{T}-\mathbb{T}^{*}\right)/2i$, and, by Kirchhoff’s law of thermal radiation, $\Phi$ is the object’s angle integrated absorption Greffet and Nieto-Vesperinas (1998). (A synopsis of scattering formalism, along with a derivation of (2), is provided in Supplemental Material Molesky et al. (2019a).) For a passive object, scattered power must be positive for any incident field. As such, (1) simultaneously dictates that all singular values of the $\mathbb{T}$-operator must be smaller than the material figure of merit $\zeta$

[TABLE]

which was similarly derived in Ref. Miller et al. (2016) for polarization fields, and that $\operatorname{Im}\left[\mathbb{T}\right]$ is positive-definite.

As $\operatorname{Im}\left[\mathbb{G}^{\mathrm{vac}}\right]$ is real-symmetric positive-definite, it can be expressed via a singular value decomposition as

[TABLE]

where, as supported by our later analysis, each $\rho_{i}$ (eigenvalue) can be equated to the outgoing radiative flux of the $i$th mode—the $i$th radiative efficacy of the domain. (The set $\left\{\rho_{i}\right\}$ plays an analogous role to the coupling coefficients $\left\{g_{ij}\right\}$ used by D.A.B. Miller in setting limits on far-field optical communication Miller (2000).) Consider (2) using this expansion, $\Phi=\left(2/\pi\right)\sum_{i}\rho_{i}\operatorname{Im}\left[\langle\textbf{q}_{i}\left|\mathbb{T}\right|\textbf{q}_{i}\rangle\right]-\rho_{i}^{2}\left|\langle\textbf{q}_{i}\left|\mathbb{T}\right|\textbf{q}_{i}\rangle\right|^{2}-$ $\left(2/\pi\right)\sum_{\left\{\left(i,j\right)|i\neq j\right\}}\rho_{i}\rho_{j}\left|\langle\textbf{q}_{i}\left|\mathbb{T}\right|\textbf{q}_{j}\rangle\right|^{2}.$ Now, take $\mathbb{T}_{\text{opt}}$ to be a general operator described by the properties $(\mathbb{T}_{\text{opt}})^{\text{T}}=\mathbb{T}_{\text{opt}}$ (reciprocity), $\left\lVert\mathbb{T}_{\text{opt}}\right\rVert\leq\zeta$ (passivity), and $\operatorname{Im}\left[\mathbb{T}_{\text{opt}}\right]$ positive-definite (passivity), ignoring all other physical constraints that any true $\mathbb{T}$-operator must satisfy. In this context, two characteristics of any maxima of $\Phi$ are clear. First, as $\left(\forall~{}i,j\right)~{}\rho_{i}\rho_{j}\left|\langle\textbf{q}_{i}\left|\mathbb{T}_{\text{opt}}\right|\textbf{q}_{j}\rangle\right|^{2}\geq 0,$ the appearance of any cross-terms ($\langle\textbf{q}_{i}\left|\mathbb{T}_{\text{opt}}\right|\textbf{q}_{j}\rangle$) will always decrease $\Phi$. Therefore, to maximize $\Phi$, a general operator $\mathbb{T}_{\text{opt}}$ must be diagonalized in the basis of $\operatorname{Im}\left[\mathbb{G}^{\mathrm{vac}}\right]$, (4). Second, the complex phase of $\left<\textbf{q}_{i}\left|\mathbb{T}_{\text{opt}}\right|\textbf{q}_{i}\right>$ only influences the first (positive) piece of the sum, and so the value of $\Phi$ peaks when $\left(\forall i\right)~{}\text{atan}\left(\operatorname{Im}\left[\langle\textbf{q}_{i}\left|\mathbb{T}\right|\textbf{q}_{i}\rangle\right]/\operatorname{Re}\left[\langle\textbf{q}_{i}\left|\mathbb{T}\right|\textbf{q}_{i}\rangle\right]\right)=\pi/2.$ Together, these two considerations show that achievable values of $\Phi$ are bounded by taking $\mathbb{T}_{\text{opt}}$ to be diagonalized by (4) with purely imaginary eigenvalues: $\mathbb{T}_{\text{opt}}=\sum_{i}i\tau_{i}\left|\textbf{q}_{i}\rangle\langle\textbf{q}_{i}\right|$ with $\left(\forall i\right)~{}\tau_{i}\in[0,\zeta]$. As such,

[TABLE]

and maximizing the contribution of each $\tau_{i}$ yields

[TABLE]

That is, based on the criterion $\zeta\rho_{i}\geq 1/2$, each channel in (6) produces either the Landauer limited contribution of $1/4$ Molesky et al. (2019b), or the material limited $\zeta\rho_{i}-\left(\zeta\rho_{i}\right)^{2}$.

**Interpretation—**In terms of the $\mathbb{T}$ operator, the total power extracted from any incident field $\left|\textbf{E}_{\text{inc}}\right>$ by an object is $P_{\text{ext}}=k_{o}\left<\textbf{E}_{\text{inc}}\right|\operatorname{Im}\left[\mathbb{T}\right]\left|\textbf{E}_{\text{inc}}\right>/\left(2Z\right)$. Comparing with (1) and (2), $\Phi$ thus amounts to the difference of the extracted ($\operatorname{Im}\left[\mathbb{T}\right]$) and scattered ($\mathbb{T}^{\dagger}\operatorname{Im}\left[\mathbb{G}^{\mathrm{vac}}\right]\mathbb{T}$) power for free-space states. The separation of these two forms persists throughout the derivation of the bounds, representing the linear and quadratic terms of (5). $\Phi_{\text{opt}}$ results from their connected physics.

In real space, $\operatorname{Tr}\left\{\operatorname{Im}\left[\mathbb{G}^{\mathrm{vac}}\right]\right\}=\sum_{i}\rho_{i}$ is the integral of the local density of free-space states over the domain of the object. Following (6), the total power that can be extracted by an object, the first term of (2), is thus bounded by its ability to interact with radiative modes, $\operatorname{Tr}\left\{\operatorname{Im}\left[\mathbb{T}\right]\operatorname{Im}\left[\mathbb{G}^{\mathrm{vac}}\right]\right\}=\sum_{i}\tau_{i}\rho_{i}$, which is maximized (independently) under complete saturation of material response, $\left(\forall i\right)~{}\tau_{i}=\zeta$. Relatedly, this form is also the result of applying the per-volume (shape independent) optical response limit of Ref. Miller et al. (2016) to integrated absorption, and is similar to the light trapping bound of Ref. Callahan et al. (2012). Due to these connections with prior work,

[TABLE]

serves as a useful comparison for $\Phi_{\text{opt}}$, and is subsequently referred to as the quasi-static bound. This name is chosen as (7) follows from the assumption that the interaction of the object with any incident field is identically material limited, which can occur in quasi-static settings. This does not mean that $\Phi_{\text{qs}}$ is valid only under the quasi-static approximation. Like $\Phi_{\text{opt}}$, $\Phi_{\text{qs}}$ is a mathematical bound derived from Maxwell’s equations, albeit for any selection of parameters $\Phi_{\text{opt}}\leq\Phi_{\text{qs}}.$

In (5) and (6) this extracted power contribution is suppressed by scattering (radiative) losses, which are captured in the quadratic term in (2) as the coupling of the polarization currents generated within an object back to free-space modes: originating through the operator $\mathbb{T}^{\dagger}\operatorname{Im}\left[\mathbb{G}^{\mathrm{vac}}\right]\mathbb{T}$, each $\tau_{i}$ represents the ability of the object to convert a given field into a current, and each radiative efficacy $\rho_{i}$ the conversion of a current into outgoing radiative flux. Equivalently, the presence of strong polarization currents, necessary for strong per-volume absorption, leads to radiative loses, and these loses limit possible absorption. If $\zeta\rho_{i}>1/2$, mirroring the observed dependencies of absorption ($\propto V$) and scattering ($\propto V^{2}$) seen in highly subwavelength metallic antennas Hamam et al. (2007); Ruan and Fan (2010), the growth of radiative losses with increasing $\tau_{i}$ can potentially surpass the growth of the extracted power, inducing saturation. As both processes are rooted in the same conversion between radiative fields and polarization currents, this critical coupling occurs at the compelling value of $\tau_{i}\rho_{i}=1/2$ Pendry (1983, 1999), the probability of a maximally entropic Bernoulli process, resulting in the Landauer limit value of $\Phi_{\text{opt}}=1/4$.

**Analysis—**The practical usefulness of (6) stems from its favorable mathematical properties. Namely, (6) monotonically increases with $\zeta$ or any $\rho_{i}$, and, as proved in Supplementary Material Molesky et al. (2019a), each $\rho_{i}$ increases if the object grows (domain monotonicity). This allows us to freely decouple any true object from an imagined encompassing region of space (bounding domain). A mismatch between the domain of the object and the domain of $\operatorname{Im}\left[\mathbb{G}^{\mathrm{vac}}\right]$ must technically reduce $\left\lVert\mathbb{T}\right\rVert$ below $\zeta$, but without any modification (6) remains an upper bound on $\Phi$. That is, the result of (6) for any particular bounding domain is applicable to any object that can be enclosed (as well as any sub-domain).

The procedure for calculating $\Phi_{\text{opt}}$ is straightforward for any bounding geometry (e.g. wires, disks, spheres, extended films, stars, disconnected patches, etc.). Precisely, the set of singular values $\left\{\rho_{i}\right\}$ of the domain can always be computed by forming a real space matrix representation of $\operatorname{Im}\left[\mathbb{G}^{\mathrm{vac}}\right]$,

[TABLE]

with every $r$ multiplied by a hidden $k_{o}$, and then performing a singular value decomposition of the result Polimeridis et al. (2013, 2015). Here, to facilitate further investigation, we will focus on the high symmetry case of a ball where semi-analytic evaluation is manageable (expressions for films, as well as minor additional details, are given in Supplemental Material Molesky et al. (2019a)). Nevertheless, we stress that determining $\Phi_{\text{opt}}$ for domains lacking symmetry does not raise any meaningful computationally difficulties.

For this geometry two types of singular values arise

[TABLE]

where $J_{\ell}\left(-\right)$ is the $\ell$th Bessel function of the first kind with an additional factor of $2\pi$ included in its argument, each $\ell$ (spherical harmonic) index has a multiplicity of $\left(2\ell+1\right)$, and $R$ is the radius of the ball normalized by the wavelength. Using standard properties of Bessel functions, it can be shown that for values of $R\gg\ell$, each of these singular values tends to the asymptote $2\pi^{2}R$, and that for any combination of arguments $\rho^{\left(1\right)}_{\ell}<\pi\left(\ell+1\right)\left(\pi R\right)^{2\ell+1}/\left(2\Gamma^{2}\left(\ell+3/2\right)\right)+2\pi\ell\left(\pi R\right)^{2\ell+5}/\left(\left(2\ell+5\right)\left(2\ell+3\right)\Gamma^{2}\left(\ell+5/2\right)\right)$ and $\rho^{\left(2\right)}_{\ell}<2\pi\left(\pi R\right)^{2\ell+3}/\left(\left(2\ell+3\right)\Gamma^{2}\left(\ell+3/2\right)\right)$ (asymptotically approached for small values of $R$). These forms reveal two prescient general features. First, in the limit of small domains ($R\ll 1$), with “small” being determined by the value of $\zeta$, only the first singular value of the first type contributes, and this triply degenerate (dipole) mode is responsible for the initial volume scaling necessitated by the physical meaning of the bounds. Second, the radial growth of the singular values shows that the saturation condition (impact of radiative losses) plays a major role in limiting radiative thermal emission and integrated-absorption in wavelength scale volumes. (For $\zeta=10^{6}$, Fig. 1 (a), radiative losses lead to order of magnitude deviations of $\Phi_{\text{opt}}$ from $\Phi_{\text{qs}}$ beyond $R\approx 0.003\lambda$.) As visually confirmed by Fig. 1 panel (a), as the domain grows an increasing number of channels (multipoles) saturate causing “steps” to appear in $\Phi_{\text{opt}}$, and these steps lead to successively larger deviation with $\Phi_{\text{qs}}$ that ultimately regularize the initial volumetric scaling. Results for films, Fig. 1 (b), are qualitatively similar. However, since the domain is infinite, the steps associated with saturation are now blended into a continuum, and the large characteristic size limit is approached from below rather than above. From a practical perspective, the fact that $\Phi_{\text{opt}}$ can achieve near ideal absorptivity for very small film thickness and moderate values of $\zeta$ is quite remarkable, a finding that is tacitly supported by a number of recent studies in 2D materials and meta-surfaces Thongrattanasiri et al. (2012); Akselrod et al. (2015); Kim et al. (2018); Nong et al. (2018). Crucially, in either case, for any value of $\zeta$, $\Phi_{\text{opt}}$ asymptotes to a geometric perfect absorber (the blackbody limit).

The asymptotic behavior of the singular values also reveals general characteristics of the dependence of $\Phi_{\text{opt}}$ on the material figure of merit $\zeta$. Applying Sterling’s approximation to the bounding expressions given above, for $\left(\ell\gg e\pi R\right)$ we have $\rho^{\left(2\right)}_{\ell}\approx\left(e\pi R/\ell\right)^{2\ell+1}/4$ and $\rho^{\left(1\right)}_{\ell}\approx\left(e\pi R/\ell\right)^{2\ell+3}/2$, to arbitrary accuracy as $\ell$ becomes large. Fix $R$, and suppose that $\zeta=\rho^{\left(2\right)}_{k}$ ($\rho^{\left(1\right)}_{k}$ is analogous). Using the fact that $e\pi R/\left(k+\ell\right)<e\pi R/k$ the remaining (unsaturated) linear contribution of $\Phi_{\text{opt}}$ is then bounded by $9\left(e\pi R\right)^{3}/\left(4\left(k^{2}-\left(e\pi R\right)^{2}\right)\right)$. Hence, as $\zeta$ saturates increasingly higher spherical harmonics, the contribution of the remaining unsaturated harmonics becomes increasingly small compared to the contribution of the newly saturated harmonic, $\approx\left(2k+1\right)/4$. But, saturation of the $\ell$th singular value (in the large $\ell$ limit) requires

[TABLE]

which has a sub-logarithmic dependence between $\ell$ and $\zeta$. Due to domain monotonicity, the above material scaling result for a ball is applicable to all compact (finite sized) objects.

This bound on material quality scaling is well matched to the features of the $\Phi_{\text{opt}}$ curves in Fig. 1 panel (a). Once the radius has surpassed $\approx\lambda$, geometric increases in $\zeta$ ($\times 10^{2}$) produce relatively minute changes in the bounds. This behavior also appears for smaller radii at larger values of $\zeta$, but this range is not of great practical relevance since materials with $\zeta$ surpassing $\approx 10^{8}$ are quite rare. For instance, in the optical to infrared, $\omega\in(0.5\text{--}15)\mu\text{m}$, $\zeta\left(\omega\right)$ has a peak value of approximately $1.7\times 10^{3}$ for gold, $2.4\times 10^{3}$ for tungsten, $2.2\times 10^{4}$ for silicon carbide, $6.8\times 10^{9}$ for silicon, $3.3\times 10^{7}$ for gallium arsenide, and $5.9\times 10^{7}$ for gallium phosphide Palik (1998).

**Optimizations—**Case evidence for the tightness of (6) is presented in Fig. 2. Using a gradient topology optimization algorithm Jensen and Sigmund (2011); Molesky et al. (2018), see Supplemental Material for details Molesky et al. (2019a), structures nearly achieving $\Phi_{\text{opt}}$ have been discovered for two widely different domain sizes ($R=0.05\lambda$ and $R=0.5\lambda$) and a variety of metallic and dielectric susceptibilities. In Fig. 2, these media are grouped by imaginary susceptibility, corresponding to four different values of $\operatorname{Im}\left[\chi\right]$, $\left\{0.5,1,2,4\right\}$, with the remaining variation in $\zeta$ occurring due to $\operatorname{Re}\left[\chi\right]$. Explicit values of $\operatorname{Re}\left[\chi\right]$ are given for circled points, providing a sense of the range considered. As was previously remarked by O. Miller et al. Miller et al. (2016), $\Phi_{\text{qs}}$ is attained for a plane wave polarized along the axis of an ellipsoidal metallic nanoparticle, given a properly chosen aspect ratio. For small values of $\zeta$ this ratio is near unity and resonant metallic structures ($\operatorname{Re}\left[\chi\right]\approx-3$) matching both bounds are easily discovered. As $\zeta$ moves to moderate values, the aspect ratio required for an ellipsoidal particle to match $\Phi_{\text{qs}}$ becomes increasingly extreme. Due to our chosen spherical boundary, discovered structures begin to deviate considerably from $\Phi_{\text{qs}}$, but continue to come within a factor of $2$ of $\Phi_{\text{opt}}$ up to $\zeta=10^{3}$. Past this point, numerical issues impede our present algorithms and it remains to be seen how much of the roughly order of magnitude headroom allowed by $\Phi_{\text{opt}}$ is accessible.

Results for the larger domain, Fig. 2 (b), show similarly good agreement. An example structure is depicted in the right inset (full view and planar cut), corresponding to the rightmost green square in the plot. Comparing with the assumptions made in deriving (6), the $\mathbb{T}$ operator for this structure ($\chi=20+4i$, $\Phi=0.60~{}\Phi_{\text{opt}}$) is indeed found to be nearly diagonal in the basis of $\operatorname{Im}\left[\mathbb{G}^{\text{vac}}\right]$ and have almost completely imaginary eigenvalues (for supporting data see Supplemental Material Molesky et al. (2019a)).

Remarks—There are a few points that should be considered when using (6), or comparing to prior literature. First, $\Phi_{\text{opt}}$ is a bound on thermal emission and integrated absorption for a given domain and $\zeta$ factor. By choosing different geometries and material parameters, (6) can be applied to any desired context, but the confining volume is an essential feature. Second, there is no universal guarantee of tightness. Beyond the demonstrated agreement of the bounds with known quasi-static and ray optics asymptotics, the only a priori guarantee is domain monotonicity; there are likely volumes and material parameters where the value of $\Phi_{\text{opt}}$ will be larger than the true $\Phi$ of any practical structure. Next, while we have only considered single wavelengths, there is no reason the bounds can not be applied to finite frequency ranges. The derivation of $\Phi_{\text{opt}}$ presented above does not incorporate any spectral sum rules (derived from causality), such as the fact that $\mathbb{T}_{\text{opt}}$ should obey Kramers-Kronig dispersion relations, but for resonant absorption or thermal emission simply multiplying the bound by the width of the resonance should not produce a substantially looser bound than $\Phi_{\text{opt}}$ at the peak wavelength. (As an expedient, taking $\Phi_{\text{opt}}$ to be the peak value of a Lorentzian function of width $\Delta\omega=\omega~{}\mathrm{Im}\chi/|\chi|$ is likely a fair approximation.) Finally, as suggested in the introduction, $\Phi_{\text{opt}}$ can be interpreted as the extension of prior multipole analysis McLean (1996); Hamam et al. (2007); Kwon and Pozar (2009); Yu et al. (2010a, b); Ruan and Fan (2010); Hugonin et al. (2015); Jia et al. (2015); Yang et al. (2017), or communication limits Miller (2000, 2007), to general domains with the crucial addition that an upper bound is set on the number modes which may contribute through the pseudo-rank of the imaginary part of the vacuum Green function ($\operatorname{Im}\left[\mathbb{G}^{\mathrm{vac}}\right]$) and the material figure of merit ($\zeta$) (3). We foresee this rank revealing capability potentially providing a number of benefits for future practical design and optimization. We also note that much of what has been developed in this manuscript is applicable not only to generalized electromagnetic scattering (for incident planewaves or dipolar emitters with applications to solar cells, light-emitting diodes, and single-photon emitters), but also to quantum mechanics, acoustics, and other wave physics.

I Acknowledgments

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. DMR1719875), the Defense Advanced Research Projects Agency (DARPA) under agreement HR00111820046, and the National Science and Engineering Research Council of Canada under PDF-502958-2017. The views, opinions and/or findings expressed herein are those of the authors and should not be interpreted as representing the official views or policies of any institution. We thank Jason Necaise for performing instructive calculations of the bounds in cylindrical coordinates.

II Supplementary Information

$\mathbb{T}$**-operator definition—**Following the Lippmann-Schwinger approach to scattering Lippmann and Schwinger (1950), the $\mathbb{T}$-operator is formally defined through the self-consistent equation

[TABLE]

as

[TABLE]

with the “inc” subscript denoting an initial (freely incident) current (field), $\mathbb{G}^{\mathrm{vac}}$ the vacuum Green function, and $\mathbb{V}$ a generalized electromagnetic susceptibility. A lack of any subscript indicates total field quantities.

**Heat transfer to thermal emission—**In two recent (related) articles Jin et al. (2019); Molesky et al. (2019b) we have established that the heat transfer between any two bodies, labeled $a$ and $b$, can be written in terms of the scattering $\mathbb{T}$-operator as

[TABLE]

where $\mathbb{M}_{a}$ and $\mathbb{M}_{b}$ are mutual scattering operators defined as

[TABLE]

To produce an expression for thermal emission, we will evaluate (13) in the limit that one of the objects, here chosen to be $b$, tends to a blackbody. Abstractly, a blackbody is an encompassing, infinitely large, region capable of perfectly absorbing any incident field: physically, a spherical material shell of inner radius $r_{b}$ and outer radius $R_{b}$, in the simultaneous limit $r_{b}\rightarrow\infty$ and $\text{Im}~{}\chi_{b}\rightarrow 0$ with $(R_{b}-r_{b})\rightarrow\infty$ and $\text{Im}~{}\chi_{b}\left(R_{b}-r_{b}\right)\rightarrow 1$. We begin by breaking (13) into

[TABLE]

where $\mathbb{W}_{b}^{-1}=\mathbb{I}_{b}-\mathbb{V}_{b}\mathbb{G}_{b}$ is the current dressing operator (producing a total current from an initial current). As $b$ tends to the blackbody limit, the above definition shows that $\mathbb{W}_{b}\rightarrow\mathbb{I}_{b}$, i.e. that the total electric current density tends to the free current density. Hence, $\mathbb{T}_{b}=\mathbb{V}_{b}\mathbb{W}_{b}$ tends to zero, and the mutual scattering operators $\mathbb{M}_{a}$ and $\mathbb{M}_{b}$ become $\mathbb{I}_{a}$ and $\mathbb{I}_{b}$ respectively. Therefore,

[TABLE]

with $\operatorname{Im}\left[\mathbb{V}_{b}\right]\rightarrow 0$ and the volume of $b$ becoming arbitrarily large.

Now, for any collection of objects, the fluctuation dissipation theorem Eckhardt (1984) states that

[TABLE]

which in our notation becomes

[TABLE]

To apply this result to (13), we imagine a fictitious addition to body $b$ that perfectly overlaps with $a$. Since $\operatorname{Im}\left[\chi_{b}\right]\rightarrow 0$ and body $a$ has a finite extent, this addition will have no material effect on the value produced by (13). However, its inclusion shows that consistent application of the blackbody limit must result in

[TABLE]

Applying this result to (13), using the simplification of (16), we find that

[TABLE]

where in the middle expression $\mathbb{V}_{a}$ has been used for the blackbody limit inside body $a$.

Domain monotonicity– Let $U$ be a subdomain of $D$, and

[TABLE]

be the vacuum Green functions for each of these volumes. Assume that these sums are finite. The first singular vector, corresponding to the largest singular value, is equivalent to the complex vector field $\textbf{f}\left(\textbf{r}\right)$ maximizing the integral

[TABLE]

subject to the constraint, $\int\limits_{V}d\textbf{r}~{}\textbf{f}^{*}\left(\textbf{r}\right)\textbf{f}\left(\textbf{r}\right)=1$. It is clear that moving to a larger volume is always favourable in this context. In the worst case $\textbf{f}\left(\textbf{r}\right)$ simply does not change. Therefore, the largest singular value of $\operatorname{Im}\left[\mathbb{G}^{\mathrm{vac}}\right]$ is domain monotonic. Now, suppose that all singular values up to $i=N$ have been shown to be domain monotonic, so that $\left(\forall i\leq N\right)u_{i}\leq d_{i}$, and consider $i=N+1$. Recall that the $N+1$ singular vector $|\textbf{d}_{N+1}\rangle$ is defined by the property of maximizing $\langle\textbf{d}_{N+1}\left|\operatorname{Im}\left[\mathbb{G}^{\mathrm{vac}}\right]\right|\textbf{d}_{N+1}\rangle$ subject to the constraints $\left(\forall i\leq N\right)\langle\textbf{d}_{i}|\textbf{d}_{N+1}\rangle=0$ and $\langle\textbf{d}_{N+1}|\textbf{d}_{N+1}\rangle=1$. Take $\pi_{{}_{U}}|\textbf{d}_{i}\rangle$ to be the projection of $|\textbf{d}_{i}\rangle$ onto the subdomain $U$. If $\sum_{i=1}^{N}\pi_{{}_{U}}|\textbf{d}_{i}\rangle$ spans $\sum_{i=1}^{N}|\textbf{u}_{i}\rangle$, then we essentially return to the case of the first singular vector. $|\textbf{u}_{N+1}\rangle$ is orthogonal to all $\pi_{{}_{U}}|\textbf{d}_{i}\rangle\ni i\leq N$, and zero outside $U$. Hence, it is orthogonal to all $|\textbf{d}_{i}\rangle\ni i\leq N$. As $\langle\textbf{u}_{N+1}\left|\operatorname{Im}\left[\mathbb{G}^{\mathrm{vac}}_{{}_{D}}\right]\right|\textbf{u}_{N+1}\rangle=\langle\textbf{u}_{N+1}\left|\operatorname{Im}\left[\mathbb{G}^{\mathrm{vac}}_{{}_{U}}\right]\right|\textbf{u}_{N+1}\rangle$, with the added freedom of the additional volume we must have $d_{N+1}\geq u_{N+1}$. If $\sum_{i=1}^{N}\pi_{{}_{U}}|\textbf{d}_{i}\rangle$ does not span $\sum_{i=1}^{N}|\textbf{u}_{i}\rangle$, then some elements of $\left\{|\textbf{u}_{i}\rangle\ni i\leq N\right\}$ are orthogonal to all $\pi_{{}_{U}}|\textbf{d}_{i}\rangle$. Select one such vector and denote it as $|\textbf{u}_{\alpha}\rangle$. Since $u_{\alpha}\geq u_{N+1}$, by the preceding argument $d_{N+1}\geq u_{N+1}$.

**Singular values for balls—**To derive the singular values of $\operatorname{Im}\left[\mathbb{G}^{\mathrm{vac}}\right]$ for a ball we have followed the formulation in Tsang et al. Tsang et al. (2004). In this work, it is shown that

[TABLE]

where

[TABLE]

and $j_{\ell}\left(R\right)$ is the $\ell$th spherical with additional factor of $2\pi$ included in its argument. In these definition, $\textbf{V}_{\ell m}^{\left(\alpha\right)}$ are the vector spherical harmonics, obeying the orthogonality conditions

[TABLE]

(As the spherical harmonics are unaffected by projection into a connected volume, the program given below is also valid for shells and other homeomorphic domains.) Comparing (23) with a standard singular value decomposition, we equate the value of the inner product of the vector pairs forming the above outer products with the singular values of $\operatorname{Im}\left[\mathbb{G}^{\mathrm{vac}}\right]$. Performing the angular integrals over the vector spherical harmonics, results in two types of singular values quoted in the main text

[TABLE]

**Singular values for films—**Calculation of the singular values for films of thickness $h$ follows a similar procedure to that of a ball. Using results from Tsang et al. Tsang et al. (2004) and Krüger et al. Krüger et al. (2012), the imaginary part of the Green function can generally be decomposed in terms of “regularized” spectral basis as $\operatorname{Im}\left[\mathbb{G}^{\mathrm{vac}}\right]=\sum_{j\in\mathrm{prop}}|\mathbf{E}^{\mathrm{reg}}_{j}\rangle\langle\mathbf{E}^{\mathrm{reg}}_{j}|$. The “propagating” basis functions of this set are orthogonal to one other (though not necessarily self-normalizable), and form a complete set both at the origin and at infinity. In particular, considering the in-plane wavevector $\mathbf{k}=k_{x}\mathbf{e}_{x}+k_{y}\mathbf{e}_{y}\in\mathbb{R}^{2}$, the basis functions can be written as plane waves $\mathbf{E}^{\mathrm{reg}}_{s,p}(\mathbf{k},\mathbf{r})$, with $s\in\{-1,1\}$ denoting the parity and $p\in\{M,N\}$ the polarization. The vacuum Green function in this basis then takes on the form

[TABLE]

with singular vector functions given by

[TABLE]

in position space, where $k_{z}=\sqrt{k_{0}^{2}-|\mathbf{k}|^{2}}$ is real and nonnegative by virtue of the restriction to propagating waves. The inner products may then all be written as

[TABLE]

due to the orthogonality of these basis functions. This immediately yields the desired singular values corresponding to each $\mathbf{k}$:

[TABLE]

**Inverse design—**To explore the largest possible design space, the optimizations shown in the text result from the “topology” (density) approach Molesky et al. (2018), in which each pixel (permittivity value) within the chosen bounding domain is considered as an independent design parameter. The key to the tractability of such large-scale optimizations is the use of gradient-based optimization algorithms (the method of moving asymptotes is employed in our algorithm Svanberg (2002)). To make use of these approaches, each pixel is initially treated as continuous. That is, at each position in the domain, we begin by assign a parameter $\lambda_{i}\in\left[0,1\right]$ and state that the susceptibility of that pixel is

[TABLE]

where $\chi_{\text{mat}}$ is the material susceptibility. An initial optimization is then carried out using NLOPT Johnson (2014), producing a “gray” structure in which many $\lambda_{i}$ take on intermediate ($\neq\left\{0,1\right\}$) values. In subsequent optimizations (for the same structure) these values are then binarized (i.e. forcing $\lambda_{i}=0$ or $\lambda_{i}=1$) by enclosing $\lambda_{i}$ in function, $\chi_{i}=f_{n}\left(\lambda_{i}\right)\chi_{\text{mat}}$, which is slowly changed from the linear relation given above to a smooth approximation of the step function Jensen and Sigmund (2011).

The primary challenge of this type of inverse design approach for the problem we have considered lies in the cost of computing $\Phi$ for inhomogeneous medium, which typically needs to be evaluated thousands of times in a single optimization iteration. To surmount this difficulty, we have exploited our previously discussed fluctuating volume current formulation Polimeridis et al. (2015); Jin et al. (2016). This framework allows for two major simplifications. First, it removes the necessity of simulating space outside the bounding domain, as would be necessary to accurately calculate emitted thermal radiation or angle-integrated absorption using a finite difference approach. Second, it allows the central matrix-vector multiplication that is used by the iterative inversion solver Polimeridis et al. (2015), $\mathbb{G}^{\mathrm{vac}}U_{i}$ in (33), to be computed via fast-Fourier transforms. To further reduce computational cost, we also use the fact that $\operatorname{Im}\left[\mathbb{G}^{\mathrm{vac}}\right]$ has low pseudo-rank. (In a similar spirit to the derivation of the optimal $\mathbb{T}$-operator given the main text, we use foreknowledge of the number of modes that will possibly contribute to select an appropriate algorithm and pre-allocate computational resources.) This allows us to formulate the scattering inversion problem, (33), directly in terms of a singular value decomposition (SVD) of $\operatorname{Im}\left[\mathbb{G}^{\mathrm{vac}}\right]=Q\Sigma Q^{\dagger}$, which can be approximated to any accuracy with efficient randomized methods Hochman et al. (2014). (For a full discussion of this procedure see Refs. Halko et al. (2011); Martinsson et al. (2011); Polimeridis et al. (2015).) Specifically, starting from an equivalent trace form of heat transfer to the one given in the main text, the singular value decomposition of $\operatorname{Im}\left[\mathbb{G}^{\mathrm{vac}}\right]$ allows us to recast $\Phi$ as

[TABLE]

where the $F$ subscript denotes the Frobenius norm. Evaluation of (32) requires only one inverse solve for each column $Q_{i}$ of the low rank decomposition:

[TABLE]

(32) is also well suited to calculation of the gradient of $\Phi$. Starting from this result, using (31), direct application of matrix calculus shows that

[TABLE]

and the total gradient is then

[TABLE]

As before, only a number solves equal to the rank of the singular value decomposition of the imaginary part of the Green function are required to evaluate (35) for all design parameters.

**Singular values of $\mathbb{T}$ from inverse design—**As a point of comparison with the assumptions made in deriving $\Phi_{\text{opt}}$, we have explored the $\mathbb{T}$-operator for the dielectric structure of the rightmost green square in Fig. 2 ($\chi=20+4i$, $\Phi=0.6~{}\Phi_{\text{opt}}$) in the basis of $\operatorname{Im}\left[\mathbb{G}^{\text{vac}}\right]$. This physical realization proves to nearly satisfy both the assumption of asymmetry ($\text{atan}\left(\text{Im}\left(\langle\textbf{q}_{i}\left|\mathbb{T}\right|\textbf{q}_{i}\rangle\right)/\text{Re}\left\{\langle\textbf{q}_{i}\left|\mathbb{T}\right|\textbf{q}_{i}\rangle\right\}\right)=\pi/2$) and simultaneous diagonalizability. Specifically, we find that

[TABLE]

The ratio of the true values of $\langle\textbf{q}_{i}\left|\mathbb{T}\right|\textbf{q}_{i}\rangle$ in comparison to the ideal values determined by $\Phi_{\text{opt}}$ are given in Tab.1 below. These values are naturally grouped by a singular value type index ($\left\{1,2\right\}$) and the spherical harmonic $\ell$ index. For clarity, an average is taken over the (ideally constant) harmonic $m$ subindex. The final column gives the weight that each entry makes to $\Phi_{\text{opt}}$.