#  Ab-initio Simulations of Coherent Phonon-Induced Pumping of Carriers in
Zirconium Pentatelluride

Tao Jiang  Ames National Laboratory, U.S. Department of Energy, Ames, Iowa
50011, USA  Peter P. Orth  Ames National Laboratory, U.S. Department of
Energy, Ames, Iowa 50011, USA  Department of Physics and Astronomy, Iowa State
University, Ames, Iowa 50011, USA  Department of Physics, Saarland University,
66123 Saarbrücken, Germany  Liang Luo  Ames National Laboratory, U.S.
Department of Energy, Ames, Iowa 50011, USA  Lin-Lin Wang  Ames National
Laboratory, U.S. Department of Energy, Ames, Iowa 50011, USA  Feng Zhang  Ames
National Laboratory, U.S. Department of Energy, Ames, Iowa 50011, USA
Department of Physics and Astronomy, Iowa State University, Ames, Iowa 50011,
USA  Cai-Zhuang Wang  Ames National Laboratory, U.S. Department of Energy,
Ames, Iowa 50011, USA  Department of Physics and Astronomy, Iowa State
University, Ames, Iowa 50011, USA  Jin Zhao  ICQD/Hefei National Laboratory
for Physical Sciences at Microscale, and Key Laboratory of Strongly-Coupled
Quantum Matter Physics, Chinese Academy of Sciences, and Department of
Physics, University of Science and Technology of China, Hefei, China
Department of Physics and Astronomy, University of Pittsburgh, Pittsburgh,
Pennsylvania, USA  Synergetic Innovation Center of Quantum Information &
Quantum Physics, University of Science and Technology of China, Hefei, China
Kai-Ming Ho  Ames National Laboratory, U.S. Department of Energy, Ames, Iowa
50011, USA  Department of Physics and Astronomy, Iowa State University, Ames,
Iowa 50011, USA  Jigang Wang  Ames National Laboratory, U.S. Department of
Energy, Ames, Iowa 50011, USA  Department of Physics and Astronomy, Iowa State
University, Ames, Iowa 50011, USA  Yong-Xin Yao<EMAIL_ADDRESS>Ames
National Laboratory, U.S. Department of Energy, Ames, Iowa 50011, USA
Department of Physics and Astronomy, Iowa State University, Ames, Iowa 50011,
USA

###### Abstract

Abstract

Laser-driven coherent phonons can act as modulated strain fields and modify
the adiabatic ground state topology of quantum materials. Here we use time-
dependent first-principles and effective model calculations to simulate the
effect of the coherent phonon induced by strong terahertz electric field on
electronic carriers in the topological insulator ZrTe5. We show that a
coherent $A_{\text{1g}}$ Raman mode modulation can effectively pump carriers
across the band gap, even though the phonon energy is about an order of
magnitude smaller than the equilibrium band gap. We reveal the microscopic
mechanism of this effect which occurs via Landau-Zener-Stückelberg tunneling
of Bloch electrons in a narrow region in the Brillouin zone center where the
transient energy gap closes when the system switches from strong to weak
topological insulator. The quantum dynamics simulation results are in
excellent agreement with recent pump-probe experiments in ZrTe5 at low
temperature.

##  I Introduction

Coherent phonons that are excited by laser pulses in the THz or mid-infrared
frequency range [1, 2, 3] can provide nonthermal pathways for the dynamical
control of quantum phases of condensed matter [4, 5, 6, 7]. Recent
experimental demonstrations include ultrafast phononic manipulation of
magnetic orders [8, 9, 10], of insulator to metal phase transitions [11, 12,
13], and a transient enhancement of martensitic phase [14] and superconducting
correlations [15, 16, 17, 18]. In topological quantum materials, coherent
phonon excitations were shown to induce dynamical switching between different
topological phases by modifying the crystal symmetry and by tuning strain
fields [19, 20, 21, 22, 23, 24, 25, 26, 27]. Particularly, recent coherent
phonon pumping work provides compelling evidence of light-induced Dirac points
[25], Weyl nodes [26] and enhanced stability of topological systems [28, 29].

The theoretical understanding and first-principles simulations of light-
excited electron-ion quantum systems are challenging, but significant progress
has been achieved in recent years [30, 31, 32, 33, 34]. Specifically to
describe the ultrafast electronic and spin dynamics associated with phonon
excitations, microscopic theories have been developed for the light-induced
insulator-to-metal structural phase transition [35], the switching of magnetic
orders [36], and the enhancement of superconducting correlations through
symmetry-allowed electron-phonon coupling [37, 38, 39, 40, 41, 42, 43, 44, 45,
46]. Here, we theoretically investigate the switching between strong and weak
topological insulators (STI and WTI) induced by THz-driven coherent Raman
phonon excitations in the model Dirac system ZrTe5 [25]. Although a
qualitative picture has previously been established using static density
functional theory (DFT) calculations [25], the ultrafast dynamics of the
laser-driven system such as the observed continuous increase of electronic
carrier density after the THz pump and the underlying mechanism, calls for a
more in-depth quantum dynamics simulations. The progress in this direction is
highly appealing to the experimental community, who have been actively
pursuing THz-driven quantum dynamics in various quantum materials recently
[47, 48, 49, 50].

In this paper, we simulate the coherent phonon-induced carrier dynamics in
ZrTe5 in the framework of time-dependent Schrödinger equation with DFT basis
functions. We complement the DFT-based dynamics simulations by an effective
model calculations that captures the essentials of the microscopic mechanism.
Our detailed numerical analysis shows that the switching between STI and WTI,
which necessarily involves the closing of the bulk gap, creates a small but
finite volume in momentum space, where effective two-level systems (TLSs)
undergo avoided level crossings. This results in Landau-Zener-Stückelberg
(LZS) tunnelling [51, 52, 53, 54, 55, 56] and leads to an increase of the
carrier concentration during several cycles of the coherent phonon modulation.
Our time-dependent Schrödinger equation simulations predict the dynamics of
the phonon-induced carrier concentration in quantitative agreement with
experiment.

##  II Results and discussion

###  II.1 Summary of previous pump-probe experimental results

To facilitate the presentation, we summarize the key results of the THz pump-
THz probe experiment on ZrTe5 at $4.2$ K, which is described in detail in Ref.
[25]. These results motivate the numerical simulations in this work. In the
experiment, an intense THz-pump pulse with an $E$-field trace plotted in Fig.
1(a) is incident normally on the ZrTe5 single crystal. The THz pump-induced
coherent phonon emission from the sample is observed after the pump pulse
between $2.5$ ps $\lesssim t\lesssim 5.8$ ps, as highlighted in Fig. 1(b). The
coherent phonon emission lasts for about five full cycles and its dominant
spectral peak at $f_{\text{ph}}=1.2$ THz after Fourier transformation matches
the $A_{1g}$ Raman mode at the Brillouin zone center. By performing THz pump
and THz probe measurements using the same pump pulse, the THz probe
differential transmission $\Delta E/E_{0}$, which is proportional to the
change of carrier density $\Delta n$, is obtained and plotted in Fig. 1(c).
The carrier density continuously increases after the pump pulse as long as the
coherent phonon emission is observed. It saturates for $t\gtrsim 5.8$ ps,
which coincides with the loss of phonon coherence. $\Delta E/E_{0}$ decays
back to zero after about $120$ ps [25].

In this paper, we focus on simulating this intriguing carrier excitation
dynamics for the time period $2.5$ ps $\lesssim t\lesssim 5.8$ ps, where the
coherent phonon excitation is present. The residual pump pulse is negligible
during this time period, and the coherent phonon excitation can be treated as
preexisting, i.e., without explicitly modelling the light-driven phonon
generation process. The strong correlation between the carrier generation and
the coherent phonon emission suggests a charge excitation mechanism assisted
by a coherent Raman vibration. Indeed, by adiabatically following the $A_{1g}$
phonon trajectory, static DFT calculations have revealed that the electronic
state of the system undergoes a topological transition between STI to critical
Dirac point (DP) to WTI [25] (see also Appendix IV). This suggests the
importance of the associated closing of the bulk band gap and potentially
further topological effects in the carrier pumping process [25]. This makes a
detailed quantum dynamics simulations of the physical process highly
desirable.

Figure 1:  THz pump-THz probe experimental results of ZrTe5 at $4.2$ K. (a)
Normalized THz pump $E$-field as a function of pump delay time $t$, with the
maximal value of $E(t)$ being $E_{\text{max}}=736$ kV/cm. (b) Phonon emission
as a function of $t$. The time-period after the pump pulse, where a coherent
phonon emission is observed, is highlighted in red. This region is the focus
of this work. (c) Normalized THz probe differential transmission $\Delta
E(t)/E_{0}$ as a function of $t$. $\Delta E(t)$ is the differential
transmission of the THz probe peak field $E_{0}$ measured with and without the
THz pump pulse. The right $y$-axis labels the corresponding change of the
carrier density $\Delta n$ in unit of $10^{16}$ cm-3.

###  II.2 Effective model description

Model setup. We first study a toy model that qualitatively captures the
dynamical carrier generation observed experimentally in ZrTe5. This model
includes the essential physics underlying this phenomenon which is the
topological phase transition (driven by coherent phonon oscillations) and a
resulting inter-band charge excitation. For simplicity, we consider a Kitaev
chain model [57], which in momentum space is represented by the following two-
orbital spinless fermion Bogoliubov-de Gennes (BdG) Hamiltonian:

$H(k)=(-2\nu\cos(k)-\mu)\tau_{z}+2\Delta\sin(k)\tau_{y}\,.$ (1)

