Near field versus far field in radiative heat transfer between two-dimensional metals

Using the standard fluctuational electrodynamics framework, we analytically calculate the radiative heat current between two thin metallic layers, separated by a vacuum gap. We analyse different contributions to the heat current (travelling or evanescent waves, transverse electric or magnetic polarization) 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 when the separation between the layers greatly exceeds the thermal photon wavelength, and the coupling is of electrostatic nature. For well-conducting metals, the evanescent contribution dominates at separations smaller than the thermal wavelength and is mainly due to magnetostatic coupling, in agreement with earlier works on bulk metals.


I. INTRODUCTION
Spatially separated objects may exchange heat via electromagnetic fluctuations [1][2][3][4][5]. This radiative heat transfer arises due to electric charge density and current fluctuations inside the constituting materials, and is usually described within the phenomenological framework of fluctuational electrodynamics (FED) [1,2,6], for which the critical inputs are the material response functions and the system geometry. It is now well known that in the near-field limit, energy may tunnel via evanescent electromagnetic waves causing a strong enhancement of the heat transfer, as has been observed experimentally (see the reviews [7][8][9][10] and references therein).
Many theoretical works have been dedicated to different material systems in the near-field regime ( [7][8][9][10] and references therein), in which various models for material response have been employed and different dominant channels for heat transfer identified. The common wisdom is that the evanescent modes dominate the heat transfer when the spatial separation d λ T ≡hc/T , the wavelength of photons at temperature T (hereh and c are the Planck constant and the speed of light, respectively, and we set the Boltzmann constant to unity). Indeed, for d >λ T the evanescent waves with the typical frequency ω ∼ T /h decay exponentially outside the material, while at d λ T the region of the wave vectors k occupied by evanescent waves, k ∼ 1/d, is larger than that of travelling states, k ∼ 1/λ T [8]. Importance of magnetic coupling in the near-field heat transfer between well-conducting metals has been emphasised [11,12]. In the extreme near-field limit, heat transfer due to the electrostatic Coulomb interaction has also been studied [13][14][15][16][17][18][19].
Here we revisit this old problem, focusing on the two-dimensional (2D) geometry, and study the radiative heat current between two thin metallic sheets in vacuum within the standard FED framework. We find two qualitatively different types of behaviour, depending on the value of the two-dimensional dc conductivity σ 2D of the * jonathan.wise@lpmmc.cnrs.fr sheets. For poor conductors characterised by the condition G ≡ 2πσ 2D /c 1 (we use CGS units throughout the paper, in SI units G = (σ 2D /2) µ 0 /ε 0 ), the heat transfer turns out to be dominated by the evanescent modes at distances d extending well beyondλ T , and the main coupling mechanism in the near field is electrostatic (Coulomb interaction between electrons in the two layers). For G 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 of parameters the near-field transfer is dominated by magnetostatic (inductive) coupling between currents in the layers.
The parameter G characterises the impedance mismatch between a two-dimensional (2D) metal and vacuum; its importance is not restricted to the heat transfer problem and is rather general. Notably, two distinct regimes in the behaviour of 2D plasmon polaritons for G < 1 and G > 1 have been identified [20][21][22][23][24][25]. In our heat transfer problem, we find no sharp distinction between G > 1 and G < 1, but rather a smooth crossover between the two limiting situations.
The rest of the paper is organised as follows. In Sec. II we specify the model and sketch the calculation; both are rather standard. In Sec. III we present various regimes of the heat transfer and the associated analytical expressions for the heat current, according to the material properties and experimental conditions. In Sec. IV we discuss the relation of our results to the well-studied case of heat transfer between bulk semi-infinite metals, and compare our theory to available experimental results. All details of calculations are given in three appendices.

