comment

# Phonon-assisted carrier cooling in $h$-BN/graphene van der Waals
heterostructures

Sangkha Borah  Okinawa Institute of Science and Technology Graduate
University, Onna-son, Okinawa 904-0495, Japan  Dinesh Yadav  Institute of
Physics, University of Augsburg, 86135 Augsburg, Germany  Okinawa Institute of
Science and Technology Graduate University, Onna-son, Okinawa 904-0495, Japan
Maxim Trushin  Institute for Functional Intelligent Materials, National
University of Singapore, 117544, Singapore  Department of Materials Science
and Engineering, National University of Singapore, 117575 Singapore  Centre
for Advanced 2D Materials, National University of Singapore, Singapore 117546
Fabian Pauly<EMAIL_ADDRESS>Institute of Physics, University of
Augsburg, 86135 Augsburg, Germany  Okinawa Institute of Science and Technology
Graduate University, Onna-son, Okinawa 904-0495, Japan

(August 27, 2024)

###### Abstract

Being used in optoelectronic devices as ultra-thin conductor-insulator
junctions, detailed investigations are needed about how exactly $h$-BN and
graphene hybridize. Here, we present a comprehensive ab initio study of hot
carrier dynamics governed by electron-phonon scattering at the $h$-BN/graphene
interface, using graphite (bulk), monolayer and bilayer graphene as benchmark
materials. In contrast to monolayer graphene, all multilayer structures
possess low-energy optical phonon modes that facilitate carrier
thermalization. We find that the $h$-BN/graphene interface represents an
exception with a comparatively weak coupling between low-energy optical
phonons and electrons. As a consequence, the thermalization bottleneck effect,
known from graphene, survives hybridization with $h$-BN, but is substantially
reduced in all other bilayer and multilayer cases considered. In addition, we
show that the quantum confinement in bilayer graphene does not have a
significant influence on the thermalization time compared to graphite and that
bilayer graphene can hence serve as a minimal model for the bulk counterpart.

Graphene, $h$-BN, interface, van der Waals heterostructures, hot carriers,
thermalization, density functional theory, electron-phonon interaction,
optoelectronics

##  I Introduction

How does an excited electron or hole lose its energy? This question is at the
heart of photovoltaic technology, as the answer offers a hint for assessing
the prospects of a given material in solar cell applications [1]. Van der
Waals (vdW) heterostructures [2] of atomically thin crystals emerge as a
flexible platform that allows for adjustment of the scattering channels and
subsequent control of optoelectronic [3, 4, 5, 6] and photovoltaic [1, 7, 8,
9] properties. The simplest realizations may be bilayer heterostructures made
of a two-dimensional (2D) semiconductor and graphene [10, 11, 12]. The
honeycomb symmetry of graphene’s carbon lattice results in a linear electron
dispersion and high optical phonon excitation energies, which all together
leads to a slow transfer of photocarrier energy to the lattice [13, 14]. On
the semiconducting side, however, the electronic dispersion at the band edge
is parabolic, and electron-phonon (EP) scattering typically turns out to be
much stronger than in graphene [15, 16, 17, 13]. This can lead to the
situation that photoexcited carriers in the semiconducting layer dissipate
their energy much faster than in the attached graphene. Such an asymmetric
behavior of photocarrier energy dissipation on the two sides of a
heterojunction can create a temperature gradient, resulting in an intrinsic
photothermoelectric effect [10]. Combining different materials in the
prototypical bilayer stacks is expected to reveal the distinct properties of
such heterojunctons. Along these lines, a recent review [18] suggests six
strategies for improving the performance of vdW optoelectronic devices. To
make a deliberate choice between the possible strategies, we need to obtain a
better understanding of the photocarrier evolution in vdW heterostructures,
which may strongly differ from that in isolated 2D layers.

A common simplifying assumption in modeling optoelectronic properties of vdW
bilayers is that the two 2D materials are in electrical contact, but that
their single-layer properties remain otherwise unaffected [19]. Following the
electron affinity rule [20] (also known as the Schottky-Mott rule for metal-
semiconductor junctions [21]), one can figure out the energy band alignment
across the interface, see Refs. [22, 23, 24, 25, 26] for recent bilayer
examples. This simple approach might not necessarily provide a complete
physical picture, especially for 2D materials, as the interfacial region in
vdW bilayers represents an integral part of the structure. Hence, vdW
heterostructures made of two monolayers should be seen as a new single
material possessing properties qualitatively different from those of each
individual monolayer.

Since photocarrier relaxation processes take place on femtosecond or
picosecond time scales, they are challenging to probe experimentally [17, 27].
In addition, atomically thin samples require a high sensitivity of optical
detection. Important recent advances on the experimental side have been
achieved by breakthroughs in time-resolved spectroscopy, notably time- and
angle-resolved photoemission spectroscopy (TR-ARPES) and pump-probe techniques
[28, 29, 30, 31]. On the theoretical side, the study of ultrafast hot carrier
relaxation processes requires state-of-the-art techniques [32, 33], which
solve the many-body physics in the time domain [34]. Recent developments in ab
initio computational modelling techniques are capable of resolving individual
scattering processes and provide insights beyond experiments [35, 36, 37, 34,
38, 39, 32, 40, 41].

Photocarrier dynamics in graphene has been studied experimentally [17, 27, 29,
42, 43] as well as theoretically using empirical models [14, 13, 42, 44] and
recently in terms of a parameter-free first principles approach [38, 45].
Graphene is typically characterized on a substrate such as silicon oxide
(SiO2) [46, 47] or silicon carbide (SiC) [48]. Rough substrate surfaces and
impurities on the surface or the graphene itself act as charge traps and lead
to an uneven charge distribution that reduces the carrier mobility in graphene
[49, 46]. On the other hand, being atomically thin, $h$-BN exhibits a smooth
surface almost free from dangling bonds and charge traps [50], a high-
temperature stability [51], and a nearly commensurate lattice structure to
graphene. These qualities have proven $h$-BN to be an ideal substrate for
graphene. Devices based on $h$-BN/graphene have shown high carrier mobilities
[50], ballistic transport [52, 53], and quantum Hall effects [54, 55]. There
are a number of theoretical studies on electronic and phononic properties as
well as electron-phonon couplings for the $h$-BN/graphene system [56, 57, 58,
59, 60]. The role of optical and acoustical phonons in carrier thermalization
is studied well in graphene [45, 38, 42], but the effects of substrate phonons
and an increasing number of graphene or insulating $h$-BN layers are yet to be
explored. These are aspects that we address in the present work.

Here, we use ab initio methods to uncover the effects of interlayer
hybridization on photocarrier dynamics in vdW heterostructures prototyping
graphene and $h$-BN as constituting monolayers. As compared to the monolayer,
new interlayer phonon modes arise in bilayer and multilayer structures. Caused
by the weak interlayer coupling, low-energy optical modes will emerge in
particular from what are acoustical or flexural modes in the decoupled
monolayers. The photocarrier thermalization bottleneck in graphene occurs for
excitation energies below the high energy optical phonon modes [38], i.e.
electronic excitation energies below around 150 meV. It is thus of particular
interest to explore, how these new optical phonon modes at low energies impact
carrier dynamics and whether they remove the thermalization bottleneck. We
will show that photocarrier dynamics nearly preserves monolayer features, if
graphene is placed on $h$-BN. In contrast the bottleneck is strongly reduced
in bilayer graphene, which resembles graphite from the photocarrier
thermalization point of view.

The rest of the paper is organized as follows. We present our theoretical
approach to study photocarrier dynamics in Sec. II. Results are discussed in
Sec. III, followed by the summary and outlook in Sec. IV.

##  II Theoretical approach

###  II.1 Electronic and phononic properties

We combine density functional theory (DFT) [61, 62] with density functional
perturbation theory (DFPT) [63] to determine the EP scattering. With the Kohn-
Sham eigenstates $|n\mathbf{k}\rangle$ and the derivative of the self-
consistent potential $\partial_{p\mathbf{q}}V$, EP coupling matrix elements
can be computed as,

$g_{mn,p}(\mathbf{k,q})=\frac{1}{\sqrt{2\omega_{p\mathbf{q}}}}\left\langle
m\mathbf{k+q}|\partial_{p\mathbf{q}}V|n\mathbf{k}\right\rangle,$ (1)

where $m,n$ denote electronic band indices, and $\hbar\omega_{p\mathbf{q}}$ is
the energy of the phonon mode $p$ at wave vector $\mathbf{q}$ [34]. The matrix
element $g_{mn,p}(\mathbf{k,q})$ defines an electronic transition from the
initial state $|n\mathbf{k}\rangle$ to the final state $|m\mathbf{k+q}\rangle$
by scattering from the phonon $p\mathbf{q}$. We determine the electron self-
energy $\Sigma_{n\mathbf{k}}(T)$ due to EP interaction within the Migdal
approximation [34, 35],

