Ultrafast Radiative Heat Transfer

Light absorption in conducting materials produces heating of their conduction electrons, followed by relaxation into phonons within picoseconds, and subsequent diffusion into the surrounding media over longer timescales. This conventional picture of optical heating is supplemented by radiative cooling, which typically takes place at an even lower pace, only becoming relevant for structures held in vacuum or under extreme conditions of thermal isolation. Here we reveal an ultrafast radiative cooling regime between neighboring plasmon-supporting graphene nanostructures in which noncontact heat transfer becomes a dominant channel. We predict that>50% of the electronic heat energy deposited on a graphene disk can be transferred to a neighboring nanoisland within a femtosecond timescale. This phenomenon is facilitated by the combination of low electronic heat capacity and large plasmonic field concentration displayed by doped graphene. Similar effects should take place in other van der Waals materials, thus opening an unexplored avenue toward efficient heat management in ultrathin nanostructures.


I. INTRODUCTION
Optical, electrical, and mechanical dissipation in nanoscale devices produces heat accumulation that can result in structural damage and poor performance. Understandably, heat management constitutes an important aspect when designing thermoelectric [1], optoelectronic [2], electromechanical [3], and photovoltaic [4] elements, as well as recently proposed thermal analogs of electronic devices [5,6]. However, the relatively slow thermal conduction in most materials [7] imposes a serious limitation. Finding new means of cooling nanostructures is therefore critical. An interesting possibility is provided by coupling to radiative degrees of freedom. Indeed, the absorption and emission of radiation by a material structure contributes to reach thermal equilibrium with other surrounding structures and the electromagnetic environment. This is the dominant cooling channel for thermally isolated structures [8], in which energy is released through the emission of photons with wavelengths ∼ λ T = 2πhc/k B T (i.e., the thermal wavelength at temperature T ). When the structures are separated by vacuum gaps of large size compared with λ T , the Planck and Kirchhoff laws determine the exchanged power [9]. In contrast, for neighboring objects separated by a small distance compared with λ T , radiative heat transfer is dominated by additional channels mediated by evanescent waves [10][11][12]. These can produce rates exceeding the black-body limit by several orders of magnitude, en-hanced by near-field coupling of resonances supported by the nanostructures, thus emerging as a potentially relevant transfer mechanism in solid state devices.
Following pioneering observations of near-field radiative energy transfer between two conducting plates [10,11], a theoretical explanation was offered [12] based on the effect of thermal fluctuations in the electrical current of the involved surfaces. Further experimental [13][14][15][16][17][18][19][20][21][22] and theoretical [5, studies have corroborated this interpretation of radiative heat transfer between structures of varied morphologies. This subject has generated fundamental insights that include important corrections due to nonlocal [30], phonon [27,42], and photonic band [47] effects, as well as magnetic polarization [34]. Additionally, retardation, radiation emission, and crossed electric-magnetic terms in the optical response have been shown to severely modify the transfer power [50]. However, the so far observed and predicted transfer rates are slow compared with dissipative transport through the surrounding media, in which heat can cause undesired effects. This situation persists even when the interaction between neighboring structures is enhanced due to strong resonant excitations, such as plasmons in noble metals.
In this context, graphene plasmons can be advantageous because their frequencies lie in the mid-infrared, which is the spectral region for thermal interactions under attainable temperatures. Indeed, plasmon energies in graphene nanostructures scale as ∼ E F /D with the Fermi energy E F and the characteristic size D (e.g., the diameter for a disk). Doping levels as high as E F ∼ 1 eV have been reported through electrostatic gating [54], and even higher values through chemical doping [55,56], manifesting themselves in the opening of a 2E F gap for vertical optical transitions [54,57]. However, plasmons are only well defined at energies below ∼ E F due to the narrowing of the gap as their momentum increases [58]. For reference, a 20 nm disk supported on silica and doped to E F = 1 eV exhibits a dipolar plasmon at ≈ 0.4 eV [58]. This explains why experiments have only explored mid-infrared plasmons, as higher energies require smaller structures, whose fabrication can be challenging.
An additional advantage of graphene lies in its large electrical tunability, which enables an active control of these phenomena. In a related context, electrical modulation of thermal emission of radiation has been accomplished in gated nanostructured graphene [59], while an optical-to-thermal converter has been proposed to be capable of efficiently transforming an optical pump into light emitted at longer mid-infrared wavelengths [60]. Electrical control of radiative heat transfer between graphene-coated surfaces or between extended graphene and other materials has been also proposed [61][62][63][64][65].
The competing mechanism (relaxation into phonons) was initially thought to be rather slow in graphene [66] (nanosecond scale), a prediction that was subsequently corrected to much shorter timescales (picoseconds) due coupling of hot charge carriers to optical phonons [67] and so-called supercollision cooling [68]. The latter is consistent with experimental observations [69,70]. Recent calculations have also identified a remarkably fast rate of radiative transfer between graphene films [62,71], graphene nanoribbons [72], and extended heterostructures of graphene and hexagonal BN [73], although all of them involve picosecond or even longer timescales. However, we need much faster transfer rates in order to prevent most of the electronic heat from being absorbed into phonons. We accomplish such a goal in this paper by resorting to graphene nanostructures capable of sustaining plasmons within an energy range that is commensurate with k B T . Incidentally, radiative energy transfer from graphene electrons to optical phonons in a silica substrate has been argued to explain the measured saturation of conductivity in the carbon layer and provide a viable way of observing quantum friction [74].
Here we exploit the extraordinary optical and thermal properties of graphene to show that ultrafast radiative heat transfer can take place between neighboring nanoislands. The commonly accepted scheme for dissipation of the thermal energy produced by electronic and optical inelastic losses (i.e., energy transfer to valence and conduction electrons of the system, followed by relaxation into phonons and subsequent heat flow into the surrounding media) is here challenged by the radiative transfer mechanism taking place between neighboring structures within femtosecond timescales, thus overcoming electron relaxation into the atomic lattice. Using attainable graphene nanostructure designs, we find that ultrafast radiative heat transfer produces thermalization of two neighboring islands that results in >50% of the electronic heat of the hot one being radiatively transferred to its colder neighbor. This extraordinary phenomenon is made possible by the large plasmonic field concentration that mediates the coupling between the neighboring graphene structures, as well as by the low specific electronic heat of this material [58]. In particular, plasmons in this material exhibit unprecedentedly large electrical tunability accompanied by strong confinement of the measured fields [75,76], which have recently enabled high mid-infrared sensitivity in the detection of proteins [77] and other organic molecules [78]. In a similar fashion, the ultrafast radiative heat transfer phenomenon here investigated can be actively switched on and off by gating the graphene structures.

A. Radiative heat transfer between graphene nanodisks
We focus on the system depicted in Fig. 1, consisting of two parallel coaxial graphene nanodisks of diameters D 1 and D 2 , separated by a distance d between carbon planes, doped to Fermi levels E F1 and E F2 , and having electronic temperatures T 1 > T 2 . For simplicity, we consider the disks to be placed in vacuum, as the conclusions of this work remain the same when the disks are surrounded by a dielectric material such as BN (e.g., ∼ 3.2, see Fig. 8). Heat is radiatively transferred from the hotter disk to the colder one as a result of thermal fluctuations in both disks, whose interaction is mediated by their self-consistent electromagnetic response. In fact, for the small size of the structures under consideration compared with the thermal wavelengths λ T (with = 1, 2), retardation and magnetic response effects can be dismissed, so we only need to deal with charge fluctuations and their Coulomb interaction.
We calculate the heat transfer power (HTP) as the net balance of the work done by the thermally fluctuating charges of the hotter disk on the colder one minus the work done on the former by the fluctuating charges of the latter. This leads to a classical electromagnetic expression involving thermal fluctuations, which are evaluated by means of the fluctuation-dissipation theorem [79,80]. A detailed self-contained derivation is offered in the Appendix, leading to a compact expression [Eq. (A6)] that is proportional to the integral over the exchanged frequency ω. The integrand consists of the difference between the Bose-Einstein occupation numbers n = [exp(hω/k B T ) − 1] −1 of the two disks at their respective temperatures T , multiplied by a loss function that is determined by the disk susceptibilities χ . The latter are dominated by plasmonic modes, which allow us to formulate a description in terms of plasmon wave functions (PWFs) [81,82]. Only the lowest-order PWFs contribute significantly to the HTP for the range of geometrical parameters under consideration. Their explicit form (see Appendix), as well as full details on the PWFbased susceptibilities, are given in the Appendix. For coaxial disks (Fig. 1), we find that modes of different az-FIG. 1: Sketch of the structure considered for ultrafast radiative heat transfer. We study heat transfer between two parallel coaxial graphene disks placed in vacuum and separated by a small distance d. Each disk = 1, 2 is characterized by its diameter D , Fermi energy E F , and electron temperature T , with T1 > T2.
imuthal number m do not mix, so we can separate their contributions to the HTP received by disk 2 as (and also P 1 = −P 2 ), where Tr[. . . ] stands for the trace, the matrix ∆ m = (I − χ m 2 · v m · χ m 1 · v m ) −1 accounts for multiple scattering between the disks, v m describes their mutual Coulomb interaction, and I is a unit matrix. The matrices v m and χ m contain elements projected on the PWFs with m azimuthal symmetry (see Appendix for detailed expressions). Incidentally, the leading (2 − δ m0 ) factor reflects the fact that m and −m modes yield the same contribution.
In this formalism, the optical response of graphene is described through its surface conductivity σ, for which we adopt the local-RPA model [58,83,84] [see Eq. (F1) in the Appendix]. We remark that, besides the explicit dependence of n on T , the temperature enters σ through the chemical potential as well (see Appendix). It should be noted that, in contrast to extended graphene, the lack of translational invariance in nanostructures prevents us from using the full nonlocal RPA conductivity [85,86]. However, a full RPA description of the optical absorption of the system under consideration based on a previous implementation for finite structures [87] reveals that nonlocal effects only play a small role (see Fig. 9). We further analyze heat transfer between closely spaced extended graphene films, and more specifically, the contribution coming from parallel wave-vector components ∼ 2π/D , for which we find that nonlocal effects are also small for the graphene parameters under consideration (see Sec. K and Fig. 11), and therefore, we also expect them to be small for disks of diameter D .
Incidentally, as the HTP of Eq. (1) is an integrated quantity, it is not too sensitive to the model used for the graphene conductivity σ. This is corroborated in Fig.  10(a,b), where we compare results obtained using either the local-RPA or the Drude model [Eq. (F1) with the E integral set to zero]. Only small discrepancies between the two models are observed at small separations d in the resulting HTP. Actually, the small d region is most sensitive to elements of the formalism such as the inclusion of multiple scattering in the optical response of the disks [∆ m matrices in Eq. (1), see Fig. 12 for a comparison with results obtained by setting ∆ m = I]. We also observe a mild dependence of the HTP on the value of the intrinsic electronic decay time (Fig. 13), which we set tohτ −1 = 10 meV throughout this work. Additionally, we find good convergence of the HTP with the number of m's and PWFs used in the calculations (Fig. 14).
We stress that the relatively high temperatures under consideration (thousands of degrees) refer to the electronic gas of the material, which can be reached by optical pumping in the ultrafast regime [88][89][90].
The disk separation dependence of the HTP is studied in Fig. 2(a) (solid curves) for 20 nm graphene disks doped to a Fermi level E F = 0.2 eV, with the hotter disk at different temperatures T 1 (see labels) and the colder one at room temperature T 2 = 300 K. In general, higher temperatures T 1 lead to larger HTP, due in part to the (n 1 − n 2 ) factor in Eq. (1). At large separations d D , only dipole-dipole interactions between the disks contribute efficiently to the transfer, leading to a 1/d 6 dependence, in agreement with the asymptotic expression of Eq. (A7) (see Appendix). A smooth convergence of the full calculation [Eq. (1)] to this limit [Eq. (A7)] is observed in the additional calculations presented in Fig. 15. The nearfield character of heat transfer is further emphasized by considering the extension of the dominant dipole plasmon away from the disks (i.e., the electric-field amplitude decays by 1/e over a distance ∼ D/2π, as estimated from the out-of-plane decay of plasmons in extended graphene for an equivalent wavelength ∼ D), which explains the low slope in the curves of Fig. 2(a) at small d's.
As a reference, we compare these results with the HTP for gold disks of the same diameter [ Fig. 2(a), broken curves], which we describe through an effective surface conductivity obtained from the measured dielectric function [92] Au as σ Au = iωt(1 − Au )/4π, where we take a thickness t = 2 nm. This approximation, which is reasonable because we are considering a small value of t compared with the diameter (20 nm), allows us to apply the same formalism as for graphene [Eq. (1)]. Despite the larger thickness of the gold disks, their HTP is much smaller than for graphene. In fact, plasmons in the graphene disks lie in the mid-infrared region for the parameters under consideration (i.e., they energies are commensurate with k B T 1 ), while those of the gold disks appear at much higher energies, and thus do not contribute efficiently to the heat transfer. This mismatch is partly alleviated at the highest temperature under consideration (T 1 = 5000 K), for which gold and graphene disks exhibit similar HTPs in the large d limit.
As an additional comparison, the left arrows in Fig. 2(a) show an estimate obtained from the Stefan-Boltzmann law [32] for radiative heat transfer between two blackbodies of an area equal to that of the present disks. As anticipated above, graphene outperforms blackbodies by several orders of magnitude.
The strength of their optical response influences the ability of the disks to transfer energy radiatively. This is examined in Fig. 2(b), where we plot the absorption cross-section of one of the graphene disks considered in Fig. 2(a). An intense plasmon feature is observed in the 0.2-0.4 eV region, whose temperature dependence is inherited from the conductivity [Eq. (F1)]. The dashed line in Fig. 2(b) shows the relation between the temperature and the photon energy according to Wien's law (i.e., the value ofhω at the maximum of ω 3 n (ω) as a function of T ). This is relevant for the analysis of Eq. (1), in which a factor ω n (ω) appears explicitly, whereas the remaining ω 2 factor comes from the low ω limit of the Im{χ m } matrices [obviously, the full ω dependence of the integrand of Eq. (1) is more complex, as shown in Fig. 10(g,h), but an analysis based on Wien's law is still informative]. Additionally, the response functions entering the trace in Eq.
(1) display maxima near the plasmons, and therefore, the overlap between the dashed line and the plasmon in Fig.  2(b) indicates that this excitation contributes efficiently to the HTP, thus providing a criterium for optimization. Incidentally, the plasmon dispersion and strength follow nonmonotonic behaviors resulting from the complex interplay between the increase in both the density of free charge carriers and the number of decay channels associated with single-electron transitions.
The electronic heat capacity provides a relation between the temperature and the amount of energy strored in the electron gas. In this respect, graphene is also advantageous relative to traditional plasmonic materials such as gold because its heat capacity is orders of magnitude smaller [ Fig. 2(c)] as a result of its conical band structure, in contrast to the parabolic dispersion of gold conduction electrons. In consequence, cooling the graphene electrons requires transferring a smaller amount of heat, thus making the process potentially faster.

