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

Alessandro Principi<EMAIL_ADDRESS>Department of
Physics and Astronomy, University of Manchester, Manchester M13 9PL, UK
Klaas-Jan Tielrooij  Catalan Institute of Nanoscience and Nanotechnology
(ICN2), BIST & CSIC, Campus UAB, 08193 Bellaterra, Barcelona, Spain

###### Abstract

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 relatively slow
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 non-retarded 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 dominate the cooling dynamics
under certain circumstances. We show that this cooling mechanism can lead to a
sub-picosecond time scale, significantly faster than the cooling dynamics
normally observed in Dirac materials.

##  I Introduction

Two-dimensional (2D) Dirac materials [1, 2, 3, 4] have been extensively
studied in the last two decades for the diverse range of intriguing properties
they harbor, which could in turn enable a wealth of novel practical
applications [5, 6, 7, 8]. This is the case of graphene [1, 9, 10, 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 [13, 5, 14, 15, 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 photoexcitation with a fluence on the order of
$100\mu{\rm J/cm}^{2}$. 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 phonons of the material [21] or of the
encapsulant [22, 23]. 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.

---
Figure 1:  Schematics of Coulomb cooling of hot surface Dirac fermions into
bulk states. Dirac fermions have a constant Fermi velocity $v_{\rm F}$, and
exhibit a hot thermal distribution characterized by temperature and chemical
potential $T_{\rm s}$ and $\mu_{\rm s}$, respectively. Similarly, bulk states
are described in terms of the effective mass $m$ and band gap $2\Delta$. The
Fermi distribution of bulk states has chemical potential $\mu_{\rm b}$ (set to
zero in later calculations) and temperature $T_{\rm b}<T_{\rm s}$.

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 mechanisms 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 non-retarded Coulomb interactions (see Fig. 1). Thus we term
this mechanism “Coulomb cooling” [24].

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 non-retarded 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 [2, 3, 12, 1, 9, 10,
11], i.e.

$\displaystyle{\cal H}_{\rm s}=\hbar v_{\rm
F}\sum_{{\bm{k}},\alpha,\beta}{\hat{\psi}}^{\dagger}_{{\bm{k}},\alpha,{\rm
s}}{\bm{k}}\cdot{\bm{\sigma}}_{\alpha,\beta}{\hat{\psi}}_{{\bm{k}},\beta,{\rm
s}}~{},$ (1)

where $v_{\rm F}$ is the Fermi velocity and ${\bm{\sigma}}$ 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),
${\hat{\psi}}^{\dagger}_{{\bm{k}},\alpha,{\rm s}}$
(${\hat{\psi}}_{{\bm{k}},\alpha,{\rm s}}$) creates (destroys) a surface
particle of two-dimensional wavevector ${\bm{k}}$ (momentum $\hbar{\bm{k}}$)
and (pseudo)spin projection $\alpha$. The band energy is $\varepsilon^{({\rm
s})}_{{\bm{k}},\lambda}=\lambda\hbar v_{\rm F}|{\bm{k}}|$, where $\lambda=+$
($\lambda=-$) 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 re-defining the
spin quantization axes.] Eq. (1) is in form identical for graphene and
topological insulators. However, their electrons exhibit different
degeneracies $N_{\rm f}$. While the surface states of topological insulators
are helical [2], thus yielding $N_{\rm f}=1$, in graphene they feature a full
spin-valley degeneracy [1]. Thus, for graphene, the number of fermion flavors
is $N_{\rm f}=4$.

We model the bulk as a two-band system confined into the half-space $z>0$ and
described by the Hamiltonian

$\displaystyle{\cal H}_{\rm b}=\sum_{{\bm{k}},\eta}\varepsilon^{({\rm
b})}_{{\bm{k}},\eta}{\hat{\psi}}^{\dagger}_{{\bm{k}},\eta,{\rm
b}}{\hat{\psi}}_{{\bm{k}},\eta,{\rm b}}~{}.$ (2)

Here, $\varepsilon^{({\rm
b})}_{{\bm{k}},\eta}=\eta\big{[}\Delta+\hbar^{2}|{\bm{k}}|^{2}/(2m)\big{]}$,
where $2\Delta$ is the bulk-band gap and $m$ is the bulk-band mass (assumed to
be the same for both valence and conduction bands), while
${\hat{\psi}}^{\dagger}_{{\bm{k}},\eta,{\rm b}}$
(${\hat{\psi}}_{{\bm{k}},\eta,{\rm b}}$) creates (destroys) a bulk particle of
three-dimensional wavevector ${\bm{k}}=({\bm{k}}_{\parallel},k_{z})$, in band
$\eta$. Here, ${\bm{k}}_{\parallel}$ and $k_{z}$ are the wavevectors parallel
and perpendicular to the surface $z=0$, $\eta=+$ denotes the conduction band,
while $\eta=-$ 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: $\Psi_{{\bm{k}},\eta,{\rm
b}}({\bm{r}},z)=\sqrt{2/V}e^{i{\bm{k}}_{\parallel}\cdot{\bm{r}}}\sin(k_{z}z)$,
where ${\bm{r}}$ is a vector along the surface.

We assume that bulk and surface electrons are coupled electrostatically by
long-range instantaneous Coulomb interactions [25], whose Hamiltonian is

$\displaystyle{\cal H}_{\rm
sb}=\frac{1}{2}\int_{0}^{\infty}dz\sum_{{\bm{q}}}V_{{\rm
sb}}({\bm{q}},z){\hat{n}}_{{\bm{q}},{\rm s}}{\hat{n}}_{-{\bm{q}},{\rm
b}}(z)~{}.$ (3)

Here, ${\bm{q}}$ is a two-dimensional wavevector along the surface of the
topological insulator, while ${\hat{n}}_{{\bm{q}},{\rm s}}$ and
${\hat{n}}_{-{\bm{q}},{\rm b}}(z)$ are the 2D-Fourier transforms of the
surface and bulk density operators. The precise form of the interaction
$V_{{\rm sb}}({\bm{q}},z)$ is determined by the solution of the associated
Poisson (electrostatic) problem, as we proceed to show.

###  II.1 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 wavevector ${\bm{q}}$ in the surface sheet. In response to this added
charge, induced bulk and surface densities are generated. These are named
$n_{\rm b}({\bm{r}},z)$ and $n_{\rm s}({\bm{r}})$, respectively. The resulting
Poisson equation (in Gaussian units) is

$\displaystyle{\bm{\nabla}}\cdot\big{[}\epsilon(z){\bm{\nabla}}\phi({\bm{r}},z)\big{]}=-4\pi
e^{2}n({\bm{r}},z)~{},$ (4)