II. THE MODEL
We consider two identical 2D metal sheets held at different temperatures T 1 and T 2 , embedded in vacuum and separated by a gap of width d. A more realistic configuration would be to place a medium with a dielectric constant ε in the half-space behind each sheet, since in experiments the layers are placed on a substrate. For the sake of simplicity, we focus on ε = 1 in most of the paper, arXiv:2101.05515v1 [cond-mat.mes-hall] 14 Jan 2021 and check for the effect of the substrate when specifically needed.
We model the metal sheets as infinitely thin layers, characterised by a local 2D Drude conductivity, with τ being the electron momentum relaxation time, assumed 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 conductivity, and (ii) field variation over the layer thickness. For atomically thin materials, such as doped graphene or transition metal dichalcogenides, condition (ii) is irrelevant, and condition (i) holds at distances d √ a 2D (a 2D and being the 2D screening radius and the electron mean free path, respectively) [18]. For thin but macroscopic layers of conventional metals, condition (ii) imposes that the thickness must be small compared both to the typical wavelength of the waves dominating the heat transfer (which may be rather short for evanescent waves) and to the skin depth at the typical frequency of these waves, while condition (i) requires the wavelength and the skin depth to be longer than the electron mean free path in the metal.
Our calculation of the heat current between the metals follows the standard FED procedure. The fluctuating inplane surface currents j (α) (r, t) in each sheet obey the fluctuation-dissipation theorem, where k is the in-plane two-dimensional wavevector and l, m = x, y label the orthogonal in-plane directions and T = T 1 or T 2 . These currents appear as sources in Maxwell's equations, whose solution in the presence of the conducting sheets determines the fluctuating electric fields E(r, t). Then, the heat current J (per unit area) from layer 1 to layer 2 is given by the average Joule loss power (per unit area) j(2) · E (2) − j(1) · E (1) , wherej (α) is the surface current in layer α, induced by the electric field E (α) in this layer, which, in turn, is produced by the fluctuating current in the other layer (see Appendix A for explicit expressions, rather standard). We emphasize that for thin layers, the Joule losses j(2) · E (2) − j(1) · E (1) are not equal to the average normal component of the Poynting vector in the gap between the layers. The reason is that some part of the radiation emitted by layer 1 may pass through layer 2 and escape to infinity, and vice versa. Whether this escaped radiation should be included in the heat current or not, depends on the precise measurement setup, which may collect this escaped radiation or not. Our calculation thus assumes that the escaped radiation is lost. Note that for two semi-infinite metals (the most studied setup), everything is collected inside the metals, so the Poynting vector and the Joule losses match exactly.
In the planar geometry considered here, the solutions of Maxwell's equations are classified by in-plane wave vector k, frequency ω, and two polarisations p = TE, TM -transverse electric and transverse magnetic, respectively, for which the electric or the magnetic field vector is parallel to the layers and perpendicular to k. The contributions to the heat current from modes with different k, ω, p add up independently, so the heat current J(T 1 , T 2 ) is given by by an integral over k and ω, and a sum over the polarisations. The integral splits in two contributions: the interior of the light cone, ω > ck hosts travelling modes, while in the region ω < ck the solutions are evanescent. The resulting heat current is comprised of four additive contributions (TM and TE, travelling and evanescent). Which contribution dominates, depends on the material conductivity, as well as the system temperature and length scales.
In the extreme near field limit, k ω/c, the TM mode field is mostly electric and longitudinal, while the magnetic field is smaller by a factor ∼ ω/(ck); these modes represent the electrostatic coupling by the Coulomb interaction between charge density fluctuations in the two layers. At the same time, for TE modes the field is mostly magnetic, while the electric field is smaller by a factor ∼ ω/(ck); these modes represent magnetostatic coupling, where the magnetic field established by transverse current fluctuations in one layer drives eddy currents in the second layer.
In Appendix B we perform analytically the k, ω integrals and derive simple asymptotic expressions for the heat current according to the separation and the temperature. For each expression, we can identify the dominant contribution (TM or TE, travelling or evanescent). Our results are approximate; one can describe the heat transfer much more precisely by solving Maxwell's equations for finite-thickness slabs with a material-specific frequency dependence of the conductivity and numerically evaluating the integrals, as routinely done in many works. However, simple approximate expressions (i) are rather useful when a quick estimate of the heat current is needed, and (ii) offer a general insight into the dominant physical mechanisms responsible for the heat transfer and enable one to characterise different possibilities.