B. Ultrafast radiative heat transfer regime
We study the heat transfer dynamics by considering the electronic heat Q deposited on each graphene disk and the evolution of these quantities according to the equationsQ where P are the transfer powers given by Eq. (1), while τ ph is a phenomenological electron relaxation time (to phonons) that we approximate as 1 ps, a value of the order of what is observed in pump-probe experiments [67,93]. We note that the electronic heat of each disk depends on the electronic temperature T as Also, the transfer powers P 1 and P 2 = −P 1 [Eq. (1)] implicitly depend on both temperatures T 1 and T 2 . In order to make this clearer, we provide equations equivalent to Eqs. (2) at the end of the Appendix with a more explicit dependence on the temperatures, along with details of the numerical solution method. It should be pointed out that, because the electronic heat capacity in graphene is much smaller than that associated with the lattice, the temperature reached by the system when electrons and phonons are in thermal equilibrium is much smaller than (c) Temperature dependence of the electronic heat capacity for one of the graphene (blue curve, see Appendix) and gold (red curve, taken from Ref. [91]) nanodisks considered in (a). (d) Illustrative example of the femtosecond dynamics of the electronic thermal energy in two graphene nanodisks under the conditions of (a) for a separation d = 1 nm, with initial temperatures T1 = 1000 K and T2 = 300 K. The electronic thermal energy is shown for both the initially hot (orange curve) and cold (cyan curve) nanodisks, as well as their sum (black curve).
the electron temperatures here considered after optical pumping. For this reason, we neglect the lattice in our analysis.
As an illustrative example, we show in Fig. 2(d) the evolution of Q according to Eqs. (2) for the two graphene disks considered in Fig. 2(a) when they are prepared at initial temperatures T 1 = 1000 K and T 1 = 300 K: the cold disk more than doubles its electronic energy after ∼ 200 fs of evolution (peak of cyan curve), when it has gained nearly the same amount of energy as the one dissipated to the atomic lattice (decay of black curve). Notably, the disks reach mutual thermal equilibrium after only ∼ 250 fs, well before full relaxation takes place.
A more detailed study of the heat transfer dynamics is presented in Fig. 3 for 20 nm graphene disks separated a distance of 1 nm and doped to a Fermi energy of 0.2 eV. The color plot of Fig. 3(a) shows the HTP as a function of the temperatures in the two disks. tial conditions at the plot axes (i.e., with one of the disks at 300 K and the other one at higher temperature). The evolution is along the direction of the arrows, with positions at specific times indicated by the dashed curves. Interestingly, the evolution toward the diagonal (thermal equilibrium) is characterized by a significant increase in the temperature of the colder disk (∆T ∼ 400 K) within the first 100-200 fs, much faster than relaxation to the atomic lattice. This evolution involves the transfer of a large fraction of electronic heat to the colder disk, as shown in Fig. 3(b): when the disks are prepared at 1000 K and 300 K initial temperatures, nearly 50% of the electronic heat of the hot disk is transferred to the cold one within the first ∼ 200 fs. We remark that fast transfers take place over a wide temperature range down to substantially smaller T 's [see Fig. 3(a)]. These conclusions are maintained when considering larger disks (40 nm diameter) or wider separations (3 nm), as shown in Fig.  18. They are also maintained when considering higher doping levels (Fig. 19), well above the dipole plasmon energy, a condition for which nonlocal effects are particularly negligible. These supplementary figures also show that the results are robust with respect to variations in the disk diameters (e.g., similar conclusions are obtained for two dissimilar disks with diameters differing by a few nanometers).
In practical implementations, optical pumping with femtosecond laser pulses grants us access into the ultrafast regime, allowing us to reach high electron temperatures such as those considered in this work [94][95][96]. Additionally, the amount of optically absorbed energy depends on the pump frequency relative to the plasmons of the system [97]. This idea can be exploited to pump neighboring graphene disks in such a way that one of them absorbs much more energy than the other, just by tuning the pump laser near the plasmon of one of the disks and away from the plasmons of the other disk. We thus need disks of either different diameters or different Fermi levels. We consider the latter possibility, which can be realized in practice through the variation in intrinsic doping produced by an asymmetric dielectric environment, or also by creating different potential landscapes through an asymmetric doping geometry. The system under investigation is depicted in the inset of Fig. 4(a): two 20 nm graphene disks, separated by 1 nm, initially placed at 300 K, and doped to Fermi energies 0.2 eV and 0.3 eV, respectively. We consider optical pumping at a photon energy of 0.17 eV with a fluence of 150 mJ m −2 . The pulse energy is closer to the lower doping disk [ Fig. 4(a)], and thus, this is the one that reaches a higher temperature. For simplicity, we assume instantaneous pumping (i.e., a δ-function temporal profile of the pulse), which rapidly elevates the electron temperatures to T 1 ∼ 1200 K and T 2 ∼ 500 K [ Fig. 4(b), left end]. Interestingly, although the plasmons in the two disks are off-resonance before irradiation, optical pumping produces a larger blue shift in the hotter disk, bringing it on resonance with the initially bluer plasmon of the colder disk. Ultrafast radiative heat transfer is again observed, leading to mutual equilibrium between the disks (T 1 ≈ T 2 ) within ∼ 500 fs, which is accompanied by nearly 60% of the electronic heat of disk 1 being transferred to disk 2. We remark that higher that 50% transferred energy fraction is made possible by the doping asymmetry, which directly affects the heat capacity (see Appendix). An interesting question for future studies relates to the maximum energy fraction that can be transferred in optimized structures.