$\displaystyle\Sigma_{n\mathbf{k}}(T)$ $\displaystyle=$
$\displaystyle\sum_{m,p}\int_{\text{BZ}}\frac{d^{3}q}{\Omega_{\text{BZ}}}\left|g_{mn,p}(\mathbf{k},\mathbf{q})\right|^{2}$
$\displaystyle\times$
$\displaystyle\Bigg{[}\frac{N_{p\mathbf{q}}(T)+f^{(0)}_{m\mathbf{k+q}}(T)}{\varepsilon_{n\mathbf{k}}-(\varepsilon_{m\mathbf{k+q}}-\varepsilon_{\text{F}})+\hbar\omega_{p\mathbf{q}}+\text{i}\eta}$
$\displaystyle+\frac{N_{p\mathbf{q}}(T)+1-f^{(0)}_{m\mathbf{k+q}}(T)}{\varepsilon_{n\mathbf{k}}-(\varepsilon_{m\mathbf{k+q}}-\varepsilon_{\text{F}})-\hbar\omega_{p\mathbf{q}}+\text{i}\eta}\Bigg{]},$

where $\varepsilon_{\text{F}}=0$ is the Fermi energy,
$f^{(0)}_{n\mathbf{k}}(T)=1/[\exp(\frac{\varepsilon_{n\textbf{k}}-\varepsilon_{\text{F}}}{k_{\text{B}}T})+1]$
is the Fermi-Dirac distribution,
$N_{p\mathbf{q}}(T)=1/[\exp(\frac{\hbar\omega_{p\mathbf{q}}}{k_{\text{B}}T})-1]$
is the Bose-Einstein distribution, $\Omega_{\text{BZ}}$ is the volume of the
Brillouin zone (BZ), and $\eta$ is a small broadening parameter. The first
term in the parentheses in Eq. (II.1) represents the self-energy contribution
from the absorption of phonons and the second one stems from their emission.
In order to accurately map the EP scattering events in the whole BZ, it is
important that the self-energy calculations are performed on fine grids over
the electronic ($\mathbf{k}$) and phononic ($\mathbf{q}$) states. For this
purpose we use an interpolation scheme based on maximally localized Wannier
functions (MLWF) [64] in the software package perturbo [65].

To determine ground state electronic and phononic properties, DFT within the
local density approximation (LDA), as implemented in quantum espresso [66], is
used. Our calculations are performed with Slater exchange and Perdew-Wang
correlation (exchange-correlation functional "SLA PW NOGX NOGC") [66]. Core
electrons are modelled through optimized norm-conserving Vanderbilt
pseudopotentials [67]. We employ plane wave (PW) basis sets with a kinetic
energy cutoff of 90 Ry and a charge density cutoff of 360 Ry. The energy
smearing is set to 0.02 Ry. We optimize unit cells and atoms therein using the
Broyden-Fletcher-Goldfarb-Shanno algorithm [68], applying tight convergence
criteria on forces ($10^{-6}$ Ry/a.u.) and total energies ($10^{-8}$ Ry). To
avoid artificial interactions, we separate periodic images through a vacuum of
16 Å (monolayer graphene) or 20 Å (bilayer graphene and $h$-BN/graphene) along
the $z$-direction for the 2D systems and truncate the Coulomb interaction in
the out-of-plane direction (flag "assume_isolated=2D" [69]). VdW corrections
are included for bilayer systems and graphite according to a scheme of Grimme
[70, 66, 71]. Calculations of electronic wavefunctions in the BZ are performed
on a $\Gamma$ centered $36\times 36\times 1$ coarse $\mathbf{k}$ mesh for the
2D systems, while we choose a $12\times 12\times 12$ $\mathbf{k}$ mesh for
bulk graphite.

We compute EP couplings in Eq. (1) and the self-energy in Eq. (II.1) with the
help of the perturbo code [65]. For that purpose the integration in Eq. (II.1)
is carried out on a $1800\times 1800\times 1$ $\mathbf{k}$ mesh for the 2D
systems, along with a $\mathbf{q}$ mesh containing $10^{7}$ points, which we
sample from a uniform distribution. In addition we truncate Coulomb
interactions for the 2D systems (flag "system2d=True"). For bulk graphite, we
select a $1800\times 1800\times 6$ $\mathbf{k}$ mesh with $10^{7}$ uniformly
distributed $\mathbf{q}$ points. Together with $\eta=10$ meV, we find this to
be adequate to achieve converged results. Beyond that we assume that EP
couplings are sufficiently weak to neglect their influence on electron
wavefunctions and phonon dynamical matrices [65]. Anharmonic effects, causing
the renormalization of phononic properties, are also disregarded for
simplicity [35].

###  II.2 Time-evolution of excited charge carriers

We describe the time evolution of the electronic occupation
$f_{n\textbf{k}}(t,T)$ in terms of the Boltzmann equation within the
relaxation time approximation (RTA) [72] via

$\frac{df_{n\textbf{k}}(t,T)}{dt}=-\frac{f_{n\textbf{k}}(t,T)-f^{\text{(0)}}_{n\textbf{k}}(T)}{\tau_{n\textbf{k}}{(T)}}.$
(3)

Assuming the excitation to take place at $t=0$, the solution is

$f_{n\textbf{k}}(t,T)=f^{(0)}_{n\textbf{k}}(T)+e^{-\frac{t}{\tau_{n\textbf{k}}{(T)}}}[f_{n\textbf{k}}(0,T)-f^{\text{(0)}}_{n\mathbf{k}}(T)],$
(4)

where $f_{n\textbf{k}}(0,T)$ denotes the initial occupation of hot carriers
and $f^{\text{(0)}}_{n\mathbf{k}}(T)$ is the previously defined Fermi-Dirac
distribution. We determine the temperature-, band- and wavevector-dependent
scattering rates $\tau_{n\mathbf{k}}^{-1}{(T)}$ that appear in Eqs. (3) and
(4) from the imaginary part of the self-energy in Eq. (II.1) as

$\tau_{n\mathbf{k}}^{-1}(T)=\frac{2}{\hbar}\text{Im}[\Sigma_{n\textbf{k}}(T)].$
(5)

In our analysis, the initial distribution $f_{n\textbf{k}}(0,T)$ of hot
carriers is generated using a Fermi-Dirac distribution at temperature $T$ and
Gaussian distributions with peak values at excitation energies
$\pm(\xi+\Delta/2)$ above and below the Fermi level for holes and electrons,
respectively. In detail,

$f_{n\textbf{k}}(0,T)=f^{(0)}_{n\textbf{k}}(T)\begin{cases}+\lambda_{\text{e}}{(\xi)}e^{-\frac{(\varepsilon_{n\textbf{k}}-\xi-\Delta/2)^{2}}{2\sigma^{2}}},&\varepsilon_{n\textbf{k}}\geq\varepsilon_{\text{F}},\\\
-\lambda_{\text{h}}{(\xi)}e^{-\frac{(\varepsilon_{n\textbf{k}}+\xi+\Delta/2)^{2}}{2\sigma^{2}}},&\varepsilon_{n\textbf{k}}<\varepsilon_{\text{F}}.\end{cases}$
(6)

In the above equation, $\Delta$ is the bandgap of the material, $\sigma=8.47$
meV is a small broadening, and $\xi$ is the excess energy. We furthermore
define the full population as

$P(E,t,T)=\sum_{n\mathbf{k}}\delta(E-\varepsilon_{n\mathbf{k}})\times\begin{cases}f_{n\mathbf{k}}(t,T),&E\geq\varepsilon_{\text{F}},\\\
[1-f_{n\mathbf{k}}(t,T)],&E<\varepsilon_{\text{F}},\end{cases}$ (7)

the time-independent equilibrium population $P_{0}(E,T)$ is obtained from Eq.
(7) by replacing $f_{n\mathbf{k}}(t,T)$ with $f^{(0)}_{n\mathbf{k}}(T)$, and
the difference yields $P_{\text{exc}}(E,t,T)=P(E,t,T)-P_{0}(E,T)$. Electron
and hole densities are related to the energy-, time- and temperature-dependent
photoexcited carrier population $P_{\text{exc}}(E,t,T)$ at $t=0$ through

$\begin{split}n_{\text{e}}\alpha&=\int_{\varepsilon_{\text{F}}}^{\varepsilon_{\text{max}}}P_{\text{exc}}(E,{0},T)dE,\\\
n_{\text{h}}\alpha&=\int_{-\varepsilon_{\text{max}}}^{\varepsilon_{\text{F}}}P_{\text{exc}}(E,{0},T)dE,\end{split}$
(8)

where $\alpha$ is the in-plane area of the unit cell for monolayer graphene,
bilayer graphene and $h$-BN/graphene, or the volume of the unit cell for
graphite. Furthermore, in the expressions we integrate from the Fermi energy
$\varepsilon_{\text{F}}=0$ to a sufficiently large
$\varepsilon_{\text{max}}=4$ eV (see also Fig. 1). For each material, we
choose initial photoexcited carrier densities of electrons and holes to be
equal, $n_{\text{h}}=n_{\text{e}}$, and constant with regard to different
temperatures $T$ and excess energies $\xi$ considered. Rearranging Eq. (8)
expresses the dimensionless functions
$\lambda_{\text{e}}{\color[rgb]{0,0,0}(\xi)}$ and
$\lambda_{\text{h}}{\color[rgb]{0,0,0}(\xi)}$ in Eq. (6) as