III. RESULTS
For temperature-independent relaxation time, the heat current naturally splits into the difference J(T 1 , T 2 ) = J(T 1 ) − J(T 2 ). The detailed analysis of different asymptotic regimes of the k, ω integrals results in several asymptotic expressions for J(T ) in different parametric ranges of d. The magnitude of the TM and TE travelling contributions are sensitive to the dimensionless conductivity parameter G = 2πσ 2D /c so we proceed to present the results sequentially for small and large values of G.
J it (T ) = 1 12 valid in the corresponding regions of the (1/d, T ) plane, schematically shown in Fig. 1 (right). The contributions given by Eqs. (5a)-(5c) are the TE evanescent contributions, while Eqs. (5d)-(5f) are the sums of travelling contributions from both polarisations which are of the same order. For G 1, the travelling channels support resonant Fabry-Perot (FP) modes. In (lt) and (it) regions, many sharp FP modes contribute significantly to the heat current. In the high-temperature case (ht) the FP modes are overdamped since the conductivity σ(ω) becomes small at high frequencies. For temperatures lower than the first mode cutoff energy, T πhc/d, the contributions from the FP modes are exponentially suppressed. However, the prefactor in front of the small thermal exponen-tial turns out to be larger than the evanescent contribution (5c) in (he2) region. Thus, the FP additive contribution is potentially significant for cτ /d < √ G, where it is dominated by the first FP mode: In Fig. 1, the areas with wavy hatching indicate the regions where the heat transfer is dominated by travelling wave contributions. For G 1, the evanescent waves dominate at separations up to d ∼ G −1/3λ T , parametrically larger than the commonly used condition for the near field, d λ T [y x in Fig. 1 (left)]; the reason for such behaviour is that the low-temperature TM evanescent contribution is determined by k 1/d for which the exponential suppression is not efficient. Moreover, for G 1 the near-field transfer is dominated by the TM evanescent contribution, basically, by electrostatic (capacitive) coupling between the two layers. This happens because in a poor conductor, the charge density response is not fast enough to dynamically screen the fluctuating Coulomb field.
For G 1, the commonly used inequality d λ T does become the accurate condition for evanescent contribution dominance, except for high temperatures where the Drude conductivity is suppressed by high frequency. In a large part of the near-field region of the parameter plane the heat current is governed by TE evanescent modes, which correspond to magnetostatic (inductive) coupling between the layers. As discussed in Refs. [11,12] for bulk metals, large conductivity leads to efficient screening of the electric fields, so the magnetostatic coupling becomes more important. The electrostatic coupling takes over only at very short distances or low temperatures, However, for very small d the finite layer thickness may become important, and/or the assumption of the local response, Eq. (1), may break down. Taking, for example, a 10 nm-thick gold film with the bulk plasma frequency ω p = 0.6 × 10 16 s −1 and relaxation time τ = 6 fs [26] gives G ≈ 3.6. Sinceλ T = 7.6 µm at T = 300 K, in such structure the crossover to electrostatics occurs at a few nanometers. We note that in the (ld) regime, the heat transfer is mainly determined by rather small wave vectors k ∼ (Gλ T d) −1/2 [18], so that even at d = 1 nm we obtain 1/k > ∼ 100 nm, and the local response assumption should still be formally valid. However, at nanometric distances other physical effects may come into play (electron or phonon tunnelling, surface roughness, etc.), so for ultrathin films of conventional metal we expect the Coulomb mechanism to be relevant mostly at low temperatures.