III. CONCLUDING REMARKS
Our prediction of ultrafast radiative heat transfer in graphene provides a fundamentally unique scenario: radiative coupling is capable of evacuating electronic heat from a nanoisland to a surrounding structure fast enough to prevent substantial relaxation into the atomic lattice. This is accomplished with attainable geometrical and material parameters: tens of nanometers in lateral size D in structures that can be patterned through state-of-the-art lithography [77,98] and bottom-up synthesis [99][100][101]; vertical separations of a few nanometers, as provided by van der Waals atomic layer spacers [102][103][104]; tenths of electronvolts Fermi energy E F , controllable through electrical gating [54,57]; and electron temperatures T of thousands of degrees reached by ultrafast optical pumping [88][89][90]105].
Although we have focused on disks for computational convenience, we expect our conclusions to be maintained for other geometries of similar lateral size because the heat transfer power is a frequency-integrated quantity that should be qualitatively independent of the actual spectral position of the plasmon modes, as long as they overlap with Wien's law [see Fig. 2(b)] and they are highly correlated with each other in the two islands. This correlation can be facilitated if the islands are nearly identical in shape and size. Actually, this is a condition that can be accomplished through lateral patterning of a stack formed by two graphene films and an atomically thin van der Waals layer spacer, using for example e-beam lithography.
In practice, the disks could have intrinsic doping due to interaction with a dielectric environment, which can change the Fermi energy by as much as ∼ 0.3 eV. Obviously, because the disks do not have electrical connectivity, their control through electrostatic gating presents a challenge. However, gating should be possible in a configuration consisting of neighboring graphene ribbons, which can be biased and exposed to distant gates. The contacts can be placed far from the ribbon region in which heat transfer takes place, while the gates can also be 100s nm away and thus should not affect the heat transfer.
Our choice of parameters leads to graphene plasmon energies [58]hω mν ∼ e E F /(−πη mν D) (as estimated from a Drude model description for the graphene conductivity, see Table 1 for values of the eigenvalue η mν associated with disk plasmons) that are commensurate with k B T (i.e., they overlap the broad spectral peak of thermal emission, see Fig. 10). As a consequence, the characteristic time interval τ RHT required to radiatively transfer a sizable fraction of the electronic heat energy is reduced to the femtosecond domain.
A simple dimensional analysis reveals that the HTP is proportional to E F /D, provided the ratios of disk diameters and temperatures, as well as d/D and the quantity E F /DT 2 , are kept constant (see also Figs. 16 and 17). The optimum temperature at which maximum transfer takes place scales as T ∝ E F /D. Additionally, we find the scaling τ RHT ∝ E F D 3 with Fermi energy and lateral size, and therefore, low doping levels and small sizes enable faster cooling. These conclusions are consistent with the detailed numerical analysis of τ RHT presented in Fig.  20.
We stress that the formalism developed in the Appendix can be readily applied to study radiative coupling assisted by fluctuations of other types of excitations  besides plasmons, such as optical phonons in 2D polar materials, whose relative characteristic transfer time deserves further analysis. Another interesting possibility consists in combining more than two structures. This could be used to accelerate the rate of heat evacuation and achieve greater control over the spatial flow of radiative heat transfer. Higher transfer rates could be also obtained through lateral shape optimization or by relying on other carbon allotropes such as carbon nanotubes. Additionally, similar fast transfers should be enabled by a wide range of existing atomic-scale materials capable of sustaining confined optical excitations [? ] (e.g., exciton polaritons in dichalcogenides). Besides the fundamental interest of this line of research, electronic cooling via radiative heat transfer constitutes a promising avenue to effectively suppress relaxation to the atomic lattice, thus preventing thermal damage in nanoscale devices.
Appendix A: Theory of radiative heat transfer We consider two structures labeled by the index = 1, 2, each of them assumed to be in internal thermal equilibrium at a temperature T . Radiative heat transfer can take place if T 1 = T 2 , mediated by electromagnetic interaction at characteristic frequencies ∼ k B T /h [50]. We further assume the corresponding light wavelengths ∼ 2πhc/k B T to be much smaller than the size of the structures. The response of the latter can be then described in the quasistatic limit through their susceptibilities χ (r, r , ω), which are defined as the induced charge density distribution at r produced by a unit potential point source oscillating with frequency ω at r . The charge density induced in the structure by a monochromatic potential φ(r) exp(−iωt) + c.c. is then given by d 3 r χ (r, r , ω)φ(r ) exp(−iωt) + c.c. Incidentally, although the emission of radiation away from the system is not accounted for within the quasistatic limit, this is a negligible contribution for the small structures under consideration, in which radiative heat transfer and relaxation to the atomic lattice occur at a much faster rate.
We express the net power received by structure 2 as the work P 2←1 done on 2 by charges fluctuating in 1 minus the work P 1←2 done on 1 by charges fluctuating in 2. It is enough to calculate the latter in detail, because the former is simply obtained by interchanging the subindices 1 and 2 in the resulting expression. We start from P 1←2 = − d 3 r j 1 (r, t) · ∇φ 2 (r, t) , which is the work exerted by the electric field −∇φ 2 (r, t) produced by fluctuations in 2, acting on the current j 1 (r, t) of 1. Here, . . . denotes the average over thermal fluctuations, the space integral extends over the entire 3D space, and the function j 1 is a distribution that vanishes outside the graphene and exhibits a singularity at the edge. Integrating the ∇ operator by parts, writing the electric potential φ 2 in terms of the charge ρ 2 via the Coulomb potential v(r, r ) (e.g., v = 1/ |r − r | in a homogeneous medium of permittivity ), and using the continuity equation ∇ · j 1 = −∂ t ρ 1 , we find P 1←2 = − d 3 rd 3 r ∂ t (ρ 1 (r, t)) v(r, r )ρ 2 (r , t) , or equivalently, where we have expressed the charges in frequency space ω and replaced ∂ t by −iω. The last line of Eq. (A1) implicitly defines a matrix notation in which r and r are used as matrix indices, while the dot indicates matrix multiplication. In this notation, ρ are column vectors, v and χ are matrices, and ρ T is the transpose of ρ .
The self-consistent charges ρ produced by the fluctuating charge ρ fl 2 are now obtained from the relations where we work in the frequency domain and use the matrix notation introduced above. We remark that ρ fl 2 (r, ω) vanishes for r outside structure 2, while χ (r, r , ω) vanishes for r or r outside . By construction, v(r, r ) only needs to be evaluated for r and r sitting at different structures. Inserting the solution of these equations into Eq. (A1), we find whereas I is the unit matrix (i.e., δ(r − r )). Now, the average over thermal fluctuations can be carried out using the fluctuation-dissipation theorem [79,80,106,107] ρ fl (r, ω)ρ fl (r , ω ) = −4πhδ δ(ω + ω ) [n (ω) + 1/2] Im {χ (r, r , ω)} , where n (ω) = [exp(hω/k B T ) − 1] −1 is the Bose-Einstein distribution at temperature T (i.e., for structure ). A detailed self-contained derivation of Eq. (A4) is offered in Sec. H. We find Eq. (A2) to reduce to Finally, the net power received by 2 is obtained from where the matrix ∆ [see Eq. (A3)] accounts for multiple scattering between the two structures. Incidentally, the latter cannot be ignored at short separations, as shown in Fig. 12. From the invariance of the expression in the square brackets of Eq. (A6) under exchange of the subindices 1 and 2 (see Sec. I), we confirm the expected result P 1 = −P 2 .
Finally, for structures separated by a large distance d compared to their sizes, in virtue of induced-charge neutrality (i.e., d 3 r χ (r, r , ω) = 0 for each ), the leading contribution to v is the dipole-dipole interaction. For parallel disks placed in vacuum, like the ones considered throughout this work, neglecting multiple scattering (i.e., taking ∆ = I), we find from Eq. (A6) is the polarizability of disk along a direction x parallel to it. An extra factor of 2 has been introduced in Eq. (A7) to account for the two equivalent orthogonal directions in the planes of the disks. The convergence of Eq. (A6) toward Eq. (A7) is illustrated by calculations presented in Fig. 15.

