of 10
PHYSICAL REVIEW RESEARCH
3
, 023072 (2021)
Toward precise simulations of the coupled ultrafast dynamics of
electrons and atomic vibrations in materials
Xiao Tong and Marco Bernardi
*
Department of Applied Physics and Materials Science, California Institute of Technology, Pasadena, California 91125, USA
(Received 13 October 2020; accepted 6 April 2021; published 26 April 2021)
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 subpicosecond 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. Combining 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 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.
DOI:
10.1103/PhysRevResearch.3.023072
I. INTRODUCTION
Time-domain spectroscopies open a window on the
dynamics of electrons, phonons, and various elementary ex-
citations, probing matter on timescales characteristic of the
interactions among its constituents. In particular, recently de-
veloped 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 theoretical
tools. The latter should ideally be able to predict or interpret
ultrafast measurements and shed light on the underlying mi-
croscopic processes and coupling between different degrees
of freedom. A common scenario is the coupled dynamics of
excited electrons and phonons, which governs a wide range of
phenomena such as carrier and lattice equilibration, photoe-
mission and the associated linewidths, structural transitions,
and superconductivity. Heuristic approaches such as two-
temperature models (TTMs) [
14
] or kinetic equations with
adjustable interactions [
15
] are routinely adopted to study
aspects of nonequilibrium electron and phonon dynamics. In
particular, TTMs can successfully model the dynamics of ma-
terials excited above the damage threshold [
16
,
17
] and treat
*
bmarco@caltech.edu
Published by the American Physical Society under the terms of the
Creative Commons Attribution 4.0 International
license. Further
distribution of this work must maintain attribution to the author(s)
and the published article’s title, journal citation, and DOI.
electron interactions in that regime in conjunction with molec-
ular dynamics [
18
,
19
]. However, at lower excitations and in
the context of modeling pump-probe experiments, TTMs and
related approaches are not geared toward quantitative pre-
dictions [
20
] and are mainly useful as tools to qualitatively
interpret experimental results.
First-principles computational methods based on density
functional theory (DFT) and related techniques have made
great strides in modeling electron and phonon interactions and
nonequilibrium dynamics [
21
26
]. An important approach is
real-time time-dependent DFT (rt-TDDFT) [
27
29
], which
propagates in time the electronic Kohn-Sham equations while
treating atomic motions using Ehrenfest forces. While rt-
TDDFT has been widely successful for studying electron
dynamics in molecules [
30
,
31
], important open challenges re-
main, including reaching simulation times longer than
1ps,
extracting phonon-related quantities when treating periodic
systems (crystals) [
32
], and better understanding, improving,
and validating the accuracy of exchange-correlation function-
als governing the interactions of electrons and nuclei [
33
35
].
Nonadiabatic molecular dynamics is also a valuable
method to model ultrafast dynamics in molecules [
36
] and
solids [
37
].
A second family of approaches—on which this work
hinges—employs perturbation theory to compute the ma-
trix elements for electron and phonon interactions [
38
], and
combines them with the real-time Boltzmann transport equa-
tion (rt-BTE) [
24
,
39
] or quantum kinetic equations [
40
]
to investigate ultrafast electron dynamics. 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 Boltzmann
transport equations (BTEs) for coupled electrons and phonons
while fully taking into account their interactions has remained
2643-1564/2021/3(2)/023072(10)
023072-1
Published by the American Physical Society
XIAO TONG AND MARCO BERNARDI
PHYSICAL REVIEW RESEARCH
3
, 023072 (2021)
an elusive goal due to computational cost and the need for
complex algorithms and workflows.
Here we show explicit simulations of the coupled dynamics
of nonequilibrium electrons and phonons for timescales up to
tens of picoseconds, using accurate electron-phonon (
e
-ph)
and phonon-phonon (ph-ph) interactions validated through
transport calculations. We time-step the electron and phonon
BTEs—a large set of coupled integro-differential equations—
using a parallel numerical approach to efficiently compute the
relevant collision integrals and obtain time-dependent electron
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, compute spectroscopic signals
such as transient absorption and diffuse x-ray scattering, and
analyze the dominant scattering channels, identifying the slow
rise of flexural 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.
II. RESULTS
A. Numerical approach
We solve the coupled electron and phonon BTEs [
41
]in
a homogeneously excited region of a material using a fourth-
order Runge-Kutta algorithm. At each time
t
, we propagate
the electronic populations
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 wavevec-
tor. Using an in-house modified version of our
PERTURBO
code [
42
], we solve the BTEs as a set of integro-differential
equations:
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 electron-phonon (
e
-
ph) interactions in the electron BTE, while
I
ph-
e
and
I
ph-ph
are, respectively, the collision integrals for phonon-electron
(ph-
e
) and phonon-phonon (ph-ph) interactions in the phonon
BTE. The collision integrals are obtained using
e
-ph matrix
elements computed in equilibrium at 0 K and ph-ph matrix
elements computed in equilibrium at 300 K. Detailed expres-
sions for these integrals are given in Appendix
B
.
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). A typical calculation with our scheme involves
time-stepping
10
6
coupled integro-differential equations,
each with 10
7
–10
9
scattering processes in the collision inte-
grals, a formidable computational challenge addressed with
an efficient algorithm combining message-passing interface
(MPI) and multi-threading (
O
pen
MP
) parallelizations. As
memory effects are not included, our BTEs amount to a
Markovian dynamics of the density matrix in which off-
diagonal coherences, which are typically short lived (
<
5 fs),
are neglected [
15
]. Using this scheme, we can reach simula-
tion times up to
100 ps, with a femtosecond time step.
B. Relaxation of excited carriers
In a first set of simulations, we model photoexcited elec-
tron 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 ex-
cited carriers within 10–50 fs [
43
,
44
]. 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 [
45
47
]. We model
the carrier dynamics in graphene after this initial thermaliza-
tion by setting the electron and hole populations at time zero
to Fermi-Dirac distributions at 4000 K temperature, while set-
ting the initial phonon populations to their equilibrium value
at 300 K.
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 Appendix
B
). In the first 100 fs, the hot electron and
hole distributions become narrower in energy as the carriers
cool to the Dirac point in about 300 fs. This rapid cooling
process is achieved by emitting optical phonons [
48
] through
intravalley and intervalley
e
-ph scattering. As a result, the op-
tical phonons generated within 1 ps possess wavevectors close
to the BZ