where $n({\bm{r}},z)=n_{\rm b}({\bm{r}},z)\Theta(z)+n_{\rm
s}({\bm{r}})\delta(z)+e^{-i{\bm{q}}\cdot{\bm{r}}}\delta(z)$ is the total
electron density, and $\epsilon(z)=\epsilon_{\rm b}\Theta(z)+\epsilon_{\rm
vac}\Theta(-z)$ is the dielectric function. Here, $\epsilon_{\rm b}$ is the
relative dielectric constant of the topological insulator, which accounts for
the screening due to filled bands, while $\epsilon_{\rm vac}=1$ [$\Theta(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

$\displaystyle\partial_{z}\big{[}\epsilon(z)\partial_{z}\phi_{\bm{q}}(z)\big{]}-q^{2}\epsilon(z)\phi_{\bm{q}}(z)=-4\pi
e^{2}n_{\bm{q}}(z)~{},$ (5)

where $n_{\bm{q}}(z)=n_{{\bm{q}},{\rm b}}(z)\Theta(z)+n_{{\bm{q}},{\rm
s}}\delta(z)+\delta(z)$. Within linear response, we write $n_{{\bm{q}},{\rm
s}}=\chi_{\rm s}({\bm{q}},\omega)\phi_{\bm{q}}(0)$, where $\chi_{\rm
s}({\bm{q}},\omega)$ is the density-density response function (polarizability)
of the surface states, and

$\displaystyle n_{{\bm{q}},{\rm b}}(z)=\int_{0}^{\infty}dz^{\prime}\chi_{\rm
b}({\bm{q}},\omega,z,z^{\prime})\phi_{\bm{q}}(z^{\prime})\simeq{\bar{\chi}}_{\rm
b}({\bm{q}},\omega)\phi_{\bm{q}}(z)~{}.$

Here, $\chi_{\rm b}({\bm{q}},\omega,z,z^{\prime})\simeq{\bar{\chi}}_{\rm
b}({\bm{q}},\omega)\delta(z-z^{\prime})$ is the density-density response
function [26] (polarizability) of bulk states treated within a semi-local
approximation. Equations of this section contain the response of the doped
bulk, embodied by $\chi_{\rm b}({\bm{q}},\omega,z,z^{\prime})$ and
${\bar{\chi}}_{\rm b}({\bm{q}},\omega)$, 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 ($k_{\rm B}T_{\rm b}$) is much
smaller than the bulk-band gap ($2\Delta$). Using these expressions, the
Poisson equation (5) can be split into the two half-spaces as

$\displaystyle\left\\{\begin{array}[]{ll}\epsilon_{\rm
b}(\partial_{z}^{2}-q^{2})\phi_{\bm{q}}(z)=-4\pi e^{2}{\bar{\chi}}_{\rm
b}({\bm{q}},\omega)\phi_{\bm{q}}(z)&{\rm if}~{}z>0\vspace{0.2cm}\\\
(\partial_{z}^{2}-q^{2})\phi_{\bm{q}}(z)=0&{\rm if}~{}z<0\end{array}\right..$
(9)

The boundary conditions, at the surface $z=0$ and at $z\to\pm\infty$, are
$\phi_{\bm{q}}(0^{+})=\phi_{\bm{q}}(0^{-})$, $\phi_{\bm{q}}(z\to\pm\infty)=0$,
and

$\displaystyle\epsilon_{\rm b}\partial_{z}\phi_{\bm{q}}(z)\Big{|}_{z\to
0^{+}}-\partial_{z}\phi_{\bm{q}}(z)\Big{|}_{z\to 0^{-}}=-4\pi
e^{2}(n_{{\bm{q}},{\rm s}}+1)~{}.$

Introducing the bulk Thomas-Fermi wavevector [26] $q_{\rm
TF}^{2}(q,\omega)=-4\pi e^{2}{\bar{\chi}}_{\rm
b}({\bm{q}},\omega)/\epsilon_{\rm b}$, the solution of Eqs. (9) with the
boundary conditions above is

$\displaystyle\left\\{\begin{array}[]{ll}\phi_{\bm{q}}(z)={\bar{\phi}}_{{\bm{q}}}e^{-\sqrt{q^{2}+q_{\rm
TF}^{2}(q,\omega)}z}&{\rm if}~{}z>0\vspace{0.2cm}\\\
\phi_{\bm{q}}(z)={\bar{\phi}}_{{\bm{q}}}e^{qz}&{\rm
if}~{}z<0\end{array}\right.~{},$ (14)

where, using Eq. (II.1), we get

$\displaystyle{\bar{\phi}}_{{\bm{q}}}=\frac{4\pi
e^{2}}{\epsilon\sqrt{q^{2}+q_{\rm TF}^{2}(q,\omega)}+q-4\pi
e^{2}\chi_{s}({\bm{q}},\omega)}~{}.$ (15)

In the limit of zero frequency and small wavevectors, $\chi_{\rm
s}({\bm{q}},\omega)$ becomes the densities of states of surface electrons.

##  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^{({\rm s})}_{{\bm{k}},\lambda}$ and $f^{({\rm
b})}_{{\tilde{\bm{k}}},\eta}$ at the temperatures $T_{\rm s}$ and $T_{\rm b}$
and chemical potentials $\mu_{\rm s}$ and $\mu_{\rm b}$, respectively. As in
Sect. II, ${\bm{k}}$ and ${\tilde{\bm{k}}}$ are their two- and three-
dimensional wavevectors, respectively. We recall that $\lambda=\pm 1$ is used
to denote the two surface bands (together forming a Dirac cone), while
$\eta=\pm 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^{({\rm
s})}_{{\bm{k}},\lambda}$ is [22]

$\displaystyle\partial_{t}f^{({\rm s})}_{{\bm{k}},\lambda}=-{\cal
I}_{{\bm{k}},\lambda}^{(\rm sb)}~{},$ (16)

where ${\cal I}_{{\bm{k}},\lambda}^{(\rm sb)}$ 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. (16) by the energy of the surface
state, $\varepsilon^{({\rm s})}_{{\bm{k}},\lambda}$, and sum over all
wavevectors ${\bm{k}}$ and all values of the surface-band index $\lambda$. The
left-hand side of the so-obtained equation yields the time derivative of the
energy stored in surface states, $\partial_{t}E_{\rm s}$. The latter is
rewritten as $\partial_{t}E^{\rm s}=C_{\rm s}\partial_{t}T_{s}$, where the
heat capacity of surface states is defined as [22]

$\displaystyle C_{\rm
s}=\sum_{\lambda}\int\frac{d^{2}{\bm{k}}}{(2\pi)^{2}}\xi^{({\rm
s})}_{{\bm{k}},\lambda}\left(-\frac{\partial f^{({\rm
s})}_{{\bm{k}},\lambda}}{\partial\xi^{({\rm
s})}_{{\bm{k}},\lambda}}\right)\left[\frac{\xi^{({\rm
s})}_{{\bm{k}},\lambda}}{T_{\rm s}}+\frac{\partial\mu_{\rm s}}{\partial T_{\rm
s}}\right].$ (17)

Here, $\xi^{({\rm s})}_{{\bm{k}},\lambda}=\varepsilon^{({\rm
s})}_{{\bm{k}},\lambda}-\mu_{s}$. The derivative of the surface chemical
potential $\mu_{\rm s}\equiv\mu_{\rm s}(T_{\rm s})$ with respect to
temperature is obtained by imposing the conservation of the surface electron
density [22]. Defining the power dissipated into bulk states as ${\cal
Q}=\sum_{{\bm{k}},\lambda}{\cal I}_{{\bm{k}},\lambda}^{(\rm
sb)}\varepsilon^{({\rm s})}_{{\bm{k}},\lambda}\equiv\gamma C_{\rm s}(T_{\rm
s}-T_{\rm b})$, Eq. (16) yields the cooling rate [22]

$\displaystyle\gamma=\frac{{\cal Q}}{(T_{\rm s}-T_{\rm b})C_{\rm s}}~{}.$ (18)

In the following we first define the collision integral and then calculate
${\cal Q}$, and thus $\gamma$.

##  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. (16) reads

$\displaystyle{\cal I}_{{\bm{k}},\lambda}^{(\rm sb)}$ $\displaystyle=$
$\displaystyle 2\frac{2\pi}{\hbar{\cal
A}^{2}L_{z}}\sum_{{\bm{k}}^{\prime},{\tilde{\bm{k}}}^{\prime},{\bm{q}}}\sum_{\lambda^{\prime}}\sum_{\eta,\eta^{\prime}}\int_{-\infty}^{\infty}d\omega
V_{{\bm{q}},{\bm{k}}^{\prime},{\tilde{\bm{k}}}^{\prime}}^{2}F_{{\bm{k}},\lambda;{\bm{k}}+{\bm{q}},\lambda^{\prime}}\delta(\varepsilon_{{\bm{k}},\lambda}^{({\rm
s})}-\varepsilon_{{\bm{k}}+{\bm{q}},\lambda^{\prime}}^{({\rm
s})}+\omega)\delta(\varepsilon_{{\bm{k}}^{\prime},\eta}^{({\rm
b})}-\varepsilon_{{\tilde{\bm{k}}}^{\prime},\eta^{\prime}}^{({\rm
b})}-\omega)\delta({\bm{k}}^{\prime}_{\parallel}-{\tilde{\bm{k}}}^{\prime}_{\parallel}-{\bm{q}})$
(19) $\displaystyle\times$ $\displaystyle\big{[}f_{{\bm{k}},\lambda}^{({\rm
s})}f_{{\bm{k}}^{\prime},\eta}^{({\rm
b})}(1-f_{{\bm{k}}+{\bm{q}},\lambda^{\prime}}^{({\rm
s})})(1-f_{{\tilde{\bm{k}}}^{\prime},\eta^{\prime}}^{({\rm
b})})-(1-f_{{\bm{k}},\lambda}^{({\rm s})})(1-f_{{\bm{k}}^{\prime},\eta}^{({\rm
b})})f_{{\bm{k}}+{\bm{q}},\lambda^{\prime}}^{({\rm
s})}f_{{\tilde{\bm{k}}}^{\prime},\eta^{\prime}}^{({\rm b})}\big{]}~{},$

where the factor $2$ upfront accounts for the spin degeneracy of bulk states,
$\omega$ and ${\bm{q}}$ are the transferred energy and momentum parallel to
the surface, ${\bm{k}}^{\prime}_{\parallel}$ and
${\tilde{\bm{k}}}^{\prime}_{\parallel}$ are the components of the three-
dimensional momenta ${\bm{k}}^{\prime}$ and ${\tilde{\bm{k}}}^{\prime}$
parallel to the surface, ${\cal A}$ is the surface area of the topological
insulator and $L_{z}$ its extension in the third dimension. In this equation,
$F_{{\bm{k}},\lambda;{\bm{k}}+{\bm{q}},\lambda^{\prime}}$ is the squared
matrix element of the surface-electron density operator between incoming and
outgoing scattering states. Since bulk states are standing waves, the
$z$-components of the three-dimensional momenta ${\bm{k}}^{\prime}$ and
${\tilde{\bm{k}}}^{\prime}$ are taken to be positive. The matrix element of
the screened Coulomb interaction
$V_{{\bm{q}},{\bm{k}}^{\prime},{\tilde{\bm{k}}}^{\prime}}$ is obtained by
integrating it over the incoming and outgoing scattering states, i.e.

$\displaystyle V_{{\bm{q}},{\bm{k}}^{\prime},{\tilde{\bm{k}}}^{\prime}}^{2}$
$\displaystyle=$ $\displaystyle
4\left|\int_{0}^{\infty}dz\sin(k_{z}^{\prime}z)\sin({\tilde{k}}_{z}^{\prime}z)e^{-\sqrt{q^{2}+q_{\rm
TF}^{2}(q,\omega)}z}\right|^{2}$ (20) $\displaystyle\times$
$\displaystyle\big{|}{\bar{\phi}}_{{\bm{q}}}\big{|}^{2}~{}.$

This matrix element is manipulated in App. A to give

$\displaystyle V_{{\bm{q}},{\bm{k}}^{\prime},{\tilde{\bm{k}}}^{\prime}}^{2}$
$\displaystyle\simeq$
$\displaystyle\big{|}{\bar{\phi}}_{{\bm{q}}}\big{|}^{2}\frac{\pi^{2}}{4\sqrt{q^{2}+q_{\rm
TF}^{2}(q,\omega)}}\delta(k_{z}^{\prime}-{\tilde{k}}_{z}^{\prime})$ (21)
$\displaystyle\equiv$ $\displaystyle
V_{\bm{q}}^{2}\delta(k_{z}^{\prime}-{\tilde{k}}_{z}^{\prime})~{}.$

To obtain this expression, we have evaluated the integral on the right-hand
side of Eq. (20) and approximated it in the limit of small $q$ and $q_{\rm
TF}(q,\omega)$ [the latter is justified when the bulk is undoped and $k_{\rm
B}T_{\rm b}\ll 2\Delta$]. In this case, the weight of the integral is located
around the lines ${\tilde{k}}_{z}^{\prime}=\pm k_{z}^{\prime}$. Thus, we
replaced Lorenzians of width $\sqrt{q^{2}+q_{\rm TF}^{2}(q,\omega)}$ with
delta-functions. Finally, we used that ${\tilde{k}}_{z}^{\prime}$ and
$k_{z}^{\prime}$ must be taken to be positive. Putting Eq. (21) back into Eq.
(19), after some lengthy but straightforward algebra we get

$\displaystyle{\cal Q}$ $\displaystyle=$
$\displaystyle\frac{4}{\pi\hbar}\int\frac{d^{2}{\bm{q}}}{(2\pi)^{2}}V_{\bm{q}}^{2}\int_{0}^{\infty}d\omega\omega\big{[}n^{({\rm
s})}(\omega)-n^{({\rm b})}(\omega)\big{]}$ (22) $\displaystyle\times$
$\displaystyle\Im m\chi_{\rm b}({\bm{q}},\omega)\Im m\chi_{\rm
s}({\bm{q}},\omega)~{},$

where $n^{({\rm s/b})}(\omega)=\big{[}e^{\hbar\omega/(k_{\rm B}T_{\rm
s/b})}-1\big{]}^{-1}$. In Eq. (22),

$\displaystyle\Im m\chi_{\rm s}({\bm{q}},\omega)$ $\displaystyle=$
$\displaystyle-\frac{\pi}{{\cal
A}}\sum_{{\bm{k}},\lambda,\lambda^{\prime}}(f^{({\rm
s})}_{{\bm{k}},\lambda}-f^{({\rm
s})}_{{\bm{k}}+{\bm{q}},\lambda^{\prime}})F_{{\bm{k}},\lambda;{\bm{k}}+{\bm{q}},\lambda^{\prime}}$
(23) $\displaystyle\times$
$\displaystyle\delta(\varepsilon_{{\bm{k}},\lambda}^{({\rm
s})}-\varepsilon_{{\bm{k}}+{\bm{q}},\lambda^{\prime}}^{({\rm s})}+\omega)~{},$

and

$\displaystyle\Im m\chi_{\rm b}({\bm{q}},\omega)$ $\displaystyle=$
$\displaystyle-\frac{\pi}{{\cal
A}L_{z}}\sum_{{\bm{k}}^{\prime},\eta,\eta^{\prime}}(f^{({\rm
b})}_{{\bm{k}}^{\prime},\eta}-f^{({\rm
b})}_{{\bm{k}}^{\prime}-{\bm{q}},\eta^{\prime}})$ (24) $\displaystyle\times$
$\displaystyle\delta(\varepsilon_{{\bm{k}}^{\prime},\eta}^{({\rm
b})}-\varepsilon_{{\bm{k}}^{\prime}-{\bm{q}},\eta^{\prime}}^{({\rm
b})}-\omega)~{},$

are the imaginary parts of the density-density response functions of surface
[27, 28, 29] and bulk [26] states, respectively. The function $\Im m\chi_{\rm
b}({\bm{q}},\omega)$, not to be confused with the function ${\bar{\chi}}_{\rm
b}({\bm{q}},\omega)$ introduced in Sect. II.1, is calculated in App. B.

The cooling rate $\gamma$ is obtained by inserting ${\cal Q}$ in Eq. (22) back
into the definition (18), and reads

$\displaystyle\gamma$ $\displaystyle=$ $\displaystyle-\frac{4}{\pi\hbar
C_{s}}\int\frac{d^{2}{\bm{q}}}{(2\pi)^{2}}V_{\bm{q}}^{2}\int_{0}^{\infty}d\omega\omega\frac{n^{({\rm
s})}(\omega)-n^{({\rm b})}(\omega)}{T_{s}-T_{\rm b}}$ (25)
$\displaystyle\times$ $\displaystyle\Im m\chi_{\rm b}({\bm{q}},\omega)\Im
m\chi_{\rm s}({\bm{q}},\omega)~{}.$

This equation is the key result of this paper.

\begin{overpic}[width=216.81pt]{fig2a} \put(35.0,35.0){(a)} \end{overpic}
---
\begin{overpic}[width=216.81pt]{fig2b} \put(35.0,35.0){(b)} \end{overpic}
Figure 2:  Panel (a): the cooling rate for Dirac-like surface states ($v_{\rm
F}=0.5\times 10^{6}~{}{\rm m}/{\rm s}$) in proximity to an undoped 3D bulk
($n_{b}=0$, which translates into $\mu_{\rm b}=0$) kept at room temperature
($T_{\rm b}=300~{}{\rm K}$, much smaller than the gap energy: $k_{\rm B}T_{\rm
b}\ll 2\Delta$), calculated from Eq. (18) and plotted as a function of
temperature. Curves exhibit a maximum at a surface electron temperature
corresponding to half the bulk-band gap, $\Delta=100~{}{\rm meV}$ (dashed
line). Panel (b): the cooling rate is very weakly dependent on surface carrier
density. The parameter used in these calculations are given at the beginning
of Sect. V and are recalled here for convenience. The bulk electron mass is
set to [30, 31] $m=0.21m_{e}$, where $m_{e}=9.1\times 10^{-31}~{}{\rm kg}$ is
the bare electron mass, while the dielectric constant of the (undoped)
topological insulator is taken to be [32] $\epsilon_{\rm b}=10$.

##  V Results

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 sub-picosecond cooling times.

We start with the surface states of the topological insulator. We numerically
evaluate Eq. (18), with the heat capacity of surface states and power lost to
bulk states given by Eqs. (17) and (22), respectively. For a sake of
definiteness, in the numerical calculations, we set the surface Fermi velocity
[33] $v_{\rm F}=0.5\times 10^{6}~{}{\rm m}/{\rm s}$, while the bulk is left
undoped, i.e. its electron density is $n_{b}=0$ (which translates into
$\mu_{\rm b}=0$). The bulk temperature is set to $T_{\rm b}=300~{}{\rm K}$,
which is much smaller than the bulk-band gap [34] $2\Delta=200~{}{\rm meV}$.
Under these conditions, the bulk bands are nearly unpopulated and we can thus
safely take $q_{\rm TF}(q,\omega)=0$ in Eqs. (15) and (20). Finally, the bulk
electron mass is set to [30, 31] $m=0.21m_{e}$, where $m_{e}=9.1\times
10^{-31}~{}{\rm kg}$ is the bare electron mass, while the dielectric constant
of the (undoped) topological insulator is taken to be [32] $\epsilon_{\rm
b}=10$. This value reflects the fact that the topological insulator we
describe has no charge carriers in the bulk at low temperature (contrary to,
e.g., common topological insulators such as ${\rm Bi}_{2}{\rm Se}_{3}$ which
are metallic). 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 $k_{\rm B}T\sim 10-100~{}{\rm 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. 2 we show the cooling rate for surface-to-bulk Coulomb cooling in
topological insulators (with Fermi energy in the bulk band gap). calculated
from Eq. (18). In Panel (a) we present our numerical results for three
different values of the surface electronic density and as a function of
temperature $T_{\rm s}$. The three chosen densities as $n_{\rm s}=0~{}{\rm
cm}^{-2}$, i.e. an undoped system with Fermi energy at the surface Dirac
crossing, $n_{\rm s}=5\times 10^{11}~{}{\rm cm}^{-2}$ and $n_{\rm
s}=10^{12}~{}{\rm 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. At this temperature, excitations of surface states
occur with typical energies equal to $k_{\rm B}T_{\rm s}$, and become resonant
with the bulk particle-hole excitations, whose weight grows sharply for
energies larger than $2\Delta$. The latter fact is seen in the plots of the
bulk density-density response (the bulk absorption spectrum) given in App. B.

\begin{overpic}[width=216.81pt]{fig3a} \put(35.0,35.0){(a)} \end{overpic}
---
\begin{overpic}[width=216.81pt]{fig3b} \put(35.0,35.0){(b)} \end{overpic}
Figure 3:  The cooling rate for Dirac-like states in graphene ($v_{\rm
F}=10^{6}~{}{\rm m}/{\rm s}$) in proximity to an undoped 3D bulk identical to
that of Fig. 2, calculated from Eq. (18). Panel (a): as a function of
temperature, it exhibits a maximum at a surface electron temperature
corresponding to half the bulk-band gap $\Delta$ (dashed line). Panel (b): the
cooling rate is very weakly dependent on surface carrier density. The
parameters used in these plots are given in Sect. V and Fig. 2.

In Panel (b) we show the cooling rate as a function of the surface density
$n_{\rm s}$ and for three values of the temperature, $T_{\rm s}=500,1000$ and
$1500~{}{\rm 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. 3 we show the cooling rate for a graphene sheet in proximity of a
small-gap bulk material as a function of temperature [in panel (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] $N_{\rm f}=4$, while the Fermi velocity is
doubled, i.e. we use [1] $v_{\rm F}=10^{6}~{}{\rm m/s}$. Comparing Figs. 2 and
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 $\omega\simeq\hbar v_{\rm F}q$. This in
turn implies that, for a given value of momentum $\hbar q$, 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.

##  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 non-retarded 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, at least as a matter of principle, near-field radiative coupling
could out-compete other cooling mechanisms. However, a 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.

##  VII Acknowledgement

A.P. acknowledges support from the European Commission under the EU Horizon
2020 MSCA-RISE-2019 programme (project 873028 HYDROTRONICS) and of the
Leverhulme Trust under the grant 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

The integration of Eq. (20) yields

$\displaystyle
V_{{\bm{q}},{\bm{k}}^{\prime},{\tilde{\bm{k}}}^{\prime}}^{2}=\big{|}{\bar{\phi}}_{{\bm{q}}}\big{|}^{2}\left(\frac{\sqrt{q^{2}+q_{\rm
TF}^{2}(q,\omega)}}{\big{[}q^{2}+q_{\rm
TF}^{2}(q,\omega)\big{]}+(k_{z}^{\prime}-{\tilde{k}}_{z}^{\prime})^{2}}-\frac{\sqrt{q^{2}+q_{\rm
TF}^{2}(q,\omega)}}{\big{[}q^{2}+q_{\rm
TF}^{2}(q,\omega)\big{]}+(k_{z}^{\prime}+{\tilde{k}}_{z}^{\prime})^{2}}\right)^{2}$
(26)

We now observe that, for small $q$ and $q_{\rm TF}^{2}(q,\omega)$, the regime
of interest for the processes we are describing, the two terms on the right-
hand side of Eq. (26) are sharply peaked around
${\tilde{k}}_{z}^{\prime}=k_{z}^{\prime}$ and
${\tilde{k}}_{z}^{\prime}=-k_{z}^{\prime}$, respectively. Thus, we approximate
them with two delta function, making sure that their total integral (over the
variable ${\tilde{k}}_{z}^{\prime}$) remains unchanged. We thus obtain

$\displaystyle V_{{\bm{q}},{\bm{k}}^{\prime},{\tilde{\bm{k}}}^{\prime}}^{2}$
$\displaystyle\simeq$ $\displaystyle\left(\frac{4\pi
e^{2}}{(\epsilon_{b}+1)q-4\pi
e^{2}\chi_{s}({\bm{q}},\omega)}\right)^{2}\frac{\pi^{2}}{4\sqrt{q^{2}+q_{\rm
TF}^{2}(q,\omega)}}\big{[}\delta(k_{z}^{\prime}-{\tilde{k}}_{z}^{\prime})+\delta(k_{z}^{\prime}+{\tilde{k}}_{z}^{\prime})\big{]}$
(27) $\displaystyle\simeq$ $\displaystyle\left(\frac{4\pi
e^{2}}{(\epsilon_{b}+1)q-4\pi
e^{2}\chi_{s}({\bm{q}},\omega)}\right)^{2}\frac{\pi^{2}}{4\sqrt{q^{2}+q_{\rm
TF}^{2}(q,\omega)}}\delta(k_{z}^{\prime}-{\tilde{k}}_{z}^{\prime})~{}.$

In the last line we noticed that $k_{z}^{\prime}$ and
${\tilde{k}}_{z}^{\prime}$ must be taken as positive, since the bulk
wavefunctions describe standing waves (negative wavevectors correspond to the
same wavefunction).

\begin{overpic}[width=106.23808pt]{fig4a} \put(30.0,150.0){(a)} \end{overpic} | \begin{overpic}[width=106.23808pt]{fig4b} \put(30.0,150.0){(b)} \end{overpic}
---|---
\begin{overpic}[width=106.23808pt]{fig4c} \put(30.0,150.0){(c)} \end{overpic} | \begin{overpic}[width=106.23808pt]{fig4d} \put(30.0,150.0){(d)} \end{overpic}
Figure A1:  The density-density response function of bulk electrons in units
of $\Delta^{2}/(\hbar v_{\rm F})^{3}$ ($v_{\rm F}=0.5\times 10^{6}~{}{\rm
m/s}$), for a fixed value of the wavevector $q$ and as a function of energy
$\omega$ (in units of the half-gap $\Delta$). In each panel we show three
curves, one for each temperature ($T=10,300$ and $600~{}{\rm K}$,
respectively). Panel (a): Here $q=0.5q_{\Delta}$ and the chemical potential is
$\mu_{\rm b}=0$ (undoped system). Panel (b): Same as in panel (a), but for
$q=1.5q_{\Delta}$. Panel (c): Here $q=0.5q_{\Delta}$ and the chemical
potential is $\mu_{\rm b}=1.4\Delta$ (n-doped system). Panel (d): Same as in
panel (c), but for $q=1.5q_{\Delta}$. The parameters used in these plots are
the same as those used in Sect. V, except that $\Delta=50~{}{\rm meV}$. We
also defined $q_{\Delta}=\Delta/(\hbar v_{\rm F})$.

##  Appendix B The density-density response of bulk electrons

The imaginary part of the density-density function in Eq. (24) is given by

$\displaystyle\Im
m\chi_{b}({\bm{q}},\omega)=-\pi\sum_{\eta,\eta^{\prime}}\int\frac{d^{3}{\bm{k}}}{(2\pi)^{3}}\big{[}f(\eta\varepsilon_{{\bm{k}}}+\eta\Delta-\mu_{\rm
b})-f(\eta\varepsilon_{{\bm{k}}}+\eta\Delta+\omega-\mu_{\rm
b})\big{]}\delta(\omega+\eta\varepsilon_{{\bm{k}}}+\eta\Delta-\eta^{\prime}\varepsilon_{{\bm{k}}+{\bm{q}}}-\eta^{\prime}\Delta)~{}.$

Here we rewrote
$\varepsilon_{{\bm{k}},\eta}=\eta\Delta+\eta\varepsilon_{\bm{k}}$, introducing
$\varepsilon_{\bm{k}}=k^{2}/(2m)$. We also defined
$f(\xi)=\big{[}e^{\xi/(k_{\rm B}T_{\rm b})}+1\big{]}^{-1}$, and introduced the
chemical potential of bulk bands $\mu_{\rm b}$. In what follows, we assume
$\omega>0$ and analyze separately the intraband ($\eta^{\prime}=\eta$) and
interband ($\eta^{\prime}\neq\eta$) contributions to Eq. (B).

###  B.1 Intraband term

In this case, $\eta^{\prime}=\eta$. Hence,

$\displaystyle\Im m\chi_{b}^{(\rm
intra)}({\bm{q}},\omega)=-\frac{1}{4\pi}\sum_{\eta}\int_{0}^{\infty}dkk^{2}\int_{0}^{\pi}d\theta\sin\theta\big{[}f(\eta\varepsilon_{{\bm{k}}}+\eta\Delta-\mu_{\rm
b})-f(\eta\varepsilon_{{\bm{k}}}+\eta\Delta+\omega-\mu_{\rm
b})\big{]}\delta(\omega+\eta\varepsilon_{{\bm{k}}}-\eta\varepsilon_{{\bm{k}}+{\bm{q}}})~{}.$

The delta function implies that

$\displaystyle\omega+\eta\frac{k^{2}}{2m}-\eta\frac{k^{2}+q^{2}+2kq\cos(\theta_{0})}{2m}=0\quad\Rightarrow\quad\cos(\theta_{0})=\frac{m}{kq}\left(\frac{q^{2}}{2m}-\eta\omega\right)~{}.$
(30)

Solutions exist for

$\displaystyle-\frac{kq}{m}\leq\frac{q^{2}}{2m}-\eta\omega\leq\frac{kq}{m}\quad\Rightarrow\quad
k\geq\frac{m}{q}\left|\frac{q^{2}}{2m}-\eta\omega\right|\equiv k_{0}$ (31)

Therefore,

$\displaystyle\Im m\chi_{b}^{(\rm intra)}({\bm{q}},\omega)$ $\displaystyle=$
$\displaystyle-\frac{1}{4\pi}\sum_{\eta}\int_{0}^{\infty}dkk^{2}\big{[}f(\eta\varepsilon_{{\bm{k}}}+\eta\Delta-\mu_{\rm
b})-f(\eta\varepsilon_{{\bm{k}}}+\eta\Delta+\omega-\mu_{\rm
b})\big{]}\frac{m}{kq}\int_{0}^{\pi}d\theta\delta\big{(}\theta-\theta_{0}\big{)}$
$\displaystyle=$ $\displaystyle-\frac{m}{4\pi
q}\sum_{\eta}\int_{k_{0}}^{\infty}dkk\big{[}f(\eta\varepsilon_{{\bm{k}}}+\eta\Delta-\mu_{\rm
b})-f(\eta\varepsilon_{{\bm{k}}}+\eta\Delta+\omega-\mu_{\rm b})\big{]}$
$\displaystyle=$ $\displaystyle-\frac{m^{2}}{4\pi
q}\sum_{\eta}\int_{k_{0}^{2}/(2m)}^{\infty}d\varepsilon\big{[}f(\eta\varepsilon+\eta\Delta-\mu_{\rm
b})-f(\eta\varepsilon+\eta\Delta+\omega-\mu_{\rm b})\big{]}$ $\displaystyle=$
$\displaystyle-\frac{m^{2}}{4\pi
q}\sum_{\eta}\eta\int_{k_{0}^{2}/(2m)}^{\infty}d\varepsilon\big{[}f(\varepsilon+\Delta-\eta\mu_{\rm
b})-f(\varepsilon+\Delta+\eta\omega-\eta\mu_{\rm b})\big{]}$ $\displaystyle=$
$\displaystyle-\frac{m^{2}k_{\rm B}T}{4\pi
q}\sum_{\eta}\eta\left\\{\ln\left[1+\exp\left(\frac{\eta\mu_{\rm
b}-\Delta}{k_{\rm B}T_{\rm b}}-\frac{k_{0}^{2}}{2mk_{\rm B}T_{\rm
b}}\right)\right]-\ln\left[1+\exp\left(\frac{\eta\mu_{\rm
b}-\eta\omega-\Delta}{k_{\rm B}T_{\rm b}}-\frac{k_{0}^{2}}{2mk_{\rm B}T_{\rm
b}}\right)\right]\right\\}~{}.$

###  B.2 Interband term

In this case, $\eta^{\prime}=-\eta$. Hence,

$\displaystyle\Im m\chi_{b}^{(\rm inter)}({\bm{q}},\omega)$ $\displaystyle=$
$\displaystyle-\frac{1}{4\pi}\sum_{\eta}\int_{0}^{\infty}dkk^{2}\int_{0}^{\pi}d\theta\sin\theta\big{[}f(\eta\varepsilon_{{\bm{k}}}+\eta\Delta-\mu_{\rm
b})-f(\eta\varepsilon_{{\bm{k}}}+\eta\Delta+\omega-\mu_{\rm b})\big{]}$
$\displaystyle\times$
$\displaystyle\delta(\omega+2\eta\Delta+\eta\varepsilon_{{\bm{k}}}+\eta\varepsilon_{{\bm{k}}+{\bm{q}}})~{}.$

It is clear that, for the delta function not to vanish, it must be $\eta=-1$.
This in turn implies that

$\displaystyle\omega-2\Delta-\frac{k^{2}}{2m}-\frac{k^{2}+q^{2}+2kq\cos(\theta_{0})}{2m}=0\quad\Rightarrow\quad\cos(\theta_{0})=\frac{m}{kq}\left(\omega-2\Delta-\frac{2k^{2}+q^{2}}{2m}\right)~{}.$
(34)

Solutions exist for

$\displaystyle-\frac{kq}{m}\leq\omega-2\Delta-\frac{2k^{2}+q^{2}}{2m}\leq\frac{kq}{m}\quad\Rightarrow\quad\left\\{\begin{array}[]{l}{\displaystyle
k^{2}-kq+2m\Delta-m\omega+\frac{q^{2}}{2}\leq 0\vspace{0.2cm}}\\\
{\displaystyle k^{2}+kq+2m\Delta-m\omega+\frac{q^{2}}{2}\geq
0}\end{array}\right.\quad\Rightarrow\quad|k_{-}|<k<k_{+}~{},$ (36)

where

$\displaystyle
k_{\pm}=\frac{q\pm\sqrt{4(m\omega-2m\Delta)-2q^{2}}}{2}~{},~{}~{}\omega>2\Delta+\frac{q^{2}}{2m}~{}.$
(37)

Therefore,

$\displaystyle\Im m\chi_{b}^{(\rm inter)}({\bm{q}},\omega)$ $\displaystyle=$
$\displaystyle-\frac{1}{4\pi}\int_{0}^{\infty}dkk^{2}\big{[}f(-\varepsilon_{{\bm{k}}}-\Delta-\mu_{\rm
b})-f(-\varepsilon_{{\bm{k}}}-\Delta+\omega-\mu_{\rm
b})\big{]}\frac{m}{kq}\int_{0}^{\pi}d\theta\delta\big{(}\theta-\theta_{0}\big{)}$
(38) $\displaystyle=$ $\displaystyle-\frac{m}{4\pi
q}\int_{|k_{-}|}^{k_{+}}dkk\big{[}f(-\varepsilon_{{\bm{k}}}-\Delta-\mu_{\rm
b})-f(-\varepsilon_{{\bm{k}}}-\Delta+\omega-\mu_{\rm b})\big{]}$
$\displaystyle=$ $\displaystyle\frac{m^{2}}{4\pi
q}\int_{k_{-}^{2}/(2m)}^{k_{+}^{2}/(2m)}d\varepsilon\big{[}f(\varepsilon+\Delta+\mu_{\rm
b})-f(\varepsilon+\Delta-\omega+\mu_{\rm b})\big{]}$ $\displaystyle=$
$\displaystyle\frac{m^{2}k_{\rm B}T}{4\pi
q}\left\\{\ln\left[1+\exp\left(-\frac{\mu_{\rm b}+\Delta}{k_{\rm B}T_{\rm
b}}-\frac{k_{-}^{2}}{2mk_{\rm B}T_{\rm
b}}\right)\right]-\ln\left[1+\exp\left(-\frac{\mu_{\rm b}+\Delta}{k_{\rm
B}T_{\rm b}}-\frac{k_{+}^{2}}{2mk_{\rm B}T_{\rm b}}\right)\right]\right.$
$\displaystyle-$ $\displaystyle\left.\ln\left[1+\exp\left(-\frac{\mu_{\rm
b}-\omega+\Delta}{k_{\rm B}T_{\rm b}}-\frac{k_{-}^{2}}{2mk_{\rm B}T_{\rm
b}}\right)\right]+\ln\left[1+\exp\left(-\frac{\mu_{\rm
b}-\omega+\Delta}{k_{\rm B}T_{\rm b}}-\frac{k_{+}^{2}}{2mk_{\rm B}T_{\rm
b}}\right)\right]\right\\}~{}.$

###  B.3 Results

Fig. A1 shows results for the imaginary part of the bulk density-density
response function, $\Im m\chi_{b}({\bm{q}},\omega)$, where its intra- and
inter-band parts are given in Eqs. (B.1) and (38), respectively. The
parameters used in these plots are those given in Sect. V, with the exception
of the half-gap which has been set to $\Delta=50~{}{\rm meV}$ for convenience
reasons. Different panels correspond to different values of the chemical
potential and of the wavevector $q$ [in units of $q_{\Delta}=\Delta/(\hbar
v_{\rm F})$] used. In each plot the density-density response function of bulk
electrons is shows in units of $\Delta^{2}/(\hbar v_{\rm F})^{3}$, where
$v_{\rm F}=0.5\times 10^{6}~{}{\rm m/s}$, as a function of energy $\omega$.
The latter is in units of the half-gap $\Delta$. In each panel we show three
curves, one for each temperature, $T=10,300$ and $600~{}{\rm K}$,
respectively.

## References

  * Castro Neto et al. [2009] A. H. Castro Neto, F. Guinea, N. M. R. Peres, K. S. Novoselov, and A. K. Geim, Rev. Mod. Phys. 81, 109 (2009), URL https://link.aps.org/doi/10.1103/RevModPhys.81.109.
  * Hasan and Kane [2010] M. Z. Hasan and C. L. Kane, Rev. Mod. Phys. 82, 3045 (2010).
  * Hasan and Moore [2011] M. Z. Hasan and J. E. Moore, Annual Review of Condensed Matter Physics 2, 55 (2011), URL https://doi.org/10.1146/annurev-conmatphys-062910-140432.
  * Wehling et al. [2014] T. Wehling, A. Black-Schaffer, and A. Balatsky, Advances in Physics 63, 1 (2014), URL https://doi.org/10.1080/00018732.2014.927109.
  * Fu et al. [2007] L. Fu, C. L. Kane, and E. J. Mele, Phys. Rev. Lett. 98, 106803 (2007), URL https://link.aps.org/doi/10.1103/PhysRevLett.98.106803.
  * Burkov and Hawthorn [2010] A. A. Burkov and D. G. Hawthorn, Phys. Rev. Lett. 105, 066802 (2010), URL https://link.aps.org/doi/10.1103/PhysRevLett.105.066802.
  * Pesin and MacDonald [2012] D. Pesin and A. H. MacDonald, Nature Materials 11, 409 (2012), URL https://doi.org/10.1038/nmat3305.
  * Lupi and Molle [2020] S. Lupi and A. Molle, Applied Materials Today 20, 100732 (2020), ISSN 2352-9407, URL https://www.sciencedirect.com/science/article/pii/S2352940720301785.
  * Kotov et al. [2012] V. N. Kotov, B. Uchoa, V. M. Pereira, F. Guinea, and A. H. Castro Neto, Rev. Mod. Phys. 84, 1067 (2012), URL https://link.aps.org/doi/10.1103/RevModPhys.84.1067.
  * Das Sarma et al. [2011] S. Das Sarma, S. Adam, E. H. Hwang, and E. Rossi, Rev. Mod. Phys. 83, 407 (2011), URL https://link.aps.org/doi/10.1103/RevModPhys.83.407.
  * Massicotte et al. [2021] M. Massicotte, G. Soavi, A. Principi, and K.-J. Tielrooij, Nanoscale 13, 8376 (2021), URL http://dx.doi.org/10.1039/D0NR09166A.
  * Qi and Zhang [2011] X.-L. Qi and S.-C. Zhang, Rev. Mod. Phys. 83, 1057 (2011), URL https://link.aps.org/doi/10.1103/RevModPhys.83.1057.
  * Kane and Mele [2005] C. L. Kane and E. J. Mele, Phys. Rev. Lett. 95, 146802 (2005), URL https://link.aps.org/doi/10.1103/PhysRevLett.95.146802.
  * Fu and Kane [2007] L. Fu and C. L. Kane, Phys. Rev. B 76, 045302 (2007), URL https://link.aps.org/doi/10.1103/PhysRevB.76.045302.
  * Teo et al. [2008] J. C. Y. Teo, L. Fu, and C. L. Kane, Phys. Rev. B 78, 045426 (2008), URL https://link.aps.org/doi/10.1103/PhysRevB.78.045426.
  * Yan and Zhang [2012] B. Yan and S.-C. Zhang, Reports on Progress in Physics 75, 096501 (2012), URL https://doi.org/10.1088/0034-4885/75/9/096501.
  * Politano et al. [2017] A. Politano, L. Viti, and M. S. Vitiello, APL Materials 5, 035504 (2017), URL https://doi.org/10.1063/1.4977782.
  * Pandey et al. [2021] A. Pandey, R. Yadav, M. Kaur, P. Singh, A. Gupta, and S. Husale, Scientific Reports 11, 832 (2021), URL https://doi.org/10.1038/s41598-020-80738-8.
  * Kovalev et al. [2021] S. Kovalev, K. J. Tielrooij, J. C. Deinert, I. Ilyakov, N. Awari, M. Chen, A. Ponomaryov, M. Bawatna, T. V. A. G. de Oliveira, L. M. Eng, et al., npj Quantum Materials 6, 84 (2021), URL https://doi.org/10.1038/s41535-021-00384-9.
  * Deinert et al. [2021] J.-C. Deinert, D. Alcaraz Iranzo, R. Pérez, X. Jia, H. A. Hafez, I. Ilyakov, N. Awari, M. Chen, M. Bawatna, A. N. Ponomaryov, et al., ACS Nano 15, 1145 (2021), URL https://doi.org/10.1021/acsnano.0c08106.
  * Pogna et al. [2021] E. A. A. Pogna, X. Jia, A. Principi, A. Block, L. Banszerus, J. Zhang, X. Liu, T. Sohier, S. Forti, K. Soundarapandian, et al., ACS Nano 15, 11285 (2021), URL https://doi.org/10.1021/acsnano.0c10864.
  * Principi et al. [2017] A. Principi, M. B. Lundeberg, N. C. H. Hesp, K.-J. Tielrooij, F. H. L. Koppens, and M. Polini, Phys. Rev. Lett. 118, 126804 (2017), URL https://link.aps.org/doi/10.1103/PhysRevLett.118.126804.
  * Tielrooij et al. [2018] K.-J. Tielrooij, N. C. H. Hesp, A. Principi, M. B. Lundeberg, E. A. A. Pogna, L. Banszerus, Z. Mics, M. Massicotte, P. Schmidt, D. Davydovskaya, et al., Nature Nanotechnology 13, 41 (2018), URL https://doi.org/10.1038/s41565-017-0008-8.
  * Mihnev et al. [2015] M. T. Mihnev, J. R. Tolsma, C. J. Divin, D. Sun, R. Asgari, M. Polini, C. Berger, W. A. de Heer, A. H. MacDonald, and T. B. Norris, Nature Communications 6, 8105 (2015), URL https://doi.org/10.1038/ncomms9105.
  * Principi et al. [2018] A. Principi, E. van Loon, M. Polini, and M. I. Katsnelson, Phys. Rev. B 98, 035427 (2018), URL https://link.aps.org/doi/10.1103/PhysRevB.98.035427.
  * Giuliani and Vignale [2005] G. F. Giuliani and G. Vignale, _Quantum Theory of the Electron Liquid_ (Cambridge University Press, 2005).
  * Wunsch et al. [2006] B. Wunsch, T. Stauber, F. Sols, and F. Guinea, New Journal of Physics 8, 318 (2006), URL https://doi.org/10.1088/1367-2630/8/12/318.
  * Barlas et al. [2007] Y. Barlas, T. Pereg-Barnea, M. Polini, R. Asgari, and A. H. MacDonald, Phys. Rev. Lett. 98, 236601 (2007), URL https://link.aps.org/doi/10.1103/PhysRevLett.98.236601.
  * Principi et al. [2009] A. Principi, M. Polini, and G. Vignale, Phys. Rev. B 80, 075418 (2009), URL https://link.aps.org/doi/10.1103/PhysRevB.80.075418.
  * Orlita et al. [2015] M. Orlita, B. A. Piot, G. Martinez, N. K. S. Kumar, C. Faugeras, M. Potemski, C. Michel, E. M. Hankiewicz, T. Brauner, i. c. v. Drašar, et al., Phys. Rev. Lett. 114, 186401 (2015), URL https://link.aps.org/doi/10.1103/PhysRevLett.114.186401.
  * Lang et al. [2012] M. Lang, L. He, F. Xiu, X. Yu, J. Tang, Y. Wang, X. Kou, W. Jiang, A. V. Fedorov, and K. L. Wang, ACS Nano 6, 295 (2012), URL https://doi.org/10.1021/nn204239d.
  * Kuznetsov et al. [2021] K. Kuznetsov, P. Kuznetsov, A. Frolov, S. Kovalev, I. Ilyakov, A. Ezhov, and G. Kitaeva, Optical Engineering 60, 1 (2021), URL https://doi.org/10.1117/1.OE.60.8.082012.
  * Tang et al. [2013] C. S. Tang, B. Xia, X. Zou, S. Chen, H.-W. Ou, L. Wang, A. Rusydi, J.-X. Zhu, and E. E. M. Chia, Scientific Reports 3, 3513 (2013), URL https://doi.org/10.1038/srep03513.
  * Arakane et al. [2012] T. Arakane, T. Sato, S. Souma, K. Kosaka, K. Nakayama, M. Komatsu, T. Takahashi, Z. Ren, K. Segawa, and Y. Ando, Nature Communications 3, 636 (2012), URL https://doi.org/10.1038/ncomms1639.