A. Comparison to the bulk case
The results presented in the previous section show two qualitatively different pictures of the near-field heat transfer between two metallic layers, depending on the value of their dimensionless 2D dc conductivity: for G 1, the heat transfer is mostly due to electrostatic coupling between the layers, up to distances significanly exceedingλ T , while for G 1 the near-field magnetostatic coupling dominates up to distances d ∼λ T , in close analogy with earlier results on bulk metals. This picture is consistent with the results of Ref. [27] where the very same problem of radiative heat transfer between parallel 2D layers was studied numerically. There, a distinction was made between thin and thick metallic films. In this formulation G is proportional to the layer thickness h (in the local approximation, the 2D conductivity is simply σ 2D = σ 3D h, where σ 3D is the bulk conductivity). In Ref. [27], the heat transfer between two theoretically imagined atomic monolayers of silver, described by a 2D Drude model with G ≈ 2, is found to be driven by TM evanescent waves, while for thicker films it is TE evanescent waves.
The peculiarity of the 2D geometry is that the 2D conductivity can be compared to two universal scales. One is the speed of light, hence the dimensionless parameter G = 2πσ 2D /c we introduced earlier. The other universal scale is the conductance quantum, e 2 /(2πh). For σ 2D < ∼ e 2 /(2πh), or G < ∼ e 2 /(hc) ≈ 1/137, the disorder is too strong, so the metallic conduction is destroyed by localization effects [28,29]. Thus, the poor conductor regime discussed above, can be realised in the interval 1/137 < ∼ G 1. The situation is quite different for the bulk metal case. The 3D conductivity σ 3D has the dimensionality of the inverse time, so that 1/(4πσ 3D ) (in CGS units, while in SI it is ε 0 /σ 3D ) has a meaning of the RC time needed to dissolve a charge density perturbation. In conventional metals this time scale is extremely short (in the attosecond range). Still, one can compare 4πσ 3D to other scales. One is the electron relaxation time τ ; typically, 4πσ 3D τ = ω 2 p τ 2 1 (ω p being the bulk plasma frequency). Moreover, at T h/τ the relaxation time drops out of the problem, so one cannot construct a dimensionless parameter out of σ 3D , which could produce different "asymptotic maps" of the kind shown in Fig. 1. The bulk case turns out to be somewhat similar to the 2D case with G 1. To see the reason for this similarity, let us recall the asymptotic expressions for the heat current between semi-infinite bulk metals, assuming T τ h (the derivation can be found in Ref. [2], we also give it in Appendix C): where the parametric intervals of d are conveniently defined in terms of two length scales: the thermal wavelengthλ T =hc/T and the normal skin depth at the thermal frequency, δ T = c/ 2πσ 3D T /h λ T . The shortestdistance expression (7a) is determined by the TM evanescent contribution and corresponds to the Coulomb limit (indeed, it does not contain the speed of light); however, the length scale δ 3 T /λ 2 T is extremely short: for 4πσ 3D = 10 17 s −1 at T = 300 K, we have δ T = 0.22 µm andλ T = 7.6 µm, so δ 3 T /λ 2 T ∼ 2Å and becomes even smaller at lower temperatures, invalidating the local approximation and making Eq. (7a) irrelevant for conventional metals. Equations (7b) and (7c) originate from the TE evanescent contribution and correspond to magnetostatic coupling [11]. Equation (7d) contains both TM evanescent and TM travelling contributions which are of the same order at such distances (only the evanescent one was evaluated in Ref. [2]); in fact, for both contributions the integral is dominated by wave vectors k very close to ω/c, and the fields vary weakly across the gap so there is no sharp physical distinction between travelling and evanescent waves. Finally, Eq. (7e) comes from the TE and TM travelling waves and is contributed by many Fabry-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 from Eqs. (5a), (5c), (3d) and (5d), respectively, by replac- 3D ω corresponds to the typical frequency scale 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 logarithmic, with the lower cutoff corresponding to δ ω ∼ d (see Appendix C for details). This replacement roughly corresponds to modelling the semi-infinite metal as an effective metallic layer whose thickness corresponds to the field penetration depth. Such effective layer is characterised by the dimensionless G eff ∼ 2πσ 3D /ω, so that G eff 1 for conventional metals and reasonable temperatures. We note that this effective layer analogy should be used with caution, since the frequency dependence of δ ω sometimes makes the convergence scale of the frequency integral different from the case of fixed layer thickness.