Appendix B: Description of graphene islands through plasmon wave functions (PWFs)
We now apply the above formalism to two parallel graphene islands placed in a homogeneous medium of permittivity and separated by a vertical distance d = |z − z | along their normal direction z. It is then convenient to use an eigenmode expansion for the response of each island [81,82]. This allows us to define a complete set of PWFs ρ j and real eigenvalues η j , where j is a mode index. More precisely, the susceptibility of the island, taking to be in the z = z plane, admits the rigorous exact expansion [81] where j runs over eigenmodes, we use the notation r = (D θ, z), θ is an in-plane coordinate vector normalized to a characteristic length of the structure D (we use the diameter for disks), and incorporates the response of the graphene through its local conductivity σ (ω). It should be noted that the latter depends on via the level of doping and the temperature (see below). The PWFs and their eigenvalues satisfy the orthogonality relation [81] For islands with the same geometrical shape (e.g., disks), the PWFs and eigenvalues are independent of size D , even if D 1 = D 2 . We can readily use Eq. (B1) to evaluate the heat transfer rate according to Eq. (A6). With some straightforward redefinitions, these equations remain the same, but now the coefficients of the matrices that they contain are labeled by eigenmode indices j instead of spatial coordinates r. More precisely, χ becomes a diagonal matrix of coefficients while the matrix elements of the Coulomb interaction reduce to when the operators to the left and right of v are referred to islands and , respectively. Incidentally, in this work we focus on disk dimers that share the same axis of symmetry; an eventual lateral displacement b between the islands is however easy to implement by adding it to D θ − D θ in the above expression. In this PWF formalism, inserting Eq. (B1) into Eq. (A8), we find that the polarizability of a graphene island along a given in-plane symmetry direction x is given by where ζ j = θ x d 2 θρ j ( θ) is a normalized plasmon dipole moment.