$\begin{split}\lambda_{\text{e}}{(\xi)}&=\frac{n_{\text{e}}\alpha}{\int_{\varepsilon_{\text{F}}}^{\varepsilon_{\text{max}}}dE{\sum\limits_{n\mathbf{k}}}\delta(E-\varepsilon_{n\textbf{k}})\exp\left(-\frac{[\varepsilon_{n\mathbf{k}}-\xi-\Delta/2]^{2}}{2\sigma^{2}}\right)},\\\
\lambda_{\text{h}}{(\xi)}&=\frac{n_{\text{h}}\alpha}{\int_{-\varepsilon_{\text{max}}}^{\varepsilon_{\text{F}}}dE{\sum\limits_{n\mathbf{k}}}\delta(E-\varepsilon_{n\textbf{k}})\exp\left(-\frac{[\varepsilon_{n\mathbf{k}}+\xi+\Delta/2]^{2}}{2\sigma^{2}}\right)}.\end{split}$
(9)

Finally, we determine the thermalization time $\tau_{\text{th}}$ through the
relation

$\frac{P_{\text{exc}}(\xi,\tau_{\mathrm{th}},T)}{P_{\text{exc}}(\xi,0,T)}=\frac{1}{e}.$
(10)

This equation expresses $\tau_{\text{th}}$ as a function of $\xi$ and $T$.

Let us note that for the materials studied in Refs. [73, 74] with gap sizes
larger than 1 eV and for temperatures below 1000 K, the thermal equilibrium
population can be ignored, and $P_{\text{exc}}(E,t,T)\approx P(E,t,T)$. Since
we study here graphene and related few-layer crystals with vanishing or low
gap sizes on the order of 100 meV, the thermal contribution needs to be taken
into account, and we thus use the excited population $P_{\text{exc}}(E,t,T)$
in Eq. (10) as well as in the following discussions.

Regarding our method, the RTA is valid for low excitation intensities, leading
to low photoexcited carrier densities $n_{\text{e}}$ and $n_{\text{h}}$. It
describes the photoexcited carrier relaxation directly to the thermalized
state, omitting any intermediate scattering events. In our formalism,
temperatures appear in electronic and phononic occupations. We keep electron
and phonon baths at the same temperature $T$ at all times, see for instance
Eq. (II.1). The approach thus corresponds to a single-temperature model.

##  III Results and Discussion

Table 1: Optimized cell parameters ($a=b$ and $c$), interlayer spacing ($d$) as well as C-C and B-N bond lengths for the studied systems.  System | Distance (Å)
---|---
$a$ |  $c$ |  $d$ |  C-C |  B-N
(a) $h$-BN/graphene |  2.47 |  20 |  3.21 |  1.43 |  1.43
(b) Monolayer graphene |  2.45 |  16 |  — |  1.41 |  —
(c) Bilayer graphene |  2.45 |  20 |  3.31 |  1.41 |  —
(d) Graphite |  2.45 |  6.61 |  3.31 |  1.41 |  —

Optimized geometrical parameters are shown in Table 1. Note that we assume
Bernal stacking for $h$-BN/graphene, bilayer graphene and graphite. The
calculated lattice constants for graphite, monolayer and bilayer graphene are
in good agreement with previous ab initio and experimental values [75, 76, 77,
78, 79, 80, 81]. The difference between the lattice constants of graphene and
$h$-BN is less than 3% [82], and hence we assume the same hexagonal unit cell
for the coupled $h$-BN/graphene system. The optimized lattice parameters for
$h$-BN/graphene match well with literature [57, 78, 83].

Figure 1: Electronic (left column) and phononic (right column) band structures
for (a) $h$-BN/graphene, (b) monolayer graphene, (c) bilayer graphene, and (d)
graphite along $\mathrm{\Gamma-M-K-\Gamma}$ directions in the irreducible
wedge of the BZ. The electronic band structures are shown for PWs (solid black
lines) and MLWFs (red dotted lines) in each case. For $h$-BN/graphene we find
a small band gap of around $150$ meV. Note that we show magnified views of
electronic dispersions, and $\Gamma$ points are thus not visible.

In Fig. 1 we display the electronic and phononic band structures of
$h$-BN/graphene, monolayer graphene, bilayer graphene, and graphite. The good
agreement of electronic band structures obtained from DFT with PWs and those
interpolated using MLWFs demonstrates the quality of the MLWF construction.
The electronic band structure of monolayer graphene features the well-known
linear energy dispersion around the Dirac point, i.e., around the K point in
the BZ. For hypothetical bilayer graphene without interlayer coupling, the
band structure would be doubled. Due to interlayer hybridization the bands
split up. Two pairs of parabolic valence and conduction bands emerge, one pair
being energetically well-separated, while the other pair forms two bands,
which touch each other at the K point and at the Fermi energy [84, 75]. When
graphene is placed on a single $h$-BN layer, the electronic band structure in
the vicinity of the Fermi energy is similar to those of graphene. Since $h$-BN
is a insulator, its valence and conduction bands are energetically well
separated from the linear band crossing at the Dirac point. However, a small
bandgap opens up at the Dirac point [85, 83], because the weak interlayer
coupling to $h$-BN distinguishes the carbon atoms in the graphene lattice.
This induces an effective electronic mass term, which opens the gap [86, 85].
Based on our DFT calculations, the electronic band gap opening amounts to
around $150$ meV, similar to literature values [56, 57]. The low-energy band
structure of graphite resembles those of bilayer graphene. Analogously,
phononic band structures of bilayer graphene and graphite are similar to each
other, and their primitive cells contain the same number and type of atoms. In
comparison to monolayer graphene with 2 atoms in the unit cell, leading to 3
acoustical (one of them turning out to be a flexural mode due to the 2D
character) and 3 optical modes, we expect 6 further optical modes in these two
materials. In particular, we observe three additional low-energy optical modes
at $\Gamma$ below 25 meV, owing to the weak interlayer coupling and a
corresponding splitting of what would be acoustical modes in fully decoupled
monolayers [80]. In the case of the vdW heterostructure $h$-BN/graphene band
degeneracies for acoustical and optical modes are typically removed due to the
difference in masses of the constituting atoms B, N and C, in contrast to
bilayer graphene and graphite [60].

Figure 2: Imaginary part of the EP self-energy Im$[\Sigma_{n\mathbf{k}}(T)]$
as a function of energy for (a) $h$-BN/graphene, (b) monolayer graphene, (c)
bilayer graphene, and (d) graphite at $T=0$ K, and the corresponding
electronic DOS.

Having discussed electronic and phononic properties, we analyze
$\text{Im}[\Sigma_{n\mathbf{k}}(T)]$ from Eq. (II.1). It is shown in Fig. 2 at
$T=0$ K within an energy interval of $\pm 3$ eV around the Fermi energy and
compared to the density of states (DOS) for each of the systems under study.
We observe that the imaginary part of the self-energy follows the electronic
DOS, as the latter essentially represents the phase space for EP scattering.
The rapid increase of $\text{Im}[\Sigma_{n\mathbf{k}}(T)]$ with electronic
energy away from valence and conduction band edges indicates an efficient
energy transfer to phonons for highly excited photocarriers. Electrons or
holes with a high energy can scatter repeatedly by phonons and quickly relax
to lower energy states, until they reach the band edges. Prominent peaks
around $\pm 2$ eV in Fig. 2 can be attributed to the slowly varying energy
bands around the M point in the BZ.

Similar to our previous study on graphene [38], we find that photocarriers
with an energy below the optical phonon energies of about 150 meV relax at a
much slower rate than above. In order to focus on this interesting region, we
will therefore perform our calculations for the remainder of this work within
the 300 meV energy window of conduction and valence band edges, i.e. within
about $\pm 300$ meV around $\varepsilon_{\text{F}}$.

Figure 3: EP scattering rates $\tau_{n\mathbf{k}}^{-1}(T)$ as a function of
energy for (a) $h$-BN/graphene, (b) monolayer graphene, (c) bilayer graphene,
and (d) graphite, evaluated at different temperatures ranging from $T=0$ to
$900$ K.

The EP scattering rates $\tau^{-1}_{n\mathbf{k}}(T)$ of Eq. (5) are shown in
Fig. 3 as a function of energy for different temperatures in the range of $0$
to $900$ K. These temperatures correspond to a realistic range, given that
melting points for graphene and $h$-BN have been reported to be as large as
4900 K [87] and 3300 K [51], respectively. For all the analyzed materials, we
observe that $\tau^{-1}_{n\mathbf{k}}(T)$ varies by orders of magnitude with
regard to energy and temperature.