B. Comparison to experiments
Values G < ∼ 1 are characteristic of atomically thin 2D materials. This is illustrated by a recent experiment [30], where two doped monolayer graphene sheets were placed on insulating silicon (ε = 11.7) and separated by a 400 nm wide vacuum gap. The Fermi energy 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 −2 K −1 was measured around room temperature. These conditions correspond to the hightemperature plasmon regime, Eq. (3b), where the substrate dielectric constant ε enters only inside the logarithmic function L [18]. Setting L = 1 in Eq. (3b) gives dJ/dT = 11 W m −2 K −1 , which agrees by order of magnitude with the experimental value.
Thin layers of conventional metals are typically characterised by G 1. Several experiments have been reported in the literature. In each case, it is important to compare the layer thickness h to the the skin depth δ ω at the relevant frequency, to ensure the layers should correspond to the 2D limit, rather than the bulk one (the latter being the case of Refs. [31,32]).
Heat transfer in a wide range of interlayer separations and temperatures was studied in Ref. [33] for two 150 nm thick tungsten layers on alumina substrates. The measured dc conductivity of the material 4πσ 3D = 0.6 × 10 18 s −1 (constant in the temperature range of the experiment) corresponds to a value of the dimensionless conductivity parameter G ≈ 150. The skin depth at T = 40K is δ T = 240 nm, and even longer at lower temperatures, so the layers are close to the 2D limit. The separation between the layers was varied over d = 1 − 300 µm, while the temperatures were T 1 = 5 K and T 2 = 10 − 40 K, corresponding to regions (he2) and (lt) in Fig. 1 (right). It can be easily checked that in these regions, the dielectric substrate plays no role as long as √ ε G, which clearly holds here. Although the numerical calculation accounting for the finite layer thickness does better in closely matching the experimental points (see Fig. 2 of Ref. [33]), our simple expressions (5c) and (5d) (i) agree with the observed values within a factor of 3 without any fitting parameters, (ii) give the correct distance dependence throughout the experiment, (iii) capture the observed approximate collapse of the rescaled data for J(T )/T 4 on a function of a single variable T d, and (iv) correctly predict the separation d ≈ 0.5λ T , at which the crossover between the near-field and the far-field regimes occurs, J he2 (T ) = J le (T ). A recent publication [34] presents measurements of heat transfer between two aluminium films of varying thicknesses h = 13 − 79 nm, separated by a fixed vacuum gap d = 215 nm and attached to silicon substrates. The experiment was performed around room temperature with one film being heated such that ∆T = 25−65 K.
An intriguing feature of the results reported in Ref. [34] is the independence of dJ/dT of the layer thickness. This agrees neither with our 2D expressions, nor with the more precise simulations done in Ref. [34]. All theoretical results point to a non-monotonic dependence of the heat current on the layer thickness or dc conductivity [the latter is also true for the bulk limit expressions (7)]. Further experimental investigations of this dependence would be interesting.

V. CONCLUSIONS
In this paper, we have performed an analytical calculation of the radiative heat current between two thin metallic layers, using the standard framework of fluctuational electrodynamics and a local 2D Drude model for the electromagnetic response of each layer. We have identified two different classes of such structures, distinguished by the dimensionless 2D dc conductivity G = 2πσ 2D /c. For poor conductors with G 1, typically represented by atomically thin 2D materials, the heat transfer is dominated by evanescent modes at distances d extending well beyondλ T , and the main coupling mechanism in this near-field regime is the Coulomb interaction between electrons in the two layers. Good conductors with G 1, such as thin films of conventional metals, behave more similarly to the bulk limit, studied in earlier works: the crossover from near to far field occurs at d ∼λ T at not too high temperatures, and the near-field transfer is dominated by magnetostatic (inductive) coupling between the layers in a wide range of parameters.
We have derived several simple approximate asymptotic expressions for the heat current valid in differ-ent parametric ranges of interlayer separation distance and temperature. Comparing these expressions with the available experimental data, we saw that they give valid order-of-magnitude estimates of the heat current and correctly capture its dependence on the distance and temperature. Better agrreement with the experimental results can be reached by a more detailed modelling of each system geometry and the dielectric response, which is strongly system-specific and lies beyond the scope of our work. Still, our approximate results offer a useful insight into the main physical mechanisms responsible for the heat transfer.