point for intravalley and K point for intervalley
scattering.
Although most of the energy stored in the excited carriers
is dissipated in the first 300 fs, surprisingly it takes more than
5 ps for the carriers to fully equilibrate, with negligible popu-
lation 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 allowing
optical phonons generated during the rapid carrier cooling
to linger for
>
10 ps, heating back the carriers via phonon
absorption.
C. Phonon equilibration
Accessing the 10–100 ps timescale characteristic of
phonon dynamics is an open challenge for existing first-
principles simulations of ultrafast dynamics. In Fig.
2
,we
analyze the long-lived equilibration of optical phonons gener-
ated by the excited carriers. The optical phonons first decay
within
5 ps to transverse acoustic (TA) and longitudi-
nal acoustic (LA) phonons, which then emit lower-energy
acoustic modes through ph-ph scattering processes, ultimately
equilibrating on a 20–100 ps timescale (not shown).
Our calculations can directly identify the main phonon
modes and relaxation pathways from the time-dependent
023072-2
TOWARD PRECISE SIMULATIONS OF THE COUPLED ...
PHYSICAL REVIEW RESEARCH
3
, 023072 (2021)
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 4000 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 Brillouin-zone-averaged carrier or phonon concentrations as a function of energy (in arbitrary units).
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 and 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 pop-
ulations and effective temperatures, taking into account mode
crossing as discussed in Appendix
B
.InFig.
3(a)
,weplot
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.
3(a)
are the mode-resolved
e
-ph coupling strengths,
g
2
ν
(see
Appendix
B
). 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, high-
lighting 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 subpi-
cosecond carrier cooling.
Although at short times the phonon distributions are still
nonthermal, we find that after 2–5 ps all phonon populations
are thermal and well approximated by hot Bose-Einstein dis-
tributions. The mode-resolved effective phonon temperatures
are computed as BZ averages of state-dependent temperatures
(see Appendix
A
) and shown in Fig.
3(b)
. We find that the
effective phonon temperatures are strongly mode dependent
throughout the simulation, their trend mirroring the respective
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- and 20-ps delay times after the initial electronic excitation.
023072-3
XIAO TONG AND MARCO BERNARDI
PHYSICAL REVIEW RESEARCH
3
, 023072 (2021)
FIG. 3. Mode-resolved phonon populations and effective tem-
peratures. (a) Time-dependent excess phonon populations