Here $\tau_{i}$ are the Pauli matrices and the Hamiltonian parameters include
onsite energy $\mu$, nearest-neighbor hopping $\nu$ and a superconducting
pairing amplitude $\Delta$. The momentum lies in the range $k\in[-\pi,\pi)$.
The model obeys particle hole symmetry $\tau_{x}H^{*}(-k)\tau_{x}=-H(k)$.

With the chemical potential fixed at zero, the model exhibits a topological
phase transition from a gapped superconductor that is trivial (NI) to one that
is topological (TI) by tuning $\mu$. The corresponding BdG band structures
together with the orbital ($\tau_{z}$) projections of the wavefunction are
shown in Fig. 2(a-e). This evolution of the band energies is qualitatively
similar to the phonon-induced topological phase transition in ZrTe5 obtained
from DFT calculations [25, 58]. The band structure in panels (a) to (d) are
obtained for $\mu/\nu$ equal to $-2.02$, $-2.00$, $-1.97$ to $-1.92$ (a-d).
The evolution of the band gap as a function of $\mu/\nu$ is shown in panel
(e). The bands in panels (a-d) are plotted together with weight of the
projection of the Bloch wave function onto the first basis orbital
$(\tau_{z}=1)$, as given by the size of the red circles. Clearly, a band
inversion occurs when the system transforms from the trivial phase
($\mu/\nu<-2$) to the topological phase ($\mu/\nu>-2$). At $\mu/\nu=-2$, the
band gap closes and a nodal point forms at $\Gamma$ point. Consistently, the
topological index, which can be defined as
$Q\equiv\text{sign}(\mu^{2}-4\nu^{2})$ [57], is $1$ in the trivial and $-1$ in
the topological phase, as labelled in panels (a,c,d). For the numerical
simulations, we set $\nu=1$ eV, and $\Delta/\nu=0.1$, which results in a mode
speed $\partial\varepsilon_{k}/\partial k=\pm 2\Delta$ at the Dirac point
($\mu/\nu=-2$) that is in qualitative agreement with the Fermi velocity in
ZrTe5.

To model the effect of the coherent phonon excitation in ZrTe5, we consider a
periodic modulation of the onsite energy which resembles the treatment of
electron-phonon coupling in the Holstein model [59]:

$\mu(t)=\mu_{0}+\mu_{1}\sin[\omega(t-t_{0})]$ (2)

resulting in the time-dependent Hamiltonian

$H(k,t)=[-2\nu\cos(k)-\mu(t)]\tau_{z}+2\Delta\sin(k)\tau_{y}\,.$ (3)

We choose $\omega/2\pi=1.2$ THz ($T=833$ fs, $\hbar\omega=4.96$ meV) to match
the experimental value of the $A_{\text{1g}}$ Raman mode frequency. We set
$\mu_{0}/\nu=-1.97$ and $\mu_{1}/\nu=0.05$, such that the gap variation at the
$\Gamma$ point is approximately the same as in ZrTe5 [25], as shown in Fig.
2(e). Because the zone-center phonon carries zero momentum ($q=0$) it does not
mix different Bloch momenta and the Hamiltonian $H(k,t)$ thus remains block
diagonal in momentum space. The quantum dynamics simulation can therefore be
performed by solving the time-dependent Schrödinger equation separately at
each $k$-point:

$i\hbar\frac{\partial}{\partial t}\ket{\psi(k,t)}=H(k,t)\ket{\psi(k,t)}\,.$
(4)

We adopt a discrete-time propagator based on a Trotter decomposition of the
state evolution [60, 61]

$\ket{\psi(k,t+dt)}=e^{-i\frac{H(k,t)}{\hbar}dt}\ket{\psi(k,t)},$ (5)

where the time step $dt\ll T$ is chosen sufficiently small compared to the
variation of the onsite energy.

Figure 2:  Toy model results of phonon-induced topological phase transition
and carrier excitation dynamics. (a-d) Band structure of the BdG Kitaev chain
model with $\mu/\nu=-2.02$, $-2.00$, $-1.97$, and $-1.92$. The red color
encodes the projection weight of the band wavefunction on the first basis
orbital ($\tau_{z}=1$). The topological index $Q=\pm 1$ is also shown. (e)
Band gap at the zone center $\Gamma$ point, $E_{G}(\Gamma)$, as a function of
$\mu/\nu\in[-2.02,-1.92]$. This behavior qualitatively mirrors the behavior of
ZrTe5 system under the $A_{\text{1g}}$ Raman phonon modulation [25]. Red
circles indicate $\mu/\nu$ values in panels (a-d). The topological region is
highlighted in blue. The vertical dashed line indicates the equilibrium value
$\mu_{0}$. (f) Time evolution of the excited state population $n_{e}(k,t)$ for
a periodic modulation $\mu(t)/\nu\in[-2.02,-1.92]$ starting at
$\mu(t=0)=\mu_{0}$ for five full cycles. The black line is for $k=0.0033\pi$
and the orange line for $k=0.033\pi$. (g) Time evolution of the energy gap
$E_{G}(k,t)$ at $k=0.0033\pi$ (black) and $0.033\pi$ (orange). (h) Excited
state population $n_{e}(t)$. The black curve is for $t_{0}=0$ and the blue one
is averaged over $10$ runs with $t_{0}$ taken from $10$ uniformly spaced
points in the interval $[0,T)$. (i) The $k$-point dependence of the offset
$\varepsilon_{0}(k)$ (9). The bottom horizontal dotted line in (e, i)
indicates the energy of the $f_{\text{ph}}=1.2$ THz $A_{\text{1g}}$ phonon
$E_{\text{ph}}$. The upper three horizontal dotted lines in (i) correspond to
$E=mE_{\text{ph}}$ with $m=10$, $11$, and $12$. The vertical dashed lines in
panels (i,j) indicate the $k$-points where
$\varepsilon_{0}(k)=mE_{\text{ph}}$. The minimal energy gap
$E_{G}^{\text{min}}(k)\equiv\min_{t\in[0,T]}E_{G}(k,t)=\absolutevalue{\Delta(k)}$
over the simulation period is also plotted as orange line for reference. (j)
The $k$-point resolved number density in the excited band, which is defined as
$W_{n}(k)\equiv\int_{0}^{5T}n_{e}(k,t)dt/\int_{0}^{5T}n_{e}(t)dt$. (k) The
$k$-dependent driving speed indicators: $\delta(k)$ in black line, LZ
transition probability $P_{\text{LZ}}(k)\equiv e^{-2\pi\delta(k)}$ in orange
line, and $1-P_{\text{LZ}}(k)$ in sky blue line.

Dynamics simulation results. The simulation starts at time $t=t_{0}$ and
evolves until time $t=t_{0}+5T$ to agree with the experimental situation (see
Fig. 1). We monitor the time-dependent excited state population $n_{e}(t)$,
which corresponds to the excited carrier density in the experiment and is tied
to the differential emission $\Delta E/E$ that is measured experimentally. In
our simulations we obtain
$n_{e}(t)=\frac{1}{2\pi}\int_{-\pi}^{\pi}dkn_{e}(k,t)\approx\sum_{k}w_{k}n_{e}(k,t)$
as a weighted sum of contributions $n_{e}(k,t)$ at each $k$-point in the
Brillouin zone $[-\pi,\pi)$. The weight $w_{k}$ is obtained as $1/N_{k}$ with
$N_{k}$ the total number of $k$-points uniformly sampled in the Brillouin
zone. Here we define

$n_{e}(k,t)=\absolutevalue{\langle\psi(k,t)|\,\psi_{c}(k,t)\rangle}^{2},$ (6)

which is the size of the projection of the one-electron wavefunction
$\psi(k,t)$ on the adiabatic conduction band $\psi_{c}(k,t)$, an eigenstate of
$H(k,t)$ with eigenvalue $\epsilon_{c}(k,t)$. We use a step size
$dt=T/1000=0.833$ fs and a uniform $k$-mesh with $N_{k}=560$ points, and find
that $n_{e}(t)$ converges to a precision of $10^{-6}$.

In Fig. 2(f) we plot the excited state population $n_{e}(k,t)$ as a function
of time $t$. We consider five modulation cycles of the onsite energy $\mu(t)$
starting with $t_{0}=0$. The black line represents the data at a momentum
point close to the zone center $k=0.0033\pi$, and the orange line is for
momentum $k=0.033\pi$. The corresponding (instantaneous) energy gap
$E_{G}(k,t)$ at the respective $k$-points versus $t$ is shown in Fig. 2(g).
Sharp variations of $n_{e}(k,t)$ are observed whenever the energy gap is
minimal. In total we observe a substantial increase of $n_{e}(k,t)$ from zero
to a finite value at the two representative $k$-points during the dynamical
process. We note that the definition of $n_{e}(k,t)$ (6) may not be unique. In
Supplementary Note 1, we compare the numerical result with that based on an
alternative definition, and show that the definition (6) gives more physically
reasonable results.

The excited state population per unit cell, $n_{e}(t)$, is plotted as a black
line in Fig. 2(h). It shows a similar behavior as $n_{e}(k,t)$ at the two
individual $k$-points shown in Fig. 2(f). By the end of the simulation,
$n_{e}(t=t_{0}+5T)$ has increased from zero to about $0.01$. To account for
the experimentally unknown initial phase of the coherent phonon oscillation,
we also provide results that are averaged over $t_{0}$ that is uniformly
sampled within the interval $[0,T]$. The averaged results
$\bar{n}_{e}(t)=\frac{1}{10}\sum_{i=0}^{9}n_{e}(t)|_{t_{0}=\frac{i}{10}T}$ are
shown in orange in Fig. 2(h) and also increase from zero to about $0.01$
during the simulations.

