Dynamic mode decomposition of nonequilibrium electron-phonon dynamics: accelerating the first-principles real-time Boltzmann equation

npj |

computationa

l

materia

l

Article

Published in partnership with the Shanghai Institute of Ceramics of the Chinese Academy of Sciences

https://doi.org/10.1038/s41524-024-01308-4

Dynamic mode decomposition of

nonequilibrium electron-phonon

dynamics: accelerating the

fi

rst-

principles real-time Boltzmann equation

Check for updates

Ivan Maliyov

, Jia Yin

, Jia Yao

1,3

,ChaoYang

& Marco Bernardi

1,3

Nonequilibrium dynamics governed by electron

–

phonon (

-ph) interactions plays a key role in

electronic devices and spectroscopies and is central to understanding electronic excitations in

materials. The real-time Boltzmann transport equation (rt-BTE) with collision processes computed

from

fi

rst principles can describe the coupled dynamics of electrons and atomic vibrations (phonons).

Yet, a bottleneck of these simulations is the calculation of

–

ph scattering integrals on dense

momentum grids at each time step. Here we show a data-driven approach based on dynamic mode

decomposition (DMD) that can accelerate the time propagation of the rt-BTE and identify dominant

electronic processes. We apply this approach to two case studies, high-

fi

eld charge transport and

ultrafast excited electron relaxation. In both cases, simulating only a short time window of ~10% of the

dynamics suf

fi

ces to predict the dynamics from initial excitation to steady state using DMD

extrapolation. Analysis of the momentum-space modes extracted from DMD sheds light on the

microscopic mechanisms governing electron relaxation to a steady state or equilibrium. The

combination of accuracy and ef

fi

ciency makes our DMD-based method a valuable tool for

investigating ultrafast dynamics in a wide range of materials.

First-principles calculations are

widely employed for modeling and

designing materials, with applications ranging from energy

to (opto)

electronic devices

–

to materials discovery

–

. Starting from the crystal

structure and atomic positions as the main inputs, these methods can

predict material properties, includin

g mechanical, electrical, magnetic, and

optical. While computing ground-state

and linear-response properties with

density functional theory (D

FT) is a decades-long effort

–

,recentworkhas

focused on modeling ultrafast dynami

cs in materials and simulating time-

domain spectroscopies from

fi

rst principles

–

. These methods focused on

nonequilibrium dynamics are a more recent research frontier with both

theoretical and computational challenges.

First-principles calculations in the time domain provide a microscopic

description of nonequilibrium dynamics in materials. These methods

propagate in time quantities characterizing the quantum dynamics, such as

the time-dependent electron wave function, density

, density matrix

,or

Green

’

sfunction

, and can also access the time-dependent atomic positions

and lattice vibrations

. Different schemes are successful in different

regimes. For example, coherent elect

ron dynamics on the attosecond time

scale can be modeled effectively using time-dependent DFT

–

,butthat

approach is not ideal for modeling phonon dynamics, which occurs on a

picosecond time scale

The real-time Boltzmann transport e

quation (rt-BTE) has emerged as

an effective tool for exploring the coupled electron and phonon dynamics

from femtosecond to nanosecond timescales

.Inthert-BTE,thetime-

dependent electron populations are o

btained by solving a set of integro-

differential equations accounting for the

–

ph scattering processes on dense

momentum grids. Following an initial excitation, the rt-BTE is propagated

in time to reach thermal equilibrium or steady state in an external

fi

eld. This

scheme employs a femtosecond time step to capture the

–

ph scattering

processes. However, evaluating the s

cattering integral at each time step

makes the rt-BTE approach computatio

nally demanding, even for materials

with a handful of atoms in the unit cell.

Department of Applied Physics and Materials Science, California Institute of Technology, Pasadena, CA 91125, USA.

Applied Mathematics and Computational

Research Division, Lawrence Berkeley National Laboratory, Berkeley, CA 94720, USA.

Department of Physics, California Institute of Technology, Pasadena, CA

91125, USA.

e-mail:

bmarco@caltech.edu

npj Computational Materials

| (2024) 10:123

1234567890():,;

Data-driven techniques are increa

singly employed in materials mod-

eling, both for accelerating computational work

fl

ows and to gain physical

insight using learning algorithms

. In particular, dynamic mode

decomposition (DMD), which was developed in the last decade to study

fl

uid dynamics, is a valuable tool to linearize dynamical problems and

reduce their dimensionality

.InDMD,explicitsimulationofashort

initial time window allows one to learn the dominant modes governing the

dynamics and extrapolate the simulation to future times at low computa-

