Toward precise simulations of the coupled ultrafast dynamics
of electrons and atomic vibrations in materials
Xiao Tong
1
and Marco Bernardi
1,
∗
1
Department of Applied Physics and Materials Science,
California Institute of Technology, Pasadena, CA 91125, USA.
Ultrafast spectroscopies can access the dynamics of electrons and nuclei at short timescales, shedding
light on nonequilibrium phenomena in materials. However, development of accurate calculations to
interpret these experiments has lagged behind as widely adopted simulation schemes are limited to
sub-picosecond timescales or employ simplified interactions lacking quantitative accuracy. Here we
show a precise approach to obtain the time-dependent populations of nonequilibrium electrons and
atomic vibrations (phonons) up to tens of picoseconds, with a femtosecond time resolution. Com-
bining first-principles electron-phonon and phonon-phonon interactions with a parallel numerical
scheme to time-step the coupled electron and phonon Boltzmann equations, our method provides
unprecedented microscopic insight into scattering mechanisms in excited materials. Focusing on
graphene as a case study, we demonstrate calculations of ultrafast electron and phonon dynamics,
transient optical absorption, structural snapshots and diffuse X-ray scattering. Our first-principles
approach paves the way for quantitative atomistic simulations of ultrafast dynamics in materials.
Time-domain spectroscopies open a window on the dy-
namics of electrons, phonons and various elementary ex-
citations, probing matter on timescales characteristic of
the interactions among its constituents. In particular,
recently developed experimental techniques enable the
characterization of the coupled dynamics of electronic
and nuclear degrees of freedom with high spatial and
time resolutions [1–4], as well as manipulation of ultrafast
structural dynamics [5–10] and phase transitions [11–13].
The vast information encoded in ultrafast spectroscopy
signals underscores a critical need for quantitative theo-
retical tools. The latter should ideally be able to pre-
dict or interpret ultrafast measurements and shed light
on the underlying microscopic processes and coupling
between different degrees of freedom. A common sce-
nario is the coupled dynamics of excited electrons and
phonons, which governs a wide range of phenomena
such as carrier and lattice equilibration, photoemission
and the associated linewidths, structural transitions and
superconductivity. Heuristic approaches such as two-
temperature models [14] or kinetic equations with ad-
justable interactions [15] are routinely adopted to study
aspects of nonequilibrium electron and phonon dynam-
ics. As they leave out important atomistic details, these
models are only occasionally geared toward quantitative
predictions [16] and are mainly useful as tools to quali-
tatively interpret experimental results.
First-principles computational methods based on den-
sity functional theory (DFT) and related techniques have
made great strides in modeling electron and phonon in-
teractions and nonequilibrium dynamics [17–22]. An im-
portant approach is real-time time-dependent DFT (rt-
TDDFT) [23–25], which propagates in time the elec-
tronic Kohn-Sham equations while treating atomic mo-
tions using Ehrenfest forces.
While rt-TDDFT has
been widely successful for studying electron dynamics
in molecules [26, 27], important open challenges remain,
including reaching simulation times longer than
∼
1 ps,
treating periodic systems (crystals) [28] and better un-
derstanding, improving and validating the accuracy of
exchange-correlation functionals governing the interac-
tions of electrons and nuclei [29–31].
Non-adiabatic
molecular dynamics is also a valuable method to model
ultrafast dynamics in molecules [32], although it is not
commonly employed for solids.
A second family of approaches
−
on which this
work hinges
−
employs perturbation theory to compute
the matrix elements for electron and phonon interac-
tions [33], and combines them with the semiclassical
Boltzmann transport equation (BTE) [20] or quantum
kinetic equations [34] to investigate ultrafast electron dy-
namics. Owing to its momentum-space formulation, this
framework is ideally suited for studies of crystals and
periodic systems. However, in spite of recent progress,
propagating in time the BTEs for coupled electrons and
phonons while fully taking into account their interactions
has remained an elusive goal due to computational cost
and the need for complex algorithms and workflows.
Here we show explicit simulations of the coupled dy-
namics of nonequilibrium electrons and phonons for
timescales up to tens of picoseconds, using accurate
electron-phonon (e-ph) and phonon-phonon (ph-ph) in-
teractions validated through transport calculations. We
time-step the electron and phonon BTEs
−
a large set
of coupled integro-differential equations
−
using a par-
allel numerical approach to efficiently compute the rele-
vant collision integrals and obtain time-dependent elec-
tron and phonon populations. Focusing on graphene,
a material with unconventional physical properties and
unique promise for ultrafast devices, we investigate the
coupled dynamics of excited electrons and phonons, com-
pute spectroscopic signals such as transient absorption
and diffuse X-ray scattering, and analyze the dominant
scattering channels, identifying the slow rise of flexural
arXiv:2009.07958v1 [cond-mat.mtrl-sci] 16 Sep 2020
2
FIG. 1: Simulated carrier and phonon dynamics in graphene.
Excited electrons (red) and holes (purple) in
graphene are initialized in a Fermi-Dirac distribution at 4,000 K electronic temperature at time zero. The
subsequent nonequilibrium dynamics is analyzed by mapping the populations
f
n
k
(
t
) for electrons and 1
−
f
n
k
(
t
) for
holes on the electronic band structure near the Dirac cone (top) and the phonon populations
N
ν
q
(
t
) on the phonon
dispersions (bottom), using a point size proportional to the population values. Shown left of each population panel
are the corresponding BZ-averaged carrier or phonon concentrations as a function of energy (in arbitrary units).
phonons as a bottleneck to equilibration. Taken together,
our work demonstrates a leap forward in simulations of
ultrafast dynamics in solids and provides a quantitative
tool to predict and interpret time-domain spectroscopies.
RESULTS
Numerical approach.
We solve the coupled electron
and phonon BTEs [35] in a homogeneously excited re-
gion of a material using a fourth-order Runge-Kutta algo-
rithm. At each time
t
, we propagate the electronic popu-
lations
f
n
k
(
t
), where
n
is the band index and
k
the crystal
momentum, and the phonon populations
N
ν
q
(
t
), where
ν
is the mode index and
q
the phonon wavevector. Using
an in-house modified version of our
Perturbo
code [36],
we solve the BTEs as a set of integro-differential equa-
tions:
∂f
n
k
(
t
)
∂t
=
I
e-ph
[
{
f
n
k
(
t
)
}
,
{
N
ν
q
(
t
)
}
]
∂N
ν
q
(
t
)
∂t
=
I
ph-e
[
{
f
n
k
(
t
)
}
,
{
N
ν
q
(
t
)
}
]
+
I
ph-ph
[
{
N
ν
q
(
t
)
}
]
,
(1)
where
I
e-ph
is the collision integral for e-ph interactions
in the electron BTE, while
I
ph-e
and
I
ph-ph
are, respec-
tively, the collision integrals for phonon-electron (ph-e)
and ph-ph interactions in the phonon BTE. Detailed ex-
pressions for these integrals are given in Methods.
The collision integrals depend explicitly on the electron
and phonon populations at the current time step, and
thus Eq. (1) is a system of integro-differential equations,
with size
N
e
=
N
b
×N
k
for electrons and
N
ph
=
N
ν
×N
q
for phonons (
N
b
and
N
k
are the number of electronic
bands and
k
-points, while
N
ν
and
N
q
are the number
of phonon modes and
q
-points). In a typical calcula-
tion, our scheme involves time-stepping
∼
10
6
coupled
integro-differential equations, each with 10
7
−
10
9
scat-
tering processes in the collision integrals, a formidable
computational challenge addressed with an efficient algo-
rithm combining MPI and
OpenMP
parallelization. As
memory effects are not included, our BTEs amount to a
Markovian dynamics of the density matrix in which the
off-diagonal coherences, which are typically short-lived
(
<
5 fs), are neglected [15]. Using this scheme, we are
able to reach simulation times of up to
∼
100 ps, with a
fs time step.
Relaxation of excited carriers.
In a first set of
simulations, we model photoexcited electron and hole
carriers in graphene by choosing appropriate carrier
populations at time zero as the initial condition of
the BTEs.
Pump-probe studies on graphene have
shown that fast electron-electron interactions rapidly
thermalize the excited carriers within 10
−
50 fs [37, 38].
Therefore, shortly after excitation, electrons and holes
are well described by a hot Fermi-Dirac distribution
with a temperature of a few thousand degrees for typical
experimental settings [39–41].
We model the carrier
dynamics in graphene after this initial thermalization by
setting the electron and hole populations at time zero
to Fermi-Dirac distributions at 4,000 K temperature,
while setting the initial phonon populations to their
equilibrium value at 300 K.
3
FIG. 2: Phonon scattering channels.
Optical phonons at Γ and K are rapidly populated in the first 300 fs due
to their strong coupling with the excited carriers, but later decay into lower-energy optical and acoustic modes
through ph-ph processes. The dominant energy and momentum-conserving phonon decay channels are shown with
red arrows at 5 ps and 20 ps delay times after the initial electronic excitation.
Figure 1 shows the time evolution of the electron,
hole and phonon populations in the first picosecond
after the initial excitation and carrier thermalization.
For representative times, we analyze both the carrier
and phonon populations (
f
n
k
for electrons, 1
−
f
n
k
for holes, and
N
ν
q
for phonons) along high-symmetry
Brillouin zone (BZ) lines, as well as BZ-averaged carrier
and phonon concentrations as a function of energy (see
Methods).
In the first 100 fs, the hot electron and
hole distributions become narrower in energy as the
carriers cool and settle near the Dirac point in about 300
fs. This rapid cooling process is achieved by emitting
optical phonons through intravalley and intervalley e-ph
scattering. As a result, the optical phonons generated
within 1 ps possess wavevectors close to the BZ Γ-point
for intravalley, and K-point for intervalley scattering.
Although most of the excess energy stored in the
excited carriers is dissipated during the first 300 fs,
surprisingly it takes more than 5 ps for the carriers
to fully reach equilibrium with negligible population
fluctuations. This slow and lingering cooling process is
due to two factors
−
one is the slow rise of out-of-plane
flexural phonons due to their weak e-ph coupling (see
below), and the other is the weak anharmonic ph-ph
coupling, which allows the optical phonons generated
during the rapid carrier cooling to linger for
>
10 ps,
heating back the carriers through phonon absorption
processes.
Phonon equilibration.
Accessing the 10
−
100 ps time
scale characteristic of phonon dynamics is an open chal-
lenge for existing first-principles simulations of ultrafast
dynamics. In Fig. 2, we analyze the long-lived equili-
bration of optical phonons generated by the excited car-
riers. The optical phonons first decay within
∼
5 ps to
transverse acoustic (TA) and longitudinal acoustic (LA)
phonons, which then emit lower-energy acoustic modes
through ph-ph scattering processes, ultimately equili-
brating on a 20
−
100 ps time scale (not shown).
Our calculations can directly identify the main phonon
modes and relaxation pathways from the time-dependent
phonon populations (see Fig. 2). In the first 300 fs, the
excess electronic energy rapidly generates modes with
strong e-ph coupling, including the A
′
1
and E
2
g
optical
modes, respectively with wavevectors near the K and Γ
points of the BZ. A slower change in phonon populations
occurs between 1
−
5 ps, when the optical phonons decay
into LA and TA modes with wavevectors halfway between
Γ and the K or M points of the BZ. Long-wavelength
(small-
q
) acoustic phonons are also generated in the same
time window. A second stage of phonon equilibration
sets in after
∼
20 ps, when energy redistribution occurs
primarily within the LA and TA branches (see Fig. 2).
We carry out a mode-by-mode analysis of the phonon
populations and effective temperatures.
In Fig. 3a,
we plot for each mode the excess population relative
to equilibrium, ∆
N
ν
(
t
) =
N
ν
(
t
)
−
N
ν
(
−∞
), where
N
ν
(
t
) are BZ-averaged phonon populations for each
mode. Also shown in Fig. 3a are the mode-resolved
e-ph coupling strengths,
g
2
ν
(see Methods).
Due to
their strong e-ph coupling, the in-plane longitudinal
optical (LO) and transverse optical (TO) modes are
rapidly emitted during the initial fast carrier cooling,
highlighting the key role of e-ph interactions at early
times [4]. The LO, TO and LA modes are populated
extensively before 1 ps, while the TA and ZO modes
are excited more gradually through ph-ph interactions
over timescales longer than 1 ps. Due to their weak
e-ph coupling, the out-of-plane flexural phonon modes
(ZA and ZO) are generated at a significantly slower rate
than in-plane phonons, resulting in a slow rise of their
populations in the first 10 ps. The ZA mode, with the
weakest e-ph scattering strength, exhibits the slowest
rise in population among all modes. In turn, the weak
e-ph coupling and slow generation of flexural phonons is
responsible for a carrier cooling bottleneck: As electrons
and holes interact with hot ZA and ZO populations,
they gain energy from the flexural modes for over 10 ps,
well beyond the initial sub-ps carrier cooling.
4
FIG. 3: Mode-resolved phonon populations and
effective temperatures. a
Time-dependent excess
phonon populations ∆
N
ν
(
t
), with e-ph coupling
strengths
g
2
ν
(in arb. units) given in the inset;
b
,
Time-dependent effective temperatures
T
ν
. All
quantities are shown separately for each of the six
phonon modes in graphene.
Although at short times the phonon distributions are
still non-thermal, we find that after 2
−
5 ps all phonon
populations are thermal and well approximated by hot
Bose-Einstein distributions. The mode-resolved effective
phonon temperatures are computed as BZ averages
of state-dependent temperatures (see Methods) and
shown in Fig. 3b. We find that the effective phonon
temperatures are strongly mode-dependent throughout
the simulation, their trend mirroring the respective
excess phonon populations in Fig. 3a. This result shows
that a two-temperature model would fail dramatically
in graphene, both because the phonon temperatures
are not well defined before 2 ps and because they are
mode-dependent at longer times.
Our calculations shed light on the entire equilibra-
tion cascade, with timescales spanning several orders
of magnitudes, from fs for carrier cooling through
e-ph interactions to ps for phonon downconversion to
equilibration via ph-ph processes over tens of ps. The
microscopic details of the scattering processes, accessed
easily in our approach due to its formulation in momen-
tum space, are out of reach for existing first-principles
simulations. As shown next, the time-dependent electron
and phonon populations further allow us to reliably
simulate various ultrafast spectroscopies employed as
probes of electron and nuclear dynamics.
Pump-probe transient absorption.
Time-resolved
differential transmission measurements are widely em-
ployed to shed light on nonequilibrium carrier dynamics
and investigate the timescales associated with electronic
processes. We build on the calculations shown above to
compute transient absorption at optical wavelengths.
In graphene, at low energies within 1
−
2 eV of the
Dirac cone, the differential transmissivity ∆
T/T
0
is pro-
portional to the transient carrier populations; in a single
particle picture and neglecting the state-dependence of
the optical transition dipoles [42],
∆
T
(
E
)
T
0
≈−
a
0
[
∆
f
e
(
E
2
)
+ ∆
f
h
(
−
E
2
)]
(2)
where
E
is carrier energy,
a
0
is the absorption coeffi-
cient of single-layer graphene,
f
e
,
h
are BZ-averaged elec-
tron or hole populations, and ∆
f
e
,
h
(
E,t
) =
f
e
,
h
(
E,t
)
−
f
e
,
h
(
E,
−∞
) the corresponding excess carrier popula-
tions relative to equilibrium.
Figure 4a compares the experimental [43] and cal-
culated ∆
T/T
0
as a function of pump-probe delay
time at 900 nm wavelength. The experimental result
exhibits two distinct temporal regions.
At sub-10
fs time delay, ∆
T/T
0
increases due to the coupled
effects of photoexcitation and electron-electron (e-e)
interactions, both of which are not treated explicitly
in our calculations but are taken into account in our
choice of an initial hot Fermi-Dirac distribution. The
second time window at
t>
10 fs, modeled here explicitly,
shows a slow time decay of ∆
T/T
0
with a time constant
of a few ps.
Our simulated differential transmission
spectrum agrees well with experiment [43] in this time
window
−
following a faster decay of ∆
T/T
0
due to
carrier cooling in the first 300 fs, we correctly predict
the onset of a slower cooling of the excited carriers
associated with the bottleneck from slowly rising flex-
ural phonons discussed above.
While previous work
focused on the role of e-e processes [43], our results
show the importance of e-ph interactions at sub-ps times.
Diffuse scattering and structural dynamics.
In-
tense experimental efforts are focused on investigating
atomic motions and structural dynamics in the time
domain [44]. Among other techniques, ultrafast elec-
tron [45] and X-ray diffraction [46] have made great
5
FIG. 4: Simulated ultrafast spectroscopies and
structural snapshots in graphene.
a
, Normalized differential transmission ∆
T/T
0
at 900
nm wavelength as a function of pump-probe delay time.
Our calculated result is compared with experiments
from Ref. [43].
b
, Structural snapshots at 0.5 ps and 5
ps after pump.
c
, Simulated ultrafast diffuse scattering
pattern, ∆
I
(
q
,t
)
/I
, at 0.5 ps and 5 ps after pump,
shown together with the reciprocal lattice.
strides and are currently able to infer atomic displace-
ment patterns and dominant phonon modes governing
the nonequilibrium structural dynamics [44, 47]. Us-
ing the time-resolved phonon populations for each vibra-
tional mode and wavevector, we can simulate and recon-
struct the structural dynamics fully from first principles.
In a classical picture, the time-dependent displacement
of atom
a
can be decomposed into lattice waves with
wavevector
q
and mode
ν
[46],
u
a
(
t
) =
∑
ν
q
A
ν
q
(
t
)
a
ν
q
e
−
i
(
ω
ν
q
t
−
q
·
r
a
−
δ
ν
q
)
(3)
where
A
ν
q
is the wave amplitude,
a
ν
q
its polarization
vector projected on atom
a
, and
δ
ν
q
a phase factor. For
quantized lattice vibrations (phonons),
ν
q
is the phonon
eigenvector obtained by diagonalizing the dynamical ma-
trix and the amplitude
A
ν
q
can be expressed in terms of
phonon populations
N
ν
q
as
A
ν
q
(
t
) =
√
[
2
N
ν
q
(
t
) + 1
]
~
m
a
ω
ν
q
,
(4)
where
m
a
is the atomic mass. The phase factors
δ
ν
q
are chosen randomly at time zero and do not affect the
dynamics. The computed classical atomic displacements
u
a
at time delays of 0.5 ps and 5 ps are visualized in
Fig. 4b. At 0.5 ps, the excitation of optical phonons from
carrier cooling drives a dominant unidirectional vibration
of atoms. At longer times, as shown in the 5 ps panel,
a larger number of phonon modes are incoherently ex-
cited through ph-ph processes and the atomic vibrations
become randomized, consistent with the thermal phonon
distributions we find at 5 ps. The quantity measured ex-
perimentally in X-ray diffraction is the averaged square
atomic displacement,
〈
u
2
a
〉
=
∑
ν
q
1
2
〈
A
2
ν
q
〉
a
, which we
show in Supplementary Video 1, where the radius of the
blue sphere around each atom represents
〈
u
2
a
〉
. The aver-
aged atomic displacements increase monotonically with
time, reaching greater values as the system approaches
thermal equilibrium.
Thermal diffuse X-ray scattering has proven success-
ful for extracting phonon dispersions and investigating
momentum-dependent phonon dynamics [48, 49]. How-
ever, interpreting measurements of the diffuse scattering
intensity
I
(
q
,t
) is challenging due to its rich informa-
tion content. For example, mode-resolved phonon con-
tributions cannot be recovered straightforwardly from
the measured signal, hindering the identification of the
dominant phonon modes governing ultrafast structural
dynamics and anharmonic phonon processes. Here we
demonstrate the inverse process of reconstructing the ex-
perimental diffuse scattering signal from the computed
mode-resolved phonon populations.
We write the scattering intensity at wavevector
q
as [46]
I
(
q
,t
)
∝
e
−
2
M
∑
ν
2
π
2
(2
N
ν
q
(
t
) + 1)
~
(
q
·
ν
q
)
2
m
a
ω
ν
q
,
(5)
where
M
is the Debye-Waller factor, and compute
the ultrafast diffuse scattering signal as the difference
between the time-dependent and equilibrium scattering
intensities, ∆
I
(
q
,t
)
/I
=
I
(
q
,t
)
−
I
(
q
,
−∞
)
I
(
q
,
−∞
)
. Our computed
diffuse scattering results for graphene are plotted in
reciprocal space in Fig. 4c at 0.5 ps and 5 ps time delays.
At 0.5 ps, we find peaks of ∆
I/I
near the BZ corners
due to optical phonons generated from intervalley
scattering during the initial sub-ps carrier cooling.
Small-
q
optical phonons are also generated extensively,
as discussed above, but are not visible in Fig. 4c as the
denominator in the differential diffuse scattering ∆
I/I
6
FIG. 5:
Time- and energy-dependent excess phonon
populations, ∆
N
(
E,t
), after coherent excitation of LO
phonons at time zero.
(the thermal equilibrium intensity
I
(
q
,
−∞
)) is large
at small-
q
due to the large equilibrium population of
acoustic phonons. The diffuse scattering result at 5 ps
captures the dissipation of the optical phonons excess
energy to the lower acoustic branches, whereby phonons
are generated throughout the BZ. While the measured
diffuse scattering signal cannot provide information on
the mode-resolved dynamics, our simulations can fill this
gap, providing detailed information on mode-dependent
scattering processes in momentum-space and their
contributions to diffuse scattering.
Ultrafast dynamics following phonon excitation.
Coherent excitation of selected phonon modes has be-
come an important approach for investigating and ma-
nipulating materials properties [50–52]. To demonstrate
the flexibility of our numerical framework, we carry out
a simulated experiment in which phonons are excited at
time zero while electrons and holes are initially kept in
their equilibrium room temperature distributions. Af-
ter populating a large excess of LO phonons near the
BZ center at time zero, we simulate the subsequent
phonon dynamics by solving the coupled electron and
phonon BTEs, as above but for different initial condi-
tions. The computed time-domain excess phonon pop-
ulations ∆
N
(
E,t
) are plotted in Fig. 5. After 1 ps,
LA phonon modes with
∼
100 meV energy are gener-
ated through ph-ph processes, and within 2 ps both LA
and TA modes with 80
−
110 meV energy absorb most of
the excess energy from the optical modes. The flexural
ZA phonons come into play only after 3 ps. Through-
out the simulation, low-energy electron-hole pairs are ex-
cited near the Dirac cone via ph-e processes, resulting
in a modest ultrasonic attenuation of the phonon dy-
namics. Overall, the excess energy initially imparted to
the excited phonons mainly dissipates through slow ph-
ph processes, the entire phonon relaxation taking tens of
ps. Our approach is uniquely able to shed light on these
longer timescales characteristic of phonon dynamics.
DISCUSSION
A key advantage of our approach is the possibility,
rather unique among existing first-principles methods for
ultrafast dynamics, to validate the interactions against
experiments, thus guaranteeing the quantitative accu-
racy of our simulated dynamics. We validate the e-ph
and ph-ph interactions employed in our calculations, re-
spectively, by computing electrical and heat transport
properties of graphene (see Methods). We obtain a room
temperature electron mobility of 186
,
000 cm
2
/Vs, in ex-
cellent agreement with experimental results in suspended
graphene [53, 54]. We also compute the thermal conduc-
tivity in the single-mode relaxation time approximation
(RTA) of the BTE, obtaining a room temperature value
of 482 W/mK in excellent agreement with previous RTA
calculations [55]; this result implies that the full solution
of the BTE (beyond the RTA) would give a thermal con-
ductivity consistent with experiment [55]. These results
show that the e-ph and ph-ph interactions employed in
our ultrafast dynamics are precise, thus our computed
timescales are expected to be quantitatively accurate.
CONCLUSIONS
In summary, we have shown a versatile numerical
framework for modeling the ultrafast coupled dynamics
of electrons and phonons. We demonstrated its accuracy
and broad applicability through simulations of pump-
probe spectroscopy, X-ray diffuse scattering, and struc-
tural and phonon dynamics. Our results shed light on
e-ph and ph-ph scattering processes governing ultrafast
dynamics in graphene, and provide valuable information
for the design of ultrafast electronic and optical devices.
We plan to make the numerical method available in a
future release of
Perturbo
[36] to equip the commu-
nity with a novel tool for simulating and interpreting
ultrafast time-domain experiments. Future extensions
will aim to treat explicitly the light excitation pulse and
the electron spin to unravel the intertwined nonequilib-
rium dynamics of electronic, structural and spin degrees
of freedom. Taken together, our first-principles approach
demonstrates a paradigm shift in computing the ultrafast
electron and atomic vibrational dynamics, bridging the
gap between theory and experiment and enabling quan-
titative predictions of ultrafast phenomena in materials.
7
METHODS
Density functional theory.
We carry out first-
principles density functional theory (DFT) calculations
on graphene with a relaxed in-plane lattice constant of
2.45
̊
A. The graphene sheet is separated from its periodic
replicas by a 9
̊
A vacuum. The ground-state electronic
structure is computed using the
Quantum ESPRESSO
code in the local density approximation (LDA) of DFT.
We use a norm-conserving pseudopotential, a 90 Ry
plane-wave kinetic energy cutoff, and a Methfessel-
Paxton smearing of 0.02 Ry. The ground-state charge
density is obtained using a 36
×
36
×
1
k
-point grid,
following which a non-self-consistent calculation is em-
ployed to obtain the Kohn-Sham eigenvalues and wave
functions on a 12
×
12
×
1
k
-point grid. To construct
maximally localized Wannier functions (WFs) with the
Wannier90
code [56], the Kohn-Sham wave functions
are first projected onto atomic
p
z
orbitals on each atom
and
sp
2
orbitals on every other atom [57, 58], for a total
of 5 wannierized bands. The WF spread is then min-
imized, and the relevant energy windows are adjusted
until the interpolated bandstructure can smoothly repro-
duce the LDA result within
∼
10 meV throughout the BZ.
Electron-phonon scattering and correction to the
phonon dispersions.
We use density functional per-
turbation theory [59] (DFPT) to compute lattice dynam-
ical properties and the e-ph perturbation potentials, and
then form the e-ph matrix elements
g
nn
′
ν
(
k
,
q
) on coarse
12
×
12
×
1
k
-point and
q
-point BZ grids using
Per-
turbo
[36]. Here and below, the e-ph matrix elements
g
nn
′
ν
(
k
,
q
) quantify the probability amplitude for an elec-
tron in Bloch state
|
n
k
〉
with energy
E
n
k
to scatter into a
final state
|
n
′
k
+
q
〉
with energy
E
n
′
k
+
q
due to emission
or absorption of a phonon with branch index
ν
, wavevec-
tor
q
, and energy
~
ω
ν
q
. The electron and phonon en-
ergies and the e-ph matrix elements are interpolated on
fine grids using WFs with
Perturbo
[36]. The average
e-ph coupling strengths
g
2
ν
used in Fig. 3a are obtained
as
g
2
ν
=
∑
nn
′
kq
|
g
nn
′
ν
(
k
,
q
)
|
2
, where the summation in-
cludes electronic states near the Dirac point and phonon
wavevectors in the entire BZ.
To account for the Kohn anomaly near
q
= K [60], we
interpolate
E
n
k
and
g
nn
′
ν
(
k
,
q
) on a random BZ grid of
10
4
points, and obtain graphene phonon dispersions
ω
ν
q
that include a previously proposed GW correction [61],
ω
ν
q
=
√
√
√
√
B
GW
q
m
+
4
γ
GW
q
N
′
k
∑
k
|
g
nn
′
ν
(
k
,
q
)
|
2
E
π
k
−
E
π
∗
k
+
q
(6)
where
π
(
π
∗
) labels the occupied (emtpy)
π
band,
N
′
k
=
10
4
is the total number of
k
-points in the random grid,
and
B
GW
q
is a parameter employed to converge the LO
and TO phonon energies at K. The rescaling factor
γ
GW
q
is computed using [62]
γ
GW
q
= 1 + (
γ
GW
−
1)
1
2
erfc
(
|
q
−
K
n
|
a
0
2
π
−
0
.
2
0
.
05
)
(7)
with
a
0
the graphene lattice constant,
γ
GW
a constant
equal to 1.61, and K
n
the nearest vector to
q
among
those equivalent to K.
To verify the convergence of the
q
-point grid and com-
pute the electrical mobility, we use
Perturbo
[36] to
calculate the state-dependent e-ph scattering rates Γ
n
k
within lowest-order perturbation theory [21],
Γ
n
k
=
2
π
~
1
N
q
∑
n
′
q
ν
|
g
nn
′
ν
(
k
,
q
)
|
2
×
[(
N
ν
q
+ 1
−
f
n
′
k
+
q
)
δ
(
E
n
k
−
~
ω
ν
q
−
E
n
′
k
+
q
)
+ (
N
ν
q
+
f
n
′
k
+
q
)
δ
(
E
n
k
+
~
ω
ν
q
−
E
n
′
k
+
q
)]
,
(8)
where
f
n
k
and
N
ν
q
are the electron and phonon equi-
librium occupations at 300 K. Convergence is achieved
for a grid of 420
×
420
×
1
q
-points (see Supplementary
Figure S1) using a Gaussian broadening of 20 meV to
approximate the
δ
functions in Eq. (8). The same grid
is employed for all ultrafast dynamics calculations.
Phonon-phonon scattering.
We use the temperature-
dependent effective potential (TDEP) method [63–65] to
obtain the interatomic force constants and ph-ph ma-
trix elements Φ
νν
′
ν
′′
(
q
,
q
′
). The latter are the probabil-
ity amplitudes for three-phonon scattering processes in
which the phonon state
|
ν
q
〉
with energy
~
ω
ν
q
scatters
to the states
|
ν
′′
q
+
q
′
〉
and
|
ν
′
q
′
〉
, or the inverse process
in which the last two phonons combine to generate
|
ν
q
〉
.
For the TDEP calculations, we prepare a number of 12
×
12
×
1 (288 atom) supercells with random thermal dis-
placements corresponding to a canonical ensemble at 300
K. To validate the ph-ph interactions and compute the
thermal conductivity, we compute the phonon scattering
rates Λ
ν
q
as [66]
Λ
ν
q
=
18
π
~
2
1
N
2
q
∑
ν
′
ν
′′
q
′
q
′′
|
Φ
νν
′
ν
′′
(
q
′
,
q
′′
)
|
2
×
[(
N
ν
′
q
′
+
N
ν
′′
q
′′
+ 1)
δ
(
ω
ν
q
−
ω
ν
′
q
′
−
ω
ν
q
′′
)
+ 2(
N
ν
′
q
′
−
N
ν
′′
q
′′
)
δ
(
ω
ν
q
−
ω
ν
′
q
′
+
ω
ν
q
′′
)]
(9)
using phonon equilibrium occupations at 300 K. We
compute and converge the ph-ph scattering rates (see
Supplementary Figure S2) with an in-house version of
Perturbo
[36]. The ph-ph scattering rates converge
for grids equal to or finer than 210
×
210
×
1
q
-points
using a 1 meV Gaussian broadening to approximate the
δ
functions in Eq. (9).
8
Coupled electron and phonon ultrafast dynamics.
We simulate the ultrafast dynamics of coupled electrons
and phonons due to e-ph and ph-ph scattering processes
in graphene in the absence of external fields. The time
evolution of the carrier and phonon distributions is ob-
tained by solving the coupled carrier and phonon BTEs:
∂f
n
k
(
t
)
∂t
=
−
2
π
~
1
N
q
∑
n
′
q
ν
|
g
nn
′
ν
(
k
,
q
)
|
2
{
δ
(
E
n
k
−
E
n
′
k
+
q
−
~
ω
ν
q
)[
f
n
k
(1
−
f
n
′
k
+
q
)(
N
ν
q
+ 1)
−
f
n
′
k
+
q
(1
−
f
n
k
)
N
ν
q
]
+
δ
(
E
n
k
−
E
n
′
k
+
q
+
~
ω
ν
q
)[
f
n
k
(1
−
f
n
′
k
+
q
)
N
ν
q
−
f
n
′
k
+
q
(1
−
f
n
k
)(
N
ν
q
+ 1)]
}
∂N
ν
q
(
t
)
∂t
=
−
4
π
~
1
N
k
∑
n
′
k
ν
|
g
nn
′
ν
(
k
,
q
)
|
2
δ
(
E
n
k
−
E
n
′
k
+
q
+
~
ω
ν
q
[
f
n
k
(1
−
f
n
′
k
+
q
)
N
ν
q
−
f
n
′
k
+
q
(1
−
f
n
k
)(
N
ν
q
+ 1)]
−
18
π
~
2
1
N
2
q
∑
ν
′
ν
′′
q
′
q
′′
|
Φ
νν
′
ν
′′
(
q
′
,
q
′′
)
|
2
{
δ
(
ω
ν
q
−
ω
ν
′
q
′
−
ω
ν
′′
q
′′
)[
N
ν
q
(
N
ν
′
q
′
+ 1)(
N
ν
′′
q
′′
+ 1)
−
N
ν
′
q
′
N
ν
′′
q
′′
(
N
ν
q
+ 1)]
+
δ
(
ω
ν
q
+
ω
ν
′
q
′
−
ω
ν
′′
q
′′
)[
N
ν
q
N
ν
′
q
′
(
N
ν
′′
q
′′
+ 1)
−
N
ν
′
q
′
(
N
ν
′′
q
′′
+ 1)(
N
ν
q
+ 1)]
}
,
(10)
which account for e-ph scattering for electrons, e-ph
scattering for phonons (denoted above as ph-e scatter-
ing) and ph-ph scattering for phonons. In Eq. (10),
f
n
k
(
t
) are time-dependent electron populations and
N
ν
q
(
t
) are time-dependent phonon populations. The
e-ph matrix elements
g
nn
′
ν
(
k
,
q
) and ph-ph matrix
elements Φ
νν
′
ν
′′
(
q
′
,
q
′′
) are computed using DFPT in
the ground state plus WF interpolation and TDEP,
respectively, as specified above.
We numerically solve Eq. (10) using the fourth-order
Runge-Kutta method with a time step of 2 fs, using
uniform fine BZ grids of 420
×
420
×
1
k
- and
q
-points
for all e-ph and ph-ph scattering processes. The ph-ph
scattering matrix elements are computed on a uniform
210
×
210
×
1
q
-point BZ grid and then Fourier
interpolated to a 420
×
420
×
1
q
-point grid. These
grids are sufficient to converge the scattering processes
and ultrafast dynamics.
To reduce the number of
scattering processes entering the BTEs, we use only
electronic states in the energy window of relevance for
our calculations, which restricts the number of
k
-points
to a few thousands, while the
q
-points span the entire
BZ. We additionally select the relevant ph-ph processes
by imposing energy conservation (within a few times
the broadening value) thereby dramatically reduc-
ing the number of ph-ph scattering processes, which
would otherwise be unmanageable for the fine
q
-point
grid employed.
The BZ-averaged energy-dependent
carrier populations
f
(
E,t
) and phonon populations
N
(
E,t
) at energy
E
, first used in Fig. 1, are obtained
respectively as
f
(
E,t
) =
∑
n
k
f
n
k
(
t
)
δ
(
E
n
k
E
) and
N
(
E,t
) =
∑
ν
q
N
ν
q
(
t
)
δ
(
~
ω
ν
q
E
) via tetrahedron inte-
gration. The average temperature of each phonon mode
at time
t
(see Fig. 3) is computed as
T
ν
(
t
) =
∑
q
T
ν
q
(
t
),
where the state-dependent temperature
T
ν
q
(
t
) is
obtained by inverting the Bose-Einstein occupation
formula,
N
ν
q
(
t
) = (
e
~
ω
ν
q
/k
B
T
ν
q
(
t
)
−
1)
−
1
, using the
phonon populations at each time step.
The phonon
populations follow a Bose-Einstein thermal distribution
at the computed temperature
T
ν
(
t
) for times greater
than 2 ps for all acoustic modes and the ZO mode, and
at times greater than 4
−
5 ps for the LO and TO optical
modes.
We develop a scheme combining MPI and
openMP
parallelization to efficiently compute the e-ph, ph-e and
ph-ph collision integrals at each time step. Briefly, (
k
,
q
)
pairs for e-ph scattering and (
q
,
q
′
,
q
′′
) triplets for ph-ph
scattering are distributed among MPI processes; using
openMP
parallelization, rapid on-node operations are
employed to compute the collision integrals within each
MPI process; the results are then added together to
form the collision integrals. This scheme leverages an
algorithm we previously employed to time-step the elec-
tron BTE [20, 36]; however, including ph-ph scattering
is dramatically more difficult due to the need to take
into account all
q
-point triplets involved in the ph-ph
processes. To accomplish this task, we developed a new
parallel implementation of the phonon BTE in this work.
Transport properties.
The electron mobility is com-
puted from a full solution of the linearized electron BTE
with the
Perturbo
code [36]. The thermal conductiv-
ity is computed using ph-ph matrix elements from TDEP
with an in-house implementation of the single-mode RTA
formula [55]. Both calculations use the converged BZ
grids given above.
DATA AVAILABILITY
The datasets generated and/or analyzed in the current
study are available from the corresponding author upon
reasonable request.
9
References
∗
E-mail: bmarco@caltech.edu
[1] Buzzi, M., F ̈orst, M., Mankowsky, R. & Cavalleri, A.
Probing dynamics in quantum materials with femtosec-
ond X-rays.
Nat. Rev. Mater.
3
, 299 (2018). URL
https:
//www.nature.com/articles/s41578-018-0024-9
.
[2] Najafi, E., Ivanov, V., Zewail, A. & Bernardi, M. Super-
diffusion of excited carriers in semiconductors.
Nat. Com-
mun.
8
, 15177 (2017). URL
https://www.nature.com/
articles/ncomms15177
.
[3] Z ̈urch, M.
et al.
Ultrafast carrier thermalization and trap-
ping in silicon-germanium alloy probed by extreme ultra-
violet transient absorption spectroscopy.
Struct. Dyn.
4
,
044029 (2017). URL
https://aca.scitation.org/doi/
abs/10.1063/1.4985056
.
[4] Na, M. X.
et al.
Direct determination of mode-projected
electron-phonon coupling in the time domain.
Sci-
ence
366
, 1231–1236 (2019). URL
https://science.
sciencemag.org/content/366/6470/1231
.
[5] Young, L.
et al.
Femtosecond electronic response of atoms
to ultra-intense X-rays.
Nature
466
, 56 (2010). URL
https://www.nature.com/articles/nature09177
.
[6] Fritz, D. M.
et al.
Ultrafast bond softening in bis-
muth: mapping a solid’s interatomic potential with X-
rays.
Science
315
, 633–636 (2007). URL
https://
science.sciencemag.org/content/315/5812/633
.
[7] Reis, D. A. & Lindenberg, A. M. Ultrafast X-ray scat-
tering in solids. In
Light Scattering in Solid IX
, 371–422
(Springer, 2006).
[8] Trigo, M.
et al.
Fourier-transform inelastic X-ray scat-
tering from time-and momentum-dependent phonon–
phonon correlations.
Nat. Phys.
9
, 790 (2013). URL
https://www.nature.com/articles/nphys2788
.
[9] Lindenberg, A.
et al.
Time-resolved X-ray diffraction
from coherent phonons during a laser-induced phase
transition.
Phys. Rev. Lett.
84
, 111 (2000).
URL
https://journals.aps.org/prl/abstract/10.1103/
PhysRevLett.84.111
.
[10] Zalden, P.
et al.
Femtosecond X-ray diffraction reveals a
liquid–liquid phase transition in phase-change materials.
Science
364
, 1062–1067 (2019). URL
https://science.
sciencemag.org/content/364/6445/1062
.
[11] De Jong, S.
et al.
Speed limit of the insulator–metal
transition in magnetite.
Nat. Mater.
12
, 882 (2013). URL
https://www.nature.com/articles/nmat3718
.
[12] Mankowsky, R., von Hoegen, A., F ̈orst, M. & Cavalleri,
A.
Ultrafast reversal of the ferroelectric polariza-
tion.
Phys. Rev. Lett.
118
, 197601 (2017).
URL
https://journals.aps.org/prl/abstract/10.1103/
PhysRevLett.118.197601
.
[13] Mitrano, M.
et al.
Possible light-induced supercon-
ductivity in K
3
C
60
at high temperature.
Nature
530
,
461 (2016). URL
https://www.nature.com/articles/
nature16522
.
[14] Wilson, R., Feser, J. P., Hohensee, G. T. & Cahill,
D. G. Two-channel model for nonequilibrium thermal
transport in pump-probe experiments.
Phys. Rev. B
88
,
144305 (2013). URL
https://journals.aps.org/prb/
abstract/10.1103/PhysRevB.88.144305
.
[15] Haug, H. & Jauho, A.-P.
Quantum kinetics in transport
and optics of semiconductors
, vol. 2 (Springer Science &
Business Media, Boston, 2008).
[16] Waldecker, L., Bertoni, R., Ernstorfer, R. & Vorberger,
J. Electron-phonon coupling and energy flow in a simple
metal beyond the two-temperature approximation.
Phys.
Rev. X
6
, 021003 (2016). URL
https://journals.aps.
org/prx/abstract/10.1103/PhysRevX.6.021003
.
[17] Murray, E.
et al.
Phonon dispersion relations and soften-
ing in photoexcited bismuth from first principles.
Phys.
Rev. B
75
, 184301 (2007). URL
https://journals.aps.
org/prb/abstract/10.1103/PhysRevB.75.184301
.
[18] Bernardi, M., Vigil-Fowler, D., Lischner, J., Neaton, J. B.
& Louie, S. G. Ab initio study of hot carriers in the first
picosecond after sunlight absorption in silicon.
Phys. Rev.
Lett.
112
, 257402 (2014). URL
http://link.aps.org/
doi/10.1103/PhysRevLett.112.257402
.
[19] Bernardi, M., Mustafa, J., Neaton, J. B. & Louie, S. G.
Theory and computation of hot carriers generated by sur-
face plasmon polaritons in noble metals.
Nat. Commun.
6
(2015). URL
http://dx.doi.org/10.1038/ncomms8044
.
[20] Jhalani, V. A., Zhou, J.-J. & Bernardi, M.
Ul-
trafast hot carrier dynamics in GaN and its impact
on the efficiency droop.
Nano Lett.
17
, 5012–5019
(2017). URL
https://pubs.acs.org/doi/abs/10.1021/
acs.nanolett.7b02212
.
[21] Bernardi, M. First-principles dynamics of electrons and
phonons.
Eur. Phys. J. B
89
, 239 (2016). URL
https:
//doi.org/10.1140/epjb/e2016-70399-4
.
[22] Sadasivam, S., Chan, M. K. & Darancet, P. Theory
of thermal relaxation of electrons in semiconduc-
tors.
Phys. Rev. Lett.
119
, 136602 (2017).
URL
https://journals.aps.org/prl/abstract/10.1103/
PhysRevLett.119.136602
.
[23] Yabana, K. & Bertsch, G.
Time-dependent local-
density approximation in real time.
Phys. Rev. B
54
,
4484 (1996).
URL
https://journals.aps.org/prb/
abstract/10.1103/PhysRevB.54.4484
.
[24] Castro, A., Marques, M. A. & Rubio, A. Propaga-
tors for the time-dependent Kohn–Sham equations.
J.
Chem. Phys.
121
, 3425–3433 (2004). URL
https://aip.
scitation.org/doi/10.1063/1.1774980
.
[25] Marques, M. A. & Gross, E. K. Time-dependent density
functional theory.
Annu. Rev. Phys. Chem.
55
, 427–455
(2004).
URL
https://www.annualreviews.org/doi/
abs/10.1146/annurev.physchem.55.091602.094449
.
[26] Nazeeruddin, M. K.
et al.
Combined experimental and
DFT-TDDFT computational study of photoelectrochem-
ical cell ruthenium sensitizers.
J. Am. Chem. Soc.
127
,
16835–16847 (2005). URL
https://pubs.acs.org/doi/
abs/10.1021/ja052467l
.
[27] Rozzi, C. A.
et al.
Quantum coherence controls
the charge separation in a prototypical artificial light-
harvesting system.
Nat. Commun.
4
, 1602 (2013). URL
https://www.nature.com/articles/ncomms2603
.
[28] Sottile, F.
et al.
TDDFT from molecules to solids: The
role of long-range interactions.
Int. J. Quantum Chem.
102
, 684–701 (2005). URL
https://onlinelibrary.
wiley.com/doi/full/10.1002/qua.20486
.
[29] Marques, M. A., Maitra, N. T., Nogueira, F. M., Gross,
E. K. & Rubio, A.
Fundamentals of Time-dependent
Density Functional Theory
, vol. 837 (Springer Science &
Business Media, Boston, 2012).
[30] Dreuw, A., Weisman, J. L. & Head-Gordon, M. Long-
range charge-transfer excited states in time-dependent
density functional theory require non-local exchange.
J.
10
Chem. Phys.
119
, 2943–2946 (2003). URL
https://aip.
scitation.org/doi/10.1063/1.1590951
.
[31] Elliott, P., Fuks, J. I., Rubio, A. & Maitra, N. T.
Universal dynamical steps in the exact time-dependent
exchange-correlation potential.
Phys. Rev. Lett.
109
,
266404 (2012). URL
https://link.aps.org/doi/10.
1103/PhysRevLett.109.266404
.
[32] Yang, J.
et al.
Simultaneous observation of nuclear
and electronic dynamics by ultrafast electron diffraction.
Science
368
, 885–889 (2020). URL
https://science.
sciencemag.org/content/368/6493/885
.
[33] Mahan, G. D.
Many-Particle Physics
(Springer Science
& Business Media, Boston, 2000).
[34] Sangalli, D. & Marini, A. Ultra-fast carriers relaxation
in bulk silicon following photo-excitation with a short
and polarized laser pulse.
Eur. Phys. Lett.
110
, 47004
(2015). URL
http://dx.doi.org/10.1209/0295-5075/
110/47004
.
[35] Mahan, G. D.
Condensed Matter in a Nutshell
(Prince-
ton University Press, Princeton, 2010).
[36] Zhou, J.-J.
et al.
Perturbo: a software package for ab
initio electron-phonon interactions, charge transport and
ultrafast dynamics.
arXiv: 2002.02045
(2020). URL
https://arxiv.org/abs/2002.02045
.
[37] Dawlaty, J. M., Shivaraman, S., Chandrashekhar, M.,
Rana, F. & Spencer, M. G. Measurement of ultrafast
carrier dynamics in epitaxial graphene.
Appl. Phys. Lett.
92
, 042116 (2008). URL
https://aip.scitation.org/
doi/abs/10.1063/1.2837539
.
[38] Winzer, T. & Mali ́c, E.
Impact of auger processes
on carrier dynamics in graphene.
Phys. Rev. B
85
,
241404 (2012). URL
https://journals.aps.org/prb/
abstract/10.1103/PhysRevB.85.241404
.
[39] Wang, H.
et al.
Ultrafast relaxation dynamics of hot op-
tical phonons in graphene.
Appl. Phys. Lett.
96
, 081917
(2010). URL
https://aip.scitation.org/doi/figure/
10.1063/1.3291615
.
[40] Sun, B., Zhou, Y. & Wu, M.
Dynamics of pho-
toexcited carriers in graphene.
Phys. Rev. B
85
,
125413 (2012). URL
https://journals.aps.org/prb/
abstract/10.1103/PhysRevB.85.125413
.
[41] Gierz, I.
et al.
Snapshots of non-equilibrium Dirac carrier
distributions in graphene.
Nat. Mater.
12
, 1119 (2013).
URL
https://www.nature.com/articles/nmat3757
.
[42] Breusing, M.
et al.
Ultrafast nonequilibrium carrier dy-
namics in a single graphene layer.
Phys. Rev. B
83
,
153410 (2011). URL
https://journals.aps.org/prb/
abstract/10.1103/PhysRevB.83.153410
.
[43] Brida, D.
et al.
Ultrafast collinear scattering and
carrier multiplication in graphene.
Nat. Commun.
4
,
1987 (2013). URL
https://www.nature.com/articles/
ncomms2987
.
[44] Lindenberg, A. M., Johnson, S. L. & Reis, D. A.
Visualization of atomic-scale motions in materi-
als via femtosecond X-ray scattering techniques.
Annu. Rev. Mater. Res.
47
, 425–449 (2017).
URL
https://www.annualreviews.org/doi/abs/10.1146/
annurev-matsci-070616-124152
.
[45] Mohammed, O. F., Yang, D.-S., Pal, S. K. & Zewail,
A. H. 4D scanning ultrafast electron microscopy: Visu-
alization of materials surface dynamics.
J. Am. Chem.
Soc.
133
, 7708–7711 (2011). URL
http://dx.doi.org/
10.1021/ja2031322
.
[46] Warren, B. E.
X-ray Diffraction
(Dover Publications,
New York, 1990).
[47] Trigo, M.
et al.
Imaging nonequilibrium atomic vibra-
tions with X-ray diffuse scattering.
Phys. Rev. B
82
,
235205 (2010). URL
https://journals.aps.org/prb/
abstract/10.1103/PhysRevB.82.235205
.
[48] Holt, M.
et al.
Determination of phonon dispersions
from X-ray transmission scattering:
The example
of silicon.
Phys. Rev. Lett.
83
, 3317 (1999).
URL
https://journals.aps.org/prl/abstract/10.1103/
PhysRevLett.83.3317
.
[49] Stern, M. J.
et al.
Mapping momentum-dependent
electron-phonon coupling and nonequilibrium phonon
dynamics with ultrafast electron diffuse scattering.
Phys.
Rev. B
97
, 165416 (2018). URL
https://journals.aps.
org/prb/abstract/10.1103/PhysRevB.97.165416
.
[50] Rini, M.
et al.
Control of the electronic phase of a
manganite by mode-selective vibrational excitation.
Na-
ture
449
, 72 (2007). URL
https://www.nature.com/
articles/nature06119
.
[51] Graves, C.
et al.
Nanoscale spin reversal by non-local
angular momentum transfer following ultrafast laser ex-
citation in ferrimagnetic GdFeCo.
Nat. Mater.
12
,
293 (2013). URL
https://www.nature.com/articles/
nmat3597
.
[52] Fausti, D.
et al.
Light-induced superconductivity in a
stripe-ordered cuprate.
Science
331
, 189–191 (2011).
URL
https://science.sciencemag.org/content/331/
6014/189
.
[53] Du, X., Skachko, I., Barker, A. & Andrei, E. Y. Ap-
proaching ballistic transport in suspended graphene.
Nat. Nanotech.
3
, 491–495 (2008). URL
https://www.
nature.com/articles/nnano.2008.199
.
[54] Bolotin, K. I.
et al.
Ultrahigh electron mobility in
suspended graphene.
Solid State Commun.
146
, 351–
355 (2008).
URL
https://www.sciencedirect.com/
science/article/pii/S0038109808001178
.
[55] Fugallo, G.
et al.
Thermal conductivity of graphene
and graphite: collective excitations and mean free paths.
Nano Lett.
14
, 6109–6114 (2014). URL
https://pubs.
acs.org/doi/10.1021/nl502059f
.
[56] Mostofi, A. A.
et al.
wannier90: A tool for obtain-
ing maximally-localised wannier functions.
Com-
put. Phys. Commun.
178
, 685–699 (2008).
URL
https://www.sciencedirect.com/science/article/
pii/S0010465507004936
.
[57] Marzari, N., Mostofi, A. A., Yates, J. R., Souza, I.
& Vanderbilt, D. Maximally localized wannier func-
tions: Theory and applications.
Rev. Mod. Phy.
84
,
1419 (2012).
URL
https://journals.aps.org/rmp/
abstract/10.1103/RevModPhys.84.1419
.
[58] Jung, J. & MacDonald, A. H.
Tight-binding
model for graphene
π
-bands from maximally lo-
calized wannier functions.
Phys. Rev. B
87
,
195450 (2013). URL
https://journals.aps.org/prb/
abstract/10.1103/PhysRevB.87.195450
.
[59] Baroni, S., De Gironcoli, S., Dal Corso, A. & Gian-
nozzi, P. Phonons and related crystal properties from
density-functional perturbation theory.
Rev. Mod. Phy.
73
, 515 (2001). URL
https://journals.aps.org/rmp/
abstract/10.1103/RevModPhys.73.515
.
[60] Piscanec, S., Lazzeri, M., Mauri, F., Ferrari, A.
& Robertson, J.
Kohn anomalies and electron-
phonon interactions in graphite.
Phys. Rev. Lett.
93
,
11
185503 (2004). URL
https://journals.aps.org/prl/
abstract/10.1103/PhysRevLett.93.185503
.
[61] Lazzeri, M., Attaccalite, C., Wirtz, L. & Mauri, F.
Impact of the electron-electron correlation on phonon
dispersion: Failure of LDA and GGA DFT function-
als in graphene and graphite.
Phys. Rev. B
78
,
081406 (2008). URL
https://journals.aps.org/prb/
abstract/10.1103/PhysRevB.78.081406
.
[62] Venezuela, P., Lazzeri, M. & Mauri, F. Theory of double-
resonant raman spectra in graphene: Intensity and line
shape of defect-induced and two-phonon bands.
Phys.
Rev. B
84
, 035433 (2011). URL
https://journals.aps.
org/prb/abstract/10.1103/PhysRevB.84.035433
.
[63] Hellman, O., Abrikosov, I. & Simak, S. Lattice dynamics
of anharmonic solids from first principles.
Phys. Rev. B
84
, 180301 (2011). URL
https://journals.aps.org/
prb/abstract/10.1103/PhysRevB.84.180301
.
[64] Hellman, O., Steneteg, P., Abrikosov, I. A. & Simak,
S. I. Temperature dependent effective potential method
for accurate free energy calculations of solids.
Phys. Rev.
B
87
, 104111 (2013). URL
https://journals.aps.org/
prb/abstract/10.1103/PhysRevB.87.104111
.
[65] Hellman, O. & Abrikosov, I. A.
Temperature-
dependent effective third-order interatomic force con-
stants from first principles.
Phys. Rev. B
88
,
144301 (2013). URL
https://journals.aps.org/prb/
abstract/10.1103/PhysRevB.88.144301
.
[66] Debernardi, A., Baroni, S. & Molinari, E.
Anhar-
monic phonon lifetimes in semiconductors from density-
functional perturbation theory.
Phys. Rev. Lett.
75
,
1819 (1995).
URL
https://journals.aps.org/prl/
abstract/10.1103/PhysRevLett.75.1819
.
Acknowledgments
The authors thank Jin-Jian Zhou for fruitful discussions.
X.T. thanks the Resnick Sustainability Institute at the
California Institute of Technology for fellowship support.
This work was partially supported by the National
Science Foundation under Grant No. DMR-1750613,
which provided for theory development, and by the
Department of Energy under Grant No. DE-SC0019166,
which provided for numerical calculations and code de-
velopment. This research used resources of the National
Energy Research Scientific Computing Center, a DOE
Office of Science User Facility supported by the Office
of Science of the U.S. Department of Energy under
Contract No. DE-AC02-05CH11231.
Author Contributions
M.B. conceived the research. X.T. developed the com-
putational codes and carried out the calculations. All
authors analyzed the results and wrote the manuscript.
Additional Information
Supplementary Information
accompany this paper
(at URL to be added by the editorial office).
Competing financial interests:
The authors declare
no competing financial interests.