Appendix A: Explicit general expression for the heat current between two thin metallic sheets
We solve Maxwell's equations for the monochromatic components of the electric and magnetic field, E kω (z) e ikr−iωt and B kω (z) e ikr−iωt in the planar geometry with the two metallic sheets placed at z = z 1 , z 2 with z 2 − z 1 = d, while the position-dependent dielectric constant ε(z) accounts for whatever (non-magnetic, isotropic) dielectric medium surrounds the layers: where e z is the unit vector in the z direction, perpendicular to the layers. The surface current in each layer α = 1, 2 consists of two contributions:j is the induced current due to the electric field, while the fluctuating currents j Because of δ lm on the right-hand side of this equation, current fluctuations are independent for any two orthogonal directions, so it is convenient to pass to the longitu-dinal and transverse basis (p and s polarisations, respectively): In this basis the solutions of Maxwell's equations decouple into transverse magnetic (TM) and transverse electric (TE) modes, whose contribution to the heat current is simply additive.
To model different metal sheets mounted on identical dielectric substrates separated by vacuum, we take ε(z 1 < z < z 2 ) = 1, ε(z < z 1 ) = ε(z > z 2 ) = ε > 1. This leads to the spatial dependence of the electric and magnetic fields ∝ e ikr±iqzz for z 1 < z < z 2 , and ∝ e ikr+iq z z , e ikr−iq z z for z > z 2 and z < z 1 , respectively. Here we defined At |ω| > ck, the metallic layers are coupled by travelling waves, while for |ω| < ck the solutions in the gap are evanescent waves, where the fields' strength decays away from the layers. The solutions are matched at z = z 1 and z = z 2 using the standard boundary conditions: continuity of the in-plane component of the electric field E , and a jump in the magnetic field in-plane component, determined by the total surface current (the fluctuatinng sources as well as the induced current σE ). The heat current from, say, sheet 1 to the sheet 2 is given by the average Joule loss power per unit area, J(T 1 , T 2 ) = j(2) · E (z 2 ) − j(1) · E (z 1 ) , determined unambiguously due to the continuity of E (z). For a temperature independent relaxation time this heat current splits into J(T 1 , is expressed in terms of reflectivities r αj and emissivities a αj for the p and s polarisations: The emissivities can also be written as where θ(x) is the Heaviside step function, and t αj are the transmittivities: Note the difference between Eqs. (A8) and Eq. (2) of Ref. [27], where the third term is absent in both polarisations. Without the third term, Eq. (A6) gives the average value of the Poynting vector in the gap between the two layers, and also counts the heat flux which is not absorbed by the metal, but irradiated to infinity behind it, due to the finite transmission. Eqs. (A8) without the third term originally appeared in Ref. [36] for the problem of heat transfer between two semi-infinite materials.
In that geometry, all heat flux transmitted through the surface is eventually absorbed by the material. In the thin layer geometry, whether the transmitted flux is detected or not, depends on the specific experimental measurement scheme. In our calculation, we assume that the transmitted radiation is lost, and thus use the full Eqs. (A8).
Appendix B: Derivation of asymptotic expressions for the heat current between two thin metallic sheets Here we derive asymptotic expressions for J(T ) in the specific case of identical sheets embedded in vacuum [σ 1 (ω) = σ 2 (ω) and ε = 1] and compute separately the travelling and evanescent wave contributions for each of the two polarisations. We quantify the contribution made by each wave type and polarisation in each region of the (1/d, T ) parameter plane, before comparing the size of the additive contributions and identifying which are dominant. It is convenient to introduce the dimensionless parameters x ≡ cτ /d and y ≡ T τ /h, as well as dimensionless integration variables: ξ = |q z |cτ instead of k [noting that k dk = ξ dξ/ (cτ ) 2 ], and η = ωτ . For the travelling waves, the integration is over the region 0 < ξ < η < ∞, while for the evanescent waves it is 0 < ξ, η < ∞.