Discussion of toy model results. The dynamics of the wavefunction at different
$k$-points is completely independent from each other [see Eq. (4)]. The above
calculation is thus composed of a collection of independent and periodically-
driven two-level systems (TLSs), which resembles the well-known problem of LZS
tunneling of a driven TLS in the presence of an avoided crossing [55, 56]. The
behavior of the TLS is largely determined by the potential ramp speed (i.e.
the oscillation frequency of the drive) and the minimal energy gap in the
avoided crossing. While the ramp speed is set by the phonon frequency ($1.2$
THz), the minimal energy gap
$E_{G}^{\text{min}}(k)\equiv\min_{t\in[0,T]}E_{G}(k,t)$ and the nature of the
diabatic level crossing is strongly $k$-point dependent. In Fig. 2 (i), we
show the minimal gap $E_{G}^{\text{min}}$ versus $k$ in the range
$k\in[-\pi/20,\pi/20]$. The minimal gap changes almost linearly from to zero
in the zone center ($k=0$) to about $60$ meV at $k=\pm\pi/20$. For reference,
we also plot the $1.2$ THz phonon mode energy ($E_{\text{ph}}=4.96$ meV) as a
dotted horizontal line. These findings suggest that the contribution to the
excited carrier density $n_{e}(t)$ arises from a small part in momentum space
around the zone center.

To obtain a more quantitative analysis, we define the following time-averaged
and $k$-resolved excited state density:

$W_{n}(k)\equiv\frac{\int_{0}^{5T}n_{e}(k,t)dt}{\int_{0}^{5T}n_{e}(t)dt}\,.$
(7)

We plot $W_{n}(k)$ in Fig. 2 (j) for a simulation with $t_{0}=0$, and find
that it peaks in a $\Gamma$-centered narrow $k$-region. Interestingly, peaks
occur in a wider range than naively expected by the condition that
$E_{G}^{\text{min}}(k)\leq E_{\text{ph}}$. This shows that higher order LZS
resonances are important, where excitations occur across a minimal band gap
that is a multiple of the driving frequency. Following the analysis of LZS
tunneling, we define diabatic states
$\ket{\tilde{\psi}(k)}\equiv\ket{\psi(k,t=0)}$. We find that the $k$-range,
where a diabatic level crossing occurs,
$\langle\tilde{\psi}_{0}(k)|\,H(k,t^{\prime})\,|\tilde{\psi}_{0}(k)\rangle=\langle\tilde{\psi}_{1}(k)|\,H(k,t^{\prime})\,|\tilde{\psi}_{1}(k)\rangle$,
at some time $0<t^{\prime}<5T$, matches well with the region of significant
$W_{n}(k)$ as highlighted in blue. We observe several distinct peaks of
$W_{n}(k)$ in this region: the location of the peak closest to $\Gamma$-point
($k=0.004\pi$) and that of the peak furthest away ($k=0.033\pi$) are close to
the $k$-points chosen for the presentation of $n_{e}(k,t)$ in (f, g). The
peaks can be understood by considering resonance condition of the multi-cycle
LZS problem as detailed in the following paragraph.

Although there is no closed form for the solution of the general LZS problem,
some analytical understanding can corroborate the numerical results [55, 56].
To facilitate the discussion, we cast the time-dependent Hamiltonian (3) into
the following form adopted in references [55, 56]:

$H(k,t)=-\frac{\Delta(k)}{2}\tau_{x}-\frac{\varepsilon(k,t)}{2}\tau_{z}\,,$
(8)

where we apply a $-\pi/2$ rotation around $z$-axis, and define
$\Delta(k)\equiv-4\Delta\sin(k)$,
$\varepsilon(k,t)\equiv\varepsilon_{0}(k)+A\sin[\omega(t-t_{0})]$, with the
offset:

$\varepsilon_{0}(k)=4\nu\cos(k)+2\mu_{0}\,,$ (9)

and the amplitude $A=2\mu_{1}$. The level crossing of the diabatic states
(eigenstates of $\tau_{z}$) occurs at $\varepsilon(k,t)=0$, which requires
$-0.09\pi\lesssim k\lesssim 0.09\pi$ for the specific parameters of the model
we set. Since the occupation probability of the upper adiabatic state is known
to be negligibly small if no diabatic level crossing occurs [55, 56], this is
consistent with the narrow $k$-range with significant electron population
transfer as shown in Fig. 2 (j). The minimal energy gap between the two
adiabatic (i.e. instantaneous) levels reads
$E_{G}^{\text{min}}(k)=\absolutevalue{\Delta(k)}=4\Delta\absolutevalue{\sin
k}$. It is reached simultaneously with the diabatic state level crossing, and
shows a linear behavior at small $k$ as plotted in Fig. 2 (i). As detailed in
references [55, 56], one can obtain an approximate analytic solution of the
LZS problem, in the slow and the fast-passage limits. These limits are
characterized by the dimensionless parameter
$\delta(k)\equiv[\Delta(k)]^{2}/[4v(k)]$, which characterizes the ratio of the
minimal gap $\absolutevalue{\Delta(k)}$ and the velocity
$v(k)\equiv\frac{\partial\varepsilon(k,t)}{\partial
t}|_{t_{c}:\varepsilon(k,t_{c})=0}=A\omega\sqrt{1-[\varepsilon_{0}(k)/A]^{2}}$
at the minimal gap. The velocity is finite and only weakly varying in the
region of interest close the $k=0$ before it rapidly drops to zero as
$k\rightarrow 0.09\pi$. As a result, we find $\delta(k)\propto k^{2}$ at small
$k$ and $\delta=1$ at $k_{1}\approx 0.03\pi$ (see Fig. 2 (k)). In terms of
this parameter, the LZ transition probability for a single passage reads
$P_{\text{LZ}}(k)\equiv e^{-2\pi\delta(k)}$. In the fast-passage limit at
$k\ll k_{1}$, we have $\delta(k)\ll 1$ and $1-P_{\text{LZ}}(k)\ll 1$; while in
the slow-passage limit, we have $\delta(k)\gg 1$ and $P_{\text{LZ}}(k)\ll 1$.
In Fig. 2 (k) we plot these driving speed indicators in the $k$-region with
significant excited state occupation. The fast-passage regime extends from the
$\Gamma$-point ($k=0$), where the minimal gap vanishes, roughly to the
position of the first peak of $W_{n}(k)$ at $k=0.0033\pi$ [see Fig. 2 (k)],
where $\delta(k)\lesssim 0.01$ and $1-P_{\text{LZ}}(k)\lesssim 0.06$. The
resonance condition in the diabatic regime reads $\varepsilon_{0}(k)=m\omega$
with integer $m$. This condition is fulfilled at $k=0.005\pi$ for $m=12$ and
at $k=0.017\pi$ for $m=11$ as shown in by the black vertical dashed lines in
Fig. 2 (i), which explains the two dominant peaks seen in Fig. 2 (j). Note
that we have $\delta(k=0.017\pi)=0.26$ and $1-P_{\text{LZ}}(k=0.017\pi)=0.81$,
which implies that $k=0.017\pi$ is near the crossover from fast to slow
passage. In fact, the resonance condition with $m=10$ gives $k=0.023\pi$ as
shown by a gray vertical dashed line in Fig. 2 (i), which is off from the
peaks in Fig. 2 (j) due to being located close to the crossover region with
$\delta(k=0.023\pi)=0.48$ and $1-P_{\text{LZ}}(k=0.023\pi)=0.95$. At larger
$k_{1}\approx 0.03\pi$, the crossover from fast to slow passage occurs and one
needs to use a more general resonance condition (see Eq. (56) in Ref. [56]).
Finally, for $\pi/20\leq k\leq 0.09\pi$ and beyond the slow-passage regime is
reached, where we find a small adiabatic state probability. The resonance
condition in the adiabatic regime reads $\frac{2A}{\pi\omega}=m$ with integer
$m$, which is not exactly fulfilled for our choice of parameters
$\frac{2A}{\pi\omega}=12.8$. One should also take into account that the
dynamics only involves five LZS oscillations due to damping effects in
experiments, which limits the total transfer into the upper adiabatic state
(even on resonance) when $P_{\text{LZ}}\ll 1$. We conclude that the majority
of the excited state population dynamics in our model occurs in the fast-
passage region and the crossover regime between the fast and the slow-passage
limits, and that the main peaks can be understood as arising from resonances
in the diabatic regime.

To summarize, our analysis demonstrates the crucial role of the phonon-induced
topological band closing for carrier excitation. This creates a finite
momentum space volume where effective TLSs experience an avoided level
crossing with a sufficiently small band gap such that carriers can be excited
through LZS tunneling.

###  II.3 First-principles quantum dynamics simulations

