The role of disorder in plasmon-assisted near-field heat transfer between two-dimensional metals

We perform a theoretical study of the near-field heat transfer between two atomically thin metallic layers, isolated galvanically but coupled by the Coulomb interaction, within the framework of fluctuational electrodynamics in the Coulomb limit. We clarify the role of disorder and spatial dispersion, and identify several distinct regimes of the heat transfer. We find that the plasmon contribution to the heat current is suppressed both in the clean and diffusive limits, but dominates in a parametrically wide crossover region at sufficiently high temperatures. In the diffusive limit, the heat transfer can be qualitatively modelled by an effective circuit theory.


I. INTRODUCTION
In metals, heat is most efficiently transported by electrons, as manifested by the Wiedemann-Franz law. In the absence of a galvanic contact between two nearby metals, heat can be transferred by electromagnetic fluctuations [1][2][3][4][5]. In the near field, when the typical distance between the metals is small compared to the thermal photon wavelength, the heat flux can be strongly enhanced compared to the Planckian radiative one, as was observed in many experiments (see reviews [6][7][8] as well as Refs. [9][10][11][12][13] for more recent experiments). Depending on the circumstances, such contactless heat transfer can be useful or harmful.
The standard framework for the theoretical description of this heat transfer is fluctuational electrodynamics [2,14]. It provides a very intuitive picture: the transferred heat can be viewed as energy dissipated in one part of the structure by fluctuating electromagnetic fields, which are created by fluctuating charges and currents in the other part of the structure. According to the fluctuation-dissipation theorem, both dissipation and fluctuations are determined by the material response functions (such as optical conductivity), which thus serve as the input to this phenomenological theory. The resulting heat current is then entirely determined by the system geometry and the model used for the material response. For the latter, models with different levels of sophistication have been used in the literature. They range from assuming a local frequency-dependent conductivity [2], to including non-local effects via spatial dispersion [15] or via surface contributions to the response [3,16]. From the literature, it is often difficult to decide which ingredients are important in which situation.
In particular, in several works dedicated to different materials in the near-field regime, the important role played by collective plasmon excitations has been pointed out [17][18][19][20][21][22][23][24][25]. However, Ref. [26] concluded that surface plasmons were unimportant for the heat transfer between * jonathan.wise@lpmmc.cnrs.fr † Deceased 15 May 2017 two bulk semi-infinite metals, and in Ref. [27] no plasmon contribution was reported for heat current between two thin metallic layers in the clean and diffusive limits.
To clarify this issue, we revisit the problem of heat transfer between two thin parallel metallic layers, assuming the separation between them to be much smaller than the thermal photon wavelength, so that the transfer is dominated by the Coulomb interaction. Including disorder and spatial dispersion within the framework of fluctuational electrodynamics in the Coulomb limit, we identify different regimes of interlayer separation and temperature where different physical ingredients are important. We find that the plasmon contribution to the heat current is suppressed both in the clean and diffusive limits, but dominates in a parametrically wide crossover region at sufficiently high temperatures. We also show that in the diffusive regime, the heat transfer can be understood in terms of an effective RC circuit [28].
The paper is organised as follows. In Sec. II we specify our model for the system, discuss the main regimes of heat transfer, and present the resulting expressions for the heat current. We give their detailed derivation in Sec. III, and verify numerically in Sec. IV. In Sec. V, we interpret our results in terms of circuit theory. In Appendix A, we give a derivation of the density response function for a 2D metal with short-range impurities, interpolating between clean and diffusive limits.