TM travelling contribution
In the dimensionless variables, the TM travelling contribution to Eq. (A6) can be rewritten exactly as The case G 1 is very simple to handle, since for ε = 1 one can neglect the reflection coefficients in the denominator of Eq. (A6), and simply set G → 0 in Eq. (B1b), since ξ < η. This gives For G 1, each layer behaves at low frequency as a well-reflecting mirror, so the structure may host Fabry-Perot modes. The Fabry-Perot modes manifest themselves as deep minima in D p ± at specific values of ξ/x = π, 2π, 3π, . . .. These minima are important when Gξ η 1 + η 2 , which is precisely the condition of good reflection. Thus, a much more elaborate analysis is needed to evaluate the integral.
Let us focus on the contributions from the region ξ x, when many modes contribute, and even if they are overdamped, e iξ/x oscillates fast. In the general case (A6) we average over the fast oscillations in the denominator which leads to the simple replacement [37]: valid as long as a αj and r αj are smooth functions of q z on the scale q z ∼ 1/d. Applying this averaging to the contribution in Eq. (B1a) leads to Note that x dropped out, and enters only through the condition ξ x. Note also that the ξ integral is always determined by the upper limit ξ ∼ η. As for the η integral, it may converge at η ∼ y when cut off by the Bose function, or, for too large y, it may be cut off by other factors in the denominator at some η y. In this latter case, one can expand the exponential in the Bose function, which becomes just y/η. We can identify three regions in y.
(i) For y √ G, the integrals separate and converge at ξ ∼ η ∼ y, so the fast oscillation condition is x y: (B5) (ii) For √ G y G, we keep 2(Gξ) 2 in the first bracket and η 4 in the second one (again, oscillations are fast when x y): (B6) (iii) For y G, we expand the Bose function, the integral converges at η ∼ G (it is convenient to write ξ = uη); the oscillations are fast when x G: Let us now pick the contributions from ξ x. Then, e iξ/x can be expanded (we again write ξ = uη): There are three possible cutoff scales for η: y, 1 + 2Gu, and (1+Gu 2 /x) −1 . Which one of the three is effective, depends on the positioning of y with repect to other scales. Again, three cases arise.
(iv) For y 1, we can neglect η 2 in the first bracket in the denominator, so the η integral converges at η ∼ y. In the second bracket, η 2 plays a role only if Gu 2 /x 1, so the second bracket can be approximated as 1+(Gu 2 η/x) 2 for any Gu 2 /x. We also assume that Gu 1, which will be verified afterwards. Then the denominator becomes 4G 2 u 2 [1 + (Gu 2 η/x) 2 ], so the u integral converges at u ∼ min{1, x/(Gy)}, giving The u integral converges at u ∼ 1 and u ∼ x/(Gy) in the two cases. In the first case, y x/G, the assumption Gu 1, as well as the condition to expand the exponential, uy/x 1, are satisfied automatically. In the second case, y x/G, both conditions translate into y Gx. (v) For y 1 but y Gu we still have η ∼ y, so the de-nominator can be approximated as 4Gu 2 η 2 (1 + Gu 2 /x) 2 : Since the convergence occurs at u ∼ min{1, x/G}, η ∼ y, the assumption y Gu is satisfied if y min{ √ Gx, G}; if so, the condition uy/x 1 to expand the exponential is satisfied automatically. Thus, Eq. (B10) is valid when 1 y min{ √ Gx, G}. (vi) For y 1, Gu, the Bose function is y/η, so we integrate over η exactly (convergence at η ∼ 1 + 2Gu), and obtain the convergence occurring at u ∼ min{ x/G, 1}. At x G the condition to expand e iξ/x is not fulfilled, since we automatically have ξ/x = uη/x ∼ 1. At x G, we have Gu 1 automatically, while uη/x ∼ G/x, so the second expression Eq. (B11) is valid at x, y G. We schematically show the regions of validity of Eqs. (B5)-(B11) in the (x, y) plane in Fig. 2(a). In the overlapping region at y x both ξ x and ξ x contributions are valid, but the Fabry-Perot contributions from ξ x naturally dominate. At y x the Fabry-Perot contributions are suppressed as e −πx/y = e −πhc/(T d) , since the temperature is lower than the first Fabry-Perot mode energy πhc/d. Nevertheless, it turns out that the prefactor in front of the exponential is large, so the contribution from the first mode (the one with the weakest exponential), coming from the narrow region around ξ = η = πx [see Fig. 2(b)] should be included together with the contribution from ξ x, as long as x, y G (otherwise, the mode is overdamped because of low reflectivity).
To pick up the first Fabry-Perot mode contribution, we approximate the Bose function by e −η/y and set η = πx everywhere else in the integrand, which is a smooth function of η. We also set ξ = πx everywhere in the integrand except the exponential e iξ/x in D p + [Eq. (B1b)]. Then we find the minimum of D p + as a function of ξ, reached at ξ min = πx(1 − x/G) + O(x 3 /G 2 ), and approximate near the minimum Then, the integration over πx < η < ∞ and −∞ < ξ − ξ min < ∞ gives e −πhc/(T d) . (B13)