Let us first concentrate on the low-temperature limit of $T=10$ K. For
graphene, scattering rates increase strongly as a function of energy around
$\varepsilon_{\text{F}}=0$, highlighting inefficient EP scattering until the
highest phonon branches at $150-200$ meV are reached. Above around 150 meV the
scattering rates increase abruptly due to emission of optical phonons [42,
38]. For the multilayer structures we see a more complex behavior of
$\tau^{-1}_{n\mathbf{k}}(T)$ with step-like features within the energy
interval of $\pm 150$ meV around the Fermi energy. For bilayer graphene and
graphite there exist pronounced plateaus in energy from 20 to 150 meV, while
for $h$-BN/graphene we see two plateaus, which range in absolute energies from
95 to 150 meV and 175 to 200 meV. Taking into account the offset $\Delta/2=75$
meV due to the band gap in $h$-BN/graphene, we attribute this behavior in
bilayer graphene, graphite and $h$-BN/graphene to the low-lying optical
phonons with energies as low as 5 meV (see Fig. 1 and Table 2).

For all the materials studied, EP scattering rates furthermore exhibit a
strong temperature dependence below 150 meV and a weaker one above, the energy
threshold coinciding again with the high-lying optical phonon energies. In our
single-temperature model, where electrons and phonons are described by the
same $T$, a higher temperature leads to a larger thermal smearing of both
electronic and phononic bath occupations. The scattering rates, hence,
increase monotonically with temperature, as illustrated in Fig. 3, where the
rates are plotted for $T=0,10,100,300,600$ and $900$ K. As discussed
exemplarily for $10$ K before, scattering rates for temperatures of 300 K or
below show a pronounced energy dependence and promptly rise, as we move away
from the Fermi level or the band edges. At $T\leq 300$ K and energies below
150 meV all systems feature rather low scattering rates. Elevated temperatures
increasingly wash out the energy dependence of $\tau^{-1}_{n\mathbf{k}}(T)$,
see especially the cases of $T=600$ and $900$ K, which we attribute to phonon
absorption. For monolayer graphene, shown in Fig. 3(b), the scattering rate
near the Fermi level increases from below $10^{-6}$ fs-1 at $T=10$ K by some
three orders of magnitude to $10^{-3}$ fs-1 at $T=600$ K. The $h$-BN/graphene
vdW heterostructure behaves in a similar manner at the valence and conduction
band edges. The corresponding values for bilayer graphene and graphite are
increased to a lesser extent from $10^{-3}$ fs-1 at $T=10$ K by around one
order of magnitude to $10^{-2}$ fs-1 at $T=600$ K. We explain this difference
by the parabolic dispersion of low-energy electron bands in multilayer
graphene systems [84, 75], which lifts the mismatch between electronic and
acoustical phonon velocities that exists in monolayer graphene. The scattering
rates for graphite feature a comparatively large spread. Since we observe the
fuzziness only in graphite, we attribute it to the lack of 2D confinement. Let
us also note that the scattering of holes and electrons is very much
comparable in monolayer graphene, in accordance with the highly symmetric DOS
visible in Fig. 2. In the case of $h$-BN/graphene, bilayer graphene, and
graphite scattering rates are slightly higher for holes than for electrons, a
feature that is also partly reflected in the DOS plots of Fig. 2.

Figure 4: Time-dependent excited photocarrier population $P_{\rm exc}(E,t,T)$
as a function of energy, measured in eV, for electrons and holes at an excess
energy $\xi=50$ meV and temperature $T=10$ K for (a) $h$-BN/graphene, (b)
monolayer graphene, (c) bilayer graphene, and (d) graphite.

The time-dependent population $P_{\text{exc}}(E,t,T)$ of excited carriers is
plotted in Fig. 4 at $T=10$ K for an excess energy $\xi=50$ meV, i.e. electron
energy of 50 meV above the conduction band minimum (CBM) and hole energy of 50
meV below the valence band minimum (VBM). At $t=0$, we excite $10^{13}$ cm-2
carriers in the 2D monolayer and bilayer systems and ${1.5\times}10^{13}$ cm-3
in graphite [see Eqs. (8) and (9)]. While monolayer and $h$-BN/graphene show
the slowest thermalization dynamics, bilayer graphene and bulk graphite
thermalize by orders of magnitude faster. Initial photoexcited carrier
populations of electrons and holes are identical in all materials due to the
excitation assumed [see Eq. (6)]. At later times, we find holes to relax
faster than electrons for $h$-BN/graphene, bilayer graphene and graphite in
Fig. 4(a,c,d). This is clearly visible from the different peak heights that
develop with time. Only for graphene electrons and holes relax at a similar
pace. This behavior is expected from the scattering rates, displayed in Fig.
3. There is an additional asymmetry, occurring in electron and hole
populations separately in the course of time, due to the energy dependence of
the scattering rates, see for instance the slight shift of peak maxima towards
$\varepsilon_{\text{F}}$ in Fig. 4.

In order to obtain an overview of photocarrier dynamics, we summarize the
thermalization times $\tau_{\text{th}}$ [see Eq. (10)] of electrons and holes
as a function of excess energy at different $T$ in Fig. 5. For excitations
with $\xi<150$ meV, i.e. below the high-lying optical phonon energies of
monolayer graphene, the thermalization time decreases rapidly with increasing
temperature in all materials studied. Simultaneously, the thermalization time
is substantially extended for the systems containing only a single graphene
layer as compared to the multilayer graphene stacks, if temperatures are
sufficiently low. Finally, beyond the highest optical phonon energy of around
200 meV, carriers thermalize quickly in all materials, demonstrating only a
weak dependence on temperature. Thermalization times then drop to a few
hundred femtoseconds or less, sharply contrasting those at smaller excess
energies $\xi$.

In Fig. 5 graphene shows nearly identical thermalization times for electrons
and holes (with dashed and solid lines coinciding). For the other systems,
electrons thermalize more slowly than holes, as expected from Figs. 3 and 4.
The overall dependence of $\tau_{\text{th}}$ on $\xi$ is very similar for all
carbon-based systems, i.e. for graphene, bilayer graphene and graphite, and
exhibits a plateau region followed by a sharp drop starting from around 150
meV for $T<300$ K. The plateau is somewhat inclined for graphene but much
flatter for bilayer graphene and graphite. The magnitude of thermalization
times at low temperature and $\xi$ is reduced by two orders of magnitude in
bilayer graphene and graphite as compared to monolayer graphene. Maximum
values of thermalization times for the $h$-BN/graphene heterostructure at low
$T$ are reduced by around an order of magnitude as compared to those of
monolayer graphene. Apart from this quantitative difference, we observe a
complex dependence on the excitation energy, where $\tau_{\text{th}}$
decreases in two steps at around 100 and 150 meV, particularly well visible at
$T<300$ K and similar to the plateau features found for the EP scattering
rates in Fig. 3. As before, the deviations from the behavior of monolayer
graphene in the multilayer systems arise from electron and phonon band
hybridization, the resulting emergence of new low-energy optical phonon modes
and the modified EP interactions.

Figure 5: Thermalization times $\tau_{\mathrm{th}}$ of electrons (solid lines)
and holes (dashed lines) as a function of excess energy $\xi$ at different
temperatures for (a) $h$-BN/graphene, (b) monolayer graphene, (c) bilayer
graphene, and (d) graphite.

In the previous paragraphs, we have seen that the absence of low-energy
optical modes leads to the thermalization bottleneck in graphene. Such low-
energy optical modes are however present in the other bi- and multilayer
systems studied here. In order to understand the role of optical phonon modes
with energies below 50 meV in the thermalization, we have computed the
corresponding EP scattering rates for electrons at the Dirac point (K) due to
the phonons at $\Gamma$ for a temperature of 10 K. In this way, we describe
the most relevant low-energy processes, i.e. the interaction of electrons at
the K valley with phonons with vanishing momentum. Let us emphasize that due
to the crossing of phonon bands (see Fig. 1), a fully mode-resolved
description of scattering rates on the whole BZ is not possible, hence we
resort to a single $\mathbf{q}$ point. In Table 2 we present the EP scattering
rates for all of the materials and all those optical modes at $\Gamma$ with
energies below 50 meV. Note that the first three acoustical phonons are
omitted, since they yield no contribution. The reason is that EP couplings
vanish for pure translations. In Fig. 6, we furthermore display the form of
the three energetically lowest optical phonon modes at the $\Gamma$ point for
$h$-BN/graphene, i.e., excluding the acoustical modes. The modes for the other
multilayer systems bilayer graphene and graphite are qualitatively similar, in
particular modes 4 and 5 are always of in-plane character while mode 6
describes an out-of-plane motion.