II. THE MODEL AND SUMMARY OF RESULTS
We consider two thin parallel metallic layers, labeled 1 and 2, separated by a distance d, and kept at different temperatures T 1 and T 2 . We assume to be in the limit d c/T 1,2 , where c is the speed of light, and we use units where Planck and Boltzmann constants are set to one. In this limit, heat transfer is dominated by the Coulomb interaction. Our analysis is based on the standard expression for the heat current per unit area [21,22,27,29,30], often called the Caroli formula, in analogy to a similar expression for electron current across a tunnel junction [31]. In Eq. (1), N 1,2 (ω) = 1/[exp(ω/T 1,2 ) − 1] are the Bose distribution functions in the two layers. In Eq. (2), (3) is the Coulomb interaction between the layers, screened in the random-phase approximation, with v q = 2πe 2 /q being the bare 2D Coulomb potential and e < 0 the electron charge. Finally, Π 1 (q, ω) and Π 2 (q, ω) are the susceptibilities determining the linear response δρ i (r, t) = Π i (q, ω) eϕ q,ω e iqr−iωt of the 2D electron density ρ i in the corresponding layer i = 1, 2 to a perturbing electrostatic potential ϕ q,ω e iqr−iωt applied to that layer. The density response function is related to the in-plane longitudinal optical conductivity, σ(q, ω) = (iω/q 2 )e 2 Π(q, ω). Eq. (1) can be derived from the Coulomb limit of fluctuational electrodynamics [21], from non-equilibrium Green's function formalism [29], or the kinetic equation [27]. The density response functions encode all material characteristics of the two layers. We model each layer as a degenerate isotropic 2D electron gas with short-range impurities. Such a system can be characterised by three parameters: (i) ν, the electronic density of states per unit area at the Fermi level including both spin projections, whose energy dependence is neglected; (ii) v F , the Fermi velocity; (iii) τ , the elastic scattering time. The Fermi momentum p F is assumed to be the largest momentum scale, . Under these assumptions, the density response function of each layer is temperature-independent and given by [32] Π(q, ω) = −ν 1 + iωτ for an arbitrary relation between ω, v F q, and τ . Equation (4) interpolates between two well-known expressions corresponding to the clean limit (τ → ∞) and the diffusive limit (ω, v F q 1/τ ): where D = v 2 F τ /2 is the diffusion coefficient. Equation (4) corresponds to a 2D analog of Mermin's prescription [33], which was recently used to model the response in monolayer graphene [18,34]. In Appendix A we give a derivation of Eq. (4) based on Boltzmann equation for electrons scattering on impurities. In spite of several assumptions underlying Eq. (4), its simplicity enables one to obtain an important insight into the interplay between the spatial dispersion [that is, the q dependence of the conductivity σ(q, ω)] and the impurity scattering. For simplicity we also assume the two metals to be identical. . For q 1/d (to the right of the hatched area), T (q, ω) is suppressed by the factor e −qd . In the hatched area, T (q, ω) is well approximated by the clean limit expression (5a), while the shaded area q 1, ωτ 1 corresponds to the diffusive limit (5b). In the white region above the shaded and hatched areas (ω 1/τ , ω > vF q) the integrand is small except, maybe, the vicinity of the symmetric and antisymmetric plasmon dispersions (upper and lower solid lines, respectively) where |V12| 2 is peaked.
The three independent material parameters, introduced above, define two important length scales: (i) the 2D screening length 1/κ = (2πe 2 ν) −1 , and (ii) the mean free path = v F τ . Typically, the screening length is very short (on the atomic scale), so we assume κ 1 and κd 1. The mean free path can vary from several nanometers to several microns, and may be larger or smaller than the separation d between the two layers.
In the isotropic model, formulated above, Eq. (1) represents a two-dimensional integral over q and ω. This integral can be rather straightforwardly evaluated numerically, but much better insight into relevant physics is obtained by studying different asymptotics of the integral in various parameter regimes. The latter approach was adopted in Ref. [27] using the two limiting expressions (5a), (5b). It turns out, however, that these expressions miss the plasmon contribution.
We show schematically the behavior of of T (q, ω) [Eq. (2)] in the (q, ω) plane for d in Fig. 1. At large ω > ∼ max{T 1 , T 2 }, the integrand in Eq. (1) is suppressed by the Bose function at ω, and this cutoff may be positioned anywhere in Fig. 1, depending on the temperatures. At large q > ∼ 1/d, the integrand is suppressed by e −qd in the numerator of Eq. (3); Fig. 1 corresponds to the case d , but for larger d the spatial cutoff may shift to the diffusive shaded area.
The strongly coupled plasmon modes (in the case of identical layers, symmetric and antisymmetric, denoted by "±") manifest themselves as poles of V 12 (q, ω). In the clean limit, τ → ∞, the plasmon frequencies are real and given by (for q κ) At finite τ , but such that ω ± τ 1, the poles acquire a small imaginary part, so |V 12 (q, ω)| 2 is peaked around the dispersions (6). At ω ± τ < ∼ 1, when the diffusive expression (5b) applies, the plasmons are overdamped and do not produce a separate contribution to the integral. In the strictly clean limit, τ → ∞, their contribution vanishes as well, since for ω > v F q the integrand vanishes because Im Π(q, ω > v F q) = 0.
For temperature-independent Π(q, ω), Eq. (1) naturally splits into the difference J(T 1 , A detailed analysis of different asymptotic regimes for the integral in Eq. (1) (given in the next section) results in several asymptotic expressions for J(T ): In the labels "l/h" denotes low/high temperature, while "c/p/d" stands for clean/plasmonic/diffusive. Here ζ(x) is the Riemann zeta function, and L(x) is a slow logarimic function defined in Eq. (15). For moderate values of log( 2 κ/4d) < 10, it can be approximated with 10% precison as Figure 2 schematically shows the domains of validity for expressions (7) in the (1/d, T ) plane. The clean and diffusive regimes were also identified in Ref. [27]. Plasmon contributions dominate in a parametrically wide region of the parameter space. Crucially, this behaviour is captured by neither of the limiting expressions (5a) nor (5b), but by the leading term of the small q expansion of the full expression (4). This is equivalent to using the Drude expression for the optical conductivity, σ(ω) = e 2 νD/(1 − iωτ ), that is, neglecting the spatial dispersion. We find that the spatial dispersion can also be neglected to describe the diffusive contribution. It becomes important only in the clean region where J(T ) is dominated by the hatched area in Fig. 1. In all three high-temperature regions, one can approximate the Bose distribution as N (ω) ≈ T /ω, so the density fluctuations can be treated classically and the resulting J(T ) ∝ T . Expressions (7) were obtained for two identical metallic layers. If they are different, but the material parameters ν, v F , τ are of the same order, our expressions are still valid as order-of-magnitude estimates. In particular, this applies to the plasmon contribution: the plasmons remain strongly hybridised even when the bare plasmonic dispersions of each layer do not match exactly. The case when the two materials are strongly different, is beyond the scope of the present paper.
Equation (4) was derived under the assumption of short-range impurities with purely s-wave scattering. If this assumption is relaxed, Eq. (4) is not valid quantitatively. Still, expressions (7) remain valid as order-ofmagnitude estimates if τ is taken to be the momentum relaxation time (see Appendix A for discussion).

III. DERIVATION OF THE RESULTS
For identical layers, Eq. (2) can be written in a slightly simpler form: Below, we analyze the contributions from the three regions of the (q, ω) plane in Fig. 1: the vicinity of the plasmon dispersions, the diffusive region, and the clean region. In each region, we separate two temperature regimes. We will see that in the high-temperature regimes, the integral (1) is dominated by a certain frequency scale, different in each regime, but always determined by q ∼ 1/d and being much lower than the temperature. Then, the thermal cutoff plays no role, and one can approximate the Bose distribution N (ω) ≈ T /ω [that is why J(T ) ∝ T in Eqs. (7b), (7d), (7f)]. In the low-temperature regimes, it is the e −qd cutoff which is ineffective, so the frequency integral is determined by ω ∼ T , while the q integral converges on a scale much smaller than 1/d. Then one can approximate e −qd ≈ 1 in the numerator of Eq. (9) and 1 + e −qd ≈ 2, 1 − e −qd ≈ qd in the denominator.

A. Plasmon contribution
As discussed earlier, the plasmon contribution comes from the region ω > 1/τ, v F q, since otherwise the plasmons are overdamped. At such frequencies one can expand Π(q, ω) to the leading order in q and approximate which corresponds to neglecting the spatial dispersion in the conductivity which takes the Drude form, σ(ω) = e 2 νD/(1−iωτ ). (Generally, the spatial dispersion can be neglected when q max{ ω/D, ω/v F }.) The integral is then dominated by the vicinities of the two plasmon dispersions, where one of the factors in the denominator of Eq. (9) is small. The plasmon contribution exists only if the plasmon frequencies ω ± at q ∼ 1/d (when the spatial cutoff becomes effective) exceed 1/τ . This gives a condition d κ 2 (the vertical line x = η 2 in Fig. 2). Let us first consider the temperature interval 1/τ T v F κ/d (the two inequalities are consistent when d κ 2 ), where the thermal cutoff plays first leaving the spatial cutoff e −qd ineffective. Then one can expand e −qd and perform the q integration first, approximating the integrand by a Lorentzian in the vicinity of each pole. The remaining frequency integral can be calculated exactly: Here ζ(x) is the Riemann zeta function, and the first/second term in the square brackets comes from the antisymmetric/symmetric plasmons, respectively. In the considered region 1/τ T v F κ/d the symmetric contribution is always small compared to the antisymmetric one.
At high temperatures, T v F κ/d, the integral is determined by the spatial factor e −qd , while the thermal cutoff is ineffective so the Bose distribution N (ω) ≈ T /ω. Taking expression (10), Eq. (9) can also be written in the form suitable for integration over ω: (12) Note that one can not just do two separate Lorentzian integrals because the separation ω 2 + − ω 2 − becomes exponentially small at q > ∼ 1/d. Fortunately, the ω integral can be calculated exactly: Then the q integral reads as leading to Eq. (7d) with the function L(x) defined as

B. Diffusive contribution
Let us focus on the contribution from the shaded region in Fig. 1: q 1/ , ω 1/τ . Then one can use expression (5b) for Π(q, ω), and since we are interested in q < ∼ 1/d κ, we have Dq 2 Dκq(1 ± e −qd )/2, which is again equivalent to neglecting the spatial dispersion in the conductivity. Thus, we can write (16) At low frequencies, when the momentum integral should converge on some scale q 1/d, the two factors in the denominator are strongly different, and it is the second factor that determines the convergence scale q ∼ ω/(Dκd). When ω ∼ T Dκ/d this scale is indeed much smaller than 1/d, so we expand e −qd , integrate over q, then over ω, and arrive at Eq. (7e). Moreover, the convergence scale T /(Dκd) 1/ provided that T κd/τ . Thus, since we always assume κd 1, Eq. (7e) is valid for the diffusive contribution everywhere below the horizontal line y ∼ η (T ∼ 1/τ ) in Fig. 2.
For T Dκ/d, the q integral is dominated by q ∼ 1/d. Then the typical frequency scale of Eq. (16) is ω ∼ Dκ/d, so for T Dκ/d the thermal cutoff is ineffective. Then we approximate N (ω) ≈ T /ω, straightforwardly integrate over ω, then over q, and arrive at Eq. (7f). Note that Eq. (7f) is valid even at T 1/τ provided that the convergence scale ω ∼ Dκ/d 1/τ , that is, to the left of the vertical line x = η 2 in Fig. 2. As we have seen, to the right of this line the plasmon contribution becomes important.

C. Clean contribution
For T 1/τ , 1/d 1/ , one should take into account the contribution from the hatched area in Fig. 1. Here one can take the limit τ → ∞ and use the expression (5a) for Π(q, ω). Then, in the integration region ω < v F q, Re Π(q, ω)v q = −κ/q, so one can neglect unity in both factors in the denominator of Eq. (9) and write For T v F /d we approximate N (ω) ≈ T /ω, straightforwardly integrate over ω between 0 and v F q, then integrate over q, and arrive at Eq. (7b). For T v F /d, in most of the integration region we have ω ∼ T v F q, so the upper limit ω < v F q is not important except the narrow region q ∼ T /v F which determines the lower cutoff of the logarithmic q integral: leading to Eq. (7a).
In the region T 1/τ , 1/d 1/ , the clean contribution and the plasmon contribution both exist and should be added up, since they come from two distinct regions in the (q, ω) plane. Thus, to determine the dominant asymptotics, one can combine the four expressions (7a)-(7d) as J = max{min{J lp , J hp }, min{J lc , J hc }} which results in the complicated shape of the boundary between the clean and plasmonic regions in Fig. 1.

IV. NUMERICAL VERIFICATION
In order to illustrate the various behaviours and crossovers indicated in Fig. 2, and to verify our asymptotic expressions (7), we evaluate numerically the integral in Eq (1) using the full response function (4).
To describe clean and plasmonic regimes, we take parameters typical for doped graphene. If the temperature is well below the Fermi energy F (counted from the Dirac point) the electrons in graphene can be viewed as a conventional 2D electron gas with the density of states per unit area ν = 2| F |/(πv 2 F ), including the valley and spin degeneracies. In Fig. 3, we plot J(T ) for two sets of parameters corresponding to 1/(κd) ≈ η 1/3 = 0.034, and 1/(κd) ≈ η 1/2 [we remind that η ≡ 1/(κ )].
To study the diffusive crossover for realistic materials, we introduce dielectric screening. Equation (3), written for metallic layers surrounded by vacuum, can be generalised to the situation when the layers are embedded in a dielectric medium. This generalisation is particularly simple when the medium is characterised by a uniaxially anisotropic dielectric constant, ε in the plane parallel to the layers, and ε ⊥ along the z direction, perpendicular to the layers. The solution of the Poisson equation in such a medium gives the 2D Coulomb potential at a distance z from a charged layer: Thus, in all expressions (7) it is sufficient to rescale In Fig. 4 we show the crossover between low-and high-temperature diffusive asymptotics (7e) and (7f) for two hole-doped tungsten diselenide monolayers embedded in boron nitride. The valence band of WSe 2 is parabolic with the hole effective mass m h being about half of the free electron mass. The spin degeneracy is lifted by a strong spin-orbit coupling, so only valley degeneracy remains, and the density of states per unit area is ν = m h /π. We take F = 50 meV and a very short = 2 nm, still consistent with κ 1, p F 1. The taken separation d = 100 nm corresponds to 1/(κd) = 0.006, well below η 2 = 0.1, and hence to the diffusive region.
We are not showing the high-temperature clean and plasmonic regimes; for realistic material parameters, they correspond to temperatures so high, that the assumptions behind our model (degenerate Fermi gas, nearfield Coulomb regime) are no longer valid. However, we checked numerically the validity of the asymptotic expressions (7b) and (7d) for J hc (T ) and J hp (T ).

V. COMPARISON TO CIRCUIT THEORY
It is interesting to relate our results to those obtained by the circuit approach [28]. This approach can be viewed as a form of fluctuational electrodynamics in a system made of lumped elements (capacitors, inductors, and resistors). Their dissipative parts provide the energy storage and thermal fluctuations (Johnson-Nyquist noise), while the reactive parts mediate the electromagnetic interactions, resulting in energy exchange. For the two metallic layers studied here, it is natural to expect an analogy to the circuit consisting of two RC contours connected by two coupling capacitors C c , as shown in Fig. 5. Indeed, the electronic excitations in each layer constitute a dissipative bath analogous to a resistor. To mimic charge oscillations within each layer, the resistor should be shunted by a capacitor. The Coulomb interac- tion between the layers resembles that between the plates of a capacitor.
The power transferred from the first to the second resistor is given by [28] where Z 1,2 (ω) = (1/R 1,2 − iωC 1,2 ) −1 is the impedance of each RC contour. As before, for simplicity we assume the two subsystems to be identical, R 1 = R 2 = R, C 1 = C 2 = C. Writing the transmission as we obtain the following asymptotic expressions for the transferred power P (T 1 , T 2 ) = P (T 1 ) − P (T 2 ): .
The intermediate "universal" regime, where the power depends only on the temperature, but not on the circuit, is present only when C C c . To relate results (24) to those of Sec. III, it is important to realise that while Eqs. (24) give the full transferred power, Eqs. (7) give the power per unit area. To make a meaningful comparison we must therefore invoke a length scale L, such that the infinite sample can be divided into squares of size L. Then Eqs. (24) describe the power transferred in each square, and the contributions of different squares can be added up independently.
Thus, the relevant length scale should be associated with the typical convergence scale of the q integral in Eq. (1).
As any circuit made of resistors and capacitors, that of Fig. 5 has an overdamped dynamics, so it corresponds to the metallic layers in the diffusive limit. Indeed, one can notice the similarity between Eqs. (16) and (23). Then, it is natural to associate the resistance R to the resistance per square 1/σ of each metallic layer, R ∼ 1/σ = 2π/(κD). The coupling capacitance is associated to the geometric capacitance between the two layers, C c ∼ L 2 /d, where the in-plane length scale L > ∼ d must be invoked because the capacitance is proportional to the area. The capacitance C should be associated to an intrinsic property of each layer, such that RC corresponds to the characteristic relaxation time of charge density modulations. In a two-dimensional system, this time depends on the in-plane length scale L and is then given by L/σ. This gives C ∼ L. Recalling the convergence scales of the q integral in Sec. III B, we associate L ∼ Dκd/T and L ∼ d in the low-and high-temperature diffusive regimes, respectively. Then Eqs. (24c) and (24a) match Eq. (7e) and (7f) at low and high temperatures, respectively. Expression (24b) does not correspond to any parametric region because at T Dκ/d the two capacitances become of the same order.

VI. CONCLUSIONS
We have studied the problem of heat transfer between two thin parallel metallic layers, mediated by the Coulomb interaction. Using a simple model for a 2D electron gas subject to scattering on short-range impurities, we described the crossover between clean and diffusive limits and showed that strongly coupled surface plasmons dominate the heat transfer in a parametrically wide region at sufficiently high temperatures, but their contribution is suppressed both in the clean and diffusive limits. We also clarified the role of the spatial dispersion in the optical conductivity, which turns out to be important only in the clean limit. In all other regimes, the effect of disorder is correctly captured by the relaxation time in the Drude conductivity.
We have shown that in the diffusive limit, the heat transfer is similar to that in an effective RC circuit. However, for this analogy to be meaningful, one must specify a length scale. This length scale should correspond to the size of regions within the infinite 2D sample where the transfer occurs independently. In other words, each region can be described by a separate circuit, and contributions from different regions add up. This length scale must be determined from the microscopic theory.
An explicit solution of this equation can be found only in the case of short-range impurities when U (p − p ) does not depend on momentum. In this case the collision integral reduces to the relaxation time approximation: where the overbar denotes the average over the momentum directions on a constant energy surface, and the relaxation time and the density of states per unit area are given by (the factor of 2 in front of the integral takes into account two spin projections) Then, Eq. (A5) gives g p = iqv p (∂f eq /∂ p ) + (1/τ ) g p iqv p − iω + 1/τ .
Averaging both sides over the momentum directions, one obtains a closed equation for g p and readily finds where S stands for the following angular average: The last two lines were written under the assumption of an isotropic dispersion p . Finally, since the electron density is given by ρ(r, t) = 2 d 2 p (2π) 2 f p (r, t) (A12) (again, the factor of 2 takes care of the spin multiplicity), the density response function can be found as Π(q, ω) = 2 d 2 p (2π) 2 g p .
Collecting all factors and assuming that −∂f eq /∂ p is a narrow peak around the Fermi energy of width ∼ T , so that the energy dependence of ν and v F can be neglected, we arrive at Eq. (4). If U (p − p ) is momentum dependent, no closed solution can be obtained even in the simplest of isotropic scattering when the scattering amplitude depends only on the difference φ − φ of the polar angles φ, φ associated to the 2D vectors p, p . Indeed, in this case the solution can be sought as a sum over polar harmonics, g p = m g m e imφ , which are the eigenfunctions of the collision integral. Different harmonics do not separate because of the second term in Eq. (A5), which is responsible for the spatial dispersion of the conductivity σ(q, ω). Only when the spatial dispersion is neglected, the Drude conductivity σ(ω) can be written in terms of the transport relaxation time τ 1 , determined by the first eigenvalue of the collision integral, −1/τ 1 . Otherwise, the result contains all eigenvalues −1/τ m>0 . Still, qualitatively, it is τ 1 that determines the relevant time scale: in the limit of Eq. (A6) all τ m>0 = τ , while in the opposite limit of small-angle scattering τ m quickly grows with m, so high harmonics are suppressed and the result is determined by the first few τ m 's.