Model and ab initio simulation method. To gain a more material-specific
understanding of the carrier excitation dynamics of the phonon-modulated ZrTe5
system, we carry out first-principles simulations based on time-dependent
Schrödinger equation with DFT basis functions. The time-dependence of the KS
Hamiltonian $\hat{\mathcal{H}}_{0}(t)$ is encoded in the ionic trajectory
$\mathbf{R}(t)$ that is set by the $A_{1g}$ coherent phonon. In the
implementation of DFT for periodic systems one often adopts a basis set with
large dimension, such as plane waves. This renders a direct manipulation of
the DFT Hamiltonian $\hat{\mathcal{H}}_{0}(t)$ cumbersome. Note that
$\hat{\mathcal{H}}_{0}(t)$ generally covers higher-energy unoccupied states
and deeper occupied states, which are likely irrelevant for the carrier
excitation dynamics in phonon-modulated ZrTe5, which we expect to be dominated
by states close to the chemical potential. Standard tight-binding downfolding
approaches, such as maximally localized Wannier function [62] and quasi-atomic
minimal basis-set orbitals method [63, 64], can be useful; but the downfolding
calculation for many snapshots along the trajectory $\mathbf{R}(t)$ in the
simulation time period can be time-consuming, and the time-dependence of the
downfolded orbitals introduces additional complexity. Here we adopt an
alternative representation where the component of the dynamical electronic
state $\ket{\Psi(t)}=\otimes_{\mathbf{k}}\ket{\Psi(\mathbf{k},t)}$ is
approximated as a linear combination of $N_{\text{b}}$ adiabatic states
$\\{\ket{\Phi_{i}(\mathbf{k},t)}\\}$ generalized to a generic
$\mathbf{k}$-point from Ref. [65, 66]:

$\ket{\Psi(\mathbf{k},t)}=\sum_{i=1}^{N_{\text{b}}}c_{i}(\mathbf{k},t)\ket{\Phi_{i}(\mathbf{k},t)},$
(10)

where $\ket{\Phi_{i}(\mathbf{k},t)}\equiv\prod_{\mu\in
S_{i}}\phi_{\mu}^{\dagger}(\mathbf{k},t)\ket{0}$ is a noninteracting single
Slater determinant state defined by a set $S_{i}$ of occupied KS orbitals
$\phi_{\mu}(\mathbf{k},t)$, which satisfies
$\hat{\mathcal{H}}_{0}(\mathbf{k},\mathbf{R}(t))\ket{\phi_{\mu}(\mathbf{k},t)}=\epsilon_{\mu}(\mathbf{k},t)\ket{\phi_{\mu}(\mathbf{k},t)}$.
Here the crystal momentum $\mathbf{k}$ is conjugate to the position vector of
the simulation unit cell.

The propagation of $\ket{\Psi(\mathbf{k},t)}$ is encoded in the time-dependent
complex amplitudes $c_{i}(\mathbf{k},t)$ and the adiabatic states
$\ket{\Phi_{i}(\mathbf{k},t)}$. Substituting Eq. (10) into the time-dependent
Schrödinger equation leads to the equation of motion (EOM) of the amplitudes

$i\hbar\frac{\partial c_{i}(\mathbf{k},t)}{\partial
t}=\sum_{j=1}^{N_{b}}H_{ij}(\mathbf{k},t)c_{j}(\mathbf{k},t)\,.$ (11)

The vibronic Hamiltonian is given by

$H_{ij}(\mathbf{k},t)=\varepsilon_{i}(\mathbf{k},t)\delta_{ij}-i\hbar
d_{ij}(\mathbf{k},t)\,.$ (12)

Here, we define $\varepsilon_{i}(\mathbf{k},t)=\sum_{\mu\in
S_{i}}\epsilon_{\mu}(\mathbf{k},t)$. The complex nonadiabatic coupling (NAC)
coefficient between a pair of distinct states
$\\{\ket{\Phi_{i}},\ket{\Phi_{j}}\\}$ is given by
$d_{ij}=\langle\Phi_{i}|\frac{\partial}{\partial t}\,|\Phi_{j}\rangle$, which
is nonzero only if there is exactly one distinct occupied KS orbital between
$\ket{\Phi_{i}}$ and $\ket{\Phi_{j}}$ due to the single Slater determinant
nature [67, 68]. The NAC can be conveniently evaluated using the finite-
difference method [67]:

$\displaystyle d_{ij}(\mathbf{k},t)\approx\frac{1}{2dt}$
$\displaystyle\left(\langle\Phi_{i}(\mathbf{k},t)|\,\Phi_{j}(\mathbf{k},t+dt)\rangle\right.$
(13)
$\displaystyle\left.-\langle\Phi_{i}(\mathbf{k},t+dt)|\,\Phi_{j}(\mathbf{k},t)\rangle\right),~{}$

which is completely determined by the state overlap matrix between consecutive
time steps.

Figure 3:  Carrier excitation dynamics of phonon-modulated ZrTe5 system from
first-principles simulations. (a) Occupancy of the two lowest-energy
conduction bands, $n_{e}(i,\mathbf{k},t)$ with $i=2,3$, as a function of
simulation time $t$ for $5$ phonon periods starting with $t_{0}=0$ ( the
equilibrium configuration with zero phonon displacement) at
$\mathbf{k}=\mathbf{b}_{1}/720+\mathbf{b}_{2}/720+\mathbf{b}_{3}/80$. Here
$(\mathbf{b}_{1},\mathbf{b}_{2},\mathbf{b}_{3})$ are the reciprocal primitive
vectors (see Appendix IV for details). (b) Time-dependence of the complex NAC
amplitude $\hbar d_{13}(\mathbf{k},t)$ between the top valence band ($i=1$)
and lowest conduction band ($i=3$) at the same $\mathbf{k}$-point as (a), with
real part in yellow and imaginary part in red. The background dark blue to
light blue indicates a band gap,
$E_{\text{G}}=\varepsilon_{2}(\mathbf{k},t)-\varepsilon_{1}(\mathbf{k},t)$,
varying from $27$ meV to $132$ meV. (c) The total excited states population
$n_{e}(t)$ integrated over the Brillouin zone as a function of $t$ for the
simulation starting with $t_{0}=0$. (d) The time-evolution of normalized
$\bar{n}_{e}(t)/\bar{n}_{e}(5T)$(black line), where $\bar{n}_{e}(t)$ is the
total excited states population $n_{e}(t)$ averaged over simulations starting
with $t_{0}$ uniformly sampled by $10$-points in a phonon period. The
normalized experimental carrier density change $\Delta n(t)/\Delta n(5T)$ are
also plotted in orange line. We obtain $\bar{n}_{e}(5T)=1.2\times 10^{16}$cm-3
from the theoretical simulation, and $\Delta n(5T)=(0.28\pm 0.12)\times
10^{-16}$ cm-3 estimated from the experiment.

Technical details of the simulation. The computational complexity of the
simulation is tied to the number $N_{\text{b}}$ of adiabatic states that are
used in Eq. (10). The Raman $A_{1g}$ phonon in ZrTe5 has a frequency of $1.2$
THz, which is one order of magnitude smaller than the equilibrium band gap. We
can thus truncate the expansion at the level of single-electron excitations
between low energy bands. We include configurations with one electron excited
from the top valence band to one of the four lowest conduction bands. We thus
keep $N_{\text{b}}=5$ adiabatic states $\Phi_{i}(\mathbf{k},t)$ in the
expansion at each $\mathbf{k}$-point. Each $\Phi_{i}$ differs by exactly one
occupied KS orbital, which therefore can also be labelled by the same index
$i$, with $i=1$ corresponding to the top valence band, and $i=2\dots
N_{\text{b}}$ to the conduction bands in ascending order of energy. The
vibronic Hamiltonian (12) can be simplified by setting
$\varepsilon_{i}(\mathbf{k},t)=\epsilon_{i}(\mathbf{k},t)$ as a rigid
potential shift. The evaluation of the NAC (13) can also be reduced to

$\displaystyle d_{ij}(\mathbf{k},t)\approx\frac{1}{2dt}$
$\displaystyle\left(\langle\phi_{i}(\mathbf{k},t)|\,\phi_{j}(\mathbf{k},t+dt)\rangle\right.$
(14)
$\displaystyle\left.-\langle\phi_{i}(\mathbf{k},t+dt)|\,\phi_{j}(\mathbf{k},t)\rangle\right)\,,$

which depends only on the $N_{\text{b}}$ KS orbitals that are kept at each
$\mathbf{k}$-point.

The EOM (11) assumes the continuity of the time dependent basis
$\ket{\Phi_{i}(\mathbf{k},t)}$ with $t$. Therefore, it is crucial to fix the
phase degree of freedom of the KS orbitals $\phi_{i}(\mathbf{k},t)$. This can
be achieved by consecutively applying a phase factor to
$\phi_{i}(\mathbf{k},t+dt)\to e^{i\theta}\phi_{i}(\mathbf{k},t+dt)$, where the
phase $e^{i\theta}\equiv O_{i}^{*}/\absolutevalue{O_{i}}$ with
$O_{i}=\langle\phi_{i}(\mathbf{k},t)|\,\phi_{i}(\mathbf{k},t+dt)\rangle$ is
determined by the overlap with the same orbital at the previous time step.
Additional complexity of the dynamics simulations originates from the presence
of time-reversal and inversion symmetry, which renders every band doubly
degenerate, and strong spin-orbit coupling in ZrTe5. The ambiguity in the
doubly degenerate bands can be partially fixed by choosing a $S_{z}$-gauge
such that the $2\times 2$ spin $S_{z}$ matrix becomes diagonal in each doubly
degenerate manifold via a unitary transformation. Numerically, we find that
the $S_{z}$ gauge transformation is not sufficient to guarantee the orbital
continuity along the dynamical path. Therefore, we propose the following
overlap gauge correction to better address the band degeneracy problem. The
simulation starts with orbitals in the $S_{z}$-gauge, and apply unitary
transformation in each doubly degenerate manifold in all following time steps,
such that each rotated orbital $\phi_{i}(\mathbf{k},t+dt)$ has maximal overlap
with the same one at the previous step. This is achieved by diagonalizing a
series of $2\times 2$ matrices
$P_{rs}^{(j)}=\langle\phi_{r}^{(j)}(\mathbf{k},t+dt)|\hat{P_{j}}\,|\phi_{s}^{(j)}(\mathbf{k},t+dt)\rangle$
with
$\hat{P_{j}}=\ket{\phi_{1}^{(j)}(\mathbf{k},t)}\bra{\phi_{1}^{(j)}(\mathbf{k},t)}$.
Here $r,s\in[1,2]$ run through the two orbitals in the $j$th degenerate
doublet. The diagonalization gives two eigenvectors, where one has a nonzero
eigenvalue and is assigned to the first orbital of the $j$th doublet. The
other eigenvector has zero eigenvalue, and is assigned to the second orbital
in the doublet. Following this procedure, we numerically find that the self-
overlap of each wavefunction at consecutive time steps always remains above
$99.99\%$. The band index exchange between different degenerate doublets,
which can be detected by checking the overlap between wavefunctions at
consecutive time steps, is not observed in the simulations reported here.