TE travelling contribution
The TE travelling contribution to Eq. (A6) can be rewritten exactly as For G 1, we may not simply set G → 0 in the denominator, as we did in the TM case: here this leads to a logarithmic divergence at ξ → 0. To see how the divergence is cut off, we note that convergence scale of the η integral is the same as in the TM case: η ∼ y if y 1 and η ∼ 1 if y 1. This gives the small-ξ cutoff scales ξ ∼ Gη and ξ ∼ G, respectively. As a result, The overall map of behaviours in parameter space is therefore equivalent to the TM travelling case given in Eq. (B2), but the TE contribution (B15) is always dominant due to the logarithmic factors. The calculation for G 1 is very similar to that of the TM travelling wave contribution. Focusing firstly on the cases where ξ x so the exponentials e iξ/x oscillate fast, the averaged contribution from Eq. (B14a) via Eq. (B3) is given by At low frequency the system Fabry-Perot modes are indicated, as in the TM case, in the minima in D s ± , this time important when Gη ξ 1 + η 2 . The integral in η may again converge at η ∼ y due to the Bose function, or something else if y is too large. We may identify the same regions as in the TM case. (i) For y √ G, we have that ξ ∼ η ∼ y, so the fast oscillation condition is x y, and we may neglect all terms in the denominator containing ξ: (B17) (ii) For √ G y G, we keep ξη 2 in the denominator in the first line of Eq. (B16) and 2 (Gη) 2 in the second line (again, oscillations are fast when x y): (B18) (iii) For y G, we expand the Bose function to give y/η and retain ξη 2 in the first line of Eq. (B16) and (ξη) 2 + 2 (Gη) 2 in the second. The integrals converge at ξ, η ∼ G so the oscillations are fast when x G: (B19) For the contributions coming from ξ x, we expand the exponential e iξ/x : Since ξ < η and G 1 we may neglect ξ in the first bracket of the denominator in the last line. This allows the simple integration over ξ: the condition for the expansion of the exponential e iξ/x becoming η x. There are three possible cutoff scales for η: y, G, and (1 + G/x) −1 . Which one of the three is effective, depends on the positioning of y with repect to other scales. Again, three cases arise. (iv) For y (1 + G/x) −1 < 1, the logarithm is expanded for small argument and the second term in the denominator is neglected since the integral converges at η ∼ y. The condition y x for the expansion of e iξ/x is satisfied automatically: (v) For (1 + G/x) −1 y G, the integral is still determined by η ∼ y, but the second term in the denominator dominates. e iξ/x may be expanded when y x: G, the Bose function is y/η and the integral converges at η ∼ G so we retain the logarithm, and e iξ/x may be expanded as long as G x: (B24) We schematically show the regions of validity of Eqs. (B17)-(B24) for G 1 in the (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 contribution should be included together with the ξ x contributions as long as x, y G. The same procedure is performed whereby the minimum of D s + [Eq. (B14b)] near ξ = πx is found, allowing the integrand to be approximated by a Lorentzian. The minimum and therefore the eventual contribution is found to be identical to the TM case, Eq. (B13).

TM evanescent contribution
The TM evanescent contribution to Eq. (A6) can be rewritten exactly as This integral turns out to be exactly identical to that already calculated in Ref. [18] when the spatial dispersion of the conductivity is neglected (namely, Eqs. (1), (10) and (11) of Ref. [18]). That is to say that in the present system the Coulomb limit (c → ∞) amounts to taking only the exact TM evanescent contributions to the heat current, while neglecting the rest. Thus, we can simply rewrite the results of Ref. [18] in terms of G: where L(x) is a slow logarithmic function approximately given by: L(x) ≈ 4 ln 3 x 1 + (ln x)/ ln(1 + ln x) .
The domains of validity of the contributions are shown in Fig. 4. Note that expression (B26a) equals the travelling contribution (B9).

TE evanescent contribution
The TE evanescent contribution to Eq. (A6) can be rewritten exactly as Despite the apparent similarity to the corresponding TE travelling contribution Eq. (B14a), there is no longer oscillatory behaviour in the denominator, so the resulting contributions are completely different. In η there are two possible decay scales: η ∼ y from the Bose function, and η ∼ ξ/(G + ξ) fromD s +D s − . In the low temperature case y ξ/(G + ξ) < 1 where the temperature cutoff is effective, expanding e −ξ/x ≈ 1 leads to logarithmic divergence at ξ → ∞. The large ξ cutoff scale is therefore given by the decay scale of the exponential, ξ ∼ x, leading to the result [valid for y x/(G + x)]: J e TE =h G 2 π 2 c 2 τ 4 ∞ 0 η 3 dη e η/y − 1 ln x Gη = π 2 G 2 T 4 15h 3 c 2 lnh c GT d .
(B29) For high temperatures y ξ/(G + ξ) the Bose function is y/η and it is convenient to perform integration over η first keepingD s ± exact: where the integrand may decay due to the exponential or the denominator. If x G the exponential is clearly active and terms in ξ may be neglected in the denominator (the expansion of the Bose function is valid for y x/G): If x G, expansion of e −ξ/x ≈ 1 in Eq. (B30) again leads to logarithmic divergence at ξ → ∞. As in the low temperature case, the divergence is cut off by ξ ∼ x (the