Appendix C: PWFs for disks
In the disk geometry, the azimuthal number m provides a natural way of classifying the PWFs. More precisely, we can label them using a double index (mν) and separate the radial and azimuthal dependences as We insist that these PWFs are the same for both disks in a dimer, as they are independent of disk size, and therefore, we drop the disk index for them. We also note that the PWFs are doubly degenerate for m > 0 (i.e., they share the same eigenvalue η mν and radial component ρ mν (θ) for both sine and cosine azimuthal dependences). We obtain the radial component ρ mν (θ) by solving the Maxwell equations numerically using the boundary-element method [108] (BEM) for a selfstanding disk of small thickness t ∼ D/100 compared with its diameter D. The disk is described by a dielectric function = 1 + 4πiσ/ωt, where σ is the Drude graphene conductivity (the actual model used for σ is irrelevant, as the PWFs depend only on geometry and not on the specifics of the material). In the limit of small damping, the plasmons emerge as sharp, spectrally-isolated features in the local density of optical states (LDOS) [109]. We average the LDOS over a set of off-center locations in order to access different m's efficiently. The radial components of the PWFs are then retrieved from the induced charge density, while the eigenvalues are derived from the resonance condition η mν = Re{iσ/ωD} at the corresponding LDOS peak maximum.   (C2). The values of m, ν, and ν cover the ranges considered in Fig. 5 and Table 1. All diagonal entries (ν = ν ) are 1 by construction. We only show ν ≥ ν values because the results are invariant under exchange of these two indices.
Our calculated radial PWFs, already normalized according to Eq. (C2), are shown in Fig. 5 for the lowest values of (mν), while their associated eigenvalues are given in Table 1. The orthogonality for ν = ν is rather satisfactory, as illustrated in Table 2, which shows the values obtained by numerically evaluating the left-hand side of Eq. (C2).
Upon insertion of the disk PWFs in Eq. (B4), we find that v jj is diagonal by blocks (two blocks per m, corresponding to the two different azimuthal symmetries of Eqs. (C1) and each of them contributing the same to the HTP). As χ ,jj is diagonal, this allows us to write P 2 as a sum over m's, essentially reflecting the fact that only modes of the same symmetry undergo mutual Coulomb interaction. The integrand of Eq. (A6) then becomes an analytical function (see expressions for n , χ , and ∆ above), except for the integral over radial wave functions in v jj , for which we derive a computationally convenient expression in Sec. J. We finally write Eq. (1) for the HTP, where the explicit dependence of the involved matrices on m is indicated.
Only m = 1 PWFs exhibit nonzero dipole moments ζ ν contributing to the polarizability α in Eq. (B5). More precisely, ζ ν is 0.84, 0.40, 0.11, and 0.08 for ν = 1 − 4, respectively. We use these coefficients and Eq. At zero temperature, the Fermi energy E F describes a charge-carrier doping density n subject to the relation [110] where v F ≈ 10 6 m s −1 is the Fermi velocity. This expression assumes a conical electronic band structure, which provides an accurate description for electron energies E up to a couple of electronvolts away from the Dirac point [111]. For concreteness, we consider doping with electrons, as exactly the same results are obtained when doping with holes within the conical band approximation. At finite temperature T , the population of electronic states is given by the Fermi-Dirac distribution where µ is the chemical potential. The latter depends on temperature in such a way that the electron density is maintained constant. Here, A is the graphene area, the factor of 4 originates in valley and spin degeneracies, k is the parallel wave vector, E =hv F k > 0 is the electron energy in the upper Dirac cone, f T (E) is the electron population in that cone, and 1 − f T (−E) is the hole distribution in the lower cone. Recasting the sum over k into an integral (i.e., k → (A/2π) ∞ 0 k dk ), and defining x =hv F k /k B T , Eq. (D1) becomes Direct numerical integration of Eq. (D2) allows us to obtain E F /k B T as a function of µ/k B T . The result is plotted as a pink solid curve in Fig. 6. Additionally, the large and small asymptotic T limits of Eq. (D2) (see pink labels in Fig. 6) suggest the following approximate relation which is in excellent agreement with the full solution of Eq. (D2) (cf. pink-solid and dashed-orange curves in Fig. 6). Also note that approximate [112][113][114] and asymptotic [115,116] values for the Drude weight have been proposed to work well in different limits, although they lack the universal accuracy of Eq. (D3).