First-principles quantum dynamics simulation results. The analysis of the one-
dimensional (1D) toy model results shows that the main contribution to the
carriers in the excited band resulted from LZS tunneling in a narrow region of
momentum space around the zone center $k\in[-\pi/20,\pi/20]$. Expecting a
similar behavior for the realistic 3D model of ZrTe5, we use a dense
$360\times 360\times 120$ uniform $\mathbf{k}$-grid covering the full
Brillouin zone for the following dynamics simulations and use a shift of
$\Delta=\frac{0.5}{360}(\mathbf{b}_{1}+\mathbf{b}_{2})+\frac{0.5}{120}\mathbf{b}_{3}$
from the $\Gamma$-point for each $\mathbf{k}$-point. Here, $\mathbf{b}_{i}$
are the reciprocal basis vectors. The center-shifted $\mathbf{k}$-mesh
therefore excludes the $\Gamma$-point, where the band gap closes at certain
times and additional gauge correction is otherwise needed.

In Fig. 3(a) we present the excited state population $n_{e}(i,\mathbf{k},t)$
as a function of simulation time $t$ for five phonon cycles, starting at
$t_{0}=0$ with zero phonon displacement. Here we define

$n_{e}(i,\mathbf{k},t)\equiv\absolutevalue{c_{i}(\mathbf{k},t)}^{2}\,,$ (15)

which is equivalent to the definition used in the toy model analysis. At
momentum $\mathbf{k}=\mathbf{b}_{3}/120+\Delta$ adjacent to the zone center,
an electron is gradually excited from the top valence band ($i=1$) to the two
lowest conduction bands ($i=2,3$) due to nonadiabatic effects. The complex NAC
amplitude $\hbar d_{13}(\mathbf{k},t)$ between the $1$st and $3$rd bands at
the same $\mathbf{k}$-point is plotted in Fig. 3(b). The yellow line denotes
the real part and the red line the imaginary part. The line width indicates
the numerical noise, which is found to have negligible impact on the state
population dynamics. The maximal difference in the dynamical state populations
from the simulation using the (noisy) NAC amplitudes versus using smoothed
data via application of a Savitzky–Golay filter is only about $10^{-3}$.
Clearly, Fig. 3(a,b) shows that a sharp transition of the state population
$n_{e}$ occurs at the peaks of the NAC. This also coincides with a minimum of
the band gap,
$E_{\text{G}}=\varepsilon_{2}(\mathbf{k},t)-\varepsilon_{1}(\mathbf{k},t)$, as
indicated by the blue shading in the background of Fig. 3(b). We note that
$d_{12}(\mathbf{k},t)$ has a similar time-dependence as $d_{13}(\mathbf{k},t)$
[both bands have degenerate energies
$\varepsilon_{2}(\mathbf{k},t)=\varepsilon_{3}(\mathbf{k},t)]$, yet with
slightly larger amplitude. This difference in the NAC induces a larger
electron occupancy in the $i=2$ band compared to the one with $i=3$, as shown
in Fig. 3(a). In contrast, the NACs from $i=1$ to $i=4,5$ are smaller by more
than one order of magnitude, resulting in negligibly small carrier excitations
to these bands, $n_{e}(i,\mathbf{k},t)<10^{-6}$ for $i=4,5$.

The first-principles dynamics simulations allow for a direct and quantitative
comparison to experiment. First, in Fig. 3(c) we show the excited state
carrier density,
$n_{e}(t)=\sum_{\mathbf{k}}w_{\mathbf{k}}\sum_{i=2}^{3}n_{e}(i,\mathbf{k},t)$
as a function of time $t$. It increases from zero to about $n_{e}(t=5T)\approx
2.4\times 10^{16}$cm-3 at the end of the simulation $t=5T$ and exhibits
qualitatively similar sharp transitions near dynamical band gap minimum as the
momentum resolved quantity in Fig. 3(a). Within the $360\times 360\times 120$
uniform $\mathbf{k}$-mesh of the Brillouin zone, we find that the dominant
contributions come from
$\mathbf{k}=\frac{l}{360}\mathbf{b}_{1}+\frac{m}{360}\mathbf{b}_{2}+\frac{n}{120}\mathbf{b}_{3}+\Delta$
with $l,m,n=0,\pm 1,\pm 2,-3$. Next, we we account for the fact that in pump-
probe experiments, the time-trace of differential transmission is obtained as
an average over multiple runs and ZrTe5 samples exhibit some degree of
electronic heterogeneity and nanostrip junctions, as observed in the THz
nanoimaging [69]. To capture these phenomena on average, we define
$\bar{n}_{e}(t)=\frac{1}{10}\sum_{i=0}^{9}n_{e}(t)|_{t_{0}=\frac{i}{10}T}$,
which is an average over simulations at $10$ different starting times. In Fig.
3 (d) we directly compare the time-dependence of a normalized
$\bar{n}_{e}(t)/\bar{n}_{e}(t=5T)$ (black), to the experimental data $\Delta
n(t)/\Delta n(t=5T)$ (orange). Both curves exhibit a similar growth pattern of
the carrier density over time, and we also find the carrier density at the end
of the simulation $\bar{n}_{e}(t=5T)\approx 2.4\times 10^{16}$ cm-3 to be in
good agreement with the one estimated from experiment $\Delta n(5T)=(0.28\pm
0.12)\times 10^{-16}$ cm-3, considering that there is electronic heterogeneity
present in the experimental ZrTe5 sample [69]. The numerical estimation of
$n_{e}(t=5T)$ and $\bar{n}_{e}(t=5T)$ reported here also includes a factor of
$2$ to take into account the double degeneracy of the top valence band in
ZrTe5 system.

##  III Conclusions

We report detailed first-principle and effective model simulations of the
carrier excitation dynamics in coherent phonon-modulated ZrTe5. Our results
shed new light on recent pump-probe experiments [25] by providing a clear
intuitive explanation of the experimental results. Both first-principle and
effective model calculations reveal the importance of the phonon-induced
topological phase transition in ZrTe5 and the associated closing of the bulk
gap for the observed excitation of carriers. We show that the excitations
occur via Landau-Zener-Stückelberg tunneling in a series of time-dependent
avoided level crossings of Bloch states located in a narrow region of momentum
space around the zone center. Our detailed time-dependent Schrödinger equation
simulations further show that the dominant tunneling occurs between the
highest valence band and the lowest doubly degenerate conduction bands, while
excitations to the next higher bands are negligible due to small transition
matrix elements. We predict that the carrier density increases gradually with
time and reaches a final value of $2.4\times 10^{16}$cm-3 at $t=5T$ when
phonon coherence is lost in the experiment. These results are in good
quantitative agreement with experiment. Our work thus demonstrates that the
coherent charge excitation process in topological quantum materials such as
ZrTe5 can be understood and predicted quantitatively by first-principles
quantum dynamics simulations.

Figure 4:  DFT band structure calculations for the adiabatic topological phase
switching induced by the $A_{\text{1g}}$ Raman phonon mode in ZrTe5. (a-d)
Band structure along $\mathbf{k}$-path $Z-\Gamma-Y$ with distortion parameter
$\lambda=-3.0,-2.2,0.0,and3.0$, respectively. The red color encodes projection
weight of $5p$-orbitals of Ted to the wavefunctions. We define $Z=(0,0,0.5)$
and $Y=(-0.5,0.5,0)$ in the reciprocal space with basis vectors
$(\mathbf{b}_{1},\mathbf{b}_{2},\mathbf{b}_{3})$. (e) Band gap at the zone
center, $E_{\text{G}}(\Gamma)$, as a function of distortion parameter
$\lambda$. Here $\lambda=1$ represented $0.033$ Ådisplacement for Zr, $0.035$
Åfor Ted, $0.032$ Åfor Tea, and $0.017$ Åfor Tez [25]. Red circles indicate
the $\lambda$ values for panels (a-d). The STI region is highlighted in blue
and WTI region in white, with boundary line indicating the Dirac point
position. (f) The $A_{\text{1g}}$ phonon mode in the conventional cell of
ZrTe5. Green spheres represent Zr atoms, orange spheres for apical Te atoms
(Tea), silver spheres for zigzag Te atoms (Tez), and purple spheres for
dimerized Te atoms (Ted). The arrows indicate the atomic displacement vectors
of the $A_{\text{1g}}$ mode.

##  IV Methods