tional cost. Recent work has employed DMD to study electron dynamics

described by model Hamiltonians with purely electronic interactions

Yet, to date, DMD has not been applied to more computationally intensive

fi

rst-principles studies.

In this work, we combine DMD with

fi

rst-principles calculations of

nonequilibrium electron dynamics, using the framework of the rt-BTE in

the presence of

–

ph collisions and external

fi

elds. We show that DMD

provides an order-of-magnitude computational speed-up while retaining

the full accuracy of the

fi

rst-principles rt-BTE. In addition, DMD reveals key

momentum-space temporal patterns and achieves a signi

fi

cant dimen-

sionality reduction of the nonequilibr

ium physics. Our results include both

high-

fi

eld transport and transient excited-state dynamics and are accom-

panied by a careful characterization of convergence with respect to the size

of the sampling window during which DMD learns the dominant modes.

Taken together, this work provides the b

lueprint for combining data-driven

methods with

fi

rst-principles calculations to study nonequilibrium

dynamics in real materials.

Results

First-principles rt-BTE

We describe the electron distribution u

sing the time-dependent populations

(

), which quantify the occupa

tion of each electronic state

∣

,where

the electron crystal momentum, and

is the band index (from now on, we

omit the band index to simplify the notation). Starting from an initial

distribution at time zero,

(

= 0), in the rt-BTE, the populations evolve

according to

∂

_

∇

Þþ

;

where

is the collision integral accounting for

-ph scattering pro-

cesses in momentum space

and

includesany external

fi

elds applied to the

system.

The rt-BTE simulations use de

nse momentum grids to accurately

describe scattering between electroni

c states via absorption and emission of

phonons. The required grid sizes ar

e typically >100 × 100 × 100 for both

electron and phonon momenta. We time-step Eq. (

) using explicit solvers

(Euler or fourth-order Runge

–

Kutta) or more advanced Strang splitting

techniques

. The collision integral inclu

des a summation over the phonon

momentum grid and is evaluated at least once per time step using a parallel

algorithm implemented in the

PERTURBO

code

(see the

“

Methods

”

section

for details). Although here we focus on t

he dynamics of electrons interacting

with phonons, the rt-BTE formali

sm has also been extended to study

nonequilibrium phonon

and exciton dynamics

DMD learning and prediction of the dynamics

We employ DMD in combination wi

th rt-BTE simulations. The DMD

approach linearizes the dynamics by relating the states of the system at times

and

+

via a time-independent matrix

.Focusingonthe

–

dynamics, this amounts to advancing t

he electronic populations at time

using

Þ¼

;

where the populations

form a vector with size

equal to the number of

points in the electronic m

omentum grid (typically,

≈

−

). To obtain

the matrix

, we time-step the rt-BTE in a sampling window consisting of

time steps (using the

PERTURBO

code

), and then we form two matrices

and

. The populations

(

)from

−

are stacked column-wise in the

matrix

,withcolumn

corresponding to time

and containing the

populations

(

)forall

-points and bands. The populations from

are similarly stacked column-wise in the second matrix

According to equation (

), these matrices are related by

,but

computing

naively from the pseudo-inverse of

has a prohibitive cost

duetothelargesize

of the

-point grid. To circumvent this problem, in

DMD, one

fi

rst performs a truncated singular value decomposition

(SVD)

of the

matrix:

;

where

×

is a matrix with diagonal entries equal to the singular

values

arranged in decreasing order, while

×

and

×

are matrices collecting the mut

ually orthogonal singular

vectors