Appendix E: Electronic heat capacity of graphene
The heat capacity is needed to relate the electronic thermal energy Q to the electronic temperature T . By analogy to Eq. (D1), the surface density of electronic thermal energy can be calculated as where the step functions arise when subtracting the energy at T = 0 because f T =0 (E) = θ(E F − E). After some straightforward algebra, we find where the thermal coefficient explicitly depends on µ/k B T , which is in turn a function of E F /k B T [see Eq. (D2)], so we find that β is only a function of E F /k B T . Numerical evaluation of Eq. (E3) yields the results shown in Fig. 7. For E F k B T , we have β = (4/π) ∞ 0 θ 2 dθ/(1 + e θ ) ≈ 2.2958. (Incidentally, we correct this parameter here for a factor of 2 that was missing in Ref. [58].) We note that the graphene heat capacity has been widely used in previous studies [68,113,114] in the so-called degenerate limit (k B T µ).

Appendix F: Graphene conductivity
We adopt the local-RPA model for the graphene conductivity [58,83,84] where is a temperature-dependent effective Drude weight that accounts for intraband transitions and has been the object of a recent theoretical and experimental study [90]. The integral term in Eq. (F1) represents the contribution from interband transitions. Besides the explicit dependence on temperature T , we note that there is an additional dependence through the chemical potential µ. We plot the resulting µ D in Fig. 6 (red-solid curve). A reasonable approximation to this parameter is obtained by substituting E F for µ in Eq. (F2) (dashed-blue curve in Fig. 6). We assume a rather conservative value for the energy broadeninghτ −1 = 10 meV throughout this work (this corresponds to a Drude-model mobility [117] . For simplicity, we further neglect the dependence of τ on temperature and chemical potential, which could be readily incorporated following previous studies [113][114][115]. This dependence is partially absorbed in the assumed value of τ over the significant range of temperatures under consideration, although a more detailed analysis could reveal unexpected effects outside that range.

Appendix G: Time evolution
The temporal evolution of the electronic temperature is given by Eqs. (2), which we solve numerically by using a 4 th order Runge-Kutta method. It is instructive to rewrite them with the temperatures appearing in a more explicit   form. Using the Q dependence on T given by Eq. (E2), we find where C(T ) = 3 + (T /β)(dβ/dT ) is a dimensionless coefficient that varies between 3 and 4 in the large and small T limits, respectively (see β dependence on T in Fig. 7).
In the simulations of Figs. 2(d), 3, 18, and 19 we fix the initial temperatures T to prescribed values. However, in the calculation of Fig. 4 the initial temperatures are determined by the energy absorbed from a light pulse via the absorption cross-section given by Eq. (C3). Assuming a δ-function pulse of frequency ω 0 and fluence F 0 , we have Q (t = 0) = σ abs (ω 0 )F 0 . The initial temperature is then obtained by entering this value of Q in Eq. (E2). Absorption crosssection of individual (blue) and closely spaced (red) graphene disks calculated using either classical (local-RPA conductivity, solid curves) or quantum-mechanical (tight-binding combined with full RPA, as described elsewhere [87], broken curves) models. We assume a Fermi energy of 2 eV and a damping of 0.05 eV.
(H9) This is the FDT used in Eq. (6), where we drop the 'sym' subscript for clarity.
Appendix J: Computation of v jj for coaxial disks In this section, we provide a computationally efficient expression to calculate the Coulomb interaction matrix elements v jj [Eq. (14)] for coaxial disks [i.e., with the plasmon wave functions (PWFs) of Eqs. (16)]. We start by rewriting the Coulomb potential as [118] 1 where Y L are spherical harmonics, r < = min{r, r }, and r > = max{r, r }. Specifying Eq. (14) for two PWFs ρ κ mν ( θ) and ρ κ m ν ( θ ) and using Eq. (J1), we can perform the azimuthal integrals of θ and θ analytically by choosing the spatial origin at a point along the axis of revolution symmetry in between the two disks. Upon detailed examination, we find v jj to be zero unless m = m and κ = κ . Therefore, PWFs of different azimuthal symmetry do not interact. It should be also noted that only κ = c contributes to m = 0. The remaining nonzero elements are independent of κ, but they depend on m, ν, and ν as where we take r = D 1 θ and r = D 2 θ . Additionally, the spherical harmonics in this expression are evaluated at zero azimuthal angle, while the polar angles are θ 1 = tan −1 (D 1 θ/d 1 ) and θ 2 = π − tan −1 (D 2 θ /d 2 ), where d 1 and d 2 are the distances from the disks to the origin (i.e., d 1 + d 2 = d), a convenient choice being d 1 = d and d 2 = 0, so that (r < /r > ) l goes rapidly down for large l, particularly at large separations. Equation (J2) gives the (νν ) elements of the matrix v m entering Eq. (1). This expression is also useful to normalize the PWFs via Eq. (13), whose integral corresponds to v m νν with = 1, D 1 = D 2 , and d = 0.
the reflection for TE polarization vanishes in the quasistatic limit), while σ is the conductivity of the layer = 1, 2.
Putting these elements together, the transfer power per unit area becomes in agreement with the c → ∞ limit of the well-known expression for the transfer power between two planar structures [64]. We plot this quantity in Fig. 11 using the full nonlocal RPA (broken curves) and the local-RPA (solid curves) models for the conductivity. The agreement between these results indicates that nonlocal effects only play a marginal role in this study.