First-principles total energy and electronic-structure calculations for ZrTe5
are based on DFT with the exchange correlation functional in generalized
gradient approximation parametrized by Perdew, Burke, and Ernzerhof (PBE)
[70]. Van der Waals interaction is included by Grimme’s damped atom-pairwise
dispersion corrections (D2) [71]. The calculations are performed using the
Vienna Ab initio Simulation Package (VASP) [72]. We use a plane-wave cutoff
energy of $230$ eV and include spin-orbit coupling for all the calculations.
The phonon modes are calculated using the finite displacement approach as
implemented in Phonopy [73]. Specifically, we use the primitive unit cell of
experimental structure($a=3.987$ Å, $b=14.502$ Å, and $c=13.727$ Å) [74]. The
primitive vectors are $\mathbf{a}_{1}=(1.994,-7.251,0)$ Å,
$\mathbf{a}_{2}=(1.994,7.251,0)$ Å, and $\mathbf{a}_{3}=(0,0,13.727)$ Å. The
corresponding reciprocal primitive vectors are
$\mathbf{b}_{1}=(0.251,-0.069,0)2\pi$ Å-1,
$\mathbf{b}_{2}=(0.251,0.069,0)2\pi$ Å-1, and $\mathbf{b}_{3}=(0,0,0.073)2\pi$
Å-1. Highly accurate wavefunctions at specific $\mathbf{k}$-points are
generated for NAC calculations by setting an energy convergence criterion to
$10^{-9}$ eV and requiring a minimum of $20$ electronic steps.

To be self-contained, we present the key DFT band structure calculation
results for the description of the $A_{\text{1g}}$ phonon-induced adiabatic
topological phase transition observed in ZrTe5 system as reported in Ref.
[25]. Under the modulation of the $A_{\text{1g}}$ eigenmode as plotted in Fig.
4(f), the band gap $E_{\text{G}}(\Gamma)$ at zone center closes at distortion
parameter $\lambda=-2.2$ as shown in Fig. 4(e), implying a topological phase
transition along the dynamical path. This is confirmed by the band structure
analysis and topological invariant index calculation [25]. In Fig. 4(a-d), we
plot the band structure along high-symmetry $\mathbf{k}$-path $Z-\Gamma-Y$ at
phonon distortion parameter $\lambda=-3.0,-2.2,0,3.0$, decorated with red
color indicating the Ted $5p$-orbital weight. Band inversion clearly occurs
when $\lambda$ passes through $\lambda=-2.2$ the Dirac point. For
$\lambda<-2.2$, the adiabatic state of the system is in WTI, and switches to
STI for $\lambda>-2.2$.

## Data availability

All the data to generate the ﬁgures are available at figshare [75]. Data
supporting the calculations are available together with the codes at figshare
[76]. All other data are available from the corresponding authors on
reasonable request.

## Code availability

All the computer codes developed and used in this work are available open-
source at figshare [76].