N
ν
(
t
),
with
e
-ph coupling strengths
g
2
ν
(in arbitrary 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.
excess phonon populations in Fig.
3(a)
. 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 equilibration
cascade, with timescales spanning several orders of mag-
nitudes, from femtosecond for carrier cooling through
e
-ph
interactions to picosecond for phonon down-conversion to
equilibration via ph-ph processes over tens of picoseconds.
The microscopic details of the scattering processes, accessed
easily in our approach due to its formulation in momentum
space, are out of reach for existing first-principles simulations.
As shown next, the time-dependent electron and phonon pop-
ulations further allow us to reliably simulate various ultrafast
spectroscopies employed as probes of electron and nuclear
dynamics.
D. Pump-probe transient absorption
Time-resolved differential transmission measurements are
widely employed to shed light on nonequilibrium carrier
dynamics and investigate the timescales associated with elec-
tronic processes. We build on the methods shown above to
compute transient absorption at optical wavelengths.
FIG. 4. Simulated ultrafast spectroscopies and structural snap-
shots in graphene. (a) Normalized differential transmission

T
/
T
0
as a function of pump-probe delay time. Our calculated results are
compared separately with experiments from Refs. [
49
,
50
]. (b) Struc-
tural snapshots at 0.5 and 5 ps after pump. (c) Simulated ultrafast
diffuse scattering pattern,

I
(
q
,
t
)
/
I
, at 0.5 and 5 ps after pump,
shown together with the reciprocal lattice.
In graphene, at low energies within 1–2 eV of the Dirac
cone, the differential transmissivity

T
/
T
0
is proportional to
the transient carrier populations; in a single-particle picture
and neglecting the state dependence of the optical transition
dipoles [
49
],

T
(
E
)
T
0
≈−
a
0
[

f
e
(
E
2
)
+

f
h
(
E
2
)]
,
(2)
where
E
is the carrier energy,
a
0
is the absorption co-
efficient of single-layer graphene,
f
e
,
h
are BZ-averaged
electron or hole populations, and

f
e
,
h
(
E
,
t
)
=
f
e
,
h
(
E
,
t
)
f
e
,
h
(
E
,
−∞
) are the corresponding excess carrier populations
relative to equilibrium.
We carried out calculations of

T
/
T
0
and compared them
to two sets of experimental results. In the first simulation, we
set the initial electron temperature to 3600 K and probe the

T
/
T
0
over the energy window of relevance (0.6–0.9 eV) in
the experiment by Breusing
et al.
[
49
]. In the second simu-
lation, the initial electron temperature is set to 4000 K and

T
/
T
0
is computed at 900-nm wavelength to compare with
experiments by Brida
et al.
[
50
]. For both simulations, the
initial phonon temperature is 300 K and the calculated

T
/
T
0
is shown as a function of pump-probe delay time in Fig.
4(a)
.
023072-4
TOWARD PRECISE SIMULATIONS OF THE COUPLED ...
PHYSICAL REVIEW RESEARCH
3
, 023072 (2021)
In both cases, we obtain very good agreement between theory
and experiment.
The experimental result exhibits two distinct temporal re-
gions, which are more clearly seen in the experiments by
Brida
et al.
[
50
] [red data points in Fig.
4(a)
]. 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 picoseconds. Our simulated differential
transmission spectrum agrees well with experiment [
50
]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 flexural phonons dis-
cussed above. While previous work focused on the role of
e
-
e
processes [
50
], our results show the importance of
e
-ph
interactions at subpicosecond times.
E. Diffuse scattering and structural dynamics
Intense experimental efforts are focused on investigating
atomic motions and structural dynamics in the time do-
main [
51
]. Among other techniques, ultrafast electron [
52
]
and x-ray diffraction [
53
] have made great strides and are cur-
rently able to infer atomic displacement patterns and dominant
phonon modes governing the nonequilibrium structural dy-
namics [
51
,
54
]. Using the time-resolved phonon populations
for each vibrational mode and wavevector, we can simulate
and reconstruct the structural dynamics fully from first princi-
ples.
In a classical picture, the time-dependent displacement of
atom
a
can be decomposed into lattice waves with wavevector
q
and mode
ν
[
53
],
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 matrix and the am-
plitude
A
ν
q
can be expressed in terms of phonon populations
N
ν
q
as
A
ν
q
(
t
)
=
[2
N
ν
q
(
t
)
+
1] ̄
h
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 and 5 ps are visualized in Fig.
4(b)
. 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 excited through ph-ph processes and the
atomic vibrations become randomized, consistent with the
thermal phonon distributions we find at 5 ps. The quantity
measured experimentally in x-ray diffraction is the averaged
square atomic displacement,

u
2
a
=
ν
q
1
2

A
2
ν
q

a
, which we
show in the Supplemental Material [
55
], where the radius
of the blue sphere around each atom represents

u
2
a

.The
averaged atomic displacements increase monotonically with
time, reaching greater values as the system approaches ther-
mal equilibrium.
Thermal diffuse x-ray scattering has proven successful for
extracting phonon dispersions and investigating momentum-
dependent phonon dynamics [
56
,
57
]. However, interpreting
measurements of the diffuse scattering intensity
I
(
q
,
t
)is
challenging due to its rich information content. For exam-
ple, mode-resolved phonon contributions cannot be recovered
straightforwardly from the measured signal, hindering the
identification of the dominant phonon modes governing ultra-
fast structural dynamics and anharmonic phonon processes.
Here we demonstrate the inverse process of reconstructing
the experimental diffuse scattering signal from the computed
mode-resolved phonon populations.
We write the scattering intensity at wavevector
q
as [
53
]
I
(
q
,
t
)
e
2
M
ν
2
π
2
(2
N
ν
q
(
t
)
+
1) ̄
h
(
q
·

ν
q
)
2
m
a
ω
ν
q
,
(5)
where
M
is the Debye-Waller factor, and compute the ul-
trafast 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.
4(c)
at 0.5- 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 subpicosecond
carrier cooling. Small-
q
optical phonons are also generated
extensively, as discussed above, but are not visible in Fig.
4(c)
as the denominator in the differential diffuse scattering

I
/
I
[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.
F. Ultrafast dynamics following phonon excitation
Coherent excitation of selected phonon modes has become
an important approach for investigating and manipulating ma-
terial properties [
58
60
]. To demonstrate the flexibility of
our rt-BTE 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. After populating a large ex-
cess of LO phonons near the BZ center at time zero, we
simulate the subsequent phonon dynamics by solving the cou-
pled electron and phonon BTEs, as above but for different
initial conditions. The computed time-domain excess phonon
populations

N
(
E
,
t
) are plotted in Fig.
5
.After1ps,LA
phonon modes with
100 meV energy are generated 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
023072-5
XIAO TONG AND MARCO BERNARDI
PHYSICAL REVIEW RESEARCH
3
, 023072 (2021)
FIG. 5. Time- and energy-dependent excess phonon populations,

N
(
E
,
t
), after coherent excitation of LO phonons at time zero.
play only after 3 ps. Throughout the simulation, low-energy
electron-hole pairs are excited near the Dirac cone via ph-
e
processes, resulting in a modest ultrasonic attenuation of the
phonon dynamics. Overall, the excess energy initially im-
parted to the excited phonons mainly dissipates through slow
ph-ph processes, the entire phonon relaxation taking tens of
picoseconds. Our approach is uniquely able to shed light on
these longer timescales characteristic of phonon dynamics.
III. DISCUSSION
A key advantage of our rt-BTE approach is the possi-
bility, rather unique among existing first-principles methods
for ultrafast dynamics, to validate the interactions against
experiments, thus guaranteeing the quantitative accuracy of
our simulated dynamics. We validate the
e
-ph and ph-ph
interactions employed in our calculations, respectively, by
computing electrical and heat transport properties of graphene
(see Appendix
A
). We obtain a room-temperature electron
mobility of 186000 cm
2
/
V s, in excellent agreement with
experimental results in suspended graphene [
61
,
62
]. We also
compute the thermal conductivity in the single-mode relax-
ation time approximation (RTA) of the BTE, obtaining a
room-temperature value of 482 W
/
mK in excellent agreement
with previous RTA calculations [
63
]; this result implies that
the full solution of the BTE (beyond the RTA) would give a
thermal conductivity consistent with experiment [
63
]. These
results show that the interactions employed in our ultrafast dy-
namics are precise, so our computed timescales are expected
to be quantitatively accurate.
Finally, note that our method treats
e
-ph and ph-ph in-
teractions perturbatively, using values for these interactions
computed on the equilibrium structure. Therefore our ap-
proach is most suitable to model the low-excitation regime in
which the atomic positions remain fairly close to their equi-
librium values; for example, in our calculations the phonon
temperatures do not rise significantly above the initial 300 K
temperature. Regimes involving intense excitation or phase
transitions, atomic diffusion, or material damage require dif-
ferent treatments. Future work will address explicitly the
issue of updating the
e
-ph and ph-ph interactions during the
nonequilibrium dynamics.
IV. CONCLUSION
In summary, we have shown a versatile numerical frame-
work, rt-BTE, for modeling the ultrafast coupled dynamics
of electrons and phonons. We demonstrated its accuracy and
broad applicability through simulations of pump-probe spec-
troscopy, x-ray diffuse scattering, and structural and phonon
dynamics. Our results shed light on
e
-ph and ph-ph scat-
tering 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
[
42
] to equip
the community with a 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 nonequilibrium dynamics of
electronic, structural, and spin degrees of freedom. Taken to-
gether, 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 quantitative predictions of ultrafast phenomena
in materials.
ACKNOWLEDGMENTS
The authors thank Jin-Jian Zhou for fruitful discussions.
X.T. thanks the Resnick Sustainability Institute at the Califor-
nia 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 calcula-
tions and code development. This research used resources of
the National Energy Research Scientific Computing Center
(NERSC), a US Department of Energy Office of Science User
Facility located at Lawrence Berkeley National Laboratory,
operated under Contract No. DE-AC02-05CH11231.
APPENDIX A: COMPUTATIONAL DETAILS
1. Density functional theory
We carry out first-principles DFT calculations on graphene
with a relaxed in-plane lattice constant of 2.45 Å. The
graphene sheet is separated from its periodic replicas by a
9-Å vacuum. The ground-state electronic structure is com-
puted 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
employed to obtain the Kohn-Sham eigenvalues and wave
functions on a 12
×
12
×
1
k
-point grid. To construct maxi-
mally localized Wannier functions (WFs) with the
WANNIER
90
023072-6
TOWARD PRECISE SIMULATIONS OF THE COUPLED ...
PHYSICAL REVIEW RESEARCH
3
, 023072 (2021)
code [
64
], the Kohn-Sham wave functions are first projected
onto atomic
p
z
orbitals on each atom and
sp
2
orbitals on every
other atom [
65
,
66
], for a total of 5 wannierized bands. The WF
spread is then minimized, and the relevant energy windows
are adjusted until the interpolated band structure can smoothly
reproduce the LDA result within
10 meV throughout the BZ.
2. Electron-phonon scattering and correction to the
phonon dispersions
We use density functional perturbation theory (DFPT) [
67
]
to compute lattice dynamical properties and the
e
-ph pertur-
bation 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
PERTURBO
[
42
]. Here and below, the
e
-ph matrix
elements
g
nn

ν
(
k
,
q
) quantify the probability amplitude for an
electron 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
ν
, wavevector
q
,
and energy ̄
h
ω
ν
q
. The electron and phonon energies and the
e
-ph matrix elements are interpolated on fine grids using WFs
with
PERTURBO
[
42
]. The average
e
-ph coupling strengths
g
2
ν
used in Fig.
3(a)
are obtained as
g
2
ν
=
nn

kq
|
g
nn

ν
(
k
,
q
)
|
2
,
where the summation includes electronic states near the Dirac
point and phonon wavevectors in the entire BZ.
To account for the Kohn anomaly near
q
=
K[
68
], 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 in-
clude a previously proposed
GW
correction [
69
],
ω
ν
q
=
B
GW
q
m
+
4
γ
GW
q
N

k
k
|
g
nn

ν
(
k
,
q
)
|
2
E
π
k
E
π
k
+
q
,
(A1)
where
π
(
π
) labels the occupied (empty)
π
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 en-
ergies at K. The rescaling factor
γ
GW
q
is computed using [
70
]
γ
GW
q
=
1
+
(
γ
GW
1)
1
2
erfc
(
|
q
K
n
|
a
0
2
π
0
.
2
0
.
05
)
(A2)
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 compute
the electrical mobility, we use
PERTURBO
[
42
] to calculate the
state-dependent
e
-ph scattering rates

n
k
within lowest-order
perturbation theory [
25
],

n
k
=
2
π
̄
h
1
N
q
n

q
ν
|
g
nn

ν
(
k
,
q
)
|
2
×
[(
N
ν
q
+
1
f
n

k
+
q
)
δ
(
E
n
k
̄
h
ω
ν
q
E
n

k
+
q
)
+
(
N
ν
q
+
f
n

k
+
q
)
δ
(
E
n
k
+
̄
h
ω
ν
q
E
n

k
+
q
)]
,
(A3)
where
f
n
k
and
N
ν
q
are the electron and phonon equilibrium
occupations at 300 K. Convergence is achieved for a grid of
420
×
420
×
1
q
-points (see the Supplemental Material [
55
])
using a Gaussian broadening of 20 meV to approximate the
δ
functions in Eq. (
A3
). The same grid is employed for all
ultrafast dynamics calculations.
3. Phonon-phonon scattering
We use the temperature-dependent effective potential
(TDEP) method [
71
73
] to obtain the interatomic force con-
stants and ph-ph matrix elements
νν

ν

(
q
,
q

). The latter
are the probability amplitudes for three-phonon scattering
processes in which the phonon state
|
ν
q

with energy ̄
h
ω
ν
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 [
74
]
ν
q
=
18
π
̄
h
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

)] (A4)
using phonon equilibrium occupations at 300 K. We compute
and converge the ph-ph scattering rates (see the Supplemental
Material [
55
]) with an in-house version of
PERTURBO
[
42
].
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. (
A4
).
4. Transport properties
The electron mobility is computed from a full solution of
the linearized electron BTE with the
PERTURBO
code [
42
].
The thermal conductivity is computed using ph-ph matrix
elements from TDEP with an in-house implementation of
the single-mode RTA formula [
63
]. Both calculations use the
converged BZ grids given above.
APPENDIX B: 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 obtained by solving
the coupled carrier and phonon BTEs:
f
n
k
(
t
)
t
=−
2
π
̄
h
1
N
q
n

