Near field versus far field in radiative heat transfer between two-dimensional metals
NNear field versus far field in radiative heat transfer between two-dimensional metals
Jonathan L. Wise ∗ and Denis M. Basko Universit´e Grenoble Alpes and CNRS, LPMMC, 25 rue des Martyrs, 38042 Grenoble, France
Using the standard fluctuational electrodynamics framework, we analytically calculate the radia-tive heat current between two thin metallic layers, separated by a vacuum gap. We analyse differentcontributions to the heat current (travelling or evanescent waves, transverse electric or magneticpolarization) and reveal the crucial qualitative role played by the dc conductivity of the metals.For poorly conducting metals, the heat current may be dominated by evanescent waves even whenthe separation between the layers greatly exceeds the thermal photon wavelength, and the couplingis of electrostatic nature. For well-conducting metals, the evanescent contribution dominates atseparations smaller than the thermal wavelength and is mainly due to magnetostatic coupling, inagreement with earlier works on bulk metals.
I. INTRODUCTION
Spatially separated objects may exchange heat via elec-tromagnetic fluctuations [1–5]. This radiative heat trans-fer arises due to electric charge density and current fluc-tuations inside the constituting materials, and is usuallydescribed within the phenomenological framework of fluc-tuational electrodynamics (FED) [1, 2, 6], for which thecritical inputs are the material response functions andthe system geometry. It is now well known that in thenear-field limit, energy may tunnel via evanescent elec-tromagnetic waves causing a strong enhancement of theheat transfer, as has been observed experimentally (seethe reviews [7–10] and references therein).Many theoretical works have been dedicated to dif-ferent material systems in the near-field regime ([7–10]and references therein), in which various models for ma-terial response have been employed and different domi-nant channels for heat transfer identified. The commonwisdom is that the evanescent modes dominate the heattransfer when the spatial separation d (cid:28) ¯ λ T ≡ ¯ hc/T , thewavelength of photons at temperature T (here ¯ h and c are the Planck constant and the speed of light, respec-tively, and we set the Boltzmann constant to unity). In-deed, for d > ¯ λ T the evanescent waves with the typicalfrequency ω ∼ T / ¯ h decay exponentially outside the ma-terial, while at d (cid:28) ¯ λ T the region of the wave vectors k occupied by evanescent waves, k ∼ /d , is larger thanthat of travelling states, k ∼ / ¯ λ T [8]. Importance ofmagnetic coupling in the near-field heat transfer betweenwell-conducting metals has been emphasised [11, 12]. Inthe extreme near-field limit, heat transfer due to the elec-trostatic Coulomb interaction has also been studied [13–19].Here we revisit this old problem, focusing on thetwo-dimensional (2D) geometry, and study the radiativeheat current between two thin metallic sheets in vacuumwithin the standard FED framework. We find two qual-itatively different types of behaviour, depending on thevalue of the two-dimensional dc conductivity σ D of the ∗ [email protected] sheets. For poor conductors characterised by the condi-tion G ≡ πσ /c (cid:28) G = ( σ / (cid:112) µ /ε ), the heattransfer turns out to be dominated by the evanescentmodes at distances d extending well beyond ¯ λ T , and themain coupling mechanism in the near field is electrostatic(Coulomb interaction between electrons in the two lay-ers). For G (cid:29)
1, the conventional situation is recovered:the crossover from near to far field occurs at d ∼ ¯ λ T at not too high temperatures, and in a wide range ofparameters the near-field transfer is dominated by mag-netostatic (inductive) coupling between currents in thelayers.The parameter G characterises the impedance mis-match between a two-dimensional (2D) metal and vac-uum; its importance is not restricted to the heat trans-fer problem and is rather general. Notably, two distinctregimes in the behaviour of 2D plasmon polaritons for G < G > G > G <
1, but rather a smooth crossoverbetween the two limiting situations.The rest of the paper is organised as follows. In Sec. IIwe specify the model and sketch the calculation; both arerather standard. In Sec. III we present various regimesof the heat transfer and the associated analytical expres-sions for the heat current, according to the material prop-erties and experimental conditions. In Sec. IV we discussthe relation of our results to the well-studied case of heattransfer between bulk semi-infinite metals, and compareour theory to available experimental results. All detailsof calculations are given in three appendices.
II. THE MODEL
We consider two identical 2D metal sheets held at dif-ferent temperatures T and T , embedded in vacuum andseparated by a gap of width d . A more realistic config-uration would be to place a medium with a dielectricconstant ε in the half-space behind each sheet, since inexperiments the layers are placed on a substrate. For thesake of simplicity, we focus on ε = 1 in most of the paper, a r X i v : . [ c ond - m a t . m e s - h a ll ] J a n and check for the effect of the substrate when specificallyneeded.We model the metal sheets as infinitely thin layers,characterised by a local 2D Drude conductivity, σ ( ω ) = σ − iωτ , (1)with τ being the electron momentum relaxation time, as-sumed to be temperature-independent. This is the case if τ is determined by elastic scattering on static impurities.Eq. (1) neglects (i) the spatial dispersion of the conduc-tivity, and (ii) field variation over the layer thickness.For atomically thin materials, such as doped grapheneor transition metal dichalcogenides, condition (ii) is ir-relevant, and condition (i) holds at distances d (cid:29) √ a (cid:96) ( a and (cid:96) being the 2D screening radius and the elec-tron mean free path, respectively) [18]. For thin butmacroscopic layers of conventional metals, condition (ii)imposes that the thickness must be small compared bothto the typical wavelength of the waves dominating theheat transfer (which may be rather short for evanescentwaves) and to the skin depth at the typical frequency ofthese waves, while condition (i) requires the wavelengthand the skin depth to be longer than the electron meanfree path in the metal.Our calculation of the heat current between the metalsfollows the standard FED procedure. The fluctuating in-plane surface currents j ( α ) ( r , t ) in each sheet obey thefluctuation-dissipation theorem, (cid:104) j l ( r , t ) j m ( r (cid:48) , t (cid:48) ) (cid:105) = δ lm (cid:90) d k dω (2 π ) ¯ hω coth ¯ hω T Re σ ( ω ) × e i k ( r − r (cid:48) ) − iω ( t − t (cid:48) ) , (2)where k is the in-plane two-dimensional wavevector and l, m = x, y label the orthogonal in-plane directions and T = T or T . These currents appear as sources inMaxwell’s equations, whose solution in the presence ofthe conducting sheets determines the fluctuating electricfields E ( r , t ). Then, the heat current J (per unit area)from layer 1 to layer 2 is given by the average Joule losspower (per unit area) (cid:104) ˜ j (2) · E (2) (cid:105) − (cid:104) ˜ j (1) · E (1) (cid:105) , where ˜ j ( α ) is the surface current in layer α , induced by the electricfield E ( α ) in this layer, which, in turn, is produced bythe fluctuating current in the other layer (see AppendixA for explicit expressions, rather standard).We emphasize that for thin layers, the Joule losses (cid:104) ˜ j (2) · E (2) (cid:105) − (cid:104) ˜ j (1) · E (1) (cid:105) are not equal to the averagenormal component of the Poynting vector in the gap be-tween the layers. The reason is that some part of theradiation emitted by layer 1 may pass through layer 2and escape to infinity, and vice versa. Whether this es-caped radiation should be included in the heat current ornot, depends on the precise measurement setup, whichmay collect this escaped radiation or not. Our calcu-lation thus assumes that the escaped radiation is lost.Note that for two semi-infinite metals (the most studiedsetup), everything is collected inside the metals, so thePoynting vector and the Joule losses match exactly. In the planar geometry considered here, the solutionsof Maxwell’s equations are classified by in-plane wavevector k , frequency ω , and two polarisations p = TE , TM– transverse electric and transverse magnetic, respec-tively, for which the electric or the magnetic field vec-tor is parallel to the layers and perpendicular to k . Thecontributions to the heat current from modes with dif-ferent k , ω, p add up independently, so the heat current J ( T , T ) is given by by an integral over k and ω , anda sum over the polarisations. The integral splits in twocontributions: the interior of the light cone, ω > ck hoststravelling modes, while in the region ω < ck the solu-tions are evanescent. The resulting heat current is com-prised of four additive contributions (TM and TE, trav-elling and evanescent). Which contribution dominates,depends on the material conductivity, as well as the sys-tem temperature and length scales.In the extreme near field limit, k (cid:29) ω/c , the TM modefield is mostly electric and longitudinal, while the mag-netic field is smaller by a factor ∼ ω/ ( ck ); these modesrepresent the electrostatic coupling by the Coulomb in-teraction between charge density fluctuations in the twolayers. At the same time, for TE modes the field is mostlymagnetic, while the electric field is smaller by a factor ∼ ω/ ( ck ); these modes represent magnetostatic coupling,where the magnetic field established by transverse cur-rent fluctuations in one layer drives eddy currents in thesecond layer.In Appendix B we perform analytically the k , ω inte-grals and derive simple asymptotic expressions for theheat current according to the separation and the tem-perature. For each expression, we can identify the domi-nant contribution (TM or TE, travelling or evanescent).Our results are approximate; one can describe the heattransfer much more precisely by solving Maxwell’s equa-tions for finite-thickness slabs with a material-specificfrequency dependence of the conductivity and numeri-cally evaluating the integrals, as routinely done in manyworks. However, simple approximate expressions (i) arerather useful when a quick estimate of the heat current isneeded, and (ii) offer a general insight into the dominantphysical mechanisms responsible for the heat transfer andenable one to characterise different possibilities. III. RESULTS
For temperature-independent relaxation time, the heatcurrent naturally splits into the difference J ( T , T ) = J ( T ) − J ( T ). The detailed analysis of different asymp-totic regimes of the k , ω integrals results in severalasymptotic expressions for J ( T ) in different parametricranges of d . The magnitude of the TM and TE travellingcontributions are sensitive to the dimensionless conduc-tivity parameter G = 2 πσ /c so we proceed to presentthe results sequentially for small and large values of G . FIG. 1. The domains of validity for asymptotic expressions, Eqs. (3) and (5) in the parameter plane (1 /d, T ), shown inthe dimensionless variables x ≡ cτ /d , y ≡ T τ / ¯ h . The crossovers between the regimes are governed by the dimensionlessconductivity parameter G ≡ πσ /c , the left and right panels corresponding to G (cid:28)
G (cid:29)
1, respectively. The encircledlabel of each region corresponds to the subscript at J ( T ) in Eqs. (3) and (5). Solid lines indicate crossovers between differentexpressions; straight lines y/x = const are not labeled for readability (the coefficient can be deduced from the endpoints). Theblue line (solid or dashed) corresponds to d ∼ ¯ λ T . The shading color indicates the heat being predominantly carried by TMmodes (red), TE modes (green), or both (white); dominance of travelling waves is indicated by wavy hatching. For
G (cid:28)
1, the asymptotic expressions for J ( T ) are: J lp ( T ) = ζ (3)4 π T ¯ h c G d , (3a) J hp ( T ) = T πτ d L ( G cτ /d ) , (3b) J ld ( T ) = ζ (3)8 π T ¯ h c G d , (3c) J hd ( T ) = 116 π G cd T, (3d) J lte ( T ) = π G ¯ h c T ln 1 G , (3e) J hte ( T ) = 14 π G c τ T ln 1 G , (3f)valid in the corresponding regions of the (1 /d, T ) plane,schematically shown in Fig. 1 (left). The contribu-tions given by Eqs. (3a)–(3d), with labels correspondingto low-temperature plasmonic, high-temperature plas-monic, low-temperature diffusive, high-temperature dif-fusive, are the TM evanescent contributions that remainin the Coulomb limit and were calculated in Ref. [18].Equations (3e) and (3f) (with labels corresponding tolow-temperature travelling electic and high-temperaturetravelling electric) are the travelling TE contributionswhich dominate over the travelling TM contributions bythe logarithmic factor ln(1 / G ). In Eqs. (3), ζ ( x ) is theRiemann zeta function, and L ( x ) is a slow logarithmicfunction, approximately given by [18] L ( x ) ≈ x x ) / ln(1 + ln x ) . (4) For G (cid:29)
1, in addition to the expressions given inEqs. (3a)–(3c) we also have: J le ( T ) = π G ¯ h c T ln ¯ λ T G d , (5a) J he1 ( T ) = 14 π G c τ T ln cτ G d , (5b) J he2 ( T ) = ζ (3)16 π c G d T, (5c) J lt ( T ) = π T ¯ h c G , (5d) J it ( T ) = 112 T ¯ hc τ , (5e) J ht ( T ) = (cid:32) √ π + 14 π (cid:33) G c τ T, (5f)valid in the corresponding regions of the (1 /d, T ) plane,schematically shown in Fig. 1 (right). The contributionsgiven by Eqs. (5a)–(5c) are the TE evanescent contri-butions, while Eqs. (5d)–(5f) are the sums of travellingcontributions from both polarisations which are of thesame order.For G (cid:29)
1, the travelling channels support resonantFabry-Perot (FP) modes. In (lt) and (it) regions, manysharp FP modes contribute significantly to the heat cur-rent. In the high-temperature case (ht) the FP modes areoverdamped since the conductivity σ ( ω ) becomes smallat high frequencies. For temperatures lower than thefirst mode cutoff energy, T (cid:28) π ¯ hc/d , the contributionsfrom the FP modes are exponentially suppressed. How-ever, the prefactor in front of the small thermal exponen-tial turns out to be larger than the evanescent contribu-tion (5c) in (he2) region. Thus, the FP additive contri-bution is potentially significant for cτ /d < √G , where itis dominated by the first FP mode: J FP1 ( T ) = πcT [2 G + ( πcτ /d ) ] d e − π ¯ hc/ ( T d ) . (6)In Fig. 1, the areas with wavy hatching indicate the re-gions where the heat transfer is dominated by travellingwave contributions. For G (cid:28)
1, the evanescent wavesdominate at separations up to d ∼ G − / ¯ λ T , paramet-rically larger than the commonly used condition for thenear field, d (cid:28) ¯ λ T [ y (cid:28) x in Fig. 1 (left)]; the reason forsuch behaviour is that the low-temperature TM evanes-cent contribution is determined by k (cid:28) /d for whichthe exponential suppression is not efficient. Moreover,for G (cid:28)
G (cid:29)
1, the commonly used inequality d (cid:28) ¯ λ T doesbecome the accurate condition for evanescent contribu-tion dominance, except for high temperatures where theDrude conductivity is suppressed by high frequency. In alarge part of the near-field region of the parameter planethe heat current is governed by TE evanescent modes,which correspond to magnetostatic (inductive) couplingbetween the layers. As discussed in Refs. [11, 12] forbulk metals, large conductivity leads to efficient screen-ing of the electric fields, so the magnetostatic coupling be-comes more important. The electrostatic coupling takesover only at very short distances or low temperatures, d (cid:28) ¯ λ T / ( π G ) , determined by J le ( T ) ∼ J ld ( T ).However, for very small d the finite layer thickness maybecome important, and/or the assumption of the local re-sponse, Eq. (1), may break down. Taking, for example,a 10 nm-thick gold film with the bulk plasma frequency ω p = 0 . × s − and relaxation time τ = 6 fs [26]gives G ≈ .
6. Since ¯ λ T = 7 . µ m at T = 300 K, in suchstructure the crossover to electrostatics occurs at a fewnanometers. We note that in the (ld) regime, the heattransfer is mainly determined by rather small wave vec-tors k ∼ ( G ¯ λ T d ) − / [18], so that even at d = 1 nm weobtain 1 /k > ∼
100 nm, and the local response assumptionshould still be formally valid. However, at nanometricdistances other physical effects may come into play (elec-tron or phonon tunnelling, surface roughness, etc.), sofor ultrathin films of conventional metal we expect theCoulomb mechanism to be relevant mostly at low tem-peratures.
IV. DISCUSSIONA. Comparison to the bulk case
The results presented in the previous section showtwo qualitatively different pictures of the near-field heattransfer between two metallic layers, depending on thevalue of their dimensionless 2D dc conductivity: for
G (cid:28)
1, the heat transfer is mostly due to electrostaticcoupling between the layers, up to distances significanlyexceeding ¯ λ T , while for G (cid:29) d ∼ ¯ λ T , inclose analogy with earlier results on bulk metals. Thispicture is consistent with the results of Ref. [27] wherethe very same problem of radiative heat transfer betweenparallel 2D layers was studied numerically. There, a dis-tinction was made between thin and thick metallic films.In this formulation G is proportional to the layer thick-ness h (in the local approximation, the 2D conductivity issimply σ D = σ D h , where σ D is the bulk conductivity).In Ref. [27], the heat transfer between two theoreticallyimagined atomic monolayers of silver, described by a 2DDrude model with G ≈
2, is found to be driven by TMevanescent waves, while for thicker films it is TE evanes-cent waves.The peculiarity of the 2D geometry is that the 2D con-ductivity can be compared to two universal scales. Oneis the speed of light, hence the dimensionless parame-ter G = 2 πσ D /c we introduced earlier. The other uni-versal scale is the conductance quantum, e / (2 π ¯ h ). For σ D < ∼ e / (2 π ¯ h ), or G < ∼ e / (¯ hc ) ≈ / / < ∼ G (cid:28) σ D has the dimensionality of theinverse time, so that 1 / (4 πσ D ) (in CGS units, while inSI it is ε /σ D ) has a meaning of the RC time neededto dissolve a charge density perturbation. In conven-tional metals this time scale is extremely short (in theattosecond range). Still, one can compare 4 πσ D to otherscales. One is the electron relaxation time τ ; typically,4 πσ D τ = ω p τ (cid:29) ω p being the bulk plasma fre-quency). Moreover, at T (cid:28) ¯ h/τ the relaxation timedrops out of the problem, so one cannot construct a di-mensionless parameter out of σ D , which could producedifferent “asymptotic maps” of the kind shown in Fig. 1.The bulk case turns out to be somewhat similar to the2D case with G (cid:29)
T τ (cid:28) ¯ h (the deriva-tion can be found in Ref. [2], we also give it in Ap-pendix C): J a ( T ) = π
60 ¯ h ( T / ¯ h ) (2 πσ D ) d ln 2 πσ D T / ¯ h , d (cid:28) δ T ¯ λ T , (7a) J b ( T ) = ζ (3)4 π πσ ¯ h ( T / ¯ h ) c , δ T ¯ λ T (cid:28) d (cid:28) δ T , (7b) J c ( T ) = 3 ζ (3)4 π c T πσ D d ln dδ T , δ T (cid:28) d (cid:28) (¯ λ T δ T ) / , (7c) J d ( T ) = 75 ζ (7 / √ π ¯ h ( T / ¯ h ) / √ πσ cd , (¯ λ T δ T ) / (cid:28) d (cid:28) ¯ λ T , (7d) J e ( T ) = 35 ζ (9 / π / ¯ h ( T / ¯ h ) / √ πσ c , ¯ λ T (cid:28) d, (7e)where the parametric intervals of d are conveniently de-fined in terms of two length scales: the thermal wave-length ¯ λ T = ¯ hc/T and the normal skin depth at the ther-mal frequency, δ T = c/ (cid:112) πσ D T / ¯ h (cid:28) ¯ λ T . The shortest-distance expression (7a) is determined by the TM evanes-cent contribution and corresponds to the Coulomb limit(indeed, it does not contain the speed of light); how-ever, the length scale δ T / ¯ λ T is extremely short: for4 πσ D = 10 s − at T = 300 K, we have δ T = 0 . µ mand ¯ λ T = 7 . µ m, so δ T / ¯ λ T ∼ k veryclose to ω/c , and the fields vary weakly across the gap sothere is no sharp physical distinction between travellingand evanescent waves. Finally, Eq. (7e) comes from theTE and TM travelling waves and is contributed by manyFabry-Perot modes inside the gap.It is easy to see that by the order of magnitude,Eqs. (7b), (7c), (7d) and (7e) can be obtained fromEqs. (5a), (5c), (3d) and (5d), respectively, by replac-ing G = 2 πσ D /c → πσ D δ ω /c , where the skin depth δ ω = c/ √ πσ D ω corresponds to the typical frequencyscale determining the integral: it is δ T for Eqs. (7b),(7d) and (7e), determined by frequencies ω ∼ T / ¯ h , and δ ω ∼ d for Eq. (7c), where the frequency integral is log-arithmic, with the lower cutoff corresponding to δ ω ∼ d (see Appendix C for details). This replacement roughlycorresponds to modelling the semi-infinite metal as an ef-fective metallic layer whose thickness corresponds to thefield penetration depth. Such effective layer is charac-terised by the dimensionless G eff ∼ (cid:112) πσ D /ω , so that G eff (cid:29) δ ω sometimes makes the convergence scale of the frequencyintegral different from the case of fixed layer thickness. B. Comparison to experiments
Values G < ∼ ε = 11 .
7) and sepa-rated by a 400 nm wide vacuum gap. The Fermi en-ergy of 0 .
27 eV and the relaxation time τ = 100 fs give G = 0 .
6. The linear thermal conductance per unit area dJ/dT = 30 W m − K − was measured around roomtemperature. These conditions correspond to the high-temperature plasmon regime, Eq. (3b), where the sub-strate dielectric constant ε enters only inside the loga-rithmic function L [18]. Setting L = 1 in Eq. (3b) gives dJ/dT = 11 W m − K − , which agrees by order of mag-nitude with the experimental value.Thin layers of conventional metals are typically char-acterised by G (cid:29)
1. Several experiments have been re-ported in the literature. In each case, it is important tocompare the layer thickness h to the the skin depth δ ω atthe relevant frequency, to ensure the layers should cor-respond to the 2D limit, rather than the bulk one (thelatter being the case of Refs. [31, 32]).Heat transfer in a wide range of interlayer separationsand temperatures was studied in Ref. [33] for two 150 nmthick tungsten layers on alumina substrates. The mea-sured dc conductivity of the material 4 πσ = 0 . × s − (constant in the temperature range of the exper-iment) corresponds to a value of the dimensionless con-ductivity parameter G ≈ T = 40Kis δ T = 240 nm, and even longer at lower temperatures,so the layers are close to the 2D limit. The separationbetween the layers was varied over d = 1 − µ m, whilethe temperatures were T = 5 K and T = 10 −
40 K, cor-responding to regions (he2) and (lt) in Fig. 1 (right). Itcan be easily checked that in these regions, the dielectricsubstrate plays no role as long as √ ε (cid:28) G , which clearlyholds here. Although the numerical calculation account-ing for the finite layer thickness does better in closelymatching the experimental points (see Fig. 2 of Ref. [33]),our simple expressions (5c) and (5d) (i) agree with theobserved values within a factor of 3 without any fit-ting parameters, (ii) give the correct distance dependencethroughout the experiment, (iii) capture the observed ap-proximate collapse of the rescaled data for J ( T ) /T ona function of a single variable T d , and (iv) correctly pre-dict the separation d ≈ . λ T , at which the crossoverbetween the near-field and the far-field regimes occurs, J he2 ( T ) = J le ( T ).A recent publication [34] presents measurements ofheat transfer between two aluminium films of varyingthicknesses h = 13 −
79 nm, separated by a fixed vac-uum gap d = 215 nm and attached to silicon substrates.The experiment was performed around room tempera-ture with one film being heated such that ∆ T = 25 −
65 K.Taking the values ω p = 1 . × s − and τ ≈ δ T ≈
50 nm and
G ≈
40 for the thinnest layer with h = 13 nm.Then Eq. (5c) predicts dJ/dT = 250 W m − / K, whichagrees in order of magnitude with the reported value, dJ/dT = 60 W m − / K.An intriguing feature of the results reported in Ref. [34]is the independence of dJ/dT of the layer thickness. Thisagrees neither with our 2D expressions, nor with the moreprecise simulations done in Ref. [34]. All theoretical re-sults point to a non-monotonic dependence of the heatcurrent on the layer thickness or dc conductivity [the lat-ter is also true for the bulk limit expressions (7)]. Furtherexperimental investigations of this dependence would beinteresting.
V. CONCLUSIONS
In this paper, we have performed an analytical calcula-tion of the radiative heat current between two thin metal-lic layers, using the standard framework of fluctuationalelectrodynamics and a local 2D Drude model for the elec-tromagnetic response of each layer. We have identifiedtwo different classes of such structures, distinguished bythe dimensionless 2D dc conductivity G = 2 πσ D /c . Forpoor conductors with G (cid:28)
1, typically represented byatomically thin 2D materials, the heat transfer is dom-inated by evanescent modes at distances d extendingwell beyond ¯ λ T , and the main coupling mechanism inthis near-field regime is the Coulomb interaction betweenelectrons in the two layers. Good conductors with G (cid:29) d ∼ ¯ λ T at nottoo high temperatures, and the near-field transfer is dom-inated by magnetostatic (inductive) coupling between thelayers in a wide range of parameters.We have derived several simple approximate asymp-totic expressions for the heat current valid in differ- ent parametric ranges of interlayer separation distanceand temperature. Comparing these expressions with theavailable experimental data, we saw that they give validorder-of-magnitude estimates of the heat current and cor-rectly capture its dependence on the distance and tem-perature. Better agrreement with the experimental re-sults can be reached by a more detailed modelling ofeach system geometry and the dielectric response, whichis strongly system-specific and lies beyond the scope ofour work. Still, our approximate results offer a useful in-sight into the main physical mechanisms responsible forthe heat transfer. ACKNOWLEDGMENTS
We thank J.-J. Greffet, J. Pekola, B. Van Tiggelen, andC. Winkelmann for helpful and stimulating discussions.This project received funding from the European Union’sHorizon 2020 research and innovation programme un-der the Marie Sk(cid:32)lodowska-Curie Grant Agreement No.766025.
Appendix A: Explicit general expression for theheat current between two thin metallic sheets
We solve Maxwell’s equations for the monochro-matic components of the electric and magnetic field, E k ω ( z ) e i kr − iωt and B k ω ( z ) e i kr − iωt in the planar geom-etry with the two metallic sheets placed at z = z , z with z − z = d , while the position-dependent dielec-tric constant ε ( z ) accounts for whatever (non-magnetic,isotropic) dielectric medium surrounds the layers: (cid:18) i k + e z ∂∂z (cid:19) × E k ω = iωc B k ω , (A1) (cid:18) i k + e z ∂∂z (cid:19) × B k ω = 4 πc (cid:88) α =1 , δ ( z − z α ) (cid:16) j ( α ) k ω + ˜ j ( α ) k ω (cid:17) − iωc ε ( z ) E k ω , (A2)where e z is the unit vector in the z direction, perpendic-ular to the layers. The surface current in each layer α =1 , j ( α ) k ω = σ α ( ω ) E k ω ( z α )is the induced current due to the electric field, while thefluctuating currents j ( α ) k ω = ( j ( α ) − k , − ω ) ∗ are complex Gaus-sian random variables with the correlator determined bythe fluctuation-dissipation theorem (2): (cid:104) j ( α ) k ω,l j ( α (cid:48) ) k (cid:48) ω (cid:48) ,m (cid:105) = (2 π ) δ ( k + k (cid:48) ) δ ( ω + ω (cid:48) ) δ αα (cid:48) δ lm × ¯ hω coth ¯ hω T α Re σ α ( ω ) (A3)Because of δ lm on the right-hand side of this equation,current fluctuations are independent for any two orthog-onal directions, so it is convenient to pass to the longitu-dinal and transverse basis ( p and s polarisations, respec-tively): j ( α ) k ω = j ( α ) k ω k k + j ( α ) k ω e z × k k . (A4)In this basis the solutions of Maxwell’s equations decou-ple into transverse magnetic (TM) and transverse elec-tric (TE) modes, whose contribution to the heat currentis simply additive.To model different metal sheets mounted on identi-cal dielectric substrates separated by vacuum, we take ε ( z < z < z ) = 1, ε ( z < z ) = ε ( z > z ) = ε >
1. This leads to the spatial dependence of the electricand magnetic fields ∝ e i kr ± iq z z for z < z < z , and ∝ e i kr + iq (cid:48) z z , e i kr − iq (cid:48) z z for z > z and z < z , respectively.Here we defined q z = (cid:26) (cid:112) ω /c − k sign ω, | ω | > ck,i (cid:112) k − ω /c , | ω | < ck, (A5a) q (cid:48) z = (cid:26) (cid:112) εω /c − k sign ω, √ ε | ω | > ck,i (cid:112) k − εω /c , √ ε | ω | < ck. (A5b)At | ω | > ck , the metallic layers are coupled by travellingwaves, while for | ω | < ck the solutions in the gap areevanescent waves, where the fields’ strength decays awayfrom the layers. The solutions are matched at z = z and z = z using the standard boundary conditions: conti-nuity of the in-plane component of the electric field E (cid:107) ,and a jump in the magnetic field in-plane component,determined by the total surface current (the fluctuatinngsources as well as the induced current σ E (cid:107) ).The heat current from, say, sheet 1 to the sheet 2is given by the average Joule loss power per unit area, J ( T , T ) = (cid:104) ˜ j (2) · E (cid:107) ( z ) (cid:105) − (cid:104) ˜ j (1) · E (cid:107) ( z ) (cid:105) , determinedunambiguously due to the continuity of E (cid:107) ( z ). For atemperature independent relaxation time this heat cur-rent splits into J ( T , T ) = J ( T ) − J ( T ), where J ( T ) = ∞ (cid:90) dω π ¯ hωe ¯ hω/T − (cid:90) d k (2 π ) (cid:88) j = p,s a j a j | e iq z d | | − r j r j e iq z d | (A6)is expressed in terms of reflectivities r αj and emissivities a αj for the p and s polarisations: r αp = q z − q (cid:48) z /ε + 4 πσ α q z q (cid:48) z / ( εω ) q z + q (cid:48) z /ε + 4 πσ α q z q (cid:48) z / ( εω ) , (A7a) a αp = 4 | q z || q (cid:48) z /ε | (4 π Re σ α /ω ) | q z + q (cid:48) z /ε + 4 πσ α q z q (cid:48) z / ( εω ) | , (A7b) r αs = q z − q (cid:48) z − πωσ α /c q z + q (cid:48) z + 4 πωσ α /c , (A7c) a αs = 4 | q z | (4 πω Re σ α /c ) | q z + q (cid:48) z + 4 πωσ α /c | . (A7d) The emissivities can also be written as a αp = (1 − | r αp | ) θ ( | ω | − ck ) + 2 Im r αp θ ( ck − | ω | ) − (cid:12)(cid:12)(cid:12)(cid:12) q (cid:48) z εq z (cid:12)(cid:12)(cid:12)(cid:12) | t αp | θ ( √ ε | ω | − ck ) , (A8a) a αs = (1 − | r αs | ) θ ( | ω | − ck ) + 2 Im r αs θ ( ck − | ω | ) − (cid:12)(cid:12)(cid:12)(cid:12) q (cid:48) z q z (cid:12)(cid:12)(cid:12)(cid:12) | t αs | θ ( √ ε | ω | − ck ) , (A8b)where θ ( x ) is the Heaviside step function, and t αj are thetransmittivities: t αp = 2 q z q z + q (cid:48) z /ε + 4 πσ α q z q (cid:48) z / ( εω ) , (A9a) t αs = 2 q z q z + q (cid:48) z + 4 πωσ α /c . (A9b)Note the difference between Eqs. (A8) and Eq. (2) ofRef. [27], where the third term is absent in both polar-isations. Without the third term, Eq. (A6) gives theaverage value of the Poynting vector in the gap betweenthe two layers, and also counts the heat flux which is notabsorbed by the metal, but irradiated to infinity behindit, due to the finite transmission. Eqs. (A8) without thethird term originally appeared in Ref. [36] for the prob-lem of heat transfer between two semi-infinite materials.In that geometry, all heat flux transmitted through thesurface is eventually absorbed by the material. In thethin layer geometry, whether the transmitted flux is de-tected or not, depends on the specific experimental mea-surement scheme. In our calculation, we assume thatthe transmitted radiation is lost, and thus use the fullEqs. (A8). Appendix B: Derivation of asymptotic expressionsfor the heat current between two thin metallic sheets
Here we derive asymptotic expressions for J ( T ) inthe specific case of identical sheets embedded in vac-uum [ σ ( ω ) = σ ( ω ) and ε = 1] and compute separatelythe travelling and evanescent wave contributions for eachof the two polarisations. We quantify the contributionmade by each wave type and polarisation in each regionof the (1 /d, T ) parameter plane, before comparing thesize of the additive contributions and identifying whichare dominant. It is convenient to introduce the dimen-sionless parameters x ≡ cτ /d and y ≡ T τ / ¯ h , as well asdimensionless integration variables: ξ = | q z | cτ insteadof k [noting that k dk = ξ dξ/ ( cτ ) ], and η = ωτ . Forthe travelling waves, the integration is over the region0 < ξ < η < ∞ , while for the evanescent waves it is0 < ξ, η < ∞ .
1. TM travelling contribution
In the dimensionless variables, the TM travelling con-tribution to Eq. (A6) can be rewritten exactly as J tTM = ¯ h G π c τ (cid:90) ∞ η dηe η/y − (cid:90) η ξ dξD p + D p − , (B1a) D p ± ≡ | η (1 − iη ) + G ξ (1 ± e iξ/x ) | . (B1b)The case G (cid:28) ε = 1one can neglect the reflection coefficients in the denom-inator of Eq. (A6), and simply set G → ξ < η . This gives J tTM = ¯ hσ c τ (cid:90) ∞ η dη (1 + η ) ( e η/y − (cid:26) π G T / (60¯ h c ) , y (cid:28) , G T / (16 πc τ ) , y (cid:29) . (B2)For G (cid:29)
1, each layer behaves at low frequency as awell-reflecting mirror, so the structure may host Fabry-Perot modes. The Fabry-Perot modes manifest them-selves as deep minima in D p ± at specific values of ξ/x = π, π, π, . . . . These minima are important when G ξ (cid:29) η (cid:112) η , which is precisely the condition of good reflec-tion. Thus, a much more elaborate analysis is needed toevaluate the integral.Let us focus on the contributions from the region ξ (cid:29) x , when many modes contribute, and even if they areoverdamped, e iξ/x oscillates fast. In the general case (A6)we average over the fast oscillations in the denominatorwhich leads to the simple replacement [37]: a j a j | − r j r j e iq z d | → a j a j − | r j | | r j | , (B3)valid as long as a αj and r αj are smooth functions of q z on the scale q z ∼ /d .Applying this averaging to the contribution inEq. (B1a) leads to J tTM = ¯ h G π c τ (cid:90) ∞ η dηe η/y − × η (cid:90) ξ dξ [ η + η + 2 G ξη + 2( G ξ ) ][ η + η + 2 G ξη ] . (B4)Note that x dropped out, and enters only through thecondition ξ (cid:29) x . Note also that the ξ integral is alwaysdetermined by the upper limit ξ ∼ η . As for the η in-tegral, it may converge at η ∼ y when cut off by theBose function, or, for too large y , it may be cut off byother factors in the denominator at some η (cid:28) y . In thislatter case, one can expand the exponential in the Bosefunction, which becomes just y/η . We can identify threeregions in y . (i) For y (cid:28) √G , the integrals separate and converge at ξ ∼ η ∼ y , so the fast oscillation condition is x (cid:28) y : J tTM = ¯ h G π c τ (cid:90) ∞ η dηe η/y − (cid:90) η ξ dξ G ξ ) η = π T
60 ¯ h c G . (B5)(ii) For √G (cid:28) y (cid:28) G , we keep 2( G ξ ) in the first bracketand η in the second one (again, oscillations are fast when x (cid:28) y ): J tTM = ¯ h G π c τ (cid:90) ∞ η dηe η/y − (cid:90) η ξ dξ G ξ ) η = T
24 ¯ hc τ . (B6)(iii) For y (cid:29) G , we expand the Bose function, the integralconverges at η ∼ G (it is convenient to write ξ = uη ); theoscillations are fast when x (cid:28) G : J tTM = ¯ h G π c τ (cid:90) ∞ yη dη (cid:90) η u du [ η + 2( G η ) u ] η = √ G T πc τ . (B7)Let us now pick the contributions from ξ (cid:28) x . Then, e iξ/x can be expanded (we again write ξ = uη ): J tTM = ¯ h G π c τ (cid:90) ∞ η dηe η/y − ×× (cid:90) u du [(1 + 2 G u ) + η ][1 + η (1 + G u /x ) ] . (B8)There are three possible cutoff scales for η : y , 1 + 2 G u ,and (1+ G u /x ) − . Which one of the three is effective, de-pends on the positioning of y with repect to other scales.Again, three cases arise.(iv) For y (cid:28)
1, we can neglect η in the first bracket inthe denominator, so the η integral converges at η ∼ y . Inthe second bracket, η plays a role only if G u /x (cid:29)
1, sothe second bracket can be approximated as 1+( G u η/x ) for any G u /x . We also assume that G u (cid:29)
1, whichwill be verified afterwards. Then the denominator be-comes 4 G u [1 + ( G u η/x ) ], so the u integral convergesat u ∼ min { , (cid:112) x/ ( G y ) } , giving J tTM = ¯ h π G cτ d (cid:90) ∞ η dηe η/y − G ηx = (cid:26) π T / (120 ¯ h c ) , y (cid:28) x/ G ,ζ (3) T / (8 π ¯ h c G d ) , y (cid:29) x/ G . (B9)The u integral converges at u ∼ u ∼ (cid:112) x/ ( G y ) inthe two cases. In the first case, y (cid:28) x/ G , the assumption G u (cid:29)
1, as well as the condition to expand the exponen-tial, uy/x (cid:28)
1, are satisfied automatically. In the secondcase, y (cid:29) x/ G , both conditions translate into y (cid:28) G x .(v) For y (cid:29) y (cid:28) G u we still have η ∼ y , so the de-nominator can be approximated as 4 G u η (1 + G u /x ) : J tTM = ¯ h π c τ (cid:90) ∞ η dηe η/y − (cid:90) u du [1 + ( G /x ) u ] = T / ¯ h cτ ( cτ + G d ) . (B10)Since the convergence occurs at u ∼ min { , (cid:112) x/ G} , η ∼ y , the assumption y (cid:28) G u is satisfied if y (cid:28) min {√G x, G} ; if so, the condition uy/x (cid:28) (cid:28) y (cid:28) min {√G x, G} .(vi) For y (cid:29) , G u , the Bose function is y/η , so we inte-grate over η exactly (convergence at η ∼ G u ), andobtain J tTM = G T πc τ (cid:90) u du (1 + 2 G u )[1 + ( G /x ) u ] = (cid:26) T x / / (16 c τ G / ) , x (cid:28) G , G T / (12 πc τ ) , x (cid:29) G , (B11)the convergence occurring at u ∼ min { (cid:112) x/ G , } . At x (cid:28) G the condition to expand e iξ/x is not fulfilled, sincewe automatically have ξ/x = uη/x ∼
1. At x (cid:29) G , wehave G u (cid:29) uη/x ∼ G /x , so thesecond expression Eq. (B11) is valid at x, y (cid:29) G .We schematically show the regions of validity ofEqs. (B5)–(B11) in the ( x, y ) plane in Fig. 2(a). In theoverlapping region at y (cid:29) x both ξ (cid:28) x and ξ (cid:29) x contributions are valid, but the Fabry-Perot contribu-tions from ξ (cid:29) x naturally dominate. At y (cid:28) x theFabry-Perot contributions are suppressed as e − πx/y = e − π ¯ hc/ ( T d ) , since the temperature is lower than the firstFabry-Perot mode energy π ¯ hc/d . Nevertheless, it turnsout that the prefactor in front of the exponential is large,so the contribution from the first mode (the one withthe weakest exponential), coming from the narrow regionaround ξ = η = πx [see Fig. 2(b)] should be includedtogether with the contribution from ξ (cid:28) x , as long as x, y (cid:28) G (otherwise, the mode is overdamped because oflow reflectivity).To pick up the first Fabry-Perot mode contribution, weapproximate the Bose function by e − η/y and set η = πx everywhere else in the integrand, which is a smooth func-tion of η . We also set ξ = πx everywhere in the integrandexcept the exponential e iξ/x in D p + [Eq. (B1b)]. Then wefind the minimum of D p + as a function of ξ , reached at ξ min = πx (1 − x/ G ) + O ( x / G ), and approximate nearthe minimum D p + = ( π G x ) (cid:12)(cid:12)(cid:12)(cid:12) − iπx G + 1 + e iξ/x (cid:12)(cid:12)(cid:12)(cid:12) ≈ ( πx ) (cid:18) π x G (cid:19) + ( π G ) ( ξ − ξ min ) . (B12)Then, the integration over πx < η < ∞ and −∞ <ξ − ξ min < ∞ gives J tTM = πcT d [2 G + ( πcτ /d ) ] e − π ¯ hc/ ( T d ) . (B13) FIG. 2. (a) Regions of validity of Eqs. (B5)–(B11) in the ( x, y )plane. In the shaded regions there are two valid contributions.(b) The ( ξ, η ) plane with the integration domain ξ < η (theshaded area does not belong to the integration domain). Thehatched area at ξ (cid:28) x contributes to Eqs. (B9), (B10).
2. TE travelling contribution
The TE travelling contribution to Eq. (A6) can berewritten exactly as J tTE = ¯ h G π c τ (cid:90) ∞ η dηe η/y − (cid:90) η ξ dξD s + D s − , (B14a) D s ± ≡ | ξ (1 − iη ) + G η (1 ± e iξ/x ) | . (B14b)For G (cid:28)
1, we may not simply set
G → ξ →
0. To see how the diver-gence is cut off, we note that convergence scale of the η integral is the same as in the TM case: η ∼ y if y (cid:28) η ∼ y (cid:29)
1. This gives the small- ξ cutoff scales0 ξ ∼ G η and ξ ∼ G , respectively. As a result, J tTE = ¯ h G π c τ (cid:90) ∞ η dη (1 + η ) ( e η/y −
1) ln η min {G , G η } = (cid:26) [( π G T ) / (15 ¯ h c )] ln(1 / G ) , y (cid:28) , [( G T ) / (4 πc τ )] ln[( T τ ) / ( G ¯ h )] , y (cid:29) , (B15)The overall map of behaviours in parameter space istherefore equivalent to the TM travelling case given inEq. (B2), but the TE contribution (B15) is always dom-inant due to the logarithmic factors.The calculation for G (cid:29) ξ (cid:29) x so the exponentials e iξ/x oscillate fast,the averaged contribution from Eq. (B14a) via Eq. (B3)is given by J tTE = ¯ h G π c τ (cid:90) ∞ η dηe η/y − (cid:90) η ξ dξξ ( η + 1) + 2 G η × ξ ( η + 1) + 2 G ξη + 2 ( G η ) . (B16)At low frequency the system Fabry-Perot modes are in-dicated, as in the TM case, in the minima in D s ± , thistime important when G η (cid:29) ξ (cid:112) η . The integral in η may again converge at η ∼ y due to the Bose function,or something else if y is too large. We may identify thesame regions as in the TM case.(i) For y (cid:28) √G , we have that ξ ∼ η ∼ y , so the fastoscillation condition is x (cid:28) y , and we may neglect allterms in the denominator containing ξ : J tTM = ¯ h G π c τ (cid:90) ∞ η dηe η/y − (cid:90) η ξ dξ G η ) = π T
180 ¯ h c G . (B17)(ii) For √G (cid:28) y (cid:28) G , we keep ξη in the denominator inthe first line of Eq. (B16) and 2 ( G η ) in the second line(again, oscillations are fast when x (cid:28) y ): J tTE = ¯ h G π c τ (cid:90) ∞ η dηe η/y − (cid:90) η ξ dξ G ξ ) η = T
24 ¯ hc τ . (B18)(iii) For y (cid:29) G , we expand the Bose function to give y/η and retain ξη in the first line of Eq. (B16) and ( ξη ) +2 ( G η ) in the second. The integrals converge at ξ, η ∼ G so the oscillations are fast when x (cid:28) G : J tTE = ¯ h G π c τ (cid:90) ∞ yη dη (cid:90) η ξ dξ ( ξ + 2 G ) η = √ G T πc τ . (B19)For the contributions coming from ξ (cid:28) x , we expand theexponential e iξ/x : J tTE = ¯ h G π c τ (cid:90) ∞ η dηe η/y − η (1 + G /x ) × (cid:90) η ξ dξ ( ξ + 2 G η ) + ξ η . (B20) Since ξ < η and G (cid:29) ξ in the firstbracket of the denominator in the last line. This allowsthe simple integration over ξ : J tTE = ¯ h G π c τ (cid:90) ∞ η dηe η/y − η / (2 G ) ]1 + η (1 + G /x ) , (B21)the condition for the expansion of the exponential e iξ/x becoming η (cid:28) x .There are three possible cutoff scales for η : y , G , and(1 + G /x ) − . Which one of the three is effective, dependson the positioning of y with repect to other scales. Again,three cases arise.(iv) For y (cid:28) (1 + G /x ) − <
1, the logarithm is expandedfor small argument and the second term in the denomi-nator is neglected since the integral converges at η ∼ y .The condition y (cid:28) x for the expansion of e iξ/x is satisfiedautomatically: J tTE = ¯ h π c τ (cid:90) ∞ η dηe η/y − π T h c . (B22)(v) For (1 + G /x ) − (cid:28) y (cid:28) G , the integral is still deter-mined by η ∼ y , but the second term in the denominatordominates. e iξ/x may be expanded when y (cid:28) x : J tTE = ¯ h π c τ x ( x + G ) (cid:90) ∞ η dηe η/y − T / ¯ h
48 ( cτ + G d ) . (B23)(vi) For y (cid:29) G , the Bose function is y/η and the integralconverges at η ∼ G so we retain the logarithm, and e iξ/x may be expanded as long as G (cid:28) x : J tTE = G T τ π c τ x ( x + G ) ∞ (cid:90) dηη ln (cid:18) η G (cid:19) = G T πc τ . (B24)We schematically show the regions of validity ofEqs. (B17)–(B24) for G (cid:29) x, y ) plane in Fig. 3,where there is no such overlap as in the TM case Fig. 2(a).As in the TM case, the first Fabry-Perot mode con-tribution should be included together with the ξ (cid:28) x contributions as long as x, y (cid:28) G . The same procedureis performed whereby the minimum of D s + [Eq. (B14b)]near ξ = πx is found, allowing the integrand to be ap-proximated by a Lorentzian. The minimum and thereforethe eventual contribution is found to be identical to theTM case, Eq. (B13).
3. TM evanescent contribution
The TM evanescent contribution to Eq. (A6) can berewritten exactly as J eTM = ¯ h G π c τ (cid:90) ∞ η dηe η/y − (cid:90) ∞ e − ξ/x ξ dξ ˜ D p + ˜ D p − , (B25a)˜ D p ± ≡ | η (1 − iη ) + i G ξ (1 ± e − ξ/x ) | . (B25b)1 FIG. 3. Regions of validity of Eqs. (B17)–(B24) in the ( x, y )plane.
This integral turns out to be exactly identical to thatalready calculated in Ref. [18] when the spatial dispersionof the conductivity is neglected (namely, Eqs. (1), (10)and (11) of Ref. [18]). That is to say that in the presentsystem the Coulomb limit ( c → ∞ ) amounts to takingonly the exact TM evanescent contributions to the heatcurrent, while neglecting the rest. Thus, we can simplyrewrite the results of Ref. [18] in terms of G : J ld = ζ (3) T π ¯ h c G d , (B26a) J hd = c G T πd (B26b) J lp = ζ (3) T π ¯ h c G d , (B26c) J hp = T πτ d L ( G cτ /d ) , (B26d)where L ( x ) is a slow logarithmic function approximatelygiven by: L ( x ) ≈ x x ) / ln(1 + ln x ) . (B27)The domains of validity of the contributions are shown inFig. 4. Note that expression (B26a) equals the travellingcontribution (B9).
4. TE evanescent contribution
The TE evanescent contribution to Eq. (A6) can berewritten exactly as J eTE = ¯ h G π c τ (cid:90) ∞ η dηe η/y − (cid:90) ∞ e − ξ/x ξ dξ ˜ D s + ˜ D s − , (B28a)˜ D s ± ≡ | iξ (1 − iη ) + G η (1 ± e − ξ/x ) | . (B28b) FIG. 4. Regions of validity of Eqs. (B26a)–(B26d) in the ( x, y )plane.
Despite the apparent similarity to the corresponding TEtravelling contribution Eq. (B14a), there is no longer os-cillatory behaviour in the denominator, so the resultingcontributions are completely different. In η there are twopossible decay scales: η ∼ y from the Bose function, and η ∼ ξ/ ( G + ξ ) from ˜ D s + ˜ D s − .In the low temperature case y (cid:28) ξ/ ( G + ξ ) < e − ξ/x ≈ ξ → ∞ . The large ξ cutoff scale is therefore given by the decay scale of theexponential, ξ ∼ x , leading to the result [valid for y (cid:28) x/ ( G + x )]: J eTE = ¯ h G π c τ ∞ (cid:90) η dηe η/y − x G η = π G T h c ln ¯ hc G T d . (B29)For high temperatures y (cid:29) ξ/ ( G + ξ ) the Bose functionis y/η and it is convenient to perform integration over η first keeping ˜ D s ± exact: J eTE = G T πc τ ∞ (cid:90) e − ξ/x ξ dξ ( G + ξ )[( G + ξ ) − G e − ξ/x )] , (B30)where the integrand may decay due to the exponential orthe denominator. If x (cid:28) G the exponential is clearly ac-tive and terms in ξ may be neglected in the denominator(the expansion of the Bose function is valid for y (cid:29) x/ G ): J eTE = T π G c τ (cid:90) ∞ ξ dξe ξ/x − ζ (3) cT π G d . (B31)If x (cid:29) G , expansion of e − ξ/x ≈ ξ → ∞ . As in the lowtemperature case, the divergence is cut off by ξ ∼ x (the2 FIG. 5. Regions of validity of Eqs. (B29)–(B32) in the ( x, y )plane. expansion of the Bose function is valid for y (cid:29) J eTE = G T πc τ (cid:90) ∼ x ξ dξ ( ξ + G )( ξ + 2 G ) = G T πc τ ln cτ G d . (B32)The domains of validity of the TE evanescent contribu-tions are shown in Fig. 5. Appendix C: Heat current betweenthree-dimensional metallic half-spaces
In this section we give a derivation of asymptoticexpressions for the heat current between two three-dimensional semi-infinite metallic half-spaces, separatedby a vacuum gap d , essentially reproducing the results ob-tained in Ref. [2]. We take two identical metals, describedby the complex dielectric functions ε ( ω ) = 1 + 4 πiσ /ω ,where σ is the three-dimensional dc conductivity,which can be written in terms of the bulk plasma fre-quency ω p and the electron relaxation time τ as 4 πσ = ω p τ , and assumed to be temperature-independent. Forconventional metals, 4 πσ (cid:29) ω p (cid:29) /τ , and it is nat-ural to assume T (cid:28) ¯ h/τ (indeed, τ = 10 − s corre-sponds to 760 K), so that for all relevant frequencies ε ( ω ) ≈ πiσ D /ω (cid:29)
1. Focusing on the local responseregime, we assume to be in the normal skin effect regime,characterised by the frequency-dependent skin depth δ ω and its value at ω = T / ¯ h : δ ω = c √ πσ D ω , δ T ≡ c (cid:112) πσ D T / ¯ h (C1)Since the metals are semi-infinite there can be no trans-mitted radiation and therefore the Joule losses are equalunambiguously to the average Poynting vector in the gap.The heat current per unit area J ( T ) may once again be written in the form of Eq. (A6), but without the thirdterm in the emissivities in Eq. (A8) (corresponding totransmission in the two-dimensional case), and with thereflectivities being just the Fresnel coefficients [36]: r p = q z − q (cid:48) z /εq z + q (cid:48) z /ε , r s = q z − q (cid:48) z q z + q (cid:48) z , (C2)where q (cid:48) z = (cid:112) [ ε ( ω ) − ω /c ) + q z is the normal com-ponent of the complex wavevector describing the electricand magnetic fields inside the metal, while q z is the samein the vacuum gap. As in the two-dimensional case, thecontributions from travelling and evanescent waves foreach polarisation are computed separately.
1. TM travelling contribution
The contribution may be written exactly as J tTM = ¯ h π (cid:90) ∞ ω dωe ¯ hω/T − (cid:90) ω/c q z dq z (cid:16) − | r p | (cid:17) (cid:12)(cid:12) − r p e − iq z d (cid:12)(cid:12) . (C3)The Fresnel coefficient is simplified drastically by notic-ing that since cq z < ω (cid:28) πσ we may write q (cid:48) z ≈ ωc (cid:112) ε ( ω ) , r p ≈ − ωcq z (cid:112) ε ( ω ) . (C4)Focussing firstly on the case where the exponential inthe denominator is oscillating fast, we may perform thesame averaging according to Eq. (B3), valid for q z (cid:29) /d ,which translates into d (cid:29) ¯ λ T . This gives J tTM = ¯ h π c (cid:90) ∞ ¯ hω dωe ¯ hω/T − (cid:114) ω πσ = 105 ζ (9 / π / ¯ h ( T / ¯ h ) / √ πσ c . (C5)When q z (cid:28) /d , the exponential in the denominator ofEq. (C3) may be expanded as 1 − iq z d , so the denomi-nator is approximately1 − r p e − iq z d ≈ iq z (cid:20) d − (1 + i ) δ ω ω c q z (cid:21) . (C6)This results in two expressions, depending on the relationbetween d and the thermal skin depth δ T : J tTM = ¯ h π (cid:90) ∞ ω dωe ¯ hω/T − (cid:90) ω/c δ ω q z dq z [( cq z /ω ) d − δ ω ] + δ ω = 15 ζ (7 / √ π ¯ h ( T / ¯ h ) / √ πσ D cd , d (cid:29) δ T , (C7a)= π ¯ h ( T / ¯ h ) c , d (cid:28) δ T . (C7b)Note that in the first case the q z integral converges at q z ∼ ( ω/c ) (cid:112) δ ω /d (cid:28) ω/c , so the expansion of e − iq z d is3valid at d (cid:28) ¯ λ T /δ T , which is a weaker condition than d (cid:28) δ T ; this means that the small q z contribution maycoexist with that of Fabry-Perot modes, but it is sub-dominant. In the second case d (cid:28) δ T , the convergenceis at q z ∼ ω/c , so the condition q z d (cid:28) d (cid:28) δ T .
2. TE travelling contribution
The situation is quite analogous to the TM case. TheTE contribution is given by Eq. (C3), but with the re-placement r p → r s . Instead of Eq. (C4), we have q (cid:48) z ≈ ωc (cid:112) ε ( ω ) , r s ≈ − − i ) q z δ ω . (C8)In the case of fast oscillation at d (cid:29) ¯ λ T , the denom-inator is again averaged using Eq. (B3), leading to anexpression, smaller than Eq. (C5) by a factor of 3.When d (cid:28) ¯ λ T , expanding e − iq z d ≈ − iq z d , weobtain 1 − r s e − iq z d ≈ iq z [ d − (1 + i ) δ ω ], which againresults in two expressions: J tTE = ¯ h π c (cid:90) ∞ ω dωe ¯ hω/T − δ ω ( ω ) | d − (1 + i ) δ ω | = ζ (3)4 π ¯ h ( T / ¯ h ) πσ D d , d (cid:29) δ T , (C9a)= π ¯ h ( T / ¯ h ) c , d (cid:28) δ T . (C9b)In contrast to the previous TM case, the integral is alwaysdominated by q z ∼ ω/c . In this case, analogously toEq. (B13), one can also take into account the contributionof the first Fabry-Perot mode: J tTE = π cT δ d e π ¯ hc/ ( T d ) , (C10)where δ is δ ω corresponding to ω = πc/d . This expres-sion has an exponential smallness, but its prefactor isparametrically larger than Eq. (C9a).
3. TM evanescent contribution
The contribution may be written exactly as J eTM = ¯ hπ (cid:90) ∞ ω dωe ¯ hω/T − (cid:90) ∞ ˜ q z d ˜ q z (Im r p ) e − q z d (cid:12)(cid:12) − r p e − q z d (cid:12)(cid:12) , (C11)where the real integration variable ˜ q z is introducedsince q z = i ˜ q z is purely imaginary. Then q (cid:48) z = (cid:112) ε ( ω ) ( ω/c ) − ˜ q z .Let us first consider the case where the ˜ q z dominatesover ε ( ω/c ) in the square root, that is ˜ q z δ ω (cid:29)
1. Then r p ≈ − ε ( ω ) = 1 + iω πσ D , (C12) and J eTM = ¯ h π (cid:90) ∞ ω dωe ¯ hω/T − (cid:90) ∞ [ ω/ (2 πσ D )] ˜ q z d ˜ q z sinh ˜ q z d + [ ω/ ( πσ D )] = π
60 ¯ h ( T / ¯ h ) (2 πσ D ) d ln min (cid:26) πσ D T / ¯ h , δ T d (cid:27) . (C13)The ˜ q z integral is logarithmic, and is determined bya broad interval of ˜ q z from the upper cutoff ∼ /d down to the lower cutoff: for d (cid:28) c/ (2 πσ D ) it is(1 /d ) (cid:112) ω/ (2 πσ D ), while at larger distances the small ˜ q z cutoff is determined by the condition of ˜ q z δ ω (cid:29)
1. Thelogarithmic region exists at all if ε ( ω/c ) can be neglectedat ˜ q z ∼ /d , which translates into d (cid:28) δ T .In the opposite case, where we neglect ˜ q z (cid:28) /δ ω in q (cid:48) z , we still assume ˜ q z (cid:29) ( ω/c ) / √ ε ∼ ( ω/c ) δ ω , so thereflection coefficient is still close to unity: r p ≈ i ) ωc ˜ q z (cid:114) ω πσ . (C14)Then the ˜ q z integral is determined by small ˜ q z ∼ (cid:112) ω/ ( cd )[ ω/ (2 πσ D )] / (cid:28) /d , so sinh ˜ q z d ≈ ˜ q z d : J eTM = ¯ h π (cid:90) ∞ ω dωe ¯ hω/T − × (cid:90) ∞ ˜ q z d ˜ q z (cid:32) ˜ q z cdω (cid:114) πσ D ω − (cid:33) + 1 − = 45 ζ (7 / √ π ¯ h ( T / ¯ h ) / √ πσ cd . (C15)The conditions ˜ q z d (cid:28)
1, ( ω/c ) δ ω (cid:28) ˜ q z (cid:28) /δ ω result inthe requirement T / ¯ h πσ D δ T (cid:28) d (cid:28) πσ D T / ¯ h δ T . (C16)Note that the lower limit on d is smaller than δ T ,so there is an interval where Eqs. (C13) and (C15)are both valid, representing contributions from differ-ent regions of ˜ q z integration. However, when inequali-ties (C16) hold, Eq. (C15) automatically dominates overEq. (C13). Going to longer distances, where the as-sumption ( ω/c ) / √ ε (cid:28) ˜ q z is violated, is not necessarysince at such distances (well exceeding ¯ λ T ) the travellingwave contributions dominate; indeed, Eq. (C5) exceedsEq. (C15) in the common wisdom region d (cid:29) ¯ λ T .
4. TE evanescent contribution
The TE contribution is given by Eq. (C11), but withthe replacement r p → r s . If we try to proceed as in theTM case and assume first ˜ q z δ ω (cid:29)
1, we obtain an integraldiverging at small ˜ q z , invalidating the assumption.4Making the opposite assumption, q z δ ω (cid:28)
1, for thereflection coefficient we obtain the same approxima-tion (C8) with q z = i ˜ q z , which leads to J eTE = ¯ h π ∞ (cid:90) ω dωe ¯ hω/T − ∞ (cid:90) ˜ q z d ˜ q z [(sinh ˜ q z d ) /δ ω + ˜ q z ] + ˜ q z . (C17)At d (cid:29) δ ω the ˜ q z integral converges at ˜ q z ∼ /d , and theresulting logarithmic ω integral J eTE = 3 ζ (3)8 π d ∞ (cid:90) δ ω ¯ hω dωe ¯ hω/T − ζ (3)4 π c T πσ D d ln dδ T , (C18) is cut off at low frequencies by the condition d ∼ δ ω , sothat the validity condition is d (cid:29) δ T .For d (cid:28) δ T we are forced to consider ˜ q z δ ω ∼ q (cid:48) z = (cid:112) i − ˜ q z δ ω /δ ω ; however, wecan safely set d → ± r s ∼
1. Then we obtain(Im r s ) | − r s | = (Re q (cid:48) z ) | q (cid:48) z | = 1 /
24 + ˜ q z δ ω + ˜ q z δ ω (cid:112) q z δ ω , (C19)which gives J eTE = 2 πσ D ¯ h π c (cid:90) ∞ ω dωe ¯ hω/T − ζ (3)4 π πσ ¯ h ( T / ¯ h ) c , (C20)the ˜ q z integral converging at ˜ q z ∼ /δ ω , as expected. [1] S. M. Rytov, Theory of electric fluctuations and thermalradiation (Air Force Cambrige Research Center, Bedford,MA, 1953).[2] D. Polder and M. Van Hove, Phys. Rev. B , 3303 (1971).[3] M. L. Levin, V. G. Polevoi, and S. M. Rytov, Sov. Phys.JETP , 1054 (1980).[4] J. J. Loomis and H. J. Maris, Phys. Rev. B , 18517(1994).[5] J. B. Pendry, Journal of Physics: Condensed Matter ,6621 (1999).[6] S. M. Rytov, Y. A. Kravtsov, and V. I. Tatarskii, Prin-ciples of statistical radiophysics (Springer-Verlag, BerlinHeidelberg, 1989).[7] K. Joulain, J.-P. Mulet, F. Marquier, R. Carminati, andJ.-J. Greffet, Surface Science Reports , 59 (2005).[8] A. I. Volokitin and B. N. J. Persson, Rev. Mod. Phys. , 1291 (2007).[9] B. Song, A. Fiorino, E. Meyhofer, and P. Reddy, AIPAdvances , 053503 (2015).[10] S.-A. Biehs, R. Messina, P. S. Venkataram, A. W. Ro-driguez, J. C. Cuevas, and B. Ben-Abdallah, “Near-fieldradiative heat transfer in many-body systems,” (2020),arXiv:2007.05604.[11] P.-O. Chapuis, S. Volz, C. Henkel, K. Joulain, and J.-J.Greffet, Phys. Rev. B , 035431 (2008).[12] P.-O. Chapuis, M. Laroche, S. Volz, and J.-J. Greffet,Phys. Rev. B , 125402 (2008).[13] M. Prunnila and S. J. Laakso, New Journal of Physics , 033043 (2013).[14] G. D. Mahan, Phys. Rev. B , 115427 (2017).[15] Z.-Q. Zhang, J.-T. L¨u, and J.-S. Wang, Phys. Rev. B , 195450 (2018).[16] J.-S. Wang, Z.-Q. Zhang, and J.-T. L¨u, Phys. Rev. E , 012118 (2018).[17] A. Kamenev, “Near-field heat transfer between disor-dered conductors,” (2018), arXiv:1811.10187.[18] J. L. Wise, D. M. Basko, and F. W. J. Hekking, Phys.Rev. B , 205411 (2020).[19] X. Ying and A. Kamenev, Phys. Rev. B , 195426(2020).[20] A. O. Govorov and A. V. Chaplik, Sov. Phys. JETP (1989). [21] V. I. Fal’ko and D. I. Khmel’nitskii, Sov. Phys. JETP (1989).[22] V. A. Volkov and V. N. Pavlov, JETP Letters , 93(2014).[23] V. M. Muravev, P. A. Gusikhin, I. V. Andreev, and I. V.Kukushkin, Phys. Rev. Lett. , 106805 (2015).[24] P. A. Gusikhin, V. M. Muravev, A. A. Zagitova, andI. V. Kukushkin, Phys. Rev. Lett. , 176804 (2018).[25] D. O. Oriekhov and L. S. Levitov, Phys. Rev. B ,245136 (2020).[26] M. A. Ordal, R. J. Bell, R. W. Alexander, L. L. Long,and M. R. Querry, Appl. Opt. , 4493 (1985).[27] L. Wang, M. Bie, W. Cai, L. Ge, Z. Ji, Y. Jia, K. Gong,X. Zhang, J. Wang, and J. Xu, Opt. Express , 36790(2019).[28] B. Altshuler and A. Aronov, in Electron–Electron Inter-actions in Disordered Systems , Modern Problems in Con-densed Matter Sciences, Vol. 10, edited by A. Efros andM. Pollak (Elsevier, 1985) pp. 1 – 153.[29] P. A. Lee and T. V. Ramakrishnan, Rev. Mod. Phys. ,287 (1985).[30] J. Yang, W. Du, Y. Su, Y. Fu, S. Gong, S. He, andY. Ma, Nature Communications , 4033 (2018).[31] C. M. Hargreaves, Physics Letters A , 491 (1969),more precise measurements were described in the Ph. D.Thesis of C. M. Hargreaves (University of Leiden, 1973),reproduced in Ref. [9].[32] B. Song, D. Thompson, A. Fiorino, Y. Ganjeh, P. Reddy,and E. Meyhofer, Nature Nanotechnology , 509 (2016).[33] T. Kralik, P. Hanzelka, M. Zobac, V. Musilova, T. Fort,and M. Horak, Phys. Rev. Lett. , 224302 (2012).[34] P. Sabbaghi, L. Long, X. Ying, L. Lambert, S. Taylor,C. Messner, and L. Wang, Journal of Applied Physics , 025305 (2020).[35] M. F. Modest, Radiative Heat Transfer (Academic Press,San Diego, 2013).[36] A. I. Volokitin and B. N. J. Persson, Phys. Rev. B ,205404 (2001).[37] C. J. Fu and Z. M. Zhang, International Journal of Heatand Mass Transfer49