## References

  * Först  _et al._ [2011] M. Först, C. Manzoni, S. Kaiser, Y. Tomioka, Y.-n. Tokura, R. Merlin, and A. Cavalleri, Nonlinear phononics as an ultrafast route to lattice control, Nat. Phys. 7, 854 (2011).
  * Merlin [1997] R. Merlin, Generating coherent thz phonons with light pulses, Solid State Commun. 102, 207 (1997).
  * Dhar  _et al._ [1994] L. Dhar, J. A. Rogers, and K. A. Nelson, Time-resolved vibrational spectroscopy in the impulsive limit, Chem. Rev. 94, 157 (1994).
  * Basov  _et al._ [2017] D. Basov, R. Averitt, and D. Hsieh, Towards properties on demand in quantum materials, Nat. Mater. 16, 1077 (2017).
  * Mankowsky  _et al._ [2016] R. Mankowsky, M. Först, and A. Cavalleri, Non-equilibrium control of complex solids by nonlinear phononics, Rep. Prog. Phys. 79, 064503 (2016).
  * de la Torre  _et al._ [2021] A. de la Torre, D. M. Kennes, M. Claassen, S. Gerber, J. W. McIver, and M. A. Sentef, Colloquium: Nonthermal pathways to ultrafast control in quantum materials, Rev. Mod. Phys. 93, 041002 (2021).
  * Disa  _et al._ [2021] A. S. Disa, T. F. Nova, and A. Cavalleri, Engineering crystal structures with light, Nat. Phys. 17, 1087 (2021).
  * Stremoukhov  _et al._ [2022] P. Stremoukhov, A. Safin, S. Nikitov, A. Kirilyuk, _et al._ , Phononic manipulation of antiferromagnetic domains in nio, New J. Phys. 24, 023009 (2022).
  * Stupakiewicz  _et al._ [2021] A. Stupakiewicz, C. Davies, K. Szerenos, D. Afanasiev, K. Rabinovich, A. Boris, A. Caviglia, A. Kimel, and A. Kirilyuk, Ultrafast phononic switching of magnetization, Nat. Phys. 17, 489 (2021).
  * Nova  _et al._ [2017] T. F. Nova, A. Cartella, A. Cantaluppi, M. Först, D. Bossini, R. V. Mikhaylovskiy, A. V. Kimel, R. Merlin, and A. Cavalleri, An effective magnetic field from optically driven phonons, Nat. Phys. 13, 132 (2017).
  * Rini  _et al._ [2007] M. Rini, R. Tobey, N. Dean, J. Itatani, Y. Tomioka, Y. Tokura, R. W. Schoenlein, and A. Cavalleri, Control of the electronic phase of a manganite by mode-selective vibrational excitation, Nature 449, 72 (2007).
  * Caviglia  _et al._ [2012] A. D. Caviglia, R. Scherwitzl, P. Popovich, W. Hu, H. Bromberger, R. Singla, M. Mitrano, M. C. Hoffmann, S. Kaiser, P. Zubko, S. Gariglio, J.-M. Triscone, M. Först, and A. Cavalleri, Ultrafast strain engineering in complex oxide heterostructures, Phys. Rev. Lett. 108, 136801 (2012).
  * Horstmann  _et al._ [2020] J. G. Horstmann, H. Böckmann, B. Wit, F. Kurtz, G. Storeck, and C. Ropers, Coherent control of a surface structural phase transition, Nature 583, 232 (2020).
  * Song  _et al._ [2023] B. Song, X. Yang, C. Sundahl, J.-H. Kang, M. Mootz, Y. Yao, I. Perakis, L. Luo, C. Eom, and J. Wang, Ultrafast martensitic phase transition driven by intense terahertz pulses, Ultrafast Sci. 3, 0007 (2023).
  * Kaiser  _et al._ [2014] S. Kaiser, C. R. Hunt, D. Nicoletti, W. Hu, I. Gierz, H. Liu, M. Le Tacon, T. Loew, D. Haug, B. Keimer, _et al._ , Optically induced coherent transport far above t c in underdoped yba 2 cu 3 o 6+ $\delta$, Phys. Rev. B 89, 184516 (2014).
  * Hu  _et al._ [2014] W. Hu, S. Kaiser, D. Nicoletti, C. R. Hunt, I. Gierz, M. C. Hoffmann, M. Le Tacon, T. Loew, B. Keimer, and A. Cavalleri, Optically enhanced coherent transport in yba 2 cu 3 o 6.5 by ultrafast redistribution of interlayer coupling, Nat. Mater. 13, 705 (2014).
  * Mankowsky  _et al._ [2014] R. Mankowsky, A. Subedi, M. Forst, S. O. Mariager, M. Chollet, H. T. Lemke, J. S. Robinson, J. M. Glownia, M. P. Minitti, A. Frano, M. Fechner, N. A. Spaldin, T. Loew, B. Keimer, A. Georges, and A. Cavalleri, Nonlinear lattice dynamics as a basis for enhanced superconductivity in YBa2Cu3O6.5, Nature 516, 71 (2014).
  * Mitrano  _et al._ [2016] 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 k3c60 at high temperature, Nature 530, 461 (2016).
  * Garate [2013] I. Garate, Phonon-induced topological transitions and crossovers in dirac materials, Phys. Rev. Lett. 110, 046402 (2013).
  * Saha and Garate [2014] K. Saha and I. Garate, Phonon-induced topological insulation, Phys. Rev. B 89, 205103 (2014).
  * Kim and Jhi [2015] J. Kim and S.-H. Jhi, Topological phase transitions in group iv-vi semiconductors by phonons, Phys. Rev. B 92, 125142 (2015).
  * Wang  _et al._ [2017] L.-L. Wang, N. H. Jo, Y. Wu, Q. Wu, A. Kaminski, P. C. Canfield, and D. D. Johnson, Phonon-induced topological transition to a type-II Weyl semimetal, Phys. Rev. B 95, 165114 (2017).
  * Weber  _et al._ [2018] C. P. Weber, M. G. Masten, T. C. Ogloza, B. S. Berggren, M. K. L. Man, K. M. Dani, J. Liu, Z. Mao, D. D. Klug, A. A. Adeleke, and Y. Yao, Using coherent phonons for ultrafast control of the dirac node of ${\mathrm{SrMnSb}}_{2}$, Phys. Rev. B 98, 155115 (2018).
  * Sie  _et al._ [2019] E. J. Sie, C. M. Nyby, C. Pemmaraju, S. J. Park, X. Shen, J. Yang, M. C. Hoffmann, B. Ofori-Okai, R. Li, A. H. Reid, _et al._ , An ultrafast symmetry switch in a weyl semimetal, Nature 565, 61 (2019).
  * Vaswani  _et al._ [2020] C. Vaswani, L.-L. Wang, D. H. Mudiyanselage, Q. Li, P. M. Lozano, G. D. Gu, D. Cheng, B. Song, L. Luo, R. H. J. Kim, C. Huang, Z. Liu, M. Mootz, I. E. Perakis, Y. Yao, K. M. Ho, and J. Wang, Light-driven Raman coherence as a nonthermal route to ultrafast topology switching in a dirac semimetal, Phys. Rev. X 10, 021013 (2020).
  * Luo  _et al._ [2021] L. Luo, D. Cheng, B. Song, L.-L. Wang, C. Vaswani, P. Lozano, G. Gu, C. Huang, R. H. Kim, Z. Liu, _et al._ , A light-induced phononic symmetry switch and giant dissipationless topological photocurrent in ZrTe5, Nat. Mater. 20, 329 (2021).
  * Wang [2021] L.-L. Wang, Expansive open fermi arcs and connectivity changes induced by infrared phonons in $\mathrm{Zr}{\mathrm{Te}}_{5}$, Phys. Rev. B 103, 075105 (2021).
  * Yang  _et al._ [2020] X. Yang, L. Luo, C. Vaswani, X. Zhao, Y. Yao, D. Cheng, Z. Liu, R. H. Kim, X. Liu, M. Dobrowolska-Furdyna, _et al._ , Light control of surface–bulk coupling by terahertz vibrational coherence in a topological insulator, npj Quantum Mater. 5, 13 (2020).
  * Yang  _et al._ [2018] X. Yang, C. Vaswani, C. Sundahl, M. Mootz, P. Gagel, L. Luo, J. Kang, P. Orth, I. Perakis, C. Eom, _et al._ , Terahertz-light quantum tuning of a metastable emergent phase hidden by superconductivity, Nat. Mater. 17, 586 (2018).
  * Prezhdo [2021] O. V. Prezhdo, Modeling non-adiabatic dynamics in nanoscale and condensed matter systems, Acc. Chem. Res. 54, 4239 (2021).
  * Lindh and González [2020] R. Lindh and L. González,  _Quantum Chemistry and Dynamics of Excited States: Methods and Applications_ (John Wiley & Sons, 2020).
  * Nelson  _et al._ [2020] T. R. Nelson, A. J. White, J. A. Bjorgaard, A. E. Sifain, Y. Zhang, B. Nebgen, S. Fernandez-Alberti, D. Mozyrsky, A. E. Roitberg, and S. Tretiak, Non-adiabatic excited-state molecular dynamics: Theory and applications for modeling photophysics in extended molecular materials, Chem. Rev. 120, 2215 (2020).
  * Curchod and Martínez [2018] B. F. Curchod and T. J. Martínez, Ab initio nonadiabatic quantum molecular dynamics, Chem. Rev. 118, 3305 (2018).
  * Caruso and Zacharias [2023] F. Caruso and M. Zacharias, Quantum theory of light-driven coherent lattice dynamics, Phys. Rev. B 107, 054102 (2023).
  * Subedi  _et al._ [2014] A. Subedi, A. Cavalleri, and A. Georges, Theory of nonlinear phononics for coherent light control of solids, Phys. Rev. B 89, 220301 (2014).
  * Gu and Rondinelli [2018] M. Gu and J. M. Rondinelli, Nonlinear phononic control and emergent magnetism in mott insulating titanates, Phys. Rev. B 98, 024102 (2018).
  * Raines  _et al._ [2015] Z. M. Raines, V. Stanev, and V. M. Galitski, Enhancement of superconductivity via periodic modulation in a three-dimensional model of cuprates, Phys. Rev. B 91, 184506 (2015).
  * Komnik and Thorwart [2016] A. Komnik and M. Thorwart, BCS theory of driven superconductivity, Eur. Phys. J. B 89, 244 (2016).
  * Knap  _et al._ [2016] M. Knap, M. Babadi, G. Refael, I. Martin, and E. Demler, Dynamical Cooper pairing in nonequilibrium electron-phonon systems, Phys. Rev. B 94, 214504 (2016).
  * Patel and Eberlein [2016] A. A. Patel and A. Eberlein, Light-induced enhancement of superconductivity via melting of competing bond-density wave order in underdoped cuprates, Phys. Rev. B 93, 195139 (2016).
  * Kennes  _et al._ [2017] D. M. Kennes, E. Y. Wilner, D. R. Reichman, and A. J. Millis, Transient superconductivity from electronic squeezing of optically pumped phonons, Nat. Phys. 13, 479 (2017).
  * Babadi  _et al._ [2017] M. Babadi, M. Knap, I. Martin, G. Refael, and E. Demler, Theory of parametrically amplified electron-phonon superconductivity, Phys. Rev. B 96, 014512 (2017).
  * Mazza and Georges [2017] G. Mazza and A. Georges, Nonequilibrium superconductivity in driven alkali-doped fullerides, Phys. Rev. B 96, 064515 (2017).
  * Murakami  _et al._ [2017] Y. Murakami, N. Tsuji, M. Eckstein, and P. Werner, Nonequilibrium steady states and transient dynamics of conventional superconductors under phonon driving, Phys. Rev. B 96, 045125 (2017).
  * Sentef [2017] M. A. Sentef, Light-enhanced electron-phonon coupling from nonlinear electron-phonon coupling, Phys. Rev. B 95, 205111 (2017).
  * Schütt  _et al._ [2018] M. Schütt, P. P. Orth, A. Levchenko, and R. M. Fernandes, Controlling competing orders via nonequilibrium acoustic phonons: Emergence of anisotropic effective electronic temperature, Phys. Rev. B 97, 035135 (2018).
  * Luo  _et al._ [2023] L. Luo, M. Mootz, J.-H. Kang, C. Huang, K. Eom, J. Lee, C. Vaswani, Y. Collantes, E. Hellstrom, I. E. Perakis, _et al._ , Quantum coherence tomography of light-controlled superconductivity, Nat. Phys. 19, 201 (2023).
  * Kim  _et al._ [2021a] R. H. Kim, C. Huang, Y. Luan, L.-L. Wang, Z. Liu, J.-M. Park, L. Luo, P. M. Lozano, G. Gu, D. Turan, _et al._ , Terahertz nano-imaging of electronic strip heterogeneity in a dirac semimetal, ACS photonics 8, 1873 (2021a).
  * Liu  _et al._ [2020a] Z. Liu, C. Vaswani, X. Yang, X. Zhao, Y. Yao, Z. Song, D. Cheng, Y. Shi, L. Luo, D.-H. Mudiyanselage, _et al._ , Ultrafast control of excitonic rashba fine structure by phonon coherence in the metal halide perovskite ch 3 nh 3 pbi 3, Phys. Rev. Lett. 124, 157401 (2020a).
  * Liu  _et al._ [2020b] Z. Liu, C. Vaswani, L. Luo, D. Cheng, X. Yang, X. Zhao, Y. Yao, Z. Song, R. Brenes, R. J. H. Kim, J. Jean, V. Bulović, Y. Yan, K.-M. Ho, and J. Wang, Coherent band-edge oscillations and dynamic longitudinal-optical phonon mode splitting as evidence for polarons in perovskites, Phys. Rev. B 101, 115125 (2020b).
  * Landau [1932] L. D. Landau, Zur Theorie der Energieübertragung. II, Phys. Z. Sowjetunion 2, 46 (1932).
  * Zener [1932] C. Zener, Non-adiabatic Crossing of Energy Levels, Proc. R. Soc. London, Ser. A 137, 696 (1932).
  * Stückelberg [1932] E. C. G. Stückelberg, Theorie der unelastischen Stösse zwischen Atomen, Helv. Phys. Acta 5, 369 (1932).
  * Majorana [1932] E. Majorana, Atomi orientati in campo magnetico variabile, Nuovo Cim. 9, 43 (1932).
  * Shevchenko  _et al._ [2010] S. N. Shevchenko, S. Ashhab, and F. Nori, Landau–Zener–Stückelberg interferometry, Phys. Rep. 492, 1 (2010).
  * Ivakhnenko  _et al._ [2023] V. Ivakhnenko, S. N. Shevchenko, and F. Nori, Nonadiabatic Landau–Zener–Stückelberg–Majorana transitions, dynamics, and interference, Phys. Rep. 995, 1 (2023).
  * [57] A. Akhmerov, J. Sau, B. van Heck, S. Rubbert, R. Skolasiński, B. Nijholt, I. Muhammad, and T. örn Rosdahl, Topology in Condensed Matter: Tying Quantum Knots, https://topocondmat.org/ .
  * Aryal  _et al._ [2021] N. Aryal, X. Jin, Q. Li, A. M. Tsvelik, and W. Yin, Topological phase transition and phonon-space dirac topology surfaces in ${\mathrm{ZrTe}}_{5}$, Phys. Rev. Lett. 126, 016401 (2021).
  * Holstein [1959] T. Holstein, Studies of polaron motion: Part i. the molecular-crystal model, Ann. Phys. 8, 325 (1959).
  * Trotter [1959] H. F. Trotter, On the product of semi-groups of operators, Proc. Am. Math. Soc. 10, 545 (1959).
  * Nielsen and Chuang [2011] M. A. Nielsen and I. L. Chuang,  _Quantum Computation and Quantum Information: 10th Anniversary Edition_, 10th ed. (Cambridge University Press, New York, USA, 2011).
  * Marzari  _et al._ [2012] N. Marzari, A. A. Mostofi, J. R. Yates, I. Souza, and D. Vanderbilt, Maximally localized wannier functions: Theory and applications, Rev. Mod. Phys. 84, 1419 (2012).
  * Chan  _et al._ [2007] T.-L. Chan, Y. Yao, C. Wang, W. Lu, J. Li, X. Qian, S. Yip, and K. Ho, Highly localized quasiatomic minimal basis orbitals for mo from ab initio calculations, Phys. Rev. B 76, 205119 (2007).
  * Qian  _et al._ [2008] X. Qian, J. Li, L. Qi, C.-Z. Wang, T.-L. Chan, Y.-X. Yao, K.-M. Ho, and S. Yip, Quasiatomic orbitals for ab initio tight-binding analysis, Phys. Rev. B 78, 245112 (2008).
  * Granucci  _et al._ [2012] G. Granucci, M. Persico, and G. Spighi, Surface hopping trajectory simulations with spin-orbit and dynamical couplings, J. Chem. Phys. 137, 22A501 (2012).
  * Li  _et al._ [2018] W. Li, L. Zhou, O. V. Prezhdo, and A. V. Akimov, Spin–orbit interactions greatly accelerate nonradiative dynamics in lead halide perovskites, ACS Energy Lett. 3, 2159 (2018).
  * Hammes-Schiffer and Tully [1994] S. Hammes-Schiffer and J. C. Tully, Proton transfer in solution: Molecular dynamics with quantum transitions, J. Chem. Phys. 101, 4657 (1994).
  * Akimov and Prezhdo [2013] A. V. Akimov and O. V. Prezhdo, The pyxaid program for non-adiabatic molecular dynamics in condensed matter systems, J. Chem. Theory Comput. 9, 4959 (2013).
  * Kim  _et al._ [2021b] R. H. Kim, C. Huang, Y. Luan, L.-L. Wang, Z. Liu, J.-M. Park, L. Luo, P. M. Lozano, G. Gu, D. Turan, _et al._ , Terahertz nano-imaging of electronic strip heterogeneity in a Dirac semimetal, ACS photonics 8, 1873 (2021b).
  * Perdew  _et al._ [1996] J. P. Perdew, K. Burke, and M. Ernzerhof, Generalized gradient approximation made simple, Phys. Rev. Lett. 77, 3865 (1996).
  * Grimme [2006] S. Grimme, Semiempirical gga-type density functional constructed with a long-range dispersion correction, J. Comput. Chem. 27, 1787 (2006).
  * Kresse and Furthmüller [1996] G. Kresse and J. Furthmüller, Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set, Phys. Rev. B 54, 11169 (1996).
  * Togo and Tanaka [2015] A. Togo and I. Tanaka, First principles phonon calculations in materials science, Scr. Mater. 108, 1 (2015).
  * Matkovic and Matkovic [1992] T. Matkovic and P. Matkovic, Constitutional study of the titanium, zirconium and hafnium tellurides, Metalurgija (Zagreb) 31, 107 (1992).
  * Jiang  _et al._ [2023] T. Jiang, P. P. Orth, L. Luo, J. Wang, and Y.-X. Yao, Data set and scripts for the analysis and figures of quantum dynamics simulations of ZrTe5 system, https://doi.org/10.6084/m9.figshare.23786355 (2023).
  * Yao [2023] Y.-X. Yao, CyQuanDyn: coherent phonon-driven quantum dynamics simulaiton toolkit, https://doi.org/10.6084/m9.figshare.23774946 (2023).

