Ultrafast electronic heat dissipation through surface-to-bulk Coulomb coupling in quantum materials

PHYSICAL REVIEW B 106, 115422 (2022) Ultrafast electronic heat dissipation through surface-to-bulk Coulomb coupling in quantum materials Alessandro Principi1,*and Klaas-Jan Tielrooij 2 1Department of Physics and Astronomy, University of Manchester, Manchester M13 9PL, United Kingdom 2Catalan Institute of Nanoscience and Nanotechnology (ICN2), BIST & CSIC, Campus UAB, 08193 Bellaterra, Barcelona, Spain (Received 18 June 2022; revised 13 August 2022; accepted 2 September 2022; published 16 September 2022) The timescale of electronic cooling is an important parameter controlling the performance of devices based on quantum materials for optoelectronic, thermoelectric, and thermal management applications. In most conventional materials, cooling proceeds via the emission of phonons, a process that can bottleneck the carrier relaxation dynamics, thus degrading the device performance. Here we present the theory of near-field radiative heat transfer that occurs when a two-dimensional electron system is coupled via the nonretarded Coulomb interaction to a three-dimensional bulk that can behave as a very efficient electronic heat sink. We apply our theory to study the cooling dynamics of surface states of three-dimensional topological insulators and of graphene in proximity to small-gap bulk materials. The “Coulomb cooling” we introduce is alternative to the conventional phonon-mediated cooling, can be very efficient, and can dominate the cooling dynamics under certain circumstances. We show that this cooling mechanism can lead to a sub-picosecond timescale, significantly faster than the cooling dynamics normally observed in Dirac materials. DOI: 10.1103/PhysRevB.106.115422 I. INTRODUCTION Two-dimensional (2D) Dirac materials [1–4] have been extensively studied in the past two decades for the diverse range of intriguing properties they harbor, which could in turn enable a wealth of novel practical applications [5–8]. This is the case of graphene [1,9–11], arguably the most-studied Dirac material, which continues to attract significant interest because of its potential for electronic and optoelectronic applications [11]. However, graphene is just an instance of a broad family of systems which includes the 2D surface states of three-dimensional (3D) topological insulator [2,3,12]. Their inverted bulk band structure allows states to localize at the surface [5,13–16]. In ideal topological insulators, surface states do not hybridize with bulk ones, and are topologically protected against any perturbation that preserves the symmetries of the bulk [2,3,12]. Thanks to the coupling between kinetic momentum and spin, the electronic surface states of a 3D topological insulator present a physics potentially richer than graphene, and promise application to a diverse range of fields, including spintronics [6,7], optoelectronics, and photonics [17,18]. In a recent experimental work [19], the cooling dynamics of surface electrons of bismuth and antimony chalcogenides was studied with pump-probe techniques. Surprisingly, electronic heat relaxation faster than that of bulk carriers was observed. The observed bulk and surface heat-decay rates differ by about an order of magnitude, while the environment experienced by their electrons is the same. Interestingly, a decay rate of a few hundred femtoseconds was obtained under rather strong *[email protected] photoexcitation with a fluence on the order of 100 μJ/cm2.In comparison, when graphene is excited with a similar fluence, cooling is rather slow—several picoseconds—due to a phonon bottleneck effect [20]. So why does the cooling of the electrons of surface states of topological insulators occur so much faster than their bulk counterpart or than graphene? In graphene, the cooling of photoexcited electrons is ultimately limited by the emission of intrinsic optical and acoustic phonons of the material [21–23] or of the encapsulant or substrate [24–26]. In topological insulators, cooling is also typically ascribed to phonon emission on a picosecond timescale [27]. Such phonon emission mechanisms fail to account for the starkly different relaxation dynamics of surface and bulk electrons of topological insulators, unless one would postulate different electron-phonon couplings or a reduced phase space for phonon emission by surface states. In this paper, we explore a different mechanism that does not require fine-tuning material parameters, as it relies on electronic systems exclusively. This alternative mechanism is based on the notion that, contrary to electrons in conventional graphene devices, the surface electrons of topological insulators are in close proximity to a macroscopic bulk whose particle-hole excitations can occur at energies comparable to the (surface) thermal excitations. Thanks to the near-field radiative coupling between surface and bulk, the latter can act as a heat sink. Being macroscopic, the bulk can absorb large amounts of heat and therefore efficiently cool down the surface electrons. In this picture, heat is dissipated into particle-hole excitations of the bulk via nonretarded Coulomb interactions (see Fig. 1). Thus we term this mechanism “Coulomb cooling.” 2469-9950/2022/106(11)/115422(11) 115422-1 ©2022 American Physical Society ALESSANDRO PRINCIPI AND KLAAS-JAN TIELROOIJ PHYSICAL REVIEW B 106, 115422 (2022) Tb µ b Δ,m Ts µs v F FIG. 1. Schematics of Coulomb cooling of hot surface Dirac fermions into bulk states. Dirac fermions have a constant Fermi velocity vF, and they exhibit a hot thermal distribution characterized by temperature and chemical potential Tsand μs, respectively. Similarly, bulk states are described in terms of the effective mass mand band gap 2. The Fermi distribution of bulk states has chemical potential μb(set to zero in later calculations) and temperature Tb<Ts. We note that for atomic-scale heterogeneous systems with massless electrons, such as graphene interfaced with other 2D materials, e.g., hexagonal boron nitride or WS2,or3D dielectric substrates, e.g., SiO2, different relaxation processes have been described. This includes cooling of hot electrons in graphene via polar [24] and hyperbolic [25,26] substrate phonons, as well as thermionic emission [28] and direct charge transfer [29]. In bare monolayer graphene suspended over a hole, hot carrier cooling occurs via optical phonon emission, followed by coupling to acoustic phonons [23]. These mechanisms thus involve either electron-phonon coupling or charge transfer. The mechanism we describe here for 2D-3D surface-bulk TIs and graphene-TI systems does not involve electron-phonon coupling nor charge transfer. Rather it relies purely on Coulomb interaction, which so far has only been explored in multilayer graphene [30]. In the remainder of the paper, we consider a system of 2D massless Dirac fermions in proximity to a 3D gapped bulk. The two are coupled only electrostatically via nonretarded Coulomb interactions. This minimal model of a topological insulator equivalently describes a graphene sheet in proximity to a small-gap material. For both cases, using the kinetic equation for surface electrons, we derive an expression for their cooling rate. This is controlled by the convolution of surface and bulk particle-hole excitation spectra. The higher their overlap, the larger the amount of heat transferred per unit time. We show that the cooling rate reaches a maximum for surface temperatures close to half the bulk-band gap. This is interpreted as a resonance between surface electronic transitions, whose typical energy is the thermal one, and interband particle-hole bulk excitations across the band gap. We find a timescale of a few hundred femtoseconds for topological insulators, and even a few tens of femtoseconds for graphene on a small-bandgap semiconductor. II. THE MODEL We model both the surface states of the topological insulator and the electrons in graphene as a gas of massless Dirac fermions [1–3,9–12], i.e., Hs=¯hvF k,α,β ˆ ψ† k,α,sk·σα,β ˆ ψk,β,s,(1) where vFis the Fermi velocity and σis a vector of Pauli matrices. These act in the real spin space for topological insulators, and in pseudospin (sublattice) space for graphene. In Eq. (1), ˆ ψ† k,α,s(ˆ ψk,α,s) creates (destroys) a surface particle of two-dimensional wave vector k(momentum ¯hk) and (pseudo)spin projection α. The band energy is ε(s) k,λ =λ¯hvF|k|, where λ=+(λ=−) denotes the surface conduction (valence) band. [The Hamiltonian (1) is obtained from the usual one for topological insulators, which features a cross product between momentum and spin operators [2], by redefining the spin quantization axes.] Equation (1) is in form identical for graphene and topological insulators. However, their electrons exhibit different degeneracies Nf. While the surface states of topological insulators are helical [2], thus yielding Nf=1, in graphene they feature a full spin-valley degeneracy [1]. Thus, for graphene, the number of fermion flavors is Nf=4. We model the bulk as a two-band system confined into the half-space z>0 and described by the Hamiltonian Hb= k,η ε(b) k,η ˆ ψ† k,η,bˆ ψk,η,b.(2) Here, ε(b) k,η =η[+¯h2|k|2/(2m)], where 2is the bulk-band gap and mis the bulk-band mass (assumed to be the same for both valence and conduction bands), while ˆ ψ† k,η,b(ˆ ψk,η,b) creates (destroys) a bulk particle of three-dimensional wave vector k=(k,kz), in band η. Here, kand kzare the wave vectors parallel and perpendicular to the surface z=0, η=+ denotes the conduction band, while η=−stands for the valence band. For simplicity, we assume that bulk bands are spherically symmetric and have no spin structure. These approximations do not affect the temperature dynamics on a qualitative level. Assuming specular reflection at the interface z=0, the bulk eigenstates acquire the form of standing waves: k,η,b(r,z)=√2/Ve ik·rsin(kzz), where ris a vector along the surface. We assume that bulk and surface electrons are coupled electrostatically by long-range instantaneous Coulomb interactions [31], whose Hamiltonian is Hsb =1 2∞ 0 dz q Vsb(q,z)ˆnq,sˆn−q,b(z).(3) Here, qis a two-dimensional wave vector along the surface of the topological insulator, while ˆnq,sand ˆn−q,b(z)arethe2DFourier transforms of the surface and bulk density operators. The precise form of the interaction Vsb(q,z) is determined by the solution of the associated Poisson (electrostatic) problem, as we proceed to show. 115422-2 ULTRAFAST ELECTRONIC HEAT DISSIPATION THROUGH … PHYSICAL REVIEW B 106, 115422 (2022) The interaction between surface and bulk electrons In this section we derive the Coulomb interaction between surface and bulk electrons. To do so, we consider the electrostatic problem of a single charge added to a conducting 2D sheet (the surface states located at z=0) placed on top of the topological-insulator bulk which fills the half-space z>0. The half-space z<0 is instead empty. The single charge is added as a plane wave of wave vector qin the surface sheet. In response to this added charge, induced bulk and surface densities are generated. These are named nb(r,z) and ns(r), respectively. The resulting Poisson equation (in Gaussian units) is ∇·[(z)∇φ(r,z)] =−4πe2n(r,z),(4) where n(r,z)=nb(r,z)(z)+ns(r)δ(z)+e−iq·rδ(z)isthe total electron density, and (z)=b(z)+vac(−z)isthe dielectric function. Here, bis the relative dielectric constant of the topological insulator, which accounts for the screening due to filled bands, while vac =1[(z) is the Heaviside step function]. The Fourier transform of Eq. (4) in the direction parallel to the topological-insulator surface (or the graphene sheet) yields ∂z[(z)∂zφq(z)] −q2(z)φq(z)=−4πe2nq(z),(5) where nq(z)=nq,b(z)(z)+nq,sδ(z)+δ(z). Within linear response, we write nq,s=χs(q,ω)φq(0), where χs(q,ω)is the density-density response function (polarizability) of the surface states, and nq,b(z)=∞ 0 dzχb(q,ω,z,z)φq(z)≃¯χb(q,ω)φq(z).(6) Here, χb(q,ω,z,z)≃¯χb(q,ω)δ(z−z) is the densitydensity response function [32] (polarizability) of bulk states treated within a semilocal approximation. This approximation allows us to reduce the integro-differential equation to a set of differential ones. Behind it, there is the assumption that the potential and induced charge density are sufficiently localized in the zdirection. The result in Eq. (9) is consistent with our assumption and justifies a posteriori its use. We do not need to introduce any approximation in the direction parallel to the surface since we can simply take a Fourier transform and reduce equations to algebraic ones. However, even doing so we find [see below Eq. (10)] that the potential approaches a constant at small q. This means that the potential decays fast at large enough distances. Thus, the potential and induced density are localized in space. Equations of this section contain the response of the doped bulk, embodied by χb(q,ω,z,z) and ¯χb(q,ω), for completeness. However, we leave this function unspecified since, as we show below, it plays no role in the undoped regime when the bulk thermal energy (kBTb) is much smaller than the bulk-band gap (2). Using these expressions, the Poisson Eq. (5) can be split into the two half-spaces as b∂2 z−q2φq(z)=−4πe2¯χb(q,ω)φq(z),if z>0, ∂2 z−q2φq(z)=0,if z<0.(7) The boundary conditions, at the surface z=0 and at z→ ±∞,areφq(0+)=φq(0−), φq(z→±∞)=0, and b∂zφq(z)|z→0+−∂zφq(z)|z→0−=−4πe2(nq,s+1).(8) Introducing the bulk Thomas-Fermi wave vector [32] q2 TF(q,ω)=−4πe2¯χb(q,ω)/b, the solution of Eqs. (7) with the boundary conditions above is φq(z)=¯ φqe−√q2+q2 TF (q,ω)z,if z>0, φq(z)=¯ φqeqz,if z<0.(9) Equation (9) shows that both the potential and the induced density (which is linearly related to the potential) are very localized in the zdirection. This result is consistent with our initial semilocal assumption and serves as an a posteriori justification of its use. Using Eq. (8), we get ¯ φq=4πe2 q2+q2 TF(q,ω)+q−4πe2χs(q,ω) .(10) In the limit of zero frequency and small wave vectors, χs(q,ω) becomes the densities of states of surface electrons. Equation (9) is also applicable to the situation in which a graphene layer is separated by a layer of vacuum (or air) of thickness dfrom the topological insulator surface. The interaction between electrons in the two systems is therefore weighted with the exponential e−qd . III. THE COOLING RATE—GENERAL THEORY We consider the kinetic equation for electrons at the surface of the topological insulator, which interact with those in the bulk via Coulomb interactions. The surface and bulk electrons are described by the Fermi distribution functions f(s) k,λ and f(b) ˜ k,η at the temperatures Tsand Tband chemical potentials μsand μb, respectively. As in Sec. II,kand ˜ kare their twoand threedimensional wave vectors, respectively. We recall that λ=±1 is used to denote the two surface bands (together forming a Dirac cone), while η=±1 is used for the bulk conduction and valence bands. For the calculation of the cooling time, we assume that no driving field is present and that material parameters (including temperatures and chemical potentials) are isotropic. However, since the two populations are at different uniform chemical potentials and temperatures, we can study the time evolution of their distribution functions. The kinetic equation satisfied by f(s) k,λ is [25] ∂tf(s) k,λ =−I(sb) k,λ ,(11) where I(sb) k,λ is the electron-electron collision integral between surface and bulk electrons [resulting from the interaction in Eq. (3)]. This conserves their numbers separately but allows for the exchange of energy between them. We will specify the collision integral in the following subsection. First, however, we will derive the general expression for the cooling rate. To obtain the cooling rate, we multiply Eq. (11)bythe energy of the surface state, ε(s) k,λ, and sum over all wave vectors 115422-3 ALESSANDRO PRINCIPI AND KLAAS-JAN TIELROOIJ PHYSICAL REVIEW B 106, 115422 (2022) kand all values of the surface-band index λ. The left-hand side of the so-obtained equation yields the time derivative of the energy stored in surface states, ∂tEs. The latter is rewritten as ∂tEs=Cs∂tTs, where the heat capacity of surface states is defined as [25] Cs= λd2k (2π)2ξ(s) k,λ−∂f(s) k,λ ∂ξ(s) k,λ ξ(s) k,λ Ts+∂μs ∂Ts.(12) Here, ξ(s) k,λ =ε(s) k,λ −μs. The derivative of the surface chemical potential μs≡μs(Ts) with respect to temperature is obtained by imposing the conservation of the surface electron density [25]. Defining the power dissipated into bulk states as Q=k,λ I(sb) k,λ ε(s) k,λ ≡γCs(Ts−Tb), Eq. (11) yields the cooling rate [25] γ=Q (Ts−Tb)Cs .(13) In the following we first define the collision integral and then calculate Q, and thus γ. IV. THE COOLING RATE—SURFACE-TO-BULK POWER DISSIPATION To calculate the power dissipated into bulk states, we first have to obtain the collision integral due to the interaction of Eq. (3). Within the Fermi-golden-rule approximation, the collision integral on the right-hand side of Eq. (11) reads I(sb) k,λ =22π ¯hA2Lz k,˜ k,q λ η,η∞ −∞ dωV2 q,k,˜ kFk,λ;k+q,λδε(s) k,λ −ε(s) k+q,λ+ωδε(b) k,η −ε(b) ˜ k,η−ωδ(k −˜ k −q) ×f(s) k,λ f(b) k,η1−f(s) k+q,λ1−f(b) ˜ k,η−1−f(s) k,λ1−f(b) k,ηf(s) k+q,λf(b) ˜ k,η,(14) where the factor 2 upfront accounts for the spin degeneracy of bulk states, ωand qare the transferred energy and momentum parallel to the surface, k and ˜ k are the components of the three-dimensional momenta kand ˜ kparallel to the surface, Ais the surface area of the topological insulator, and Lzits extension in the third dimension. In this equation, Fk,λ;k+q,λis the squared matrix element of the surface-electron density operator between incoming and outgoing scattering states. Since bulk states are standing waves, the zcomponents of the threedimensional momenta kand ˜ kare taken to be positive. The matrix element of the screened Coulomb interaction Vq,k,˜ k is obtained by integrating it over the incoming and outgoing scattering states, i.e., V2 q,k,˜ k=4∞ 0 dzsin(k zz)sin(˜ k zz)e−√q2+q2 TF (q,ω)z 2 |¯ φq|2. (15) This matrix element is manipulated in Appendix Ato give V2 q,k,˜ k≃|¯ φq|2π2 4q2+q2 TF(q,ω) δ(k z−˜ k z) ≡V2 qδ(k z−˜ k z).(16) To obtain this expression, we have evaluated the integral on the right-hand side of Eq. (15) and approximated it in the limit of small qand qTF(q,ω) (the latter is justified when the bulk is undoped and kBTb2). In this case, the weight of the integral is located around the lines ˜ k z=±k z. Thus, we replaced Lorenzians of width √q2+q2 TF(q,ω) with δfunctions. Finally, we used that ˜ k zand k zmust be taken to be positive. Putting Eq. (16) back into Eq. (14), after some lengthy but straightforward algebra we get Q=4 π¯hd2q (2π)2V2 q∞ 0 dωω[n(s)(ω)−n(b)(ω)] ×mχb(q,ω)mχs(q,ω),(17) where n(s/b)(ω)=[e¯hω/(kBTs/b)−1]−1.InEq.(17), mχs(q,ω)=−π A k,λ,λf(s) k,λ −f(s) k+q,λFk,λ;k+q,λ ×δε(s) k,λ −ε(s) k+q,λ+ω(18) and mχb(q,ω)=− π ALz k,η,ηf(b) k,η −f(b) k−q,η ×δε(b) k,η −ε(b) k−q,η−ω(19) are the imaginary parts of the density-density response functions of surface [33–35] and bulk [32] states, respectively. The function mχb(q,ω), not to be confused with the function ¯χb(q,ω) introduced in Sec. II A, is calculated in Appendix B. The cooling rate γis obtained by inserting Qin Eq. (17) back into the definition (13), and reads γ=− 4 π¯hCsd2q (2π)2V2 q∞ 0 dωωn(s)(ω)−n(b)(ω) Ts−Tb ×mχb(q,ω)mχs(q,ω).(20) This equation is the key result of this paper. We stress that Eq. (20) describes the exchange of heat from surface to bulk electrons and vice versa. This is at the origin of the difference in Bose factors in Eq. (20). 115422-4 ULTRAFAST ELECTRONIC HEAT DISSIPATION THROUGH … PHYSICAL REVIEW B 106, 115422 (2022) V. R E S U LT S We now show our results for the cooling time of electrons in the surface state of a topological insulator and in graphene in proximity to a small-gap material. We find that the bulk behaves as an efficient heat sink for electrons, yielding subpicosecond cooling times. We start with the surface states of the topological insulator. We numerically evaluate Eq. (13), with the heat capacity of surface states and power lost to bulk states given by Eqs. (12) and (17), respectively. For the sake of definiteness, in the numerical calculations, we set the surface Fermi velocity [36] vF=0.5×106m/s, while the bulk is left undoped, i.e., its electron density is nb=0 (which translates into μb=0). The bulk temperature is set to Tb=300 K, which is much smaller than the bulk-band gap [37]2=200 meV. Under these conditions, the bulk bands are nearly unpopulated and we can thus safely take qTF(q,ω)=0inEqs.(10) and (15). Finally, the bulk electron mass is set to [38,39]m=0.21me, where me=9.1×10−31 kg is the bare electron mass, while the dielectric constant of the (undoped) topological insulator is takentobe[40]b=10. We would like to point out that the topological insulators that we describe have no charge carriers in the bulk, at low temperature, which corresponds to a Fermi energy inside the gap. Most common growth methods of binary topological insulators, such as Bi2Se3and Sn2Te3result in a Fermi energy either in the valance or conduction band. However, there are well-known methods to produce topological insulators with Fermi energy in the gap, for example by using appropriate substrate materials [41] or using mixtures with controlled stoichiometry such as Bi1.4Sb0.6Te1.51Se1.49 (BSTS) [19]. The effect of bulk thermal excitations is included in our theory via the bulk response function. Furthermore, since the typical energies of excitations here are of the order of kBTs∼10–100 meV (in particular, they can be interband transitions), the dielectric constant should be taken to be the one at intermediate frequencies, not the zero-frequency one. In Fig. 2we show the cooling rate for surface-to-bulk Coulomb cooling in topological insulators (with Fermi energy in the bulk band gap) calculated from Eq. (13). In Fig. 2(a) we present our numerical results for three different values of the surface electronic density and as a function of temperature Ts. The three chosen densities as ns=0cm −2, i.e., an undoped system with Fermi energy at the surface Dirac crossing, ns= 5×1011 cm−2and ns=1012 cm−2. the latter corresponds to a Fermi energy close to the bottom of the bulk conduction band. We find that the cooling rate depends only weakly on surface carrier concentration, while it depends quite strongly on temperature. In particular, it decreases rapidly at low temperatures (i.e., for the surface temperature approaching the bulk one). Curiously, we find that the cooling rate exhibits a maximum at a temperature approximately equal to half the bulk-band gap. This behavior can be understood by noting that surface electrons can be very hot. As shown in Fig. 2(a), their temperature can be of thousands of Kelvins. At these temperatures their Fermi distribution is very broadened. Thus, they can undergo both intraand inter-band transitions (where the bands (b) FIG. 2. (a) The cooling rate for Dirac-like surface states (vF=0.5×106m/s) in proximity to an undoped 3D bulk (nb=0, which translates into μb=0) kept at room temperature (Tb=300 K, much smaller than the gap energy: kBTb2), calculated from Eq. (13) and plotted as a function of temperature. Curves exhibit a maximum at a surface electron temperature corresponding to half the bulk-band gap, =100 meV (dashed line). (b) The cooling rate is very weakly dependent on surface carrier density. The parameter used in these calculations are given at the beginning of Sec. Vand are recalled here for convenience. The bulk electron mass is set to [38,39]m=0.21me,whereme=9.1×10−31 kg is the bare electron mass, while the dielectric constant of the (undoped) topological insulator is taken to be [40]b=10. involved here are the surface ones). At the same time bulk electrons are kept at room temperature. Thus, they experience both intraband excitations within the valence and conduction bands and interband ones from the valence to the conduction band. For a surface-electron temperature of the order of the bulk-band gap, excitations of surface states occurring with energies of the order of kBTsbecome resonant with the bulk interband excitations. The weight of such excitations grows sharply for energies larger than 2. The latter fact is seen in the plots of the bulk density-density response (the bulk absorption spectrum) given in Appendix B. However, the cooling rate must fall off rapidly at low temperatures, because interband processes become heavily suppressed. This effect is clearly visible in Fig. 2(a). We stress that, for different values of the bulk-band gap, the cooling rate is expected to follow similar functional form as shown in Fig. 2(a). The main difference being that the peak 115422-5 ALESSANDRO PRINCIPI AND KLAAS-JAN TIELROOIJ PHYSICAL REVIEW B 106, 115422 (2022) (b) FIG. 3. The cooling rate for Dirac-like states in graphene (vF=106m/s) in proximity to an undoped 3D bulk identical to that of Fig. 2, calculated from Eq. (13). (a) As a function of temperature, it exhibits a maximum at a surface electron temperature corresponding to half the bulk-band gap (dashed line). (b) The cooling rate is very weakly dependent on surface carrier density. The parameters used in these plots are given in Sec. Vand Fig. 2. occurs at different temperatures of the order of the bulk-band gap. Thus for smaller (larger) gaps, the peak temperature will be lower (higher). In Fig. 2(b) we show the cooling rate as a function of the surface density nsand for three values of the temperature, Ts=500,1000, and 1500 K. We see that the numerical results are practically independent of carrier density at the highest temperatures, and only weakly dependent at the lowest one. This can be understood as a consequence of the complete smearing of the surface occupation function at temperatures much larger than the Fermi energy. Only at the lowest temperatures and highest densities achievable in our model we start observing some deviation from perfect flatness. In Fig. 3we show the cooling rate for a graphene sheet in direct contact with a small-gap bulk material as a function of temperature [in Fig. 3(a)] and carrier density [in panel (b)]. For the bulk material, we use the same 3D topological insulator used above, i.e., the 3D parameters are taken to be the same. However, for the “surface” states (now a graphene sheet), the number of fermion flavors is set to [1]Nf=4, while the Fermi velocity is doubled, i.e., we use [1]vF= 106m/s. Comparing Figs. 2and 3, we see that graphene would exhibit dynamics approximately four-five times faster than the topological-insulator surface states. This is surprising because, thanks to the doubling of the Fermi velocity and quadrupling of the number of fermion flavors, the density of states of undoped graphene and of the surface states studied above are (accidentally) identical. This in turn implies that, in the undoped limit, they also exhibit the same heat capacity. The reason for the enhanced cooling rate is instead to be found in the typical energies of particle-hole excitations, which are different in the two systems. Due to the linear energy dispersion, the typical energy exchanged by surface and bulk states during a collision is ω≃¯hvFq. This in turn implies that, for a given value of momentum ¯hq, the energy exchanged between graphene and bulk electrons is twice the energy exchanged between surface and bulk states. Thus, the cooling dynamics proceeds at a faster pace in graphene, even if the number of interactions per unit time is the same as in a topological insulator, a fact that is reflected in the larger cooling rate. If a spacer is introduced between the graphene and bulk material, then the interaction decreases exponentially as e−qd , as shown at the end of Sec. II A. The wave vector qis of the order of the graphene Fermi wave vector. Thus, for the carrier densities and (high) temperatures explored in this paper, the interaction is expected to be suppressed for graphene-bulk distances exceeding 10 nm. We therefore expect the cooling time to grow significantly for larger spacer thicknesses. VI. DISCUSSION AND CONCLUSIONS In this paper we have developed the theory of the cooling dynamics of electrons in the surface states of a topological insulator, as well as in a graphene sheet, coupled to bulk states via nonretarded Coulomb interactions. The aim has been to explore the cooling capabilities of all-electronic surface-bulk coupling in these systems. For this reason, we have employed a simplified model: we have treated the surface states (and graphene electrons) as massless Dirac fermions characterized by a Dirac-like energy dispersion. We have neglected corrections due to, e.g., trigonal warping, which are expected to be only minor. We have considered a fully gapped bulk with particle-hole symmetric parabolic-energy bands of equal masses. We have neglected all possible couplings between surface and bulk states (such as impurity or phonon-mediated hopping) with the exception of a density-density interaction of the Coulomb type. The latter conserves the number of electrons in surface and bulk states separately, but allows the exchange of energy between them. We have thus studied how such near-field radiative coupling can efficiently transfer energy between the surface and the bulk. We have thus explored an alternative mechanism to phonon-mediated cooling, which can be very efficient and could become dominant in some circumstances. In fact, we find that Coulomb cooling can out-compete other cooling mechanisms. We find cooling rates of the order of few inverse picoseconds for the surface states of topological insulators, which agrees with a recent result reported in the literature [19]. More surprisingly, we find that the cooling time of graphene in contact with a narrow-band insulating material 115422-6 ULTRAFAST ELECTRONIC HEAT DISSIPATION THROUGH … PHYSICAL REVIEW B 106, 115422 (2022) can be as short as a tenth of picosecond. This is significantly faster than graphene cooling via optical, acoustic and substrate phonons, which typically occurs on a picosecond timescale [21–26]. A true comparison between different cooling pathways would require a significantly more detailed work to describe the impact of disorder, phonons, etc., on the cooling dynamics. All these mechanism are also potentially strongly dependent on material characteristics, and therefore it would be hard to derive general and universal trends. On the contrary, our model depends on few, experimentally available parameters (e.g., the surface Fermi velocity, the bulk-band mass, the undoped-material dielectric constant). Thus, it could be employed to study universal trends in electronic cooling of surfaces of topological insulators. We further observe that, although we have not treated phonon cooling in this paper, the interaction we describe could be thought as the result of phonon emission and re-absorption by the electrons at the surface and in the bulk, respectively. Such phonon-number-conserving processes could be easily incorporated in the theory and would result in a further enhancement of cooling rates. Finally, we stress that the theory also applies, with minor modifications, to graphene in proximity to narrow-gap materials, as we have shown above. This fact, which broadens the applicability of the present theory, has also important practical implications. Graphene is in fact one of the most studied materials for optoelectronic applications. The slowing down of cooling dynamics at high powers limits however its potential, for example, to applications such as higher harmonic generation [19]. The fact that heat dissipation could be made more efficient via near-field coupling to narrow-gap materials offers a novel way to overcome the limitation intrinsic to current graphene devices. ACKNOWLEDGMENTS A.P. acknowledges support from the European Commission under the EU Horizon 2020 MSCA-RISE-2019 programme (Project No. 873028 HYDROTRONICS) and of the Leverhulme Trust under the Grant No. RPG-2019-363. K.J.T. acknowledges funding from the European Union’s Horizon 2020 research and innovation program under Grant Agreement No. 804349 (ERC StG CUHL). APPENDIX A: MANIPULATION OF THE SURFACE-BULK INTERACTION TheintegrationofEq.(15) yields V2 q,k,˜ k=¯ φq 2⎧ ⎨ ⎩q2+q2 TF(q,ω) q2+q2 TF(q,ω)+(k z−˜ k z)2−q2+q2 TF(q,ω) q2+q2 TF(q,ω)+(k z+˜ k z)2⎫ ⎬ ⎭ 2 .(A1) We now observe that, for small qand q2 TF(q,ω), the regime of interest for the processes we are describing, the two terms on the right-hand side of Eq. (A1) are sharply peaked around ˜ k z=k zand ˜ k z=−k z, respectively. Thus, we approximate them with two δfunctions, making sure that their total integral (over the variable ˜ k z) remains unchanged. We thus obtain V2 q,k,˜ k≃4πe2 (b+1)q−4πe2χs(q,ω)2π2 4q2+q2 TF(q,ω)δ(k z−˜ k z)+δ(k z+˜ k z) ≃4πe2 (b+1)q−4πe2χs(q,ω)2π2 4q2+q2 TF(q,ω) δ(k z−˜ k z).(A2) In the last line we noticed that k zand ˜ k zmust be taken as positive, since the bulk wave functions describe standing waves (negative wave vectors correspond to the same wave function). APPENDIX B: THE DENSITY-DENSITY RESPONSE OF BULK ELECTRONS The imaginary part of the density-density function in Eq. (19) is given by mχb(q,ω)=−π η,ηd3k (2π)3[f(ηεk+η −μb)−f(ηεk+η +ω−μb)]δ(ω+ηεk+η −ηεk+q−η).(B1) Here we rewrote εk,η =η +ηεk, introducing εk=k2/(2m). We also defined f(ξ)=[eξ/(kBTb)+1]−1and introduced the chemical potential of bulk bands μb. In what follows, we assume ω>0 and analyze separately the intraband (η=η) and interband (η= η) contributions to Eq. (B1). 115422-7 ALESSANDRO PRINCIPI AND KLAAS-JAN TIELROOIJ PHYSICAL REVIEW B 106, 115422 (2022) 1. Intraband term In this case, η=η. Hence, mχ(intra) b(q,ω)=−1 4π η∞ 0 dkk2π 0 dθsin θ[f(ηεk+η −μb)−f(ηεk+η +ω−μb)]δ(ω+ηεk−ηεk+q). (B2) The δfunction implies that ω+ηk2 2m−ηk2+q2+2kq cos(θ0) 2m=0⇒cos(θ0)=m kqq2 2m−ηω.(B3) Solutions exist for −kq m⩽q2 2m−ηω ⩽kq m⇒k⩾m q q2 2m−ηω≡k0.(B4) Therefore, mχ(intra) b(q,ω)=−1 4π η∞ 0 dkk2[f(ηεk+η −μb)−f(ηεk+η +ω−μb)] m kq π 0 dθδ(θ−θ0) =− m 4πq η∞ k0 dkk[f(ηεk+η −μb)−f(ηεk+η +ω−μb)] =− m2 4πq η∞ k2 0/(2m) dε[f(ηε +η −μb)−f(ηε +η +ω−μb)] =− m2 4πq η η∞ k2 0/(2m) dε[f(ε+−ημb)−f(ε++ηω −ημb)] =−m2kBT 4πq η ηln 1+exp ημb− kBTb−k2 0 2mkBTb−ln 1+exp ημb−ηω − kBTb−k2 0 2mkBTb. (B5) 2. Interband term In this case, η=−η. Hence, mχ(inter) b(q,ω)=−1 4π η∞ 0 dkk2π 0 dθsin θ[f(ηεk+η −μb)−f(ηεk+η +ω−μb)] ×δ(ω+2η +ηεk+ηεk+q).(B6) It is clear that, for the δfunction not to vanish, it must be η=−1. This in turn implies that ω−2−k2 2m−k2+q2+2kq cos(θ0) 2m=0⇒cos(θ0)=m kqω−2−2k2+q2 2m.(B7) Solutions exist for −kq m⩽ω−2−2k2+q2 2m⩽kq m⇒⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ k2−kq +2m−mω+q2 2⩽0 k2+kq +2m−mω+q2 2⩾0⇒|k−|<k<k+, where k±=q±4(mω−2m)−2q2 2,ω>2+q2 2m.(B8) 115422-8 ULTRAFAST ELECTRONIC HEAT DISSIPATION THROUGH … PHYSICAL REVIEW B 106, 115422 (2022) (a) (b) (c) (d) FIG. 4. The density-density response function of bulk electrons in units of 2/(¯hvF)3(vF=0.5×106m/s), for a fixed value of the wave vector qand as a function of energy ω(in units of the half-gap ). In each panel we show three curves, one for each temperature (T= 10,300, and 600 K, respectively). (a) Here q=0.5qand the chemical potential is μb=0 (undoped system). (b) Same as in panel (a), but for q=1.5q. (c) Here q=0.5qand the chemical potential is μb=1.4(n-doped system). (d) Same as in panel (c), but for q=1.5q. The parameters used in these plots are the same as those used in Sec. V, except that =50 meV. We also defined q=/(¯hvF). Therefore, mχ(inter) b(q,ω)=−1 4π∞ 0 dkk2[f(−εk−−μb)−f(−εk−+ω−μb)] m kq π 0 dθδ(θ−θ0) =− m 4πqk+ |k−| dkk[f(−εk−−μb)−f(−εk−+ω−μb)] =m2 4πqk2 +/(2m) k2 −/(2m) dε[f(ε++μb)−f(ε+−ω+μb)] =m2kBT 4πqln 1+exp −μb+ kBTb−k2 − 2mkBTb−ln 1+exp −μb+ kBTb−k2 + 2mkBTb −ln 1+exp −μb−ω+ kBTb−k2 − 2mkBTb+ln 1+exp −μb−ω+ kBTb−k2 + 2mkBTb.(B9) 3. Results Figure 4shows results for the imaginary part of the bulk density-density response function, mχb(q,ω), where its intraand interband parts are given in Eqs. (B5) and (B9), respectively. The parameters used in these plots are those given in Sec. V, with the exception of the half-gap which has been set to =50 meV for convenience reasons. Different panels correspond to different values of the chemical potential and of the wave vector q[in units of q=/(¯hvF)] used. In each plot the density-density response function of bulk electrons is shows in units of 2/(¯hvF)3, where vF=0.5×106m/s, as a function of energy ω.The latter is in units of the half-gap . In each panel we show three curves, one for each temperature, T=10,300, and 600 K, respectively. 115422-9