.(Above,

†

indicates the Hermitian conjugate of

Because

contains time-dependent populations, this SVD pro-

cedure can single out the main patterns in the momentum-space

dynamics. Here, we keep only the

fi

rst

singular values (typically,

≈

10) to restrict the solution space to the leading

momentum-space

modes and then project the matrix

onto this reduced

-dimensional

space. This procedure provides the matrix

, with reduced size

which can be diagonalized straight

forwardly to obtain the dominant

DMD modes. Using this procedure, the populations at future times

are

predicted

—

that is, obtained without explicit solution of the rt-

BTE

—

using

≈

DMD

;

where

are the momentum-space DMD mod

es obtained from the matrix

,and

DMD

and

are their frequencies and amplitudes (see the

“

Methods

”

section for detailed derivations).

We summarize the main steps of

this DMD procedure, which are

illustrated in Fig.

1. Simulate the rt-BTE dynamics for the

fi

rst

steps and construct the

matrices

and

;

2. Perform SVD on

fi

nd the matrix

in the reduced

-dimensional space;

Fig. 1 | Work

fl

ow of DMD plus rt-BTE calcula-

tions.

The

fi

rst

steps of the dynamics, which make

up the sampling window for DMD learning, are

simulated by solving the rt-BTE. The resulting

populations

(

) are stacked in the

and

matrices with a relative shift of one time step. The

dynamics at later times

is predicted with DMD

using the

leading modes obtained by SVD of the

matrix

and diagonalization of the

matrix

https://doi.org/10.1038/s41524-024-01308-4

Article

npj Computational Materials

| (2024) 10:123

3. Diagonalize

fi

nd the DMD modes

and their frequencies

DMD

with

...

;

4. Obtain the mode amplitudes

from the initial condition

(

);

5. Predict the dynamics for

using Eq. (

A key parameter is the duration of the sampling window (

)required

for accurate DMD extrapolation of the dynamics beyond

.Asthecom-

putational cost of DMD is negligible, the size of the sampling window,

during which the rt-BTE is solved by explicit time-stepping, determines the

computational cost of the entire work

fl

ow.

High-

fi

eld electron dynamics

We employ our DMD-based approach to simulate time-domain electron

dynamics in an applied electric

fi

eld in the presence of

-ph collisions. We

recently demonstrated similar calculat

ions using the rt-BTE without the aid

of data-driven techniques

.Here,weusethiscasestudytoexplorethe

accuracy and ef

fi

ciency of our rt-BTE plus DMD approach, as well as

fi

optimal values for the sampling window and analyze the momentum-space

DMD modes. Our calculations focus on electrons in GaAs, where the

conduction band has three se

ts of low-energy valleys, at

and near

and

in order of increasing energy

(see the inset in Fig.

a). Upon applying an

electric

fi

eld, the electrons are accelerated to higher band energies while they

also transfer part of that excess energy to the lattice via

–

ph collisions. These

competing mechanisms lead to a steady

-state electronic distribution which

is typically reached on a picosecond to nanosecond time scale.

Our simulations begin with electrons in thermal equilibrium with the

lattice at 300 K. We apply a constant electric

fi

eld

and time-step the

electron populations until they reach the steady-state distribution,

,from

whichwecomputethemeandriftvelocity,

(

), a quantity routinely

measured in experiments

–

. Repeating this procedure for multiple

fi

eld

values allows us to construct the ful

l drift velocity versus electric

fi

eld curve

in a material, starting from linear response at low

fi

eld to velocity saturation

at high

fi

eld

Figure

a shows the time-dependent populations in four regions of the

Brillouin zone following the application of a high

fi

eld (5 kV cm

−

). Elec-

trons initially occupying the

-valley scatter to the higher-energy

-and

-valleys. As a result, the electron populations in the

-valley decrease, with

a corresponding increase in

-and

-valley populations. In regions of

momentum space between the

-and

-valleys, the populations peak at

intermediate times and then relax to lower values.

This dynamics is nontrivial because the populations evolve differently

in different momentum-space regions,

making accurate predictions chal-

lenging. Our DMD approach can learn the dominant modes governing these

intricate dynamics and extrapolate t

he time-dependent populations well

beyond the sampling window. Remarkably, we

fi

nd that a short sampling

window

−

400 fs to 2 ps out of a total simulation time of 12.5 ps

—

is suf

fi

cient to extrapolate the dynamics all th

e way to steady state, with rt-BTE and

DMD trajectories in nearly exact agre

ement outside the sampling window

(Fig.

a). This accuracy extends to the entire set of ~10

-points considered

in our simulations, providing carrier n

umber conservation within 1% error.

The ability to learn key temporal momentum-space patterns is a

consequence of the relatively rapid decay of the singular values of the

matrix used for learning the dynam

ics in the sampling window (Fig.

b).

This decay becomes slower as the sampling window increases, but it remains

signi

fi

cant even for the longest sampling window of 2 ps used here (note the

log scale in the plot). In turn, the sin

gular value decay enables a striking

dimensionality reduction, with DMD employing only

≈

10 modes to solve

the dynamicsas opposed to 10

populations

and billions of

–

ph scattering

termsinthert-BTE.

Fig. 2 | DMD simulations of electrons in GaAs in an applied electric

fi

eld. a

Time-

dependent electron populations

(

) for four electron momenta. The gray region

indicates the shortest sampling window (0.4 ps), and the red vertical line is the

longest sampling window we tested (2 ps). Solid black lines show the rt-BTE results

and orange dashed lines, the DMD predictions obtained using the longest sampling

window. The inset is a schematic of the low-energy band structure of GaAs showing

the

- and higher energy

- and

-valleys.

Singular values of the

matrix, with the

ten largest singular values used in our DMD calculations separated by a vertical line.

DMD frequencies plotted in the complex plane and shown as circles with radii

proportional to the DMD mode amplitudes

. In panels

–

, the colors indicate the

duration of the DMD sampling window according to the legend given in (

DMD

momentum-space modes

, multiplied by the corresponding amplitudes

, given

as a function of energy. The initial state

(

) is shown with a dashed line.

Convergence of the steady-state drift velocity with respect to the duration of the

DMD sampling window. The rt-BTE value is shown for reference as a dashed line.

https://doi.org/10.1038/s41524-024-01308-4

Article

npj Computational Materials

| (2024) 10:123

The choice of an ideal sampling window can rely on the appearance of

speci

fi

c DMD modes at a steady state. Figure

cshowstheDMDmode

frequencies

DMD

in the complex plane, where the imaginary part of

DMD

corresponds to the decay rate of a given mode, and the real part gives its

oscillation frequency. The populations

(

) are real-valued and are written

as a summation of complex exponentials in Eq. (

). Therefore, physically

meaningful results are possible only when Re

DMD

Þ¼

0(modes1,2)or

when

DMD

appear as conjugate pairs (modes 3

−

10). Describing the steady

state is particularly important in our simulations. In DMD, all modes with a

non-zero imaginary frequency vanish in the long time limit, with only one

mode surviving at steady state (mode 1 in Fig.

c). As the sampling window

increases, the imaginary frequency of this mode goes to zero, providing the

correct steady-state behavior. This analysis allows us to

fi

nd the minimal

sampling window required for accurate

steady-state results by monitoring

the zero-frequency mode.

The DMD eigenvector of the zero

-frequency mode (mode 1 in Fig.

determines the steady-state electron distribution

!1Þ

,while

the other modes control the transie

nt dynamics. For example, mode 2

governs electron scattering from the

-tothe

-and

-valleys, and higher

modes appearing as conjugate pairs exhibit oscillating trends in energy

(modes 3

−

10 in Fig.

d). Converging the zero-frequency mode allows us to

compute the steady-state drift velocity more ef

fi

ciently. Figure

eshowsthat

a sampling window of 1.7 ps (170 snapshots) provides a drift velocity nearly

identical to the full rt-BTE calculation, which requires much longer simu-

lation times of up to 12.5 ps (1250 snapshots). On this basis, we conclude

that DMD needs only ~10% of the dynami

cs data for accurate steady-state

predictions.

Velocity-

fi

eld curves

We also employ DMD to accelerate calculations of entire velocity-

fi

eld

curves. This requires the drift velocity for a set of electric

fi

eld values, and

thus we adopt a modi

fi

ed work

fl

ow. Following our recent work

,wegra-

dually increase the electric

fi

eld (black curve in Fig.

a) and use the steady-

statepopulations for a given

fi

eld,

, as the initial condition for the next

fi

eld

value,

+

,wherethe

fi

eld increment

is typically 100

−

200 V cm

−

As the applied

fi

eld increases, the DMD frequencies and momentum-space

modes change substantially. Therefore, for each new

fi

eld value, we repeat

DMD learning in the initial stage of t

he simulation (see the DMD sampling

regions shown as red rectangles in Fig.

a). We then predict the steady-state

populations using mode 1 from DMD,

and the drift velocity for

that

fi

eld value, and use

as the initial condition for the next

fi

eld value.

The velocity-

fi

eld curves obtained with this approach are shown in Fig.

b for GaAs and graphene and compared with rt-BTE results obtained

without DMD. Using DMD lowers signi

fi

cantly the computational cost to

obtain the full velocity-

fi

eld curves by a factor of 10.5 for GaAs and ~ 16.5 for

graphene, while fully preserving the a

ccuracy. Because the drift velocity is

computed as a weighted sum of

, the nearly exact agreement between the

DMD and full rt-BTE results demonstrates the accuracy of the DMD

populations in momentum space. The DMD ef

fi

ciency is a consequence of

its ability to capture the dominant modes in the population dynamics using

only a small number of snapshots, with a

similar accuracy regardless of the

electric

fi

eld value.Our strategy of gradually increasing the electric

fi

eld leads

to easier-to-extrapolate dynamics compared to the abrupt application of a

strong

fi

eld.

Excited electron relaxation

Next, we consider different nonequil

ibrium dynamics where the material is

initially prepared in an excited electronic state. This setting can be used, for

example, to model the effect of an optical excitation with a laser pulse

Different from the high-

fi

eld dynamics, in this case, the long-time limit is

known, and we are primarily interested in the

transient

dynamics. Following

the initial excitation, in the presence of

–

ph interactions and without any

external

fi

elds, the electrons relax to a thermal equilibrium Fermi

–

Dirac

distribution

, typically on a sub-picosecond time scale. This ultrafast

dynamics can be modeled by time-ste

pping the rt-BTE until reaching the

equilibrium Fermi-Dirac distributi

on. Using this approach, our previous

work has shown that electrons

relax to the band edge signi

fi

cantly slower

than holes in GaN semiconductors, with

implications for optoelectronic

devices

Following that work, we model an ex

cited state in GaN by placing the

electrons ~1 eV above the conduction band edge and then obtain the time-

dependent electron populations by solving the rt-BTE (see Fig.

a, b). We

employ DMD to predict this transient dynamics and

fi

nd large errors when

using a short time window of up to ~50 fs (solid red line in Fig.

c). The

correct steady state and transient dynamics are obtained by increasing the

sampling window to 200 fs (dashed orange line in Fig.

c). Our analysis of

the DMD frequencies shows that the zero-frequency mode describing

thermal equilibrium in the long-ti

me limit appears when the sampling

window reaches 100 fs (see the arrow in Fig.

d) and fully converges for

a ~200 fs sampling window. The need

for such a long sampling window

relative to the total duration of the dynamics (400 fs) makes DMD

ineffective.

To address this issue and more ef

fi

ciently study transient dynamics

with DMD, we formulate a different learning procedure that incorporates

knowledge of the equilibrium state. We focus on the difference between the

transient and equilibrium populations,

Þ¼

,asopposedto

just

(

) as we did in the high-

fi

eld example. After predicting

(t) with

DMD, we obtain the time-dependent populations

(

)byaddingbackthe

term. As

vanishes in the long-time limit (Fig.

b), the zero-frequency

DMD mode is missing when computing

(Fig.

e); all other DMD fre-

quencies associated with

are similar to those for

(

)(Fig.

d, e). We

fi

that the DMD method based on

is far more effective and requires a

signi

fi

cantly shorter sampling windo

w for accurate DMD predictions

−

using a 50 fs sampling window, we achieve results similar to DMD for

(

)

with a four times longer (200 fs) window (Fig.

c).

With this improved DMD approach, using a sampling window of

only ~12% of the total simulation time

allows us to accurately predict the

average electron relaxation rate in GaN, with a DMD computed value of

5.23 eV fs

−

in close agreement (within 0.8%) with the rt-BTE result. This

result demonstrates that our DMD approach can predict excited electron

relaxation with a high accuracy.

Discussion

The DMD approach introduced here is very ef

fi

cient: the rt-BTE is solved

explicitly on a high-performance com

puter only for a small number of initial

time steps, after which the entire dyn

amics can be computed straightfor-

wardly with DMD, using only a laptop. T

he most demanding step is carrying

out truncated SVD on the

matrix, but for comparison, this step requires

lower computational resources than even just a single rt-BTE time step.

This remarkable speed-up is achieved by reducing the dimensionality

of the rt-BTE dynamics and is

linked to the shape of the

matrix. The rt-

BTE employs a large number of

-points (about 10

−

), which equals the

number of rows of the matrix

, and a signi

fi

cantly smaller number of

snapshots in the DMD sampling window

, typically ~100 time steps, which

sets the number of columns in

. Following truncated SVD, the size of the

problem is reduced to (at most) the number of snapshots and is typically of

order 50

−

100, and thus smaller by orders of magnitude compared to the

original rt-BTE. (Note that one could use the entire set of singular values, but

here we prefer using only ~10 singular values to prevent numerical

instabilities

This ef

fi

ciency allows us to evaluate the accuracy of DMD on the

fl

halting explicit time-stepping of th

e rt-BTE when the DMD steady state or

transient dynamics are fully converged. In addition, our approach addresses

the key challenge of storing the rt-

BTE populations. This is a critical

improvement because in conventional rt-BTE simulations one needs to

store the populations

(

) on dense momentum grids for thousands of time

steps, resulting in terabytes of data. I

n contrast, after carrying out SVD in the

sampling window, DMD stores only a handful of complex frequencies and

momentum-space modes, using which the dynamics can be reconstructed

for the entire simulation.

https://doi.org/10.1038/s41524-024-01308-4

Article

npj Computational Materials

| (2024) 10:123