Table 2 reveals that the optical phonon modes at $\Gamma$ with energies below
50 meV contribute to the scattering of electrons at K for $h$-BN/graphene,
bilayer graphene, and bulk graphite. Their contributions cause a faster
relaxation of hot carriers in contrast to monolayer graphene, where such low-
energy modes are missing. For bilayer graphene and graphite, the low-energy
in-plane optical modes (modes 4 and 5) are mostly responsible for the
scattering, as indicated by the values of $\tau_{1\text{K}}^{-1}$ in Table 2.
Rates of in-plane modes are comparable for these two materials, and
thermalization in bilayer graphene and graphite therefore takes place on
similar time scales in agreement with Fig. 5. While for graphite the
contribution of the out-of-plane optical mode is very small, it is much higher
for bilayer graphene and further enhanced for the $h$-BN/graphene vdW
heterostructure. This shows that the EP scattering due to interlayer coupling
is more pronounced in 2D systems relative to their bulk counterparts. Note
also that the energy of the out-of-plane mode 6 is substantially reduced in
graphite as compared to $h$-BN/graphene and bilayer graphene. A similar
softening is also observable for the in-plane modes, although their energies
remain more similar throughout the considered multilayer materials.
Remarkably, for $h$-BN/graphene the contribution of the out-of-plane mode to
the EP scattering is comparable to the in-plane ones, but the rates
$\tau_{1\text{K}}^{-1}$ for modes 4-6 are overall around two orders of
magnitude lower than those of the important in-plane modes for bilayer
graphene and graphite. This explains why the thermalization times in
$h$-BN/graphene are longer than in these two carbon-based materials and why
the thermalization bottleneck in $h$-BN/graphene is largely preserved.

Let us now relate our outcomes to experimental results on graphene and
graphene-based heterostructures. Dawlaty _et al._ [88] have measured the
carrier relaxation times in epitaxial graphene layers on SiC wafers using
ultrafast optical pump-probe spectroscopy and observed two distinct time
scales associated with the relaxation of nonequilibrium photogenerated
carriers: an initial fast relaxation transient in the 70 to 120 fs range,
followed by a slower relaxation process in the $0.4$ to $1.7$ ps range. These
fast and slow time constants are related to carrier-carrier and carrier-phonon
scattering processes in graphene. Note that epitaxial "graphene" typically
does not represent a single layer but is in fact a stack of a few graphene
layers. Our model predicts that the phonon-related thermalization time does
not strongly depend on temperature in the optical excitation limit ($\xi>0.5$
eV in Ref. [88]) and is indeed around 0.1 ps for bilayer graphene and
graphite. In contrast, the carrier relaxation time in $h$-BN-encapsulated
graphene has been measured for near infrared excitations ($\xi\approx 0.4$ eV)
at temperatures in the range from 10 to 300 K [89] and was found to be one
order of magnitude higher, as expected from our Fig. 5.

Optically excited electrons, injected from graphene into an adjacent 2D
semiconductor, have been studied in Ref. [90] at temperatures ranging from 100
to 300 K. The electronic gap of the semiconductor filters out low-energy (i.e.
thermalized) electrons so that only non-thermalized electrons are collected,
see also Ref. [73]. The measurement shows 25 fs electron injection time from
graphene to the 2D semiconductor with up to 50% quantum yield. Our model
predicts a thermalization time of about 100 fs at the relevant excitation
energies, which matches well with the expectations. Note that the measured
quantum yield shows strong dependence on photon energy but remains nearly
constant with varying photon density (fluence). This is in good agreement with
our model, predicting a thermalization time that strongly depends on $\xi$. At
the same time, it independently justifies our assumption of linear response,
neglecting possible complex nonlinear intensity dependences.

Finally, let us discuss different possibilities to improve the theoretical
approach presented here. We note that DFT methods using exchange correlation
functionals based on LDA do not fully capture the Kohn anomalies in graphene
and graphite [91, 45]. Kohn anomalies have been argued to be responsible for a
prevalent contribution of the high-energy optical phonons to electronic
thermalization in graphite [92, 45]. Despite the shortcomings, our LDA
calculations reproduce this crucial relaxation mechanism. Nevertheless, an
improved description of Kohn anomalies through $GW$ corrections appears to be
worthwhile in future work [91, 45]. Furthermore, the use of polar corrections
[93, 94] might be advised especially for the $h$-BN/graphene system. We
checked that they do not lead to significant changes in the EP scattering
rates of Fig. 3. We did not consider polar corrections here, because for
$h$-BN/graphene they resulted in two imaginary frequencies in the phonon
dispersion at the $\Gamma$ point. Finally, we employ a single temperature for
electrons and phonons and use the RTA. Instead, coupled electron and phonon
distributions may be propagated in time, involving a reevaluation of
scattering rates in every time step [45]. In this way, the occupation of all
electronic and phononic modes in the BZ is tracked, representing a "multi-
temperature" model. In addition, methods may be utilized to better describe
the quantum behavior of electrons and phonons [95, 96].

Table 2: Phonon modes $p$ at $\mathbf{q}=\Gamma=0$ with energy $\hbar\omega_{p\Gamma}<50$ meV, and the scattering rates $\tau_{n\mathbf{k}}^{-1}(T)$ at $T=10$ K for electrons in band $n=1$ at the Dirac point $\mathbf{k}=$K due to the interaction with the corresponding phonon mode $p$ at $\mathbf{q}=\Gamma$. No data is shown for monolayer graphene, since all of its optical phonon modes are located at energies above 50 meV. In order to provide scattering rates that are consistent with the data shown in Fig. 3, we have divided the rates obtained by perturbo for $\mathbf{k}=\text{K}$ and $\mathbf{q}=\Gamma$ by the factor of $10^{7}$, which is the number of total $\mathbf{q}$ points used previously.  Mode $p$ |  $h$-BN/graphene | Monolayer graphene | Bilayer graphene | Graphite
---|---|---|---|---
|  $\hbar\omega_{p\Gamma}$ (meV) |  $\tau_{1\text{K}}^{-1}$ (fs-1) |  $\hbar\omega_{p\Gamma}$ (meV) |  $\tau_{1\text{K}}^{-1}$ (fs-1) |  $\hbar\omega_{p\Gamma}$ (meV) |  $\tau_{1\text{K}}^{-1}$ (fs-1) |  $\hbar\omega_{p\Gamma}$ (meV) |  $\tau_{1\text{K}}^{-1}$ (fs-1)
4 |  8.57 |  $2.20\times 10^{-9}$ |  — |  — |  $6.48$ |  $2.64\times 10^{-7}$ |  $4.93$ |  $5.80\times 10^{-7}$
5 |  8.57 |  $3.21\times 10^{-9}$ |  — |  — |  $6.48$ |  $2.60\times 10^{-7}$ |  $4.93$ |  $5.78\times 10^{-7}$
6 |  20.07 |  $4.49\times 10^{-9}$ |  — |  — |  $19.3$ |  $3.36\times 10^{-11}$ |  $8.48$ |  $4.59\times 10^{-16}$
Figure 6: Optical phonon modes with energies below 50 meV for the
$h$-BN/graphene system at the $\Gamma$ point. The color code of atoms and the
coordinate system are shown on the right side.

##  IV Conclusions

In conclusion, we carried out a detailed study of EP-induced scattering rates
and the related photocarrier thermalization in the $h$-BN/graphene vdW
heterostructure. In this context, monolayer and bilayer graphene as well as
bulk graphite served as reference systems. Combining DFT and DFPT with a MLWF
interpolation, we calculated the excited carrier relaxation times from first
principles, i.e. without free parameters, for a large range of temperatures
and excess energies.

We find that EP scattering rates differ significantly for the investigated
systems at low temperatures and at excitations below around $150$ meV. At
temperatures above 600 K and at higher excitation energies, the scattering
rates turn out to be quite comparable.

Monolayer graphene exhibits an extremely low scattering rate at low
temperatures and at low excitation energies, resulting from its linear band
dispersion and the absence of low-energy optical phonon modes, which leads to
a slow cooling of excited carriers. In contrast to bilayer graphene, where the
photo-carriers thermalize similarly rapidly as in graphite, replacing one
graphene layer with $h$-BN restores the hot-carrier thermalization bottleneck
found in monolayer graphene. We attribute this to a comparatively weak
coupling of low-energy optical phonon modes to the electronic states in this
particular heterostructure. Suppressed EP coupling along with a low impurity
concentration at the $h$-BN/graphene interface may be responsible for the
large electron mobility, typically observed in such structures. Understanding
interfacial coupling mechanisms, we foresee deliberate approaches to design
functional electronic devices based on vdW heterostructures.

## Acknowledgments

The authors thank Goki Eda for stimulating discussions. The computational
facilities provided by the Okinawa Institute of Science and Technology (OIST)
Graduate University are acknowledged. This research is also supported by the
Ministry of Education, Singapore, under its Research Centre of Excellence
award to the Institute for Functional Intelligent Materials (I-FIM, project
No. EDUNC-33-18-279-V12).