## Acknowledgements

This work was supported by the U.S. Department of Energy (DOE), Office of
Science, Basic Energy Sciences, Materials Science and Engineering Division,
including the grant of computer time at the National Energy Research
Scientific Computing Center (NERSC) in Berkeley, California. The research was
performed at the Ames National Laboratory, which is operated for the U.S. DOE
by Iowa State University under Contract No. DE-AC02-07CH11358.

## Author contributions

T.J. performed the DFT and quantum dynamics simulations. Y.X.Y. and T.J. wrote
the codes for dynamics simulation and analysis. P.P.O. initiated the effective
model simulation and contributed to the analysis of the results. L.L. and J.W.
performed the experimental analysis. L.L.W. and F.Z. helped with the DFT
calculations and analysis. J.Z., C.Z.W., and K.M.H. provided inputs for the
first-principle dynamics simulations. Y.X.Y., T.J., and P.P.O. wrote the
paper, with contributions from all the authors. Y.X.Y. supervised the project.

## Competing interests

The authors declare no competing interests.

Supplemental materials for:
Ab-initio Simulations of Coherent Phonon-Induced Pumping of Carriers in
Zirconium Pentatelluride

##  I Alternative way for excited state population analysis

Fig. S1:  Comparison of excited state population dynamics between two
definitions at $k=0.033\pi$. (a) Time evolution of the excited state
population $n_{e}(k,t)$ defined with respect to the adiabatic conduction band
for a periodic modulation $\mu(t)/\nu\in[-2.02,-1.92]$ starting at
$\mu(t_{0})=\mu_{0}$ for five full cycles. (b) The same as (a) but for
$n^{0}_{e}(k,t)$ defined with respect to the initial conduction band. (c) Time
evolution of the energy gap $E_{G}(k,t)$. (d) The absolute value of the
overlap between adiabatic conduction band wavefunction $\psi_{c}(k,t)$ and the
initial conduction band wavefunction $\psi_{c}(k,t=0)$ as a function of $t/T$.

For the carrier density analysis of the time-dependent BdG model in the main
text, we define the excited state population
$n_{e}(k,t)=\absolutevalue{\langle\psi(k,t)|\,\psi_{c}(k,t)\rangle}^{2}$ with
respect to the adiabatic conduction band. Here we investigate an alternative
definition of
$n^{0}_{e}(k,t)=\absolutevalue{\langle\psi(k,t)|\,\psi_{c}(k,t=0)\rangle}^{2}$,
which is defined with respect to the conduction band of the initial
unperturbed Hamiltonian.

In Fig. S1 (a,b), we plot the time evolution of $n_{e}(k,t)$ and
$n^{0}_{e}(k,t)$ at $k=0.033\pi$ starting at $\mu(t_{0})=\mu_{0}$ for five
full cycles. The major variation of $n_{e}(k,t)$ occurs reasonably when the
system approaches the minimum of the band gap $E_{G}(k,t)$, as shown in Fig.
S1 (c). In contrast, the curve of $n^{0}_{e}(k,t)$ shows strong oscillation
even away from band gap minimum, which implies $n^{0}_{e}(k,t)$ is not a
physically reasonable definition in practice. Nevertheless, we observe that if
one discards the fast oscillation, $n^{0}_{e}(k,t)$ reaches a value of about
$0.2$, in a fairly good agreement with $n_{e}(k,t)$. This can be understood by
checking the overlap between two reference conduction band wavefunctions
$\absolutevalue{\langle\psi_{c}(k,t)|\,\psi_{c}(k,t=0)\rangle}$, as plotted in
Fig. S1 (d). $\psi_{c}(k,t)$ agrees with $\psi_{c}(k,t=0)$ at $t=\frac{i}{2}T$
with $i=0,\dots 10$ where the perturbation in $\mu(t)$ vanishes.
$\psi_{c}(k,t)$ deviates more from $\psi_{c}(k,t=0)$ when the perturbation in
$\mu(t)$ increases, and the deviation grows more rapidly when the perturbation
reduces the band gap.

##  II Significance of gauge and phase corrections for the ab initio quantum
dynamics simulations

Fig. S2:  Importance of phase correction on the simulation results. (a)
Occupancy of the two lowest-energy conduction bands, $n_{e}(i,\mathbf{k},t)$
with $i=2,3$, as a function of simulation time $t$ for $5$ phonon periods
starting with $t_{0}=0$ (the equilibrium configuration with zero phonon
displacement) at
$\mathbf{k}=\mathbf{b}_{1}/720+\mathbf{b}_{2}/720+\mathbf{b}_{3}/80$. Here
$(\mathbf{b}_{1},\mathbf{b}_{2},\mathbf{b}_{3})$ are the reciprocal primitive
vectors. (b, c) The real and imaginary part of the time-dependent complex NAC
amplitude $\hbar d_{13}(\mathbf{k},t)$ between the top valence band ($i=1$)
and lowest conduction band ($i=3$) at the same $\mathbf{k}$-point as (a), with
real part in yellow and imaginary part in red. The background dark blue to
light blue indicates a band gap,
$E_{\text{G}}=\varepsilon_{2}(\mathbf{k},t)-\varepsilon_{1}(\mathbf{k},t)$,
varying from $27$ meV to $132$ meV. (d, e, f) Similar results to (a, b, c)
with the exception that the phase correction is not applied to the
calculation. (a, b, c) is the reproduction of fig.3(a, b) to facilitate the
comparison. To clearly show the random jumps in the real and imaginary part of
NAC due to the uncorrected random phases in the wavefunctions, we plot the
curves in (b,c,e,f) for every $20$ points of the time mesh.

In the main text, we have discussed the $S_{z}$ and overlap gauge corrections
and phase correction for the wavefunctions, which are crucial to produce
physically reasonable quantum dynamics simulation results by fixing the
artificial random jumps in the otherwise uncorrected time-dependent NAC
amplitudes. For a specific example, we compare the simulation results with and
without phase correction in Fig. S2 (a) and (c). The variation of the
conduction band occupations $n_{e}(i\in[2,3],\mathbf{k},t)$ are found to be
suppressed by over one order of magnitude without applying the phase
correction. This can be understood by examining the effect of phase correction
on the complex NAC amplitudes, as contrasted in Fig. S2 (b, c, e, f). While
the real and imaginary part of the time-dependent NAC amplitude is quite
smooth with phase correction, they show artificial random jumps without the
phase correction. Similar random jumps can be observed without gauge
corrections.