q
ν
|
g
nn

ν
(
k
,
q
)
|
2
{
δ
(
E
n
k
E
n

k
+
q
̄
h
ω
ν
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
+
̄
h
ω
ν
q
)[
f
n
k
(1
f
n

k
+
q
)
N
ν
q
f
n

k
+
q
(1
f
n
k
)(
N
ν
q
+
1)]
}
,
023072-7
XIAO TONG AND MARCO BERNARDI
PHYSICAL REVIEW RESEARCH
3
, 023072 (2021)
N
ν
q
(
t
)
t
=−
4
π
̄
h
1
N
k
n

k
ν
|
g
nn

ν
(
k
,
q
)
|
2
δ
(
E
n
k
E
n

k
+
q
+
̄
h
ω
ν
q
[
f
n
k
(1
f
n

k
+
q
)
N
ν
q
f
n

k
+
q
(1
f
n
k
)(
N
ν
q
+
1)]
18
π
̄
h
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)]
}
,
(B1)
which account for
e
-ph scattering for electrons,
e
-ph scattering
for phonons (denoted above as ph-
e
scattering), and ph-ph
scattering for phonons. In Eq. (
B1
),
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, respec-
tively, as specified above.
We numerically solve Eq. (
B1
) 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 thousand,
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 dramat-
ically reducing 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
,firstusedinFig.
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
)
δ
( ̄
h
ω
ν
q
E
) via tetrahedron integration. Mode crossing has been
taken into account when computing all mode-dependent quan-
tities in this work. To reliably label the phonon modes, we
ensure that the first and second derivatives of the phonon
dispersions are smooth for each phonon mode, and label the
modes according to their character at the BZ center.
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 in-
verting the Bose-Einstein occupation formula,
N
ν
q
(
t
)
=
(
e
̄
h
ω
ν
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
O
pen
MP
par-
allelization 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
O
pen
MP
paral-
lelization, 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 electron BTE [
24
,
42
]; 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 parallel
implementation of the phonon BTE in this work.
[1] M. Buzzi, M. Först, R. Mankowsky, and A. Cavalleri, Probing
dynamics in quantum materials with femtosecond x-rays,
Nat.
Rev. Mater.
3
, 299 (2018)
.
[2] E. Najafi, V. Ivanov, A. Zewail, and M. Bernardi, Super-
diffusion of excited carriers in semiconductors,
Nat. Commun.
8
, 15177 (2017)
.
[3] M. Zürch, H.-T. Chang, P. M. Kraus, S. K. Cushing, L. J.
Borja, A. Gandman, C. J. Kaplan, M. H. Oh, J. S. Prell, D.
Prendergast
et al.
, Ultrafast carrier thermalization and trapping
in silicon-germanium alloy probed by extreme ultraviolet tran-
sient absorption spectroscopy,
Struct. Dyn.
4
, 044029 (2017)
.
[4] M. Na, A. K. Mills, F. Boschini, M. Michiardi, B. Nosarzewski,
R. P. Day, E. Razzoli, A. Sheyerman, M. Schneider, G. Levy
et al.
, Direct determination of mode-projected electron-phonon
coupling in the time domain,
Science
366
, 1231 (2019)
.
[5] L. Young, E. Kanter, B. Kraessig, Y. Li, A. March, S. Pratt,
R. Santra, S. Southworth, N. Rohringer, L. DiMauro
et al.
,
Femtosecond electronic response of atoms to ultra-intense
x-rays,
Nature (London)
466
, 56 (2010)
.
[6]D.M.Fritz,D.Reis,B.Adams,R.Akre,J.Arthur,C.Blome,
P. Bucksbaum, A. L. Cavalieri, S. Engemann, S. Fahy
et al.
,Ul-
trafast bond softening in bismuth: Mapping a solid’s interatomic
potential with x-rays,
Science
315
, 633 (2007)
.
[7] D. A. Reis and A. M. Lindenberg, Ultrafast x-ray scattering in
solids, in
Light Scattering in Solid IX
, edited by M. Cardona and
R. Merlin (Springer, Berlin, 2006), pp. 371–422.
[8] M. Trigo, M. Fuchs, J. Chen, M. Jiang, M. Cammarata, S.
Fahy, D. M. Fritz, K. Gaffney, S. Ghimire, A. Higginbotham
etal.
, Fourier-transform inelastic x-ray scattering from time-and
momentum-dependent phonon–phonon correlations,
Nat. Phys.
9
, 790 (2013)
.
[9] A. M. Lindenberg, I. Kang, S. L. Johnson, T. Missalla, P. A.
Heimann, Z. Chang, J. Larsson, P. H. Bucksbaum, H. C.
Kapteyn, H. A. Padmore, R. W. Lee, J. S. Wark, and R. W.
023072-8
TOWARD PRECISE SIMULATIONS OF THE COUPLED ...
PHYSICAL REVIEW RESEARCH
3
, 023072 (2021)
Falcone, Time-Resolved X-Ray Diffraction from Coherent
Phonons During a Laser-Induced Phase Transition,
Phys. Rev.
Lett.
84
, 111 (2000)
.
[10] P. Zalden, F. Quirin, M. Schumacher, J. Siegel, S. Wei, A.
Koc, M. Nicoul, M. Trigo, P. Andreasson, H. Enquist
et al.
,
Femtosecond X-ray diffraction reveals a liquid–liquid phase
transition in phase-change materials,
Science
364
, 1062 (2019)
.
[11] S. De Jong, R. Kukreja, C. Trabant, N. Pontius, C. Chang, T.
Kachel, M. Beye, F. Sorgenfrei, C. Back, B. Bräuer
et al.
, Speed
limit of the insulator–metal transition in magnetite,
Nat. Mater.
12
, 882 (2013)
.
[12] R. Mankowsky, A. von Hoegen, M. Först, and A. Cavalleri,
Ultrafast Reversal of the Ferroelectric Polarization,
Phys. Rev.
Lett.
118
, 197601 (2017)
.
[13] M. Mitrano, A. Cantaluppi, D. Nicoletti, S. Kaiser, A. Perucchi,
S. Lupi, P. Di Pietro, D. Pontiroli, M. Riccò, S. R. Clark
et al.
, Possible light-induced superconductivity in K
3
C
60
at high
temperature,
Nature (London)
530
, 461 (2016)
.
[14] R. B. Wilson, J. P. Feser, G. T. Hohensee, and D. G. Cahill, Two-
channel model for nonequilibrium thermal transport in pump-
probe experiments,
Phys.Rev.B
88
, 144305 (2013)
.
[15] H. Haug and A.-P. Jauho,
Quantum Kinetics in Transport and
Optics of Semiconductors
, Vol. 2 (Springer, Boston, 2008).
[16] T. Ogitsu, A. Fernandez-Pañella, S. Hamel, A. A. Correa, D.
Prendergast, C. D. Pemmaraju, and Y. Ping,
Ab initio
modeling
of nonequilibrium electron-ion dynamics of iron in the warm
dense matter regime,
Phys. Rev. B
97
, 214203 (2018)
.
[17] B. I. Cho, K. Engelhorn, A. A. Correa, T. Ogitsu, C. P.
Weber, H. J. Lee, J. Feng, P. A. Ni, Y. Ping, A. J. Nelson, D.
Prendergast, R. W. Lee, R. W. Falcone, and P. A. Heimann,
Electronic Structure of Warm Dense Copper Studied by Ul-
trafast X-Ray Absorption Spectroscopy,
Phys. Rev. Lett.
106
,
167601 (2011)
.
[18] D. Duffy and A. Rutherford, Including the effects of electronic
stopping and electron-ion interactions in radiation damage sim-
ulations,
J. Phys.: Condens. Matter
19
, 016207 (2006)
.
[19] A. Tamm, M. Caro, A. Caro, G. Samolyuk, M. Klintenberg, and
A. A. Correa, Langevin Dynamics with Spatial Correlations as
a Model for Electron-Phonon Coupling,
Phys. Rev. Lett.
120
,
185501 (2018)
.
[20] L. Waldecker, R. Bertoni, R. Ernstorfer, and J. Vorberger,
Electron-Phonon Coupling and Energy Flow in a Simple Metal
Beyond the Two-Temperature Approximation,
Phys.Rev.X
6
,
021003 (2016)
.
[21] E. D. Murray, S. Fahy, D. Prendergast, T. Ogitsu, D. M. Fritz,
and D. A. Reis, Phonon dispersion relations and softening in
photoexcited bismuth from first principles,
Phys.Rev.B
75
,
184301 (2007)
.
[22] M. Bernardi, D. Vigil-Fowler, J. Lischner, J. B. Neaton, and
S. G. Louie,
Ab Initio
Study of Hot Carriers in the First Pi-
cosecond After Sunlight Absorption in Silicon,
Phys. Rev. Lett.
112
, 257402 (2014)
.
[23] M. Bernardi, J. Mustafa, J. B. Neaton, and S. G. Louie, Theory
and computation of hot carriers generated by surface plasmon
polaritons in noble metals,
Nat. Commun.
6
, 7044 (2015)
.
[24] V. A. Jhalani, J.-J. Zhou, and M. Bernardi, Ultrafast hot carrier
dynamics in GaN and its impact on the efficiency droop,
Nano
Lett.
17
, 5012 (2017)
.
[25] M. Bernardi, First-principles dynamics of electrons and
phonons,
Eur.Phys.J.B
89
, 239 (2016)
.
[26] S. Sadasivam, M. K. Y. Chan, and P. Darancet, Theory of
Thermal Relaxation of Electrons in Semiconductors,
Phys. Rev.
Lett.
119
, 136602 (2017)
.
[27] K. Yabana and G. F. Bertsch, Time-dependent local-density
approximation in real time,
Phys.Rev.B
54
, 4484 (1996)
.
[28] A. Castro, M. A. L. Marques, and A. Rubio, Propagators for
the time-dependent Kohn-Sham equations,
J. Chem. Phys.
121
,
3425 (2004)
.
[29] M. A. L. Marques and E. K. Gross, Time-dependent density
functional theory,
Annu. Rev. Phys. Chem.
55
, 427 (2004)
.
[30] M. K. Nazeeruddin, F. De Angelis, S. Fantacci, A. Selloni, G.
Viscardi, P. Liska, S. Ito, B. Takeru, and M. Grätzel, Combined
experimental and DFT-TDDFT computational study of photo-
electrochemical cell ruthenium sensitizers,
J. Am. Chem. Soc.
127
, 16835 (2005)
.
[31] C. A. Rozzi, S. M. Falke, N. Spallanzani, A. Rubio, E. Molinari,
D. Brida, M. Maiuri, G. Cerullo, H. Schramm, J. Christoffers
et al.
, Quantum coherence controls the charge separation in a
prototypical artificial light-harvesting system,
Nat. Commun.
4
,
1602 (2013)
.
[32] F. Sottile, F. Bruneval, A. G. Marinopoulos, L. K. Dash, S. Botti,
V. Olevano, N. Vast, A. Rubio, and L. Reining, TDDFT from
molecules to solids: The role of long-range interactions,
Int. J.
Quantum Chem.
102
, 684 (2005)
.
[33] M. A. L. Marques, N. T. Maitra, F. M. S. Nogueira, E. K. U.
Gross, and A. Rubio,
Fundamentals of Time-Dependent Density
Functional Theory
, Vol. 837 (Springer, Boston, 2012).
[34] A. Dreuw, J. L. Weisman, and M. Head-Gordon, Long-range
charge-transfer excited states in time-dependent density func-
tional theory require non-local exchange,
J. Chem. Phys.
119
,
2943 (2003)
.
[35] P. Elliott, J. I. Fuks, A. Rubio, and N. T. Maitra, Universal
Dynamical Steps in the Exact Time-Dependent Exchange-
Correlation Potential,
Phys. Rev. Lett.
109
, 266404 (2012)
.
[36] J. Yang, X. Zhu, J. P. F. Nunes, K. Y. Jimmy, R. M. Parrish, T. J.
Wolf, M. Centurion, M. Gühr, R. Li, Y. Liu
et al.
, Simultaneous
observation of nuclear and electronic dynamics by ultrafast
electron diffraction,
Science
368
, 885 (2020)
.
[37] A. V. Akimov and O. V. Prezhdo, The PYXAID program for
non-adiabatic molecular dynamics in condensed matter sys-
tems,
J. Chem. Theory Comput.
9
, 4959 (2013)
.
[38] G. D. Mahan,
Many-Particle Physics
(Springer, Boston, 2000).
[39] F. Caruso, Nonequilibrium lattice dynamics in monolayer
MoS
2
,
J. Phys. Chem. Lett.
12
, 1734 (2021)
.
[40] D. Sangalli and A. Marini, Ultra-fast carriers relaxation in bulk
silicon following photo-excitation with a short and polarized
laser pulse,
Europhys. Lett.
110
, 47004 (2015)
.
[41] G. D. Mahan,
Condensed Matter in a Nutshell
(Princeton Uni-
versity Press, Princeton, NJ, 2010).
[42] 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)
.
[43] 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)
.
[44] T. Winzer and E. Mali
́
c, Impact of auger processes on carrier
dynamics in graphene,
Phys.Rev.B
85
, 241404(R) (2012)
.
[45] H. Wang, J. H. Strait, P. A. George, S. Shivaraman, V. B.
Shields, M. Chandrashekhar, J. Hwang, F. Rana, M. G. Spencer,
023072-9