Fundamental limits to radiative heat transfer: theory
Sean Molesky, Prashanth S. Venkataram, Weiliang Jin, Alejandro W. Rodriguez
FFundamental limits to radiative heat transfer: theory
Sean Molesky, ∗ Prashanth S. Venkataram, ∗ Weiliang Jin, and Alejandro W. Rodriguez Department of Electrical Engineering, Princeton University, Princeton, New Jersey 08544, USA (Dated: July 9, 2019)Near-field radiative heat transfer between bodies at the nanoscale can surpass blackbody limits on thermalradiation by orders of magnitude due to contributions from evanescent electromagnetic fields, which carry noenergy to the far-field. Thus far, principles guiding explorations of larger heat transfer beyond planar structureshave assumed utility in surface nanostructuring, which can enhance the density of states, and further assumedthat such design paradigms can approach Landauer limits, in analogy to conduction. We derive fundamentalshape-independent limits to radiative heat transfer, applicable in near- through far-field regimes, that incorporatematerial and geometric constraints such as intrinsic dissipation and finite object sizes, and show that thesepreclude reaching the Landauer limits in all but a few restrictive scenarios. Additionally, we show that theinterplay of material response and electromagnetic scattering among proximate bodies means that bodies whichmaximize radiative heat transfer actually maximize scattering rather than absorption. Finally, we compare ournew bounds to existing Landauer limits, as well as limits involving bodies maximizing far-field absorption, andshow that these lead to overly optimistic predictions. Our results have ramifications for the ultimate performanceof thermophotovoltaics and nanoscale cooling, as well as related incandescent and luminescent devices.
The concept of a blackbody, derived from electromag-netic reciprocity (or detailed balance), has provided a bench-mark of the largest emission rates that can be achieved bya heated object: through nanoscale texturing, gray objectscan be designed in myriad ways to mimic the response of ablackbody at selective wavelengths [1–3], with implicationsfor a variety of technologies, including high-efficiency solarcells, selective emitters, and thermal sensors [4]. Over thepast few decades, much effort has gone toward understand-ing analogous limits to enhancements of near-field radiativeheat transfer (RHT) [5–8], supported by a rich and growingnumber of experimental [9–12] and theoretical [13–17] in-vestigations, and motivated by potential applications to ther-mophotovoltaics [18, 19], nanoscale cooling [20], and ther-mal microscopy [21, 22]. A key principle underlying furthernear-field RHT enhancements is the use of materials support-ing bound (plasmon and phonon) polaritons in the infrared,where the Planck distribution peaks at typical temperaturesprobed in experiments. This leads to strong subwavelengthresponses tied to corresponding enhancements in the den-sity of states [23–26]; consequently, the amplified near-fieldRHT spectrum exhibits a narrow lineshape, justifying focuson selective wavelengths. However, while the properties ofsuch polaritons, particularly their resonance frequencies, as-sociated densities of states, and scattering characteristics canbe modified through nanoscale texturing, only recently havecomputational methods [14–16, 27] arisen to model RHT be-tween bodies of arbitrary shapes beyond those with high sym-metry [17, 28, 29]. Furthermore, the challenge of gainingsimultaneous control over the scattering properties of largenumbers of contributing surface waves has generally pre-cluded general upper bounds on RHT.RHT between two bodies A and B in vacuum is given as P = ˆ ∞ [Π( ω, T B ) − Π( ω, T A )]Φ( ω ) d ω, (1)in terms of their local temperatures T A and T B , entering the Planck function Π( ω, T ) = (cid:126) ω/ [exp( (cid:126) ω/ ( k B T )) − (andit has been assumed, without loss of generality, that T B >T A so P > ), and the spectral function Φ( ω ) , which canbe enhanced by changing material and geometric propertiesthrough the creation of resonances and changes in the elec-tromagnetic density of states. In particular, nanostructuringmetallic surfaces or polar dielectrics makes it possible to shiftresonant frequencies from the visible or far infrared into theinfrared, such that the peak of the spectrum Φ may coincidewith the peak of the Planck distribution near room tempera-ture. It remains an open question, however, to what extent thepeak value of Φ itself may be enhanced through appropriategeometric and material choices, as well as what such optimalstructures should be.Previous attempts at deriving bounds on RHT have pri-marily focused on extended media [5–7, 30], showing thatat least for translationally invariant structures, Φ can be ex-pressed as the trace of a “transmission” matrix whose singu-lar values (corresponding to evanescent Fourier modes) eachcontribute a finite flux, bounded above by a Landauer limitin analogy with conduction [31, 32]. Aside from being re-stricted to planar geometries, these bounds turn out to be ei-ther pessimistic [7], ignoring the large densities of states thatcan arise in nanostructured and low-loss materials, or too op-timistic [5, 6], ignoring any constraints imposed by Maxwell’sequations and assuming instead that all such Fourier modes,up to an unrealistic cutoff on the order of the atomic scale,can saturate the flux [5]. From a design perspective, Lan-dauer limits present a hurdle as they rely on ad-hoc esti-mates of the number and relative contribution of radiativemodes/channels, which depend on specific material and ge-ometric features. More recent works have derived comple-mentary material limits on electromagnetic absorption in sub-wavelength regimes [33], showing that absorbed power ina medium of susceptibility χ increases in proportion to an“inverse resistivity” figure of merit, | χ | / Im χ , in principlediverging with increasing indices of refraction and decreas- a r X i v : . [ phy s i c s . c l a ss - ph ] J u l ing dissipation. Saturation of these bounds for a subwave-length absorber in the quasistatic regime can generally beachieved through the strong polarization currents arising inresonant media supporting surface plasmon or phonon po-laritons. These arguments have been extended to near-fieldRHT [8] by exploiting energy conservation and reciprocity,finding the upper bound of Φ at a polariton resonance to scalequadratically with | χ | / Im χ , corresponding to enhanced ab-sorption and emission in both objects. While near-field RHTbetween dipolar objects can attain these bounds in a dilutelimit, such a universal scaling has yet to be observed in large-area structures. This naïvely suggests room for improvementin Φ through nanostructuring via enhancements in the den-sity of states or equivalently, via saturation of modal contribu-tions, yet trial-and-error explorations and optimization proce-dures [34, 35] have failed to produce nanostructured geome-tries that bridge this gap, leading to the alternative possibilitythat 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 anal-ysis of the singular value decompositions of relevant responsequantities. In contrast to prior limits, our bounds incorporateconstraints imposed by material losses and multiple scatter-ing, and are therefore tighter; moreover, they are formulatedto be independent of object shapes while simultaneously ac-counting for finite size effects. In particular, every channel ofenergy transmission is shown to be generally prohibited fromsaturating its Landauer limit, in contrast to predictions basedon modal decompositions [5–7, 30] that neglect material prop-erties and are most applicable in the ray optics regime. Fur-thermore, the growth of RHT with decreasing material dis-sipation is shown to be strongly limited by radiative losses,in contrast to predictions based on energy-conservation limitsto material response [8] that neglect finite-size scattering ef-fects and are thus most applicable in the quasistatic regime.Finally, we discuss the difficulties in optimizing each indi-vidual energy transmission channel and show how more re-stricted bounds may be obtained by approximating each bodyto maximally absorb all incident fields in isolation. In upcom-ing papers closely related to this one, we apply these boundsto various scenarios of interest, providing predictions of themaximum RHT achievable in compact and extended geome-tries [36], and deriving related bounds on far-field thermalemission from single bodies in isolation [37]. I. HEAT TRANSFER DEFINITIONS
For two bodies A and B in vacuum [Fig. 1], the spec-tral function Φ appearing in (1) represents the average powerabsorbed in B due to fluctuating current sources in A , de-picted in and is reciprocal (invariant under interchange ofA and B). Using formal operator notation, this average ab-sorbed power can be written in terms of the susceptibilities V p , the vacuum Maxwell Green’s function G vac pq , and scat-tering T-operators T p , for p, q ∈ { A , B } . Each susceptibil- χχχχχχχ AB χχχχχχχχχχχχχχχχχχχχχ AAAAAAAAAAAAAAAAA
Figure 1. Two bodies labeled A and B exchange heat in vacuum.Each body could be compact or of infinite extent in at least one spa-tial dimension, and for given susceptibilities χ p , the optimal struc-tures may be quite complicated, but the upper bounds, which de-pend on ζ p = | χ p | / Im( χ p ) , can be evaluated in simpler boundingdomains that enclose each object while respecting any other designconstraints present. ity is written as V p = χ p I p , where each χ p is assumed tobe homogeneous, local, and isotropic. The vacuum MaxwellGreen’s function G vac solves [( c/ω ) ∇ × ( ∇× ) − I ] G vac = I in all space, and its blocks are denoted as G vac pq for sources inbody q propagating fields to body p . Finally, the T-operators T p = ( V − p − G vac pp ) − represent the total induced polarizationmoment in body p due to a localized field of unit magnitudeincident upon it. All of these quantities are reciprocal, so theyare equal to their unconjugated transposes in position space: V p = V (cid:62) p , T p = T (cid:62) p , and G vac pq = ( G vac qp ) (cid:62) . This means thatHermitian conjugation is equivalent to complex conjugation: V † p = V (cid:63)p , T † p = T (cid:63)p , and ( G vac pq ) † = G vac (cid:63)qp . Additionally, allof these quantities depend on frequency ω , though this will besuppressed for the sake of notational brevity.Given these definitions and relations (see appendices formore details), the RHT spectrum can be written as [17] Φ = 2 π Tr[ T (cid:63) B ( I B − G vac (cid:63) BA T (cid:63) A G vac (cid:63) AB T (cid:63) B ) − Im( V − (cid:63) B ) × T B ( I B − G vacBA T A G vacAB T B ) − × G vacBA T A Im( V − (cid:63) A ) T (cid:63) A G vac (cid:63) AB ] , (2)where Im( A ) = ( A − A (cid:63) ) / (2i) and Asym( A ) = ( A − A † ) / (2i) for any operator A ; if A is reciprocal, then Asym( A ) = Im( A ) . This expression is manifestly recipro-cal in A and B, and treats the T-operators of A and B on anequal footing, linked only by the Green’s function G vacBA prop-agating fields in vacuum from one body to the other. However,it is possible to write this spectrum more suggestively in termsof operator combinations that hide this reciprocity in order tomore strongly link this expression to absorbed and emittedpowers. In particular, Φ may be rewritten as Φ = 2 π Tr[ Y (cid:63) B Im( V − (cid:63) B ) Y B G vacBA T A Im( V − (cid:63) A ) T (cid:63) A G vac (cid:63) AB ] , (3)in terms of the reciprocal operator Y B = T B S AB , which isin turn written in terms of the scattering operator S AB =( I B − G vacBA T A G vacAB T B ) − : essentially, Y B is a new “dressedT-operator” describing absorption and scattering in B in thepresence of A, just as the bare T-operators T p describe absorp-tion and scattering from each body p ∈ { A , B } in isolation.Having assumed the susceptibility in each body p ∈ { A , B } to be homogeneous, uniform, and isotropic, makes it possibleto identify Im( V − (cid:63)p ) = Im( χ p ) | χ p | I p . For convenience, we de-note ζ p = | χ p | Im( χ p ) as the material response factor. Using this,we write the RHT spectrum as Φ = 2 πζ A ζ B (cid:107) Y B G vacBA T A (cid:107) (4)where (cid:107) A (cid:107) = Tr( A † A ) denotes the Frobenius norm for anyoperator A .As we show in the appendices, the RHT spectrum may al-ternatively be written as Φ = 2 π (cid:107) Q (cid:107) (5)in terms of the transmission operator Q =Im( V B ) / G BA Im( V A ) / , which in turn depends onthe total Green’s function G BA = V − T B S AB G vacBA T A V − connecting dipole sources in body A to total fields in bodyB and accounting for multiple scattering to all orders withinand between both bodies. This obeys a Landauer limit, as thesingular values of Q † Q do not exceed / , so including theprefactor /π , the contribution of each mode/channel in thetrace expression to Φ does not exceed π . We will shortlyexplain the conditions under which the Landauer bounds foreach of these contributions may be saturated. II. SINGULAR VALUE BOUNDS
Each of the operators Y B , G vacBA , and T A may be decom-posed into their respective singular values and vectors. Ingeneral, the singular vectors of each of these operators willnot bear any relation to one another. However, we prove in theappendices that an upper bound to Φ can be found by ensur-ing that the right singular vectors of Y B and the left singularvectors of T A are respectively the same as the left and rightsingular vectors of G vacBA . Intuitively, these conditions can beinterpreted as representing perfect coupling and propagationof fields due to sources in A which are perfectly absorbedin body B (while accounting for multiple scattering from the presence of body A). Formally, we write the statements Y B = (cid:88) i y i | b (cid:63)i (cid:105)(cid:104) b i | G vacBA = (cid:88) i g i | b i (cid:105)(cid:104) a i | T A = (cid:88) i τ i | a i (cid:105)(cid:104) a (cid:63)i | (6)where as a reminder, the singular values τ i , g i , and y i arenonnegative real numbers, and the vectors | a i (cid:105) and | b i (cid:105) are orthonormal among themselves so (cid:104) a i , a j (cid:105) = δ ij and (cid:104) b i , b j (cid:105) = δ ij . Our explicit conjugation of | b (cid:63)i (cid:105) in the def-inition of Y B and of (cid:104) a (cid:63)i | in the definition of T A ensure thatthey still respect reciprocity. This allows us to write a boundon the RHT spectrum as Φ = 2 πζ A ζ B (cid:88) i ( y i g i τ i ) (7)which is maximized when y i and τ i attain their maximum val-ues as functions of g i ; the singular values g i depend only onthe geometry and topology of the system as they enter G vacBA ,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 Aand B together implies that the operator Im( Y B ) − Y B (cid:16) ζ B I B + ζ A G vacBA T A T (cid:63) A G vac (cid:63) AB (cid:17) Y (cid:63) B is not only Hermitianbut also positive-semidefinite. This implies that if Y B hada nontrivial Hermitian part, any bounds on the singular val-ues y i arising from this total nonnegative far-field scatter-ing constraint must be more restrictive than if Y B werepurely anti-Hermitian, meaning Y B = i Im( Y B ) . Thus,we take Y B to be anti-Hermitian, which also means that (cid:104) b j , b (cid:63)j (cid:105) = i for every channel j . Additionally, wenote that G vacBA T A T (cid:63) A G vac (cid:63) AB = (cid:80) i g i τ i | b ( i ) (cid:105)(cid:104) b ( i ) | followsfrom our singular value decompositions, so it remains tobe proved that the operator (cid:80) j y j (cid:0) | b (cid:63)j (cid:105)(cid:104) b j | − | b j (cid:105)(cid:104) b (cid:63)j | (cid:1) − (cid:80) j y j (cid:16) ζ B + ζ A ( g j τ j ) (cid:17) | b (cid:63)j (cid:105)(cid:104) b (cid:63)j | is positive-semidefinite.To do this, we take the inner product of this operator with | b (cid:63)i (cid:105) on the right and (cid:104) b (cid:63)i | on the left for a single fixed index i (not to be confused with the imaginary unit), which yieldsthe condition y i − y i (cid:16) ζ B + ζ A ( g i τ i ) (cid:17) ≥ . This inequalityis saturated to give the largest possible y i = ζ B ζ B ζ A ( g i τ i ) interms of g i τ i , which in turn gives the largest possible flux, Φ = 2 ζ B πζ A (cid:88) i ( g i τ i ) (cid:104) ζ B ζ A ( g i τ i ) (cid:105) in terms of g i τ i .The above expression is maximized for each channel i , andfor the sum in turn, if ζ B ζ A ( g i τ i ) = 1 . This is solved to yield τ i = (cid:113) ζ B ζ A g i , and plugging this into the expression for y i yields y i = ζ B , while plugging this into the expression for Φ yields a contribution πζ A ζ B ( y i g i τ i ) = π × = π . Weinterpret this to mean that to obtain optimal heat transfer, theT-operator of body A in isolation must be engineered in a waythat depends on the presence of body B, due to both the pres-ence of the material enhancement factor ζ B and the depen-dence on the singular values g i of G vacBA propagating electro-magnetic fields in vacuum from A to B, though no higher-order scattering effects come into play. In turn, the expression y i = ζ B means that the effective T-operator of body B dressedby scattering from body A must actually exhibit maximal scat-tering , and not absorption, in the presence of body A [33],though it is more difficult to extract information about the im-plications for the T-operator of body B in isolation; we clarifythat maximal scattering includes both far-field scattering frombody B in the presence of body A, as well as absorption frombody A in the presence of body B. If these two conditions canbe met simultaneously for the given channel i , which is effec-tively a rate-matching condition relating the absorption andscattering rates of each body in the presence of the other, thenthe contribution π is exactly the Landauer upper bound onenergy transmission.However, this analysis has so far ignored constraints on τ i .In particular, nonnegativity of far-field scattering from bodyA in isolation requires that Im( T A ) − ζ A T A T (cid:63) A be a positive-semidefinite Hermitian operator [33]. Following similar stepsas above, we may see that for the purposes of bounding τ i , theloosest constraints on τ i arise if T A is purely anti-Hermitian,so taking the inner product of the above operator with respectto | a i (cid:105) on the right and (cid:104) a i | on the left yields the condition τ i ≤ ζ A . The expression τ i = (cid:113) ζ B ζ A g i is consistent withthis inequality only if the inequality g i ≥ √ ζ A ζ B also holdsfor the given index i . If not, and instead g i < √ ζ A ζ B , then τ i = ζ A must be used to maximize the contribution at index i to Φ , saturating the bound on the response of A in isola-tion. This yields y i = ζ B ζ A ζ B g i ≥ ζ B and a contribution of π ζ A ζ B g i (1+ ζ A ζ B g i ) ≤ π to Φ . We interpret this to mean that ifthe singular value g i of G vacBA falls below a threshold involv-ing the two material enhancement factors, then the optimalT-operator of body A in isolation corresponds to maximal ab-sorption, the optimal effective T-operator of body B dressedby body A evinces the effects of multiple scattering with A,and the contribution to Φ similarly shows the effects of multi-ple scattering between the two bodies and is unable to saturatethe Landauer bound for that channel.To summarize, the bound on RHT may be written as Φ opt = (cid:88) i (cid:104) π Θ( ζ A ζ B g i − π ζ A ζ B g i (1 + ζ A ζ B g i ) Θ(1 − ζ A ζ B g i ) (cid:105) (8)where Θ is the Heaviside step function. As the material re-sponse (encoded in ζ A ζ B ) increases, progressively more chan- nels may saturate the Landauer limit per channel, so that theLandauer limit (summed over all channels) is reached asymp-totically as ζ A ζ B → ∞ . However, the rate at which this diver-gence occurs depends on the general topology of the problem,as that determines how the singular values g i of G vacBA dependon the index i .We emphasize that while the singular values g i of G vacBA aretechnically restricted to the domains of the objects to give thetightest bound on heat transfer, such a restriction is less thanideal given the explicit dependence on the shapes of the ob-jects. However, as we prove in the appendices, the singularvalues g i of G vacBA are domain monotonic, meaning that theyincrease monotonically as the volumes of regions A and Bincrease; consequently, Φ opt is domain monotonic, as it ismonotonically non-decreasing with respect to g i for each i .Separately from this, the regions containing only the materialdegrees of freedom of each body can be replaced by largerregions that fully enclose each body, as the T-operators ofeach body will commute with projections into the smaller sub-spaces corresponding to the actual material degrees of free-dom. Thus, these bounds can be slightly loosened to be in-dependent of body shapes, and can then be evaluated subjectto constraints on topology and domain volumes as determinedby the desired application [Fig. 1], e.g. ellipsoids with pre-scribed aspect ratios or films of prescribed thicknesses rep-resentative of compact or extended geometries, respectively.Essentially, the effective rank of G vacBA , which determines thenumber of modes that could participate in RHT, is largely de-termined by the size and topology of the choice of boundingsurface, which represents a general and fundamental geomet-ric constraint on the bounds of RHT analogous to the generalmaterial constraints imposed by ζ p for each body p ∈ { A , B } . III. COMPARISON TO ALTERNATIVE BOUNDS
The bound for the RHT spectrum Φ in (8) may be comparedto a number of other bounds. Strictly speaking, this bound isnot necessarily the tightest general bound that could be for-mulated. In particular, using the relation T − = V − − G vacAA allows for writing (3) in terms of T A and Im( G vacAA ) withoutreference to V A . Such a procedure, in analogy with bounds onthermal emission which we detail in an upcoming paper [37],would more explicitly capture far-field radiative losses frombodies of finite size, which becomes more relevant at largeseparations where such losses may compete with RHT itself.However, as we show in the appendices, we find the result-ing bound to be intractable, requiring self-consistent solutionof systems of nonlinear equations to find the optimal singularvalues of T A . Therefore, we do not further consider such abound, and henceforth refer only to (8).With respect to prior work, the most obvious point of com-parison is the Landauer bound [5], namely Φ L = (cid:88) i π , (9)which simply depends on the number of modes participatingin RHT, without any reference to separation, geometric con-straints, multiple scattering, or even material constraints, letalone their interplay; consequently, in contrast to our bounds,there is no metric to evaluate how many participating modescan actually saturate the limit π . Even modal analyses thattechnically do not necessarily assume saturation of the Lan-dauer limits for every mode [5–7, 30] tend to neglect materialeffects, so the purely geometric arguments are valid only in theray-optical regime where blackbody limits are reproduced.An alternative bound can be derived by modifying the formof (8) such that the contributions from saturating the singularvalue bounds for τ i are used for all indices i , not only thosewhere ζ A ζ B g i < . We may write this as Φ sc = (cid:88) i π ζ A ζ B g i (1 + ζ A ζ B g i ) (10)and term this the “scalar approximation” with the subscript“sc” because (as shown in the appendices) it can equivalentlybe derived by plugging into (2) the T-operators T p = i τ p I p for each body p ∈ { A , B } , corresponding to uniform singularvalues τ p for each T p ; in the appendices, we further derivethat this bound on RHT is maximal when each body exactlysatisfies the condition for perfect isolated absorption [33] forevery incident field, meaning τ p = ζ p . As each contribution to Φ sc depends on the singular values g i in a nonmonotonic way,domain monotonicity is not obvious, though arguments basedon the embedding of the T-operators in larger spaces allowfor consideration of general bounding surfaces as geometricconstraints, making this a useful bound restricted to the classof T-operators with uniform singular values. However, as weshow in an upcoming paper, Φ sc overall can be shown to bedomain monotonic in the near-field.Another such bound neglects multiple scattering effects be-tween the two bodies, and can be derived in several ways.As we show in the appendices, this bound follows from thescalar approximation Φ sc by neglecting the scattering opera-tor S AB → I B , or alternatively by ignoring the denominators (1 + ζ A ζ B g i ) in Φ sc encoding the effects of multiple scatter-ing within each channel. Regardless of the particular deriva-tion, the end result may be written as Φ Born = (cid:88) i π ζ A ζ B g i (11)which we term the “Born bound” as its neglect of multi-ple scattering effects between the two bodies is characteris-tic of a Born approximation to the full solution of Maxwell’sequations. The domain monotonicity of this bound triv-ially follows from that of g i for each channel i , and the va-lidity of the embedding argument with respect to T p stillholds, meaning that Φ Born is also a useful bound for theRHT spectrum when considering bounding surfaces of ar-bitrary size and topology. We note that Miller et al [8]equivalently derive this bound by limiting the singular val-ues of Y B such that only far-field scattering of body B in the presence of body A needs to be nonnegative, meaning Im( Y B ) − Y (cid:63) B Im( V − (cid:63) B ) Y B should be positive-semidefinite;this leads to the bound y i ≤ ζ B , so maximization of (7) sub-ject to that constraint as well as τ i ≤ ζ A simply requires sat-uration of both of these constraints. However, we have shownthat positive-semidefiniteness of Im( Y B ) − Y (cid:63) B Im( V − (cid:63) B ) Y B ,corresponding to nonnegative scattering from body B inthe presence of body A, is a looser constraint on y i thannonnegativity of scattering from the system as a whole,corresponding to positive-semidefiniteness of Im( Y B ) − Y B (cid:0) Im( V − (cid:63) B ) + G vacBA T A Im( V − (cid:63) A ) T (cid:63) A G vac (cid:63) AB (cid:1) Y (cid:63) B , and thatonly the latter produces a bound on y i that explicitly accountsfor multiple scattering and absorption losses through body A.The Born bound not only neglects multiple scattering but alsoignores far-field radiative losses, which means that while it istechnically a looser bound, its usefulness is limited to deeplysubwavelength systems exchanging heat in the near-field.We are now in a position to relate Φ opt , Φ sc , Φ Born , and Φ L to each other via inequalities. First, the discussion of (8)immediately makes clear that because some of the contribu-tions to Φ opt equal corresponding contributions to Φ L whileall of the others equal corresponding contributions to Φ sc , thenit must be the case that Φ sc ≤ Φ opt ≤ Φ L . Additionally,because the contributions to Φ opt from each channel i neverexceed π , and in particular because the saturation of the Lan-dauer bound for each channel occurs for ζ A ζ B g i ≥ , then itmust be the case that Φ sc ≤ Φ opt ≤ Φ Born , because in theexpression for Φ Born , the contributions ζ A ζ B g i will alwaystrivially exceed / when ζ A ζ B g i ≥ , and will always ex-ceed ζ A ζ B g i (1+ ζ A ζ B g i ) when ζ A ζ B g i < . In general, it is not pos-sible to write an inequality relation between Φ Born and Φ L in all situations, because Φ Born may have some contributions π ζ A ζ B g i which fall above or below π , and the geometry de-termining g i would have to be known in order to know howmany fall above or below. Thus, we may write the overallinequalities as Φ sc ≤ Φ opt ≤ Φ Born , Φ L . (12) IV. CONCLUDING REMARKS
We have determined bounds for the RHT spectrum Φ basedpurely on algebraic arguments. In particular, we show thatthere is a tension between optimizing transmission channelsand material/geometric constraints placed on each object inisolation as well as in in the presence of the other; as a re-sult, some but not all channels can saturate previously derivedLandauer bounds, while others are restricted by the aforemen-tioned constraints. By virtue of domain monotonicity, thesebounds can be applied in a shape-independent manner, sowhile they can be evaluated analytically in highly symmet-ric bounding surfaces, they can just as easily be evaluated nu-merically in more complicated domains depending on specificdesign constraints [Fig. 1]. Similarly, the dependence on thematerial response factor ζ = | χ | / Im χ does not make ex-plicit reference to a particular frequency or material model.In comparison, the Landauer bounds yield overly optimisticpredictions, while choosing a scalar response for each ob-ject corresponding to maximal absorption of every incidentfield in isolation yields overly pessimistic predictions. Addi-tionally, we find that previous work by Miller et al [8] alsoyields overly optimistic predictions compared to our currentbounds, because those derivations neglect the geometric ef-fects of multiple scattering between the two bodies and con-sequently overestimate the optimal response of one body inthe presence of the other. We point out that while our boundsare always tighter than Landauer and Born limits for any givenbounding domain, they say nothing about which domains mayyield the tightest per-volume limits given material constraints,or whether they may in fact be attained by physically realiz-able structures.In summary, while Born and Landauer limitsare technically upper bounds on RHT, their neglect of mul-tiple scattering and material losses, respectively, render themloose compared to the bounds presented here.We emphasize that in contrast to Born limits [8], our boundsare valid from the near-field all the way through the far-fieldfor bodies of arbitrary size, as we never needed to exploitnonretarded or quasistatic approximations. We also point outthat in the near-field, G vacBA is real-valued in position space,so arguments involving the anti-Hermiticity of G vacBA T A G vacAB imply the anti-Hermiticity of T A , as Im( G vacBA T A G vacAB ) → G vacBA Im( T A ) G vacAB . In this context, Born bounds only yieldsensible results in the near-field, as they tend to diverge be-yond the near-field regime for objects of infinite extent in atleast one spatial dimension faster than expected.In a complementary paper [36], we analyze these boundsin the near-field in specific geometries of interest, particularlydipolar and extended (infinite area) bodies. We take particularnote of the scalar approximation because, though represent-ing a restricted version of our general bounds and resulting incomparably pessimistic predictions, it may prove challengingif not impossible to optimize every individual channel to sat-urate its Landauer limit when the material constraints are notsaturated. While this mere observation does not preclude thepossibility of larger RHT (as allowed by our general bounds),from a design perspective it might prove more feasible to con-sider bodies capable of maximally absorbing every incidentfield in isolation, as could occur near a polaritonic resonance.Moreover, as we show in this follow-up work, the simplic-ity of the scalar approximation in the near field lends itself tovarious mathematical manipulations which lead to the follow-ing suggestive results. First, the scalar bound is a local sta-tionary point with respect to perturbations of the T-operatorsof each body, away from the perfect isolated absorption con-dition. Second, Φ sc can be shown to be domain-monotonic,meaning that for fixed ζ p for each body p ∈ { A , B } , the near-field RHT spectrum will always increase with increasing ob-ject volumes. Together, these findings suggest that the role ofnanostructuring in enhancing near-field RHT will be limited,and this has important implications for the theoretical and ex-perimental design of devices for cooling, heat dissipation, and energy generation. ACKNOWLEDGMENTS
The authors would like to thank Riccardo Messina andPengning Chao for helpful discussions. This work was sup-ported by the National Science Foundation under Grants No.DMR-1454836, DMR 1420541, DGE 1148900, the CornellCenter for Materials Research MRSEC (award no. DMR-1719875), and the Defense Advanced Research ProjectsAgency (DARPA) under agreement HR00111820046. Theviews, opinions and/or findings expressed are those of the au-thors and should not be interpreted as representing the officialviews or policies of the Department of Defense or the U.S.Government.
Appendix A: Notation
We briefly discuss the notation used through the maintext and the appendices. A vector field v ( x ) will be de-noted as | v (cid:105) . The conjugated inner product is (cid:104) u , v (cid:105) = ´ d x u (cid:63) ( x ) · v ( x ) . An operator A ( x , x (cid:48) ) will be denotedas A , with ´ d x (cid:48) A ( x , x (cid:48) ) · v ( x (cid:48) ) denoted as A | v (cid:105) . The Her-mitian conjugate A † is defined such that (cid:104) u , A † v (cid:105) = (cid:104) A u , v (cid:105) .The anti-Hermitian part of a square operator (whose do-main and range are the same size) is defined as the operator Asym( A ) = ( A − A † ) / (2i) . Finally, the trace of an operatoris Tr( A ) = ´ d x Tr( A ( x , x )) . Through this paper, unlessstated explicitly otherwise, all quantities implicitly depend on ω , and such dependence will be notationally suppressed forbrevity. Appendix B: Proof of trace maximization via shared singularvectors
In this section, we prove the following lemma: if oper-ators A n for n ∈ { , , . . . , N } have fixed singular val-ues labeled σ ( n ) i , then the singular vectors that maximize Tr[ A A . . . A N ] are common between operators multipliedtogether. That is, the singular value decomposition of A n should follow A n = (cid:80) i σ ( n ) i | a ( n ) i (cid:105)(cid:104) a ( n +1) i | for n ∈{ , , . . . , N − } , with A N = (cid:80) i σ ( N ) i | a ( N ) i (cid:105)(cid:104) a (1) i | , wherethe vectors | a ( n ) i (cid:105) are orthonormal for each n such that (cid:104) a ( n ) i , a ( n ) j (cid:105) = δ ij . This lemma will hold even if each A n is not square, as long as A n A n +1 forms a valid nontrivialoperator product, as these can be embedded in larger spacespadded with more vanishing singular values. Thus, we restrictour consideration to square operators. Moreover, associativitymeans A n A n +1 A n +2 = ( A n A n +1 ) A n +2 , and the trace of aproduct of operators is invariant under cyclic permutations, sowe ultimately only consider maximizing the trace of a productof two operators, as maximization of the trace of products ofmore than two operators follows inductively from this.To maximize Tr[ AB ] , we start by writing A = N (cid:88) i =1 σ i | u i (cid:105)(cid:104) v i | (B1) B = N (cid:88) j =1 τ j | w j (cid:105)(cid:104) y j | (B2)where N is the size of the space; this may be larger thanthe rank of either A or B , but the point is moot because thesingular values are fixed, whether they vanish or not, andit has already been assumed that A and B are square. Wealso assume that the singular values are ordered such that σ i ≥ σ i +1 for all i ∈ { , , . . . , N − } and τ j ≥ τ j +1 for all j ∈ { , , . . . , N − } . This allows for writing Tr[ AB ] = N (cid:88) i =1 N (cid:88) j =1 σ i p ij τ j q ji (B3)in terms of p ij = (cid:104) v i , w j (cid:105) and q ji = (cid:104) y j , u i (cid:105) . As the singu-lar vectors are orthonormal, then p ij and q ji are the elementsof unitary matrices, satisfying (cid:80) Nj =1 | p ij | = (cid:80) Ni =1 | p ij | = (cid:80) Nj =1 | q ji | = (cid:80) Ni =1 | q ji | = 1 . We consider only real non-negative Tr[ AB ] , so that is maximized when σ i p ij τ j q ji are allnonnegative; this means the singular vectors can be chosenwithout loss of generality such that p ij and q ji are real andnonnegative, implying p ij and q ji are the elements of real-valued orthogonal matrices.We use induction to prove that maximizing the trace re-quires that {| u i (cid:105)} be the duals of {(cid:104) y j |} , and that {(cid:104) v i |} by the duals of {| w j (cid:105)} . The case N = 1 is trivial, asall quantities are scalars. For N = 2 , we use orthogonal-ity to note that p , = p , , q , = q , , p , = p , = (cid:113) − p , , and q , = q , = (cid:113) − q , . As a result, wemay write Tr[ AB ] = (cid:113) (1 − p , )(1 − q , )( σ τ + σ τ ) + p , q , ( σ τ + σ τ ) . As the first term in parentheses is largerthan the second term in parentheses by the nonnegative value ( σ − σ )( τ − τ ) given the ordering of singular values, hav-ing the left singular vectors of one operator not be duals of theright singular vectors of the other and vice versa could onlyincrease the trace if (cid:113) (1 − p , )(1 − q , ) + p , q , > ,but this leads to the impossible condition > ( p , − q , ) ,so we can only have p , = q , for the trace to be maximized,implying the duality result must hold.The inductive step assumes an arbitrary N − and movesfrom there to proving the statement for N . Without loss ofgenerality, we consider first the contribution of the largestsingular value τ of B , namely τ | w (cid:105)(cid:104) y | , interacting with A = (cid:80) i σ i | u i (cid:105)(cid:104) v i | in the trace. This yields the contribution (cid:80) i ( p i, − q ,i ) = 2 (1 − (cid:80) i p i, q ,i ) ≥ using the fact that (cid:80) i p i, = (cid:80) i q ,i = 1 , so this in turn gives the condition (cid:80) i p i, q ,i ≤ . The trace can be seen to be maximal when the above condition is saturated, so (cid:80) i p i, q ,i = 1 , whichimplies p i, = q ,i for every i . As this also holds when theroles of A and B are interchanged, and as this can be progres-sively carried out for each successively smaller singular valuegiven orthonormality of the singular vectors, then the dualitycondition must hold, completing the proof. Appendix C: Derivation of radiative heat transfer formulas
In this section, we derive the formula for the radiativeheat transfer spectrum between two bodies, valid even with-out assumptions about retardation, homogeneity, locality, orisotropy. The formula depends on individual T-operators andthe vacuum Green’s function, and follows a previous deriva-tion [38] which considered energy transfer by fluctuating vol-ume currents. Through further derivation, we also equate thisformula to another formula involving the susceptibilities andthe full Maxwell Green’s function, and use that to recast theheat transfer spectrum in a Landauer form, whence we provethat the singular values of the Landauer transmission operatorfor RHT do not exceed / . Finally, we prove that the for-mula for thermal emission of a single body in isolation can bederived from the formula for RHT between two bodies in vac-uum, by taking the second body to fully enclose the first andto be perfectly absorbing, thus taking on the role of a perfectlyabsorbing medium (vacuum).
1. T-operator formula
Our derivation of the heat transfer spectrum from thefluctuation–dissipation theorem for dipole sources in eachbody follows prior work [38], which we reproduce here forclarity. Consider two bodies A and B in vacuum with generalsusceptibilities V p for p ∈ { A , B } which may be inhomoge-neous, nonlocal, or anisotropic. Maxwell’s equations may bewritten in integral form as | E (cid:105) = G vac | P (cid:105) (C1) | P (cid:105) = | P (0) (cid:105) + V | E (cid:105) (C2)for the fields | E (cid:105) and total polarizations | P (cid:105) in terms of thepolarization sources | P (0) (cid:105) , after defining | E (cid:105) = (cid:20) | E A (cid:105)| E B (cid:105) (cid:21) , | P (cid:105) = (cid:20) | P A (cid:105)| P B (cid:105) (cid:21) (C3) G vac = (cid:20) G vacAA G vacAB G vacBA G vacBB (cid:21) , V − = (cid:20) V − V − (cid:21) (C4)in block form for the material degrees of freedom constitutingeach object. By defining the total T-operator via T − = (cid:20) T − − G vacAB − G vacBA T − (cid:21) (C5)where T − p = V − p − G vac pp , then Maxwell’s equations can beformally solved to yield | E (cid:105) = G vac TV − | P (0) (cid:105)| P (cid:105) = TV − | P (0) (cid:105) (C6)obtained by applying formulas for the block matrix inverse tocompute T . We also define the projection operators, I A = (cid:20) I A
00 0 (cid:21) , I B = (cid:20) I B (cid:21) , (C7)such that (abusing notation) I p is the projection onto the ma-terial degrees of freedom of body p .We consider the energy flow from fluctuating dipolesources only in body A into material degrees of freedom inbody B, noting that reciprocity would yield the same heattransfer if the roles of bodies A and B were interchanged. Thismeans | P (0) (cid:105) = (cid:20) | P (0)A (cid:105) (cid:21) defines the fluctuating sources inbody A. The heat transfer spectrum is the ensemble-averagedwork, denoted by (cid:104)· · ·(cid:105) , done by the field, Φ = 12 Re( (cid:104)(cid:104) I B E , I B J (cid:105)(cid:105) ) (C8)where | J (cid:105) = − i ω | P (cid:105) . Using the Hermiticity and idempotenceof I p yields Φ = − ω ( (cid:104)(cid:104) I B P , E (cid:105)(cid:105) − (cid:104)(cid:104) E , I B P (cid:105)(cid:105) ) , and usingthe results of (C6) gives Φ = − ω (cid:104)(cid:104) P (0)A , I A V − † T † Asym( I B G vac ) TV − I A P (0)A (cid:105)(cid:105) (C9)in terms of the fluctuating sources | P (0)A (cid:105) . As these fluctua-tions are thermal in nature, their correlations are given by thefluctuation–dissipation theorem (cid:104)| P (0)A (cid:105)(cid:104) P (0)A |(cid:105) = 4 πω Asym( V A ) (C10)(suppressing the Planck function Π as it has already been fac-tored to be separate from Φ ), yielding Φ = − π Tr(Asym( V − † A ) I A T † Asym( I B G vac ) TI A ) (C11)as the dressed radiative heat transfer spectrum.To prove equivalence of this expression for Φ to that in-volving only G vac and T p , it is useful to explicitly invokereciprocity: V (cid:62) p = V p , T (cid:62) p = T p , and ( G vac pq ) (cid:62) = G vac qp for p, q ∈ { A , B } , so this also means V † p = V (cid:63)p , T † p = T (cid:63)p ,and ( G vac pq ) † = G vac (cid:63)qp ; this means Asym( A ) = Im( A ) for A ∈ { V p , T p , G vac pp } . This allows for writing the operators Asym( I B G vac ) = (cid:20) − G vac (cid:63) AB / (2i) G vacBA / (2i) Im( G vacBB ) (cid:21) TI A = (cid:20) ( T − − G vacAB T B G vacBA ) − T B G vacBA ( T − − G vacAB T B G vacBA ) − (cid:21) I A T † = (cid:20) ( T − (cid:63) A − G vac (cid:63) AB T (cid:63) B G vac (cid:63) BA ) − ( T − (cid:63) A − G vac (cid:63) AB T (cid:63) B G vac (cid:63) BA ) − G vac (cid:63) AB T (cid:63) B (cid:21) in block matrix form, where the projection onto A al-lows for truncation to the appropriate block column orrow for notational convenience; note that I A T † shouldactually be a row vector, but has been written asa column for ease of reading. Multiplying thesematrices together, it can be noted that G vacBA ( T − − G vacAB T B G vacBA ) − = G vacBA ( I A − T A G vacAB T B G vacBA ) − T A =( I B − G vacBA T A G vacAB T B ) − G vacBA T A = S AB G vacBA T A usingthe definition of the scattering operator S AB = ( I B − G vacBA T A G vacAB T B ) − . Additionally, using the definition T − p = V − p − G vac pp , it is easy to prove that Im( T B ) − T (cid:63) B Im( G vacBB ) T B = T (cid:63) B Im( V − (cid:63) B ) T B , and likewise Im( T A ) − T A Im( G vacAA ) T (cid:63) A = T A Im( V − (cid:63) A ) T (cid:63) A . Thus, the result is (2)in the main text, as expected.
2. Derivation of Green’s function heat transfer formula
Our derivation of the bounds in the main text relies on therelationship between the heat transfer spectrum Φ written interms of the vacuum Green’s function and the T-operators ofindividual objects, to the heat transfer formula [1] Φ = 2 π Tr[Asym( V A ) G † BA Asym( V B ) G BA ] (C12)where G = ( G vac − − V A − V B ) − is the full MaxwellGreen’s function in the presence of both bodies, with the block G BA representing the fields in body B due to dipole sourcesin body A. We start with the T-operator by rewriting Φ( ω ) = 2 π Tr[Asym( V A ) V − † A T † A ( S AB G vacBA ) † T † B × V − † B Asym( V B ) V − T B S AB G vacBA T A V − ] (C13)after using Asym( V − † B ) = V − † B Asym( V B ) V − and Asym( V − † A ) = V − Asym( V A ) V − † A along with invarianceof the trace under cyclic permutations of operator products.From this, it can be seen that the two expressions for Φ( ω ) are guaranteed to be the same if the operator G BA is thesame as V − T B S AB G vacBA T A V − . We use the fact that V − σ = T − σ + G vac σσ to say that V − T B S AB G vacBA T A V − = ( I B + G vacBB T B ) S AB G vacBA ( I A + T A G vacAA ) must hold. To prove that this is equal to G BA , we use thedefinition G = G vac + G vac TG vac (C14)in conjunction with definitions of G vac and T as × blockmatrices in (C5) to write G BA = G vacBA + (cid:2) G vacBA G vacBB (cid:3) T (cid:20) G vacAA G vacBA (cid:21) (C15)for this system. Performing this matrix multiplication, rec-ognizing that G vacBA T A G vacAB ( T − − G vacBA T A G vacAB ) − = S AB − I B , using the fact definition of S AB , and collecting and can-celing terms leads to the proof of the equality G BA = V − T B S AB G vacBA T A V − .
3. Landauer bounds on heat transfer singular values
We now prove that radiative heat transfer between arbitrar-ily shaped bodies can also be expressed as the trace of a trans-mission matrix whose singular values can be bounded above,similar to previously derived bounds in planar media. This re-lation intuitively connects the finite value of the RHT boundsand approximate low rank of G BA , and can be proved as fol-lows. For this, we use the cyclic property of the trace to define Φ = 2 π Tr( Q † Q ) (C16)where Q = Im( V B ) / G BA Im( V A ) / is the heat transmis-sion operator.The definition G − = G vac − − ( V A + V B ) along with thefact that the vacuum Maxwell operator G vac − is real-valuedin position space allows for writing Asym( G ) = G † Asym( V A + V B ) G (C17)which relates dissipation in polarization currents and elec-tromagnetic fields in equilibrium. Additionally, the factthat Asym( V p ) is a Hermitian positive-definite operator foreach body p ∈ { A , B } means it has a unique square root Asym( V p ) / . Rearranging the above equation, multiplyingboth sides by V A ) / , and adding I A to both sidesgives V A ) / G † Asym( V B ) G Asym( V A ) / + 4 Asym( V A ) / G † Asym( V A ) G Asym( V A ) / + 2i(Asym( V A ) / G Asym( V A ) / − Asym( V A ) / G † Asym( V A ) / ) + I A = I A where we recognize the equality Q † Q =Asym( V A ) / G † Asym( V B ) G Asym( V A ) / . Follow-ing this substitution, this may be factored as Q † Q + (cid:16) I A + 2i Asym( V A ) / G AA Asym( V A ) / (cid:17) † × (cid:16) I A + 2i Asym( V A ) / G AA Asym( V A ) / (cid:17) = I A (C18)where G has been replaced by its blocks G AA and G AB dueto multiplications on each each side by Asym( V p ) / for p ∈ { A , B } (and likewise for G † ). This expression is thesum of two Hermitian positive-semidefinite operators equal tothe identity; though this has been done for body A, reciprocityof heat transfer yields a similar expression in terms of the op-erators for body B. Consequently, the singular values of theoperator Q † Q entering the trace expression for Φ( ω ) must allbe less than or equal to / . We emphasize that this deriva-tion is valid for compact or extended structures of arbitrary geometry, without any need to expand heat transfer in termsof incoming and outgoing plane waves specific to translation-ally symmetric systems [5].
4. Single-body thermal radiation from two-body radiative heattransfer
In this section, we prove that the formula for thermal emis-sion of a single body in isolation in vacuum can be derived bystarting from the formula for heat transfer between two bodiesin vacuum under the following conditions. We take body A tobe the thermal emitter in question, while body B is taken tofully surround body A as a shell of inner radius r B and outerradius R B with susceptibility V B = χ B I B , and take the si-multaneous limits Im( χ B ) → and ωr B /c → ∞ constrainedby ω ( R B − r B ) /c → ∞ and ω ( R B − r B ) Im( χ B ) /c → , inwhich case body B takes on the role of a perfectly absorbingmedium.To start, we note that T − = V − − G vacBB in conjunctionwith reciprocity allows for writing (2) as Φ = 2 π Tr[ T B S AB Asym( V − † B ) S A † B T † B × ( G vacAB ) † (Asym( T A ) − T † A Asym( G vacAA ) T A ) G vacAB ] . In the aforementioned size and susceptibility limits for bodyB, T B → so S AB → I B , and T B S AB Asym( V − † B ) S A † B T † B → Asym( V B ) . Using reciprocity, this yields Φ = π Tr[ G vacAB Im( V B ) G vac (cid:63) BA (Im( T A ) − T (cid:63) A Im( G vacAA ) T A ) . Fora system with a general susceptibility V and Maxwell Green’sfunction G = ( G vac − − V ) − , the relations G Im( V ) G (cid:63) = G (cid:63) Im( V ) G = Im( G ) will always hold. Considering B in iso-lation, the simultaneous constrained limits of infinite size andinfinitesimal susceptibility mean that G vacAB Im( V B ) G vac (cid:63) BA → Im( G vacAA ) . Finally, this yields the emission formula Φ = 2 π Tr[Im( G vac )(Im( T ) − T (cid:63) Im( G vac ) T )] (C19)in agreement with the formula derived by Krüger et al [17],where the subscripts A have been dropped as there is onlyone material body under consideration given that body B haseffectively vanished. Appendix D: Nonnegative far-field scattering
In this section, we derive conditions on operators describ-ing the response of bodies p ∈ { A , B } for far-field scat-tering to be nonnegative. In general, given a susceptibil-ity V and an associated T-operator T , the far-field scatteredpower from a given incident field | E inc (cid:105) is ω (cid:104) E inc , (Im( T ) − T (cid:63) Im( V − (cid:63) ) T ) E inc (cid:105) , and for this to be nonnegative, the op-erator Im( T ) − T (cid:63) Im( V − (cid:63) ) T must be positive-semidefinite.This must hold true for each body in isolation, meaning Im( T p ) − T (cid:63)p Im( V − (cid:63)p ) T p must be positive-semidefinite foreach p ∈ { A , B } . However, more conditions can be derivedwhen the two bodies are proximate to each other.As a first step, we show that the operator Y B = ( T − − G vacBA T A G vacAB ) − = T B S AB is an effective T-operator for body0B dressed by the proximity of body A. In particular, the defi-nitions (C4) and (C5) can be plugged into (C6) and rearrangedin order to write | P B (cid:105) = Y B ( G vacBA T A V − | P (0)A (cid:105) + V − | P (0)B (cid:105) ) and if only sources in A are relevant, then we may set | P (0)B (cid:105) → and define an effective incident field | E inc (cid:105) = G vacBA T A V − | P (0)A (cid:105) which depends on multiple scatteringwithin A but not on any properties of B apart from projec-tion onto its volumetric degrees of freedom. This means | P B (cid:105) = Y B | E inc (cid:105) , which is interpreted to mean that the to-tal induced polarization in B arises from the response of Bdressed in the presence of A, namely Y B , acting on the effec-tive incident field | E inc (cid:105) accounting only for body A; this isanalogous to T B which relates the total polarization inducedin B to incident fields in vacuum.Given this and the fact that | E B (cid:105) = V − | P B (cid:105) af-ter setting | P (0)B (cid:105) → , the scattered power only frombody B (in the presence of body A) may be written asthe difference between extinction and absorption powersonly from body B (in the presence of body A), namely ω (cid:0) Im( (cid:104) E inc , P B (cid:105) ) − (cid:104) E B , P B (cid:105) (cid:1) = ω (cid:104) E inc , (Im( Y B ) − Y (cid:63) B Im( V − (cid:63) B ) Y B ) E inc (cid:105) . Nonnegativity of this quantityfor any | P (0)A (cid:105) , or more generally any | E inc (cid:105) , means that Im( Y B ) − Y (cid:63) B Im( V − (cid:63) B ) Y B must be positive-semidefinite.However, nonnegativity of far-field scattering from the sys-tem in general means that upon evaluating the inverse of (C5),the operator Im( T ) − T (cid:63) Im( V − (cid:63) ) T must be positive-semidefinite, which means in turn that each of its diagonalblocks must be positive-semidefinite. Manipulating operatorsallows for showing that the bottom-right block is Im( Y B ) − Y B (cid:0) Im( V − (cid:63) B ) + G vacBA T A Im( V − (cid:63) A ) T (cid:63) A G vac (cid:63) AB (cid:1) Y (cid:63) B , whichinvolves another positive-semidefinite operator (as Im( V − (cid:63) A ) is positive-semidefinite and is multiplied on its left andright by operators which are Hermitian adjoints of eachother, leaving the positive-semidefiniteness unaffected)subtracted from the operator Im( Y B ) − Y (cid:63) B Im( V − (cid:63) B ) Y B .Therefore, the positive-semidefiniteness of Im( Y B ) − Y B (cid:0) Im( V − (cid:63) B ) + G vacBA T A Im( V − (cid:63) A ) T (cid:63) A G vac (cid:63) AB (cid:1) Y (cid:63) B yieldsstronger bounds on the singular values of Y B than thepositive-semidefiniteness of Im( Y B ) − Y (cid:63) B Im( V − (cid:63) B ) Y B alone, because the former also subtracts the absorption inbody A in the presence of body B, whereas the latter does not. Appendix E: Alternative derivations of scalar approximationand Born bound
In this section, we provide alternative derivations of thescalar approximation Φ sc and the Born bound Φ Born . We startby writing the RHT spectrum in (2) as
Φ = 2 πζ A ζ B (cid:13)(cid:13) T B S AB G vacBA T A (cid:13)(cid:13) (E1) in terms of the scattering operator S AB , after using the fact that Asym( V − † p ) = ζ p I p for each body p ∈ { A , B } . The scalarapproximation follows from choosing the T-operator for eachbody to be T p = i τ p I p for some choice of uniform singularvalues τ p ∈ [0 , ζ p ] . This gives Φ sc = τ τ πζ A ζ B (cid:80) i g i (1+ τ A τ B g i ) ,and it can be seen that the contribution from each channel i monotonically increases with increasing τ A τ B . Thus, theonly limit comes from the constraints on τ p , which there-fore must saturate at τ p = ζ p . This is equivalent to directlyplugging in T p = i ζ p I p , which immediately yields Φ sc = π ζ A ζ B (cid:13)(cid:13) S AB G vacBA (cid:13)(cid:13) where S AB = ( I B + ζ A ζ B G vacBA G vacAB ) − inthe scalar approximation. Using the singular value decom-position of G vacBA , particularly the facts that (cid:104) a ( j ) , a ( j ) (cid:63) (cid:105) = (cid:104) b ( j ) , b ( j ) (cid:63) (cid:105) i for each channel j , leads to the expression inthe main text in terms of the singular values ρ j of G vacBA .The Born bound can be derived from Φ sc = π ζ A ζ B (cid:13)(cid:13) S AB G vacBA (cid:13)(cid:13) by taking the scattering operator S AB → I B , which is exactly a Born approximation. Thisexactly gives Φ Born = π ζ A ζ B (cid:107) G vacBA (cid:107) , and substituting thesingular values of G vacBA in the Frobenius norm squared givesthe result in the main text. Alternatively, the operator S AB contributes the denominators (1 + ζ A ζ B g i ) to each channel i , so the replacement S AB → I B makes those denominatorsdisappear. In any case, the result in the main text is recoveredby neglecting multiple scattering between the two bodies. Appendix F: Proof of domain monotonicity of singular values of G vacBA In this section, we prove that the singular values g i of G vacBA are domain monotonic. The singular values of G vacBA are theeigenvalues of G vacBA ( G vacBA ) † . We consider the effects of a per-turbative addition of volume only to body A; a perturbativeeffect on body B can be considered through reciprocity, andthe proof will remain the same. This allows for writing theblock row vector of operators G vacBA = (cid:2) G vacBA G vacB∆A (cid:3) (F1)where G vacBA is the operator propagating fields in vacuum fromthe unperturbed volume A to body B, and G vacB∆A is the opera-tor propagating fields in vacuum from the perturbative volume ∆A to body B. Using reciprocity, we may then write G vacBA ( G vacBA ) † = G vacBA G vac (cid:63) A B + G vacBA G vac (cid:63) ∆AB + G vacB∆A G vac (cid:63) A B + G vacB∆A G vac (cid:63) ∆AB (F2)for which the first term G vacBA G vac (cid:63) A B is the Hermitian positive-semidefinite unperturbed operator, the terms G vacBA G vac (cid:63) ∆AB and G vacB∆A G vac (cid:63) A B vanish because the projections onto the volume A and ∆A are orthogonal to each other, and G vacB∆A G vac (cid:63) ∆AB is the Hermitian positive-semidefinite perturbation. FromRayleigh-Schrödinger perturbation theory, if ρ i is an unper-turbed singular value of G vacBA with | b i (cid:105) being the correspond-ing normalized right singular vector, then the perturbation to1 ρ i is (cid:104) b i , G vacB∆A G vac (cid:63) ∆AB b i (cid:105) , which is nonnegative by virtue ofthe positive-semidefiniteness of G vacB∆A G vac (cid:63) ∆AB . Therefore, anyincrease in the volume of a body will increase the singularvalues of G vacBA . Appendix G: Alternative bounds incorporating far-fieldradiative losses through
Im( G vacAA ) In this section, we derive an alternative bound to (8) thatinvolves
Im( G vacAA ) , thus capturing constraints on scatteringlosses purely from finite object sizes rather than through mul-tiple scattering. Starting from T − = V − − G vacAA , op-erator manipulations allow for writing T − Im( V − (cid:63) A ) T (cid:63) A =Im( T A ) − T A Im( G vacAA ) T (cid:63) A . This allows for rewriting (3) as Φ = 2 πζ B Tr [ G vac (cid:63) AB Y (cid:63) B Y B G vacBA (Im( T A ) − T A Im( G vacAA ) T (cid:63) A )] , which now hides reciprocity, as no similar transformation hasbeen made to eliminate terms giving rise to ζ B . Using the sin-gular value decomposition from (6) for Y B and G vacBA but leav-ing T A general, we find that the constraint on y i is saturatedwhen y i = (cid:0) ζ − + ζ − g i (cid:104) a i , T A T (cid:63) A a i (cid:105) (cid:1) − . Completing thesquare, one may write Im( T A ) − T A Im( G vacAA ) T (cid:63) A = Im( G vacAA ) − / × (cid:34) I A − (cid:18) Im( G vacAA ) / T A Im( G vacAA ) / − i2 I A (cid:19) × (cid:18) Im( G vacAA ) / T A Im( G vacAA ) / − i2 I A (cid:19) (cid:63) (cid:35) Im( G vacAA ) − / . This strongly suggests that the optimal T A should be diago-nalized in the same basis as Im( G vacAA ) , so if we write Im( G vacAA ) = (cid:88) i ρ i | q i (cid:105)(cid:104) q i | , then one may also write T A = i (cid:80) i τ i | q i (cid:105)(cid:104) q i | . This impliesthat the constraint on the singular values of Y B becomes y i = ζ − + ζ − g i (cid:88) j τ j |(cid:104) a i , q j (cid:105)| − so y i depends on τ j for every channel j , not just j = i . Con-sequently, we arrive at the following bound: Φ = 2 π (cid:88) i ζ B g i (cid:16)(cid:80) j τ j (1 − ρ j τ j ) |(cid:104) a i , q j (cid:105)| (cid:17)(cid:104) ζ − ζ B g i (cid:16)(cid:80) j τ j |(cid:104) a i , q j (cid:105)| (cid:17)(cid:105) (G1)for which finding the optimal values of τ i for each chan-nel i requires self-consistently solving a large set of nonlin-ear equations subject to the constraint τ i ≤ ζ A for each i .While this expression should yield tighter bounds on RHT owing to the incorporation of constraints on scattering lossesfor both objects in isolation (in addition to multiple scatter-ing), it appears to be analytically intractable and must there-fore be evaluated numerically, which we leave to future work;that said, numerical solution of the optimal values of τ i forevaluating this bound requires solving a large set of nonlin-ear polynomial equations, which is generally computationallyeasier than brute-force optimization of the RHT spectrum dueto the non-polynomial dependence of Φ on the T-operatorsin general. We point out that in the nonretarded quasistaticlimit, Im( G vacAA ) → , so all of its singular values may betaken to vanish as well; this means that its singular vectors be-come arbitrary, allowing for choosing | q i (cid:105) = | a i (cid:105) . Doing so,the above expression simplifies to Φ = π (cid:80) i ζ B g i τ i (1+ ζ − ζ B g i τ i ) ,from which it can be seen that if the material bound τ i ≤ ζ A issaturated, each contribution is identical to the correspondingcontribution from (8); if the material bound is not saturated,then the optimal τ i = √ ζ − ζ B g i leads to a larger contribu-tion than what we find in (8), yielding an overall looser bound.This corroborates the notion that our present bounds, which donot start with an explicit expression in terms of Im( G vacAA ) , donot fully account for dissipation via far-field radiative losses,while consideration of the singular values ρ i will remedy this. ∗ These two authors contributed equally.[1] W. Jin, S. Molesky, Z. Lin, and A. W. Rodriguez, “Materialscaling and frequency-selective enhancement of near-field ra-diative heat transfer for lossy metals in two dimensions via in-verse design,” Phys. Rev. B , 041403 (2019).[2] D. Costantini, A. Lefebvre, A.-L. Coutrot, I. Moldovan-Doyen,J.-P. Hugonin, S. Boutami, F. Marquier, H. Benisty, and J.-J.Greffet, “Plasmonic metasurface for directional and frequency-selective thermal emission,” Phys. Rev. Applied , 014023(2015).[3] A. W. Rodriguez, O. Ilic, P. Bermel, I. Celanovic, J. D.Joannopoulos, M. Soljaˇci´c, and S. G. Johnson, “Frequency-selective near-field radiative heat transfer between photoniccrystal slabs: A computational approach for arbitrary geome-tries and materials,” Phys. Rev. Lett. , 114302 (2011).[4] T. Asano, M. Suemitsu, K. Hashimoto, M. De Zoysa, T. Shiba-hara, T. Tsutsumi, and S. Noda, “Near-infrared–to–visiblehighly selective thermal emitters based on an intrinsic semicon-ductor,” Science Advances (2016), 10.1126/sciadv.1600499,http://advances.sciencemag.org/content/2/12/e1600499.full.pdf.[5] J. B. Pendry, “Radiative exchange of heat between nanostruc-tures,” Journal of Physics: Condensed Matter , 6621–6633(1999).[6] S.-A. Biehs, E. Rousseau, and J.-J. Greffet, “Mesoscopic de-scription of radiative heat transfer at the nanoscale,” Phys. Rev.Lett. , 234301 (2010).[7] P. Ben-Abdallah and K. Joulain, “Fundamental limits for non-contact transfers between two bodies,” Phys. Rev. B , 121419(2010).[8] O. D. Miller, S. G. Johnson, and A. W. Rodriguez, “Shape-independent limits to near-field radiative heat transfer,” Phys.Rev. Lett. , 204302 (2015). [9] S. Shen, A. Narayanaswamy, and G. Chen, “Surfacephonon polaritons mediated energy transfer betweennanoscale gaps,” Nano Letters , 2909–2913 (2009),http://dx.doi.org/10.1021/nl901208v.[10] R. St-Gelais, B. Guha, L. Zhu, S. Fan, and M. Lipson, “Demon-stration of strong near-field radiative heat transfer between in-tegrated nanostructures,” Nano Letters , 6971–6975 (2014),http://dx.doi.org/10.1021/nl503236k.[11] L. Cui, W. Jeong, V. Fernández-Hurtado, J. Feist, F. J. García-Vidal, J. C. Cuevas, E. Meyhofer, and P. Reddy, “Study ofradiative heat transfer in ångström-and nanometre-sized gaps,”Nature Communications (2017).[12] K. Kloppstech, N. Könne, S.-A. Biehs, A. W. Rodriguez,L. Worbes, D. Hellmann, and A. Kittel, “Giant heat transferin the crossover regime between conduction and radiation,” Na-ture Communications (2017).[13] C. Luo, A. Narayanaswamy, G. Chen, and J. D. Joannopoulos,“Thermal radiation from photonic crystals: A direct calcula-tion,” Phys. Rev. Lett. , 213905 (2004).[14] C. R. Otey, L. Zhu, S. Sandhu, and S. Fan, “Fluctuationalelectrodynamics calculations of near-field heat transfer in non-planar geometries: A brief overview,” Journal of QuantitativeSpectroscopy and Radiative Transfer , 3–11 (2014).[15] A. G. Polimeridis, M. T. H. Reid, W. Jin, S. G. Johnson, J. K.White, and A. W. Rodriguez, “Fluctuating volume-currentformulation of electromagnetic fluctuations in inhomogeneousmedia: Incandescence and luminescence in arbitrary geome-tries,” Phys. Rev. B , 134202 (2015).[16] A. W. Rodriguez, M. T. H. Reid, and S. G. John-son, “Fluctuating-surface-current formulation of radiative heattransfer: Theory and applications,” Phys. Rev. B , 054305(2013).[17] M. Krüger, G. Bimonte, T. Emig, and M. Kardar, “Trace for-mulas for nonequilibrium casimir interactions, heat radiation,and heat transfer for arbitrary objects,” Phys. Rev. B , 115423(2012).[18] A. Lenert, D. M. Bierman, Y. Nam, W. R. Chan, I. Celanovi´c,M. Soljaˇci´c, and E. N. Wang, “A nanophotonic solar ther-mophotovoltaic device,” Nature nanotechnology , 126–130(2014).[19] A. Karalis and J. Joannopoulos, “‘squeezing’ near-field ther-mal emission for ultra-efficient high-power thermophotovoltaicconversion,” Scientific reports , 28472 (2016).[20] B. Guha, C. Otey, C. B. Poitras, S. Fan, and M. Lip-son, “Near-field radiative cooling of nanostructures,”Nano Letters , 4546–4550 (2012), pMID: 22891815,https://doi.org/10.1021/nl301708e.[21] B. D. Boudreau, J. Raja, R. J. Hocken, S. R. Patter-son, and J. Patten, “Thermal imaging with near-field mi-croscopy,” Review of Scientific Instruments , 3096–3098(1997), https://doi.org/10.1063/1.1148248.[22] A. C. Jones and M. B. Raschke, “Thermal infrared near-fieldspectroscopy,” Nano Letters , 1475–1481 (2012), pMID:22280474, https://doi.org/10.1021/nl204201g.[23] A. I. Volokitin and B. N. J. Persson, “Radiative heat transferbetween nanostructures,” Phys. Rev. B , 205404 (2001). [24] G. Domingues, S. Volz, K. Joulain, and J.-J. Greffet, “Heattransfer between two nanoparticles through near field interac-tion,” Phys. Rev. Lett. , 085901 (2005).[25] A. I. Volokitin and B. N. J. Persson, “Near-field radiative heattransfer and noncontact friction,” Rev. Mod. Phys. , 1291–1329 (2007).[26] B. Song, Y. Ganjeh, S. Sadat, D. Thompson, A. Fiorino,V. Fernández-Hurtado, J. Feist, F. J. Garcia-Vidal, J. C. Cuevas,P. Reddy, et al. , “Enhancement of near-field radiative heat trans-fer using polar dielectric thin films,” Nature nanotechnology ,253–258 (2015).[27] M. T. H. Reid, A. W. Rodriguez, and S. G. Johnson,“Fluctuation-induced phenomena in nanoscale systems: Har-nessing the power of noise,” Proceedings of the IEEE , 531–545 (2013).[28] A. Pérez-Madrid, J. M. Rubí, and L. C. Lapas, “Heat trans-fer between nanoparticles: Thermal conductance for near-fieldinteractions,” Phys. Rev. B , 155417 (2008).[29] R. Messina, M. Tschikin, S.-A. Biehs, and P. Ben-Abdallah,“Fluctuation-electrodynamic theory and dynamics of heattransfer in systems of multiple dipoles,” Phys. Rev. B ,104307 (2013).[30] G. Bimonte, “Scattering approach to casimir forces and radia-tive heat transfer for nanostructured surfaces out of thermalequilibrium,” Phys. Rev. A , 042102 (2009).[31] S. Datta, Electronic Transport in Mesoscopic Systems , Cam-bridge Studies in Semiconductor Physics and MicroelectronicEngineering (Cambridge University Press, 1995).[32] J. C. Klöckner, M. Bürkle, J. C. Cuevas, and F. Pauly, “Lengthdependence of the thermal conductance of alkane-based single-molecule junctions: An ab initio study,” Phys. Rev. B ,205425 (2016).[33] O. D. Miller, A. G. Polimeridis, M. T. H. Reid, C. W. Hsu, B. G.DeLacy, J. D. Joannopoulos, M. Soljaˇci´c, and S. G. Johnson,“Fundamental limits to optical response in absorptive systems,”Opt. Express , 3329–3364 (2016).[34] W. Jin, R. Messina, and A. W. Rodriguez, “Overcoming lim-its to near-field radiative heat transfer in uniform planar me-dia through multilayer optimization,” Opt. Express , 14746–14759 (2017).[35] V. Fernández-Hurtado, F. J. García-Vidal, S. Fan, and J. C.Cuevas, “Enhancing near-field radiative heat transfer with si-based metasurfaces,” Phys. Rev. Lett. , 203901 (2017).[36] P. S. Venkataram, S. Molesky, W. Jin, and A. W. Ro-driguez, “Fundamental limits to radiative heat transfer: thelimited role of nanostructuring in the near-field,” (2019),arXiv:1903.07968.[37] S. Molesky, W. Jin, P. S. Venkataram, and A. W. Rodriguez,“Bounds on absorption and thermal radiation from arbitrary ob-jects,” [submitted] (2019).[38] P. S. Venkataram, J. Hermann, A. Tkatchenko, and A. W.Rodriguez, “Phonon-polariton mediated thermal radiation andheat transfer among molecules and macroscopic bodies: Nonlo-cal electromagnetic response at mesoscopic scales,” Phys. Rev.Lett.121