## References

  * Das  _et al._ [2019] S. Das, D. Pandey, J. Thomas, and T. Roy, The role of graphene and other 2D materials in solar photovoltaics, Adv. Mater. 31, 1802722 (2019).
  * Geim and Grigorieva [2013] A. K. Geim and I. V. Grigorieva, Van der Waals heterostructures, Nature 499, 419 (2013).
  * Tian  _et al._ [2016] H. Tian, M. L. Chin, S. Najmaei, Q. Guo, F. Xia, H. Wang, and M. Dubey, Optoelectronic devices based on two-dimensional transition metal dichalcogenides, Nano Res. 9, 1543 (2016).
  * Ponraj  _et al._ [2016] J. S. Ponraj, Z.-Q. Xu, S. C. Dhanabalan, H. Mu, Y. Wang, J. Yuan, P. Li, S. Thakur, M. Ashrafi, K. Mccoubrey, Y. Zhang, S. Li, H. Zhang, and Q. Bao, Photonics and optoelectronics of two-dimensional materials beyond graphene, Nanotechnology 27, 462001 (2016).
  * Castellanos-Gomez [2016] A. Castellanos-Gomez, Why all the fuss about 2D semiconductors?, Nat. Photonics 10, 202 (2016).
  * Wang  _et al._ [2012] Q. H. Wang, K. Kalantar-Zadeh, A. Kis, J. N. Coleman, and M. S. Strano, Electronics and optoelectronics of two-dimensional transition metal dichalcogenides, Nat. Nanotechnol. 7, 699 (2012).
  * Wang  _et al._ [2018] L. Wang, L. Huang, W. C. Tan, X. Feng, L. Chen, X. Huang, and K.-W. Ang, 2D photovoltaic devices: Progress and prospects, Small Methods 2, 1700294 (2018).
  * Bati  _et al._ [2019] A. S. R. Bati, M. Batmunkh, and J. G. Shapter, Emerging 2D layered materials for perovskite solar cells, Adv. Energy Mater. 10, 1902253 (2019).
  * Cheng  _et al._ [2018] K. Cheng, Y. Guo, N. Han, X. Jiang, J. Zhang, R. Ahuja, Y. Su, and J. Zhao, 2D lateral heterostructures of group-III monochalcogenide: Potential photovoltaic applications, Appl. Phys. Lett. 112, 143902 (2018).
  * Lin  _et al._ [2019] Y. Lin, Q. Ma, P.-C. Shen, B. Ilyas, Y. Bie, A. Liao, E. Ergeçen, B. Han, N. Mao, X. Zhang, X. Ji, Y. Zhang, J. Yin, S. Huang, M. Dresselhaus, N. Gedik, P. Jarillo-Herrero, X. Ling, J. Kong, and T. Palacios, Asymmetric hot-carrier thermalization and broadband photoresponse in graphene-2D semiconductor lateral heterojunctions, Sci. Adv. 5, eaav1493 (2019).
  * LaGasse  _et al._ [2019] S. W. LaGasse, P. Dhakras, K. Watanabe, T. Taniguchi, and J. U. Lee, Gate-tunable graphene - WSe2 heterojunctions at the Schottky-Mott limit, Adv. Mater. 31, 1901392 (2019).
  * Meng  _et al._ [2015] J. Meng, H.-D. Song, C.-Z. Li, Y. Jin, L. Tang, D. Liu, Z.-M. Liao, F. Xiu, and D.-P. Yu, Lateral graphene p-n junctions formed by the graphene/MoS2 hybrid interface, Nanoscale 7, 11611 (2015).
  * König-Otto  _et al._ [2016] J. C. König-Otto, M. Mittendorff, T. Winzer, F. Kadi, E. Malic, A. Knorr, C. Berger, W. A. de Heer, A. Pashkin, H. Schneider, M. Helm, and S. Winnerl, Slow Noncollinear Coulomb Scattering in the Vicinity of the Dirac Point in Graphene, Phys. Rev. Lett. 117, 087401 (2016).
  * Mihnev  _et al._ [2016] M. T. Mihnev, F. Kadi, C. J. Divin, T. Winzer, S. Lee, C.-H. Liu, Z. Zhong, C. Berger, W. A. de Heer, E. Malic, A. Knorr, and T. B. Norris, Microscopic origins of the terahertz carrier relaxation and cooling dynamics in graphene, Nat. Commun. 7, 11617 (2016).
  * Abrikosov [2017] A. A. Abrikosov,  _Fundamentals of the Theory of Metals_ (Courier Dover Publications, 2017).
  * Song and Levitov [2015] J. C. W. Song and L. S. Levitov, Energy flows in graphene: hot carrier dynamics and cooling, J. Phys.: Condens. Matter 27, 164201 (2015).
  * Brida  _et al._ [2013] D. Brida, A. Tomadin, C. Manzoni, Y. J. Kim, A. Lombardo, S. Milana, R. R. Nair, K. S. Novoselov, A. C. Ferrari, G. Cerullo, and M. Polini, Ultrafast collinear scattering and carrier multiplication in graphene, Nat. Commun. 4, 1987 (2013).
  * Lin  _et al._ [2017] S. Lin, Y. Lu, J. Xu, S. Feng, and J. Li, High performance graphene/semiconductor van der Waals heterostructure optoelectronic devices, Nano Energy 40, 122 (2017).
  * Jeong  _et al._ [2016] H. Jeong, H. M. Oh, S. Bang, H. J. Jeong, S.-J. An, G. H. Han, H. Kim, S. J. Yun, K. K. Kim, J. C. Park, Y. H. Lee, G. Lerondel, and M. S. Jeong, Metal–insulator–semiconductor diode consisting of two-dimensional nanomaterials, Nano Lett. 16, 1858 (2016).
  * Nelson [2003] J. A. Nelson,  _The physics of solar cells_ (Imperial College Press, London, 2003).
  * Bardeen [1947] J. Bardeen, Surface states and rectification at a metal semi-conductor contact, Phys. Rev. 71, 717 (1947).
  * Trushin [2018] M. Trushin, Theory of photoexcited and thermionic emission across a two-dimensional graphene-semiconductor Schottky junction, Phys. Rev. B 97, 195447 (2018).
  * Cao  _et al._ [2019] L. Cao, Y. S. Ang, Q. Wu, and L. K. Ang, Janus PtSSe and graphene heterostructure with tunable Schottky barrier, Appl. Phys. Lett. 115, 241601 (2019).
  * Cao  _et al._ [2020a] L. Cao, G. Zhou, Q. Wu, S. A. Yang, H. Y. Yang, Y. S. Ang, and L. Ang, Electrical contact between an ultrathin topological Dirac semimetal and a two-dimensional material, Phys. Rev. Appl. 13, 054030 (2020a).
  * Cao  _et al._ [2020b] L. Cao, Q. Wu, Y. S. Ang, and L. K. Ang, Tunable band alignment in boron carbon nitride and blue phosphorene van der Waals heterostructure, Nano Express 1, 020021 (2020b).
  * Cao  _et al._ [2021] L. Cao, G. Zhou, Q. Wang, L. K. Ang, and Y. S. Ang, Two-dimensional van der Waals electrical contact to monolayer MoSi2N4, Appl. Phys. Lett. 118, 013106 (2021).
  * Breusing  _et al._ [2011] M. Breusing, S. Kuehn, T. Winzer, E. Malić, F. Milde, N. Severin, J. P. Rabe, C. Ropers, A. Knorr, and T. Elsaesser, Ultrafast nonequilibrium carrier dynamics in a single graphene layer, Phys. Rev. B 83, 153410 (2011).
  * Lui  _et al._ [2010] C. H. Lui, K. F. Mak, J. Shan, and T. F. Heinz, Ultrafast Photoluminescence from Graphene, Phys. Rev. Lett. 105, 127404 (2010).
  * Shang  _et al._ [2011] J. Shang, T. Yu, J. Lin, and G. G. Gurzadyan, Ultrafast electron-optical phonon scattering and quasiparticle lifetime in CVD-grown graphene, ACS Nano 5, 3278 (2011).
  * Johannsen  _et al._ [2013] J. C. Johannsen, S. Ulstrup, F. Cilento, A. Crepaldi, M. Zacchigna, C. Cacho, I. C. E. Turcu, E. Springate, F. Fromm, C. Raidel, T. Seyller, F. Parmigiani, M. Grioni, and P. Hofmann, Direct View of Hot Carrier Dynamics in Graphene, Phys. Rev. Lett. 111, 027403 (2013).
  * Cabo  _et al._ [2015] A. G. Cabo, J. A. Miwa, S. S. Grønborg, J. M. Riley, J. C. Johannsen, C. Cacho, O. Alexander, R. T. Chapman, E. Springate, M. Grioni, J. V. Lauritsen, P. D. C. King, P. Hofmann, and S. Ulstrup, Observation of ultrafast free carrier dynamics in single layer MoS2, Nano Lett. 15, 5883 (2015).
  * Park  _et al._ [2009] C.-H. Park, F. Giustino, C. D. Spataru, M. L. Cohen, and S. G. Louie, First-Principles Study of Electron Linewidths in Graphene, Phys. Rev. Lett. 102, 076803 (2009).
  * Bernardi  _et al._ [2014] M. Bernardi, D. Vigil-Fowler, J. Lischner, J. B. Neaton, and S. G. Louie, Ab Initio Study of Hot Carriers in the First Picosecond after Sunlight Absorption in Silicon, Phys. Rev. Lett. 112, 257402 (2014).
  * Bernardi [2016] M. Bernardi, First-principles dynamics of electrons and phonons, Eur. Phys. J. B 89, 239 (2016).
  * Giustino [2017] F. Giustino, Electron-phonon interactions from first principles, Rev. Mod. Phys. 89, 15003 (2017).
  * Bernardi  _et al._ [2015] M. Bernardi, D. Vigil-Fowler, C. S. Ong, J. B. Neaton, and S. G. Louie, Ab-initio study of hot electrons in GaAs, Proc. Natl. Acad. Sci. 112, 5291 (2015).
  * Winzer  _et al._ [2010] T. Winzer, A. Knorr, and E. Malic, Carrier multiplication in graphene, Nano Lett. 10, 4839 (2010).
  * Yadav  _et al._ [2019] D. Yadav, M. Trushin, and F. Pauly, Photocarrier thermalization bottleneck in graphene, Phys. Rev. B 99, 155410 (2019).
  * Malić  _et al._ [2011] E. Malić, T. Winzer, E. Bobkin, and A. Knorr, Microscopic theory of absorption and ultrafast many-particle kinetics in graphene, Phys. Rev. B 84, 205406 (2011).
  * Calandra and Mauri [2007] M. Calandra and F. Mauri, Electron-phonon coupling and electron self-energy in electron-doped graphene: Calculation of angular-resolved photoemission spectra, Phys. Rev. B 76, 205411 (2007).
  * Winzer and Malić [2012] T. Winzer and E. Malić, Impact of Auger processes on carrier dynamics in graphene, Phys. Rev. B 85, 241404 (2012).
  * Winnerl  _et al._ [2011] S. Winnerl, M. Orlita, P. Plochocka, P. Kossacki, M. Potemski, T. Winzer, E. Malic, A. Knorr, M. Sprinkle, C. Berger, W. A. de Heer, H. Schneider, and M. Helm, Carrier Relaxation in Epitaxial Graphene Photoexcited Near the Dirac Point, Phys. Rev. Lett. 107, 237401 (2011).
  * Gierz [2017] I. Gierz, Probing carrier dynamics in photo-excited graphene with time-resolved ARPES, J. Electron Spectrosc. Relat. Phenom. 219, 53 (2017).
  * Butscher  _et al._ [2007] S. Butscher, F. Milde, M. Hirtschulz, E. Malić, and A. Knorr, Hot electron relaxation and phonon dynamics in graphene, Appl. Phys. Lett. 91, 203103 (2007).
  * Tong and Bernardi [2021] X. Tong and M. Bernardi, Toward precise simulations of the coupled ultrafast dynamics of electrons and atomic vibrations in materials, Phys. Rev. Res. 3, 023072 (2021).
  * Chen  _et al._ [2008] J.-H. Chen, C. Jang, S. Xiao, M. Ishigami, and M. S. Fuhrer, Intrinsic and extrinsic performance limits of graphene devices on SiO2, Nat. Nanotechnol. 3, 206 (2008).
  * Geringer  _et al._ [2009] V. Geringer, M. Liebmann, T. Echtermeyer, S. Runte, M. Schmidt, R. Rückamp, M. C. Lemme, and M. Morgenstern, Intrinsic and extrinsic corrugation of monolayer graphene deposited on SiO2, Phys. Rev. Lett. 102, 076102 (2009).
  * Mendes-de-Sa  _et al._ [2012] T. G. Mendes-de-Sa, A. M. B. Goncalves, M. J. S. Matos, P. M. Coelho, R. Magalhaes-Paniago, and R. G. Lacerda, Correlation between (in)commensurate domains of multilayer epitaxial graphene grown on SiC(000${\bar{1}}$) and single layer electronic behavior, Nanotechnology 23, 475602 (2012).
  * Martin  _et al._ [2008] J. Martin, N. Akerman, G. Ulbricht, T. Lohmann, J. H. Smet, K. von Klitzing, and A. Yacoby, Observation of electron–hole puddles in graphene using a scanning single-electron transistor, Nat. Phys. 4, 144 (2008).
  * Dean  _et al._ [2010] C. R. Dean, A. F. Young, I. Meric, C. Lee, L. Wang, S. Sorgenfrei, K. Watanabe, T. Taniguchi, P. Kim, K. L. Shepard, and J. Hone, Boron nitride substrates for high-quality graphene electronics, Nat. Nanotechnol. 5, 722 (2010).
  * Kim  _et al._ [2017] D.-H. Kim, H.-S. Kim, M. W. Song, S. Lee, and S. Y. Lee, Geometric and electronic structures of monolayer hexagonal boron nitride with multi-vacancy, Nano Convergence 4, 13 (2017).
  * Mayorov  _et al._ [2011] A. S. Mayorov, R. V. Gorbachev, S. V. Morozov, L. Britnell, R. Jalil, L. A. Ponomarenko, P. Blake, K. S. Novoselov, K. Watanabe, T. Taniguchi, and A. K. Geim, Micrometer-scale ballistic transport in encapsulated graphene at room temperature, Nano Lett. 11, 2396 (2011).
  * Wang  _et al._ [2013] L. Wang, I. Meric, P. Y. Huang, Q. Gao, Y. Gao, H. Tran, T. Taniguchi, K. Watanabe, L. M. Campos, D. A. Muller, J. Guo, P. Kim, J. Hone, K. L. Shepard, and C. R. Dean, One-dimensional electrical contact to a two-dimensional material, Science 342, 614 (2013).
  * Dean  _et al._ [2013] C. R. Dean, L. Wang, P. Maher, C. Forsythe, F. Ghahari, Y. Gao, J. Katoch, M. Ishigami, P. Moon, M. Koshino, T. Taniguchi, K. Watanabe, K. L. Shepard, J. Hone, and P. Kim, Hofstadter’s butterfly and the fractal quantum Hall effect in moiré superlattices, Nature 497, 598 (2013).
  * Ponomarenko  _et al._ [2013] L. A. Ponomarenko, R. V. Gorbachev, G. L. Yu, D. C. Elias, R. Jalil, A. A. Patel, A. Mishchenko, A. S. Mayorov, C. R. Woods, J. R. Wallbank, M. Mucha-Kruczynski, B. A. Piot, M. Potemski, I. V. Grigorieva, K. S. Novoselov, F. Guinea, V. I. Fal’ko, and A. K. Geim, Cloning of Dirac fermions in graphene superlattices, Nature 497, 594 (2013).
  * Jung  _et al._ [2015] J. Jung, A. M. DaSilva, A. H. MacDonald, and S. Adam, Origin of band gaps in graphene on hexagonal boron nitride, Nat. Commun. 6, 6308 (2015).
  * Giovannetti  _et al._ [2007] G. Giovannetti, P. A. Khomyakov, G. Brocks, P. J. Kelly, and J. van den Brink, Substrate-induced band gap in graphene on hexagonal boron nitride: Ab initio density functional calculations, Phys. Rev. B 76, 073103 (2007).
  * Fan  _et al._ [2011] Y. Fan, M. Zhao, Z. Wang, X. Zhang, and H. Zhang, Tunable electronic structures of graphene/boron nitride heterobilayers, Appl. Phys. Lett. 98, 083103 (2011).
  * Aggoune  _et al._ [2020] W. Aggoune, C. Cocchi, D. Nabok, K. Rezouali, M. A. Belkhir, and C. Draxl, Structural, electronic, and optical properties of periodic graphene/h-BN van der Waals heterostructures, Phys. Rev. Materials 4, 084001 (2020).
  * Slotman  _et al._ [2014] G. J. Slotman, G. A. de Wijs, A. Fasolino, and M. I. Katsnelson, Phonons and electron-phonon coupling in graphene-_h_ -BN heterostructures, Ann. Phys. 526, 381 (2014).
  * Kohn and Sham [1965] W. Kohn and L. J. Sham, Self-consistent equations including exchange and correlation effects, Phys. Rev. 140, A1133 (1965).
  * Hohenberg and Kohn [1964] P. Hohenberg and W. Kohn, Inhomogeneous electron gas, Phys. Rev. 136, B864 (1964).
  * Baroni  _et al._ [2001] S. Baroni, S. de Gironcoli, A. Dal Corso, and P. Giannozzi, Phonons and related crystal properties from density-functional perturbation theory, Rev. Mod. Phys. 73, 515 (2001).
  * Mostofi  _et al._ [2014] A. A. Mostofi, J. R. Yates, G. Pizzi, Y.-S. Lee, I. Souza, D. Vanderbilt, and N. Marzari, An updated version of wannier90: A tool for obtaining maximally-localised Wannier functions, Comput. Phys. Commun. 185, 2309 (2014).
  * Zhou  _et al._ [2021] J.-J. Zhou, J. Park, I.-T. Lu, I. Maliyov, X. Tong, and M. Bernardi, Perturbo: A software package for ab initio electron–phonon interactions, charge transport and ultrafast dynamics, Comput. Phys. Commun. 264, 107970 (2021).
  * Giannozzi  _et al._ [2017] P. Giannozzi, O. Andreussi, T. Brumme, O. Bunau, M. Buongiorno Nardelli, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, M. Cococcioni, N. Colonna, I. Carnimeo, A. Dal Corso, S. de Gironcoli, P. Delugas, R. A. DiStasio, A. Ferretti, A. Floris, G. Fratesi, G. Fugallo, R. Gebauer, U. Gerstmann, F. Giustino, T. Gorni, J. Jia, M. Kawamura, H.-Y. Ko, A. Kokalj, E. Küçükbenli, M. Lazzeri, M. Marsili, N. Marzari, F. Mauri, N. L. Nguyen, H.-V. Nguyen, A. Otero-de-la-Roza, L. Paulatto, S. Poncé, D. Rocca, R. Sabatini, B. Santra, M. Schlipf, A. P. Seitsonen, A. Smogunov, I. Timrov, T. Thonhauser, P. Umari, N. Vast, X. Wu, and S. Baroni, Advanced capabilities for materials modelling with Quantum ESPRESSO, J. Phys.: Condens. Matter 29, 465901 (2017).
  * Hamann [2013] D. R. Hamann, Optimized norm-conserving Vanderbilt pseudopotentials, Phys. Rev. B 88, 085117 (2013).
  * Fletcher [2000] R. Fletcher,  _Practical methods of optimization_, 2nd ed. (John Wiley & Sons, Chichester, 2000).
  * Sohier  _et al._ [2017] T. Sohier, M. Calandra, and F. Mauri, Density functional perturbation theory for gated two-dimensional heterostructures: Theoretical developments and application to flexural phonons in graphene, Phys. Rev. B 96, 075448 (2017).
  * Grimme [2006] S. Grimme, Semiempirical GGA-type density functional constructed with a long-range dispersion correction, J. Comp. Chem 27, 1787 (2006).
  * Barone  _et al._ [2009] V. Barone, M. Casarin, D. Forrer, M. Pavone, M. Sambi, and A. Vittadini, Role and effective treatment of dispersive forces in materials: Polyethylene and graphite crystals as test cases, J. Comp. Chem. 30, 934 (2009).
  * Lundstrom [2000] M. Lundstrom,  _Fundamentals of Carrier Transport_, 2nd ed. (Cambridge University Press, Cambridge, 2000).
  * Yadav  _et al._ [2020] D. Yadav, M. Trushin, and F. Pauly, Thermalization of photoexcited carriers in two-dimensional transition metal dichalcogenides and internal quantum efficiency of van der Waals heterostructures, Phys. Rev. Res. 2, 043051 (2020).
  * Yadav  _et al._ [2021] D. Yadav, F. Pauly, and M. Trushin, Charge-carrier thermalization in bulk and monolayer CdTe from first principles, Phys. Rev. B 103, 125113 (2021).
  * Castro Neto  _et al._ [2009] A. H. Castro Neto, F. Guinea, N. M. R. Peres, K. S. Novoselov, and A. K. Geim, The electronic properties of graphene, Rev. Mod. Phys. 81, 109 (2009).
  * Kong  _et al._ [2009] B. D. Kong, S. Paul, M. Buongiorno Nardelli, and K. W. Kim, First-principles analysis of lattice thermal conductivity in monolayer and bilayer graphene, Phys. Rev. B 80, 033406 (2009).
  * Lebégue  _et al._ [2010] S. Lebégue, J. Harl, T. Gould, J. G. Ángyán, G. Kresse, and J. F. Dobson, Cohesive properties and asymptotics of the dispersion interaction in graphite by the random phase approximation, Phys. Rev. Lett. 105, 196401 (2010).
  * Ramasubramaniam  _et al._ [2011] A. Ramasubramaniam, D. Naveh, and E. Towe, Tunable band gaps in bilayer graphene-BN heterostructures, Nano Lett. 11, 1070 (2011).
  * Liu  _et al._ [2012] Z. Liu, S.-M. Zhang, J.-R. Yang, J. Z. Liu, Y.-L. Yang, and Q.-S. Zheng, Interlayer shear strength of single crystalline graphite, Acta Mech. Sin. 28, 978 (2012).
  * Yan  _et al._ [2008] J.-A. Yan, W. Y. Ruan, and M. Y. Chou, Phonon dispersions and vibrational properties of monolayer, bilayer, and trilayer graphene: Density-functional perturbation theory, Phys. Rev. B 77, 125401 (2008).
  * Razado-Colambo  _et al._ [2018] I. Razado-Colambo, J. Avila, D. Vignaud, S. Godey, X. Wallart, D. P. Woodruff, and M. C. Asensio, Structural determination of bilayer graphene on SiC(0001) using synchrotron radiation photoelectron diffraction, Sci. Rep. 8, 10190 (2018).
  * Amorim  _et al._ [2013] R. G. Amorim, X. Zhong, S. Mukhopadhyay, R. Pandey, A. R. Rocha, and S. P. Karna, Strain- and electric field-induced band gap modulation in nitride nanomembranes, J. Phys.: Condens. Matter 25, 195801 (2013).
  * Wang  _et al._ [2017] J. Wang, F. Ma, and M. Sun, Graphene, hexagonal boron nitride, and their heterostructures: Properties and applications, RSC Adv. 7, 16801 (2017).
  * McCann and Fal’ko [2006] E. McCann and V. I. Fal’ko, Landau-Level Degeneracy and Quantum Hall Effect in a Graphite Bilayer, Phys. Rev. Lett. 96, 086805 (2006).
  * Kan  _et al._ [2012] E. Kan, H. Ren, F. Wu, Z. Li, R. Lu, C. Xiao, K. Deng, and J. Yang, Why the band gap of graphene is tunable on hexagonal boron nitride, J. Phys. Chem. C 116, 3142 (2012).
  * Semenoff [1984] G. W. Semenoff, Condensed-matter simulation of a three-dimensional anomaly, Phys. Rev. Lett. 53, 2449 (1984).
  * Zakharchenko  _et al._ [2011] K. V. Zakharchenko, A. Fasolino, J. H. Los, and M. I. Katsnelson, Melting of graphene: from two to one dimension, J. Phys.: Condens. Matter 23, 202202 (2011).
  * Dawlaty  _et al._ [2008] J. M. Dawlaty, S. Shivaraman, M. Chandrashekhar, F. Rana, and M. G. Spencer, Measurement of ultrafast carrier dynamics in epitaxial graphene, Appl. Phys. Lett. 92, 042116 (2008).
  * Jadidi  _et al._ [2016] M. M. Jadidi, R. J. Suess, C. Tan, X. Cai, K. Watanabe, T. Taniguchi, A. B. Sushkov, M. Mittendorff, J. Hone, H. D. Drew, M. S. Fuhrer, and T. E. Murphy, Tunable Ultrafast Thermal Relaxation in Graphene Measured by Continuous-Wave Photomixing, Phys. Rev. Lett. 117, 257401 (2016).
  * Chen  _et al._ [2019] Y. Chen, Y. Li, Y. Zhao, H. Zhou, and H. Zhu, Highly efficient hot electron harvesting from graphene before electron-hole thermalization, Sci. Adv. 5, eaax9958 (2019).
  * Lazzeri  _et al._ [2008] M. Lazzeri, C. Attaccalite, L. Wirtz, and F. Mauri, Impact of the electron-electron correlation on phonon dispersion: Failure of LDA and GGA DFT functionals in graphene and graphite, Phys. Rev. B 78, 081406 (2008).
  * Piscanec  _et al._ [2004] S. Piscanec, M. Lazzeri, F. Mauri, A. C. Ferrari, and J. Robertson, Kohn anomalies and electron-phonon interactions in graphite, Phys. Rev. Lett. 93, 185503 (2004).
  * Zhou and Bernardi [2016] J.-J. Zhou and M. Bernardi, Ab initio electron mobility and polar phonon scattering in GaAs, Phys. Rev. B 94, 201201 (2016).
  * Zhou  _et al._ [2018] J.-J. Zhou, O. Hellman, and M. Bernardi, Electron-Phonon Scattering in the Presence of Soft Modes and Electron Mobility in ${\mathrm{SrTiO}}_{3}$ Perovskite from First Principles, Phys. Rev. Lett. 121, 226603 (2018).
  * Kadanoff [2017] L. P. Kadanoff, Entropy is in flux V3.4, J. Stat. Phys. 167, 1039 (2017).
  * Stan  _et al._ [2009] A. Stan, N. E. Dahlen, and R. van Leeuwen, Time propagation of the Kadanoff–Baym equations for inhomogeneous systems, J. Chem. Phys. 130, 224101 (2009).