# Unexpected linear conductivity in Landau-Zener model: limitations and
improvements of the relaxation time approximation in the quantum master
equation

Ibuki Terada  Department of Physics, Ritsumeikan University, Shiga 525-8577,
Japan  Sota Kitamura  Department of Applied Physics, The University of Tokyo,
Hongo, Tokyo, 113-8656, Japan  Hiroshi Watanabe  Research Organization of
Science and Technology, Ritsumeikan University, Shiga 525-8577, Japan  Hiroaki
Ikeda  Department of Physics, Ritsumeikan University, Shiga 525-8577, Japan
Department of Physics, Ritsumeikan University, Shiga 525-8577, Japan
Department of Applied Physics, The University of Tokyo, Hongo, Tokyo,
113-8656, Japan  Research Organization of Science and Technology, Ritsumeikan
University, Shiga 525-8577, Japan  Department of Physics, Ritsumeikan
University, Shiga 525-8577, Japan

###### Abstract

The nonequilibrium steady states of quantum materials have many challenges.
Here, we highlight issues with the relaxation time approximation (RTA) for the
DC conductivity in insulating systems. The RTA to the quantum master equation
(QME) is frequently employed as a simple method, yet this phenomenological
approach is exposed as a fatal approximation, displaying metallic DC
conductivity in insulating systems within the linear response regime. We find
that the unexpected metallic behavior is caused by the fact that the density
matrix in the RTA incompletely incorporates the first order of the external
field. To solve this problem, we have derived a new calculation scheme based
on the QME that ensure correct behavior in low electric fields. Our method
reproduces well the overall features of the exact electric currents in the
whole field region. It is not time-consuming, and its application to lattice
systems is straightforward. This method will encourage progress in this
research area as a simple way to more accurately describe nonequilibrium
steady states.

††preprint: QME

Introduction— Field-induced phenomena in quantum systems under strong electric
fields have attracted much attention with the recent developments in laser
technology [1, 4, 2, 3, 5, 6]. The electromagnetic responses of condensed-
matter systems beyond the linear response regime lead to a variety of
intriguing phenomena, including bulk photovoltaic effects [7, 8, 9, 10, 11,
12] and nonreciprocal transport of quantum geometric origin [13, 14, 15, 16,
17]. Laser light with stronger intensities triggers photo-induced phase
transitions as a result of nonperturbative quantum effects, which is
experimentally demonstrated in the ultrafast time scale and investigated
theoretically [18, 19, 20, 21, 22, 23, 24, 25, 26].

When the driving field becomes strong, dissipation to environment plays a
vital role in determining the distribution function far from equilibrium,
which is generally difficult and has a long history of research [27, 28, 29,
30, 31, 32, 33, 34, 35]. Theoretical description of open dissipative systems
necessitates the nonequilibrium Green’s function [1, 36, 37, 38, 39] or the
density matrix [40, 41, 42, 43] as a fundamental quantity, rather than the
wave function. The nonequilibrium Green’s function can be calculated using the
diagrammatic approach formulated on the Schwinger-Keldysh contour, while the
density matrix is calculated on the basis of the quantum master equation (QME)
[45, 46]. Recently, nonperturbative modulation of the distribution function
due to the Landau-Zener tunneling [47, 48, 49] is formulated by combining the
nonequilibrium Green’s functions approach with analytic methods [14]. There,
it was discussed that systems with broken inversion symmetry lead to
interesting nonperturbative phenomena such as tunnel spin current and
nonreciprocal current. Such Green’s function method is a powerful tool but
requires the treatment of cumbersome dynamical phases in applications to
lattice systems. The development of QME methods as a complementary method is
therefore considered important.

The relaxation time approximation (RTA) is the simplest approximation to
describe the nonequilibrium steady state in the QME approach [50, 51, 52]. It
has been frequently used in the context of the semi-classical Boltzmann
equation and the semiconductor Bloch equation. Quite recently, however, it has
been pointed out that the RTA has a problem in its application to nonlinear
optics [43, 44]. We show here that, in addition to this problem, the RTA
treatment in insulating systems involves a fatal failure that the DC
conductivity exhibits metallic behavior in the linear response regime, despite
the absence of Fermi surfaces. The reason why this has been overlooked until
now is that transport properties are mainly studied in the perturbative regime
(i.e., the DC conductivity of metals or the AC conductivity of insulators),
whereas the tunneling current is a purely nonperturbative response seen only
under strong fields.

In this letter, we demonstrate the failure of the RTA using the Landau-Zener
model as a minimal example. We find that the fatal metallic behavior
originates from the fact that the density matrix in the RTA incorporates the
first order of the external field $E$ in an incomplete form. On the basis of
the QME, we derive a calculation scheme to ensure correct behavior in low
electric fields by sequentially incorporating the perturbation correction of
the external field. Our method reproduces well the overall features of
electric currents in the whole field region. It is not time-consuming, and its
computational time is comparable to RTA calculations.

Relaxation Time Approximation— To demonstrate a fatal failure of the RTA, let
us consider a two-band system with a finite gap for simplicity. With two
eigenenergies $\varepsilon_{k\pm}$, the band gap at each $k$ point is given by
$\Delta_{k}=\varepsilon_{k+}-\varepsilon_{k-}>0$. We introduce the DC electric
field $E$ via the Peierls substitution, $k\to k-eEt$ (hereafter $e=1$).
Following Ref. [14], we introduce the snapshot basis
$|\Phi_{k\alpha}(t)\rangle=|u_{\alpha,k-Et}\rangle e^{-i\Theta_{\alpha}(t)}$
($\alpha=\pm$) with $|u_{\alpha,k}\rangle$ being the eigenstate of the
snapshot Hamiltonian $H(k)$. The phase factor $\Theta_{\alpha}(t)$ consists of
the dynamical and Berry phases [53]. In this basis, the QME for the density
matrix $[\rho_{k}(t)]_{\alpha\beta}$ is written in the $2\times 2$ matrix form
as

$\displaystyle\frac{d\rho_{k}(t)}{dt}=-i\big{[}\mathcal{W}_{k}(t),\rho_{k}(t)\big{]}+\mathcal{D}\big{[}\rho_{k}(t)\big{]},$
(1)

where $[\mathcal{W}_{k}(t)]_{\alpha\alpha}=0$,
$[\mathcal{W}_{k}(t)]_{+-}=[\mathcal{W}_{k}(t)]^{*}_{-+}\equiv
W_{k}(t)=E\langle u_{+,k-Et}|i\partial_{k}|u_{-,k-Et}\rangle
e^{i\Theta_{+}(t)-i\Theta_{-}(t)}$ denotes the transition dipole matrix
elements. Here the dissipation to environment is described by
$\mathcal{D}\big{[}\rho_{k}(t)\big{]}$, whose form in the RTA is given by

$\displaystyle\big{[}\mathcal{D}_{k}[\rho_{k}(t)]\big{]}_{\alpha\beta}=\frac{f_{D}(\varepsilon_{k\alpha}(t))\delta_{\alpha\beta}-\big{[}\rho_{k}(t)\big{]}_{\alpha\beta}}{\tau_{\alpha\beta}},$
(2)

where $f_{D}$ is the Fermi-Dirac distribution function, and
$\varepsilon_{k\alpha}(t)=\varepsilon_{k-Et,\alpha}$.
$\tau_{++}=\tau_{--}=\tau_{1}$ and $\tau_{+-}=\tau_{-+}=\tau_{2}$ denote the
longitudinal and transverse relaxation time, respectively.

By solving Eq. (1), we obtain the density matrix of the nonequilibrium steady
state. Then the electric current $J(t)$ is calculated as
$J(t)=-e\int\frac{dk}{2\pi}\text{Tr}[\partial_{k}H(k)\rho_{k}(t)]$, which can
be decomposed into intra- and interband contributions as

$\displaystyle J(t)=J_{\rm{intra}}(t)+J_{\rm{inter}}(t),$ (3) $\displaystyle
J_{\rm{intra}}(t)=-\sum_{\alpha=\pm}\int\frac{dk}{2\pi}\frac{\partial\varepsilon_{k\alpha}(t)}{\partial
k}\big{[}\rho_{k}(t)\big{]}_{\alpha\alpha},$ (4) $\displaystyle
J_{\rm{inter}}(t)=2{\rm
Im}\left[\int\frac{dk}{2\pi}\frac{W_{k}(t)\Delta_{k}(t)}{E}\big{[}\rho_{k}(t)\big{]}_{-+}\right].$
(5)

The electric current in the nonequilibrium steady states is obtained as the
long-time limit of $J(t)$.

Let us exemplify the electric current in the RTA by a numerical calculation
for the Landau-Zener model

$H(k)=\left(\begin{matrix}vk&\delta\\\ \delta&-vk\end{matrix}\right),$ (6)

where $2\delta$ corresponds to the band gap. In the absence of the
dissipation, the tunneling probability in the $t\to\infty$ limit is exactly
given by $P_{\rm{LZ}}=\exp({-\pi E_{\rm th}/E})$ with $E_{\rm
th}=\delta^{2}/v$ [the dashed line of the main of Fig. 1(a)] [47, 54]. Figure
1(a) shows the field dependence of the RTA electric current $J^{\rm RTA}$ with
$\tau_{1}=\tau_{2}=\tau$ at the temperature $T=0.01\delta$. One can see a
remarkable increase of the current $J^{\rm RTA}$ at around $E=0.5E_{\rm th}$,
consistent with the generation of tunneling carriers. Surprisingly, an
unexpected linear $E$ dependence with the slope increasing with $\tau^{-1}$ is
observed in the low $E$ regime, indicating metallic behavior even though the
system is actually an insulator. What is responsible for this finite linear
dependence? To clarify the origin, in Fig. 1(b), we separately show the
contributions from $J^{\rm RTA}_{\rm intra}$ (dashed line) and $J^{\rm
RTA}_{\rm inter}$ (dotted line) for $\tau^{-1}=0.02\delta$. One can see that
the $E$ linear behavior comes from the interband contribution $J^{\rm
RTA}_{\rm inter}$.

Figure 1: (a) Field dependence of the DC current $J^{\rm RTA}$ in the RTA at
$T=0.01\delta$. The inset depicts $J^{\rm RTA}$ in a wide range. Solid lines
denote the DPA currents for several $\tau^{-1}$. The dotted line represents
the tunneling probability $P_{\rm{LZ}}$. One can see the exponential behavior
in the inset, but encounter an unexpected linear behavior at a low $E$ limit
in the main panel. (b) Intraband $J^{\rm RTA}_{\rm intra}$ (dashed line) and
interband $J^{\rm RTA}_{\rm inter}$ (dotted line) contributions of $J^{\rm
RTA}$ at $\tau^{-1}=0.02\delta$ as a function of $E/E_{\rm th}$. $J^{\rm
RTA}_{\rm inter}$ shows unphysical linear $E$ behavior.

In fact, the density matrix in the low $E$ limit can be calculated as

$\displaystyle\left[\rho_{k}(t)\right]_{\alpha\alpha}\sim
f_{D}(\varepsilon_{k\alpha}(t))+E\tau_{1}\frac{\partial
f_{D}(\varepsilon_{k\alpha}(t))}{\partial k},$ (7)
$\displaystyle\left[\rho_{k}(t)\right]_{+-}\sim-\frac{W_{k}(t)}{\Delta_{k}(t)-i\tau_{2}^{-1}}\delta
f_{k}(t)$ (8)

with $\delta f_{k}(t)=f_{D}(\varepsilon_{k-}(t))-f_{D}(\varepsilon_{k+}(t))$,
by solving Eq. (1) in a perturbative manner. Then, within the linear response,
the electric conductivities $\displaystyle\sigma^{\rm RTA}=\lim_{E\to 0}J^{\rm
RTA}/E$ are given by

$\displaystyle\sigma^{\rm RTA}_{\rm intra}$
$\displaystyle=\tau_{1}\sum_{\alpha=\pm}\int\frac{dk}{2\pi}\left(\frac{\partial\varepsilon_{k\alpha}}{\partial
k}\right)^{2}\left(-\frac{\partial
f_{D}}{\partial\varepsilon_{k\alpha}}\right),$ (9) $\displaystyle\sigma^{\rm
RTA}_{\rm inter}$
$\displaystyle=2\tau_{2}^{-1}\int\frac{dk}{2\pi}\frac{|\langle
u_{+,k}|i\partial_{k}|u_{-,k}\rangle|^{2}\Delta_{k}}{\tau_{2}^{-2}+\Delta_{k}^{2}}\delta
f_{k}.$ (10)

In the insulating case, the intraband conductivity $\sigma^{\rm RTA}_{\rm
intra}$ vanishes at $T=0$, which is consistent with the Boltzmann theory. On
the other hand, the interband conductivity $\sigma^{\rm RTA}_{\rm inter}$ does
not vanish even in the insulating case, which gives the slope of $J^{\rm
RTA}_{\rm inter}$ at $E\to 0$. This term increases roughly in proportion to
damping $\tau_{2}^{-1}$. The small but finite linear conductivity, which
exists despite the absence of Fermi surfaces, is a failure of the RTA in the
QME. Although the RTA is useful to describe the nonequilibrium steady state in
a simple manner, it has been recently argued that this approximation is also
problematic for the optical response [43, 44], which is another problem
different from the unphysical behavior observed here for the DC current.

As we will show later, this problem comes from the presence of
$i\tau_{2}^{-1}$ in the denominator of Eq. (8). This is due to the fact that
the RTA does not properly incorporate the first-order contribution of the
electric field $E$. In our formalism beyond the RTA, this term in the
denominator vanishes [cf. Eq. (19)]. This $i\tau_{2}^{-1}$ term also affects
the Hall current, with which the Hall conductivity does not quantize within
the RTA [50].

Beyond RTA— To resolve the failure of the RTA, here we consider to couple the
two-band insulator to a fermionic reservoir [1, 55] within the QME formalism,
and derive the dissipation term $\mathcal{D}[\rho_{k}(t)]$ microscopically. We
start with the Born-Markov master equation [45, 56]

$\frac{d\tilde{\rho}_{k}(t)}{dt}=-\int_{-\infty}^{t}\mathrm{Tr_{B}}\left[\tilde{H}_{I,k}(t),\left[\tilde{H}_{I,k}(s),\tilde{\rho}_{k}(t)\otimes\tilde{\rho}_{B}\right]\right]ds,$
(11)

where $\tilde{\rho}_{k}(t)$ is the reduced density operator and
$\tilde{\rho}_{B}$ is the thermal density operator of the bath, respectively.
Here the tilde on operators denote the interaction picture.
$\tilde{H}_{I,k}(t)$ represents the interaction term between the system and
the bath, $\tilde{H}_{I,k}(t)=\sum_{\sigma
p}V_{p}\tilde{b}^{{\dagger}}_{k\sigma
p}(t)\tilde{c}_{k\sigma}(t)+\text{H.c.}$, where $\tilde{c}_{k\sigma}(t)$ and
$\tilde{b}_{k\sigma p}(t)=\tilde{b}_{k\sigma p}(0)e^{-i\omega_{p}t}$ are
respectively the annihilation operator of an electron in the insulator and
reservoir, with momentum $k$ and pseudospin $\sigma$.
$\mathrm{Tr_{B}}\left[\cdots\right]$ means tracing out the bath degrees of
freedom. We impose the broadband condition for the spectral density of the
fermionic reservoir as

$\displaystyle\sum_{p}\pi|V_{p}|^{2}\,\delta(\omega-\omega_{p})=\Gamma~{}~{}({\rm
const.}).$ (12)

The key point of our formalism is to express Eq. (11) in the snapshot basis by
introducing the transformed field operator

$\tilde{\psi}_{k\alpha}(t)=\sum_{\sigma}\innerproduct{\Phi_{k\alpha}(t)}{\sigma}\tilde{c}_{k\sigma}(t),$
(13)

and evaluate the integral in Eq. (11) with employing the adiabatic
perturbation theory [13, 14]. Specifically, in the adiabatic limit $E\to 0$,
we can incorporate changes in the dynamical phase as

$\displaystyle W_{k}(t-s)$ $\displaystyle\sim e^{-i\Delta_{k}(t)s}W_{k}(t),$
(14) $\displaystyle|\Phi_{k\alpha}(t-s)\rangle$ $\displaystyle\sim
e^{i\varepsilon_{k\alpha}(t)s}|\Phi_{k\alpha}(t)\rangle.$ (15)

We further include the change of the parameter $\Delta k=-Es$ in a
perturbative manner, such that the resultant QME is exact up to $O(E)$ while
maintaining the unitarity of the time evolution operator up to $O(E^{2})$. We
call this approximation the dynamical phase approximation (DPA).

After all, we obtain the dissipation term in the DPA as ${\cal
D}[\rho_{k}(t)]={\cal D}_{0}+{\cal D}_{1}+{\cal D}_{2}$, where

$\displaystyle\big{[}{\cal
D}_{0}\big{]}_{\alpha\beta}=-2\Gamma\big{(}[\rho_{k}(t)]_{\alpha\beta}-f_{D}(\varepsilon_{k\alpha}(t))\delta_{\alpha\beta}\big{)},$
(16) $\displaystyle\big{[}{\cal
D}_{1}\big{]}_{\alpha\beta}=-2\Gamma\frac{[{\cal
W}_{k}(t)]_{\alpha\beta}}{\Delta_{k}(t)}\delta f_{k}(t),$ (17)
$\displaystyle\big{[}{\cal
D}_{2}\big{]}_{\alpha\beta}=2\Gamma\alpha\frac{|W_{k}(t)|^{2}}{\Delta_{k}^{2}(t)}\bigg{(}\delta
f_{k}(t)+\Delta_{k}(t)f_{D}^{\prime}(\varepsilon_{k,-\alpha}(t))\bigg{)}\delta_{\alpha\beta}.$
(18)

The subscript $n$ of ${\cal D}_{n}$ denotes the perturbation order of $W_{k}$
[59]. Note that the expression here is simplified one under the particle-hole
symmetry, $\varepsilon_{k+}(t)=-\varepsilon_{k-}(t)$. Here the zeroth order
term ${\cal D}_{0}$ corresponds to the RTA with $\tau_{1}=\tau_{2}=1/2\Gamma$,
while ${\cal D}_{1},{\cal D}_{2}$ describe the field-induced correction terms.
Specifically, in the presence of Eq. (17), the low-$E$ expression for
$[\rho_{k}(t)]_{+-}$ is replaced from Eq. (8) into

$[\rho_{k}(t)]_{+-}\sim-\dfrac{W_{k}(t)}{\Delta_{k}(t)}\delta f_{k}(t),$ (19)

with which unphysical linear term in $J^{\rm RTA}_{\rm inter}$ [Eq. (5)] is
completely canceled out.

Figure 2: Occupation number of a single electron in the upper band. Solid
lines and filled circles denote DPA results and the numerically exact results
[14], respectively. The schematic figure represents the electron dynamics in
the Landau-Zener model. The DPA results are almost consistent with the exact
results. The small deviation in due to the damping-induced excitation inherent
in open systems.

Now, let us compare the DPA calculation with the numerically exact calculation
[14]. Figures 2(a) and (b) show the occupation numbers of the upper band at
$\tau^{-1}=2\Gamma=0.02\delta$ and $T=0.01\delta$ as a function of $k(t)=-Et$.
Electrons starting from $t=t_{0}<0~{}(k>0)$ are excited by Landau-Zener
tunneling as they pass through the gap minimum at $t=0~{}(k=0)$. Then, at
$t>0~{}(k<0)$, the excited electrons decay with a relaxation time $\tau$. The
DPA results are consistent with the exact result at higher temperatures and
stronger electric fields. The small deviation in the low-$E$ regime is due to
the damping-induced excitation inherent to the present fermionic reservoir,
which is only partially included in the DPA. While this is incorrect from the
viewpoint of the rigorous treatment of the fermionic reservoir, it could be
considered as an advantage since this particular excitation should be absent
for ideal environment.

Figure 3: (a) Field dependence of the DC current $J^{\rm DPA}$ in the DPA. The
inset depicts $J^{\rm RTA}$ in the wide range. Solid lines denote the DPA
currents for several $\tau^{-1}$, the dashed line represents the RTA current
at $\tau^{-1}=0.02\delta$ and the filled circles the result of the numerically
exact result. In the inset, one can see that the DPA result reproduces the
exact result well. $J^{\rm DPA}$ almost vanishes at $E<0.2E_{\rm th}$. We can
verify in (b) that $J^{\rm DPA}$ has no linear $E$ terms with an accuracy of
less than $10^{-7}\delta$ at $E<0.06E_{\rm th}$. At intermediate region $E\sim
0.4E_{\rm th}$, $J^{\rm DPA}$ deviates from the exact result by $\sim
10^{-3}\delta$ but roughly increases with the tunneling probability $P$,
implying that the excited carriers are carrying the electric current.

Finally, we illustrate the field dependence of the electric current in the DPA
in Fig. 3. Solid lines in Fig. 3(a) denote the DPA currents for several
damping parameters, the dotted line represents the RTA current at
$\tau^{-1}=0.02\delta$, and the filled circles are the numerically exact
results. It is clear that the RTA fails to capture qualitative features of the
exact result. On the other hand, as shown in the inset of Fig. 3(a), the DPA
result well describes the exact one. The most significant improvement is the
disappearance of the linear $E$ dependence in the interband contribution. This
is due to the fact that the $i\tau^{-1}$ term in the off-diagonal term of the
density matrix in the RTA is completely canceled by correctly treating the
electric field up to the first order. This also improves the behavior of the
intraband contribution, where the exponential behavior is more pronounced. As
a result, the total current in the low-field region is also greatly improved
and shows exponential behavior. Indeed, $J^{\rm DPA}$ almost vanishes at
$E<0.2E_{\rm th}$. It can be verified in Fig. 3(b) that $J^{\rm DPA}$ has no
linear $E$ terms with an accuracy of less than $10^{-7}\delta$ at
$E<0.06E_{\rm th}$. At intermediate region $E\sim 0.4E_{\rm th}$, $J^{\rm
DPA}$ has an error of $\sim 10^{-3}\delta$. This is thought to be due to the
damping-induced excitation and/or the fact that the DPA calculation was
terminated up to $\mathcal{D}_{2}$. Although further refinements to
incorporate these effects are available in principle, the overall features of
the DC current are described well enough qualitatively and semi-quantitatively
in our DPA. In addition, it has a great advantage that the computational time
required to obtain these results is almost the same as for the RTA.

Conclusion— This letter highlights issues with the RTA for the DC current in
insulating systems and proposes an improvement based on the QME. The RTA is
frequently employed as a simple method, yet this phenomenological approach is
exposed as a fatal approximation, displaying metallic DC conductivity in
insulating systems within the linear response regime. This is because of the
incomplete inclusion of the first-order contribution of the electric field in
the density matrix. We reevaluate the QME and incorporate the dynamical phase
of the transition matrix, thereby correctly capturing the first-order terms of
the electric field. We obtained a new scheme that accurately describes the
insulating behavior. It was demonstrated that it accurately predicts the
correct off-diagonal terms in the density matrix and is in semi-quantitative
agreement with the numerically exact result. This also improves the
calculation of quantum Hall effects. While numerically exact calculation can
be performed for the present fundamental model, it is often time-consuming.
Our method is not rigorous but correctly captures the overall features in the
DC current. Furthermore, it does not require much computation time and is
straightforward to apply to lattice systems. The DPA is an alternative to the
RTA, which describes nonequilibrium steady states more correctly. We believe
that this method will encourage progress in this research field.

###### Acknowledgements.

We are grateful to Y. Michishita, K. Takasan, A. Oguri, M. Sato, N. Kawakami,
T. Morimoto and T. Oka for useful comments. This work was supported by KAKENHI
Grants No. 19H01842, No. 19H05825. SK is supported by JSPS KAKENHI Grant No.
20K14407.

## References

  * [1] H. Aoki, N. Tsuji, M. Eckstein, M. Kollar, T. Oka, and P. Werner, Rev. Mod. Phys. 86, 779 (2014).
  * [2] T. Oka and S. Kitamura: Annu. Rev. Condens. Matter Phys. 10 387 (2019).
  * [3] André Eckardt, Rev. Mod. Phys. 89, 011004 (2017).
  * [4] D. N. Basov, R. D. Averitt,and D. Hsieh, Nat. Mater. 16, 1077–1088 (2017).
  * [5] A. de. la. Torre, D. M. Kennes, M. Claassen, S. Gerber, J. W. McIver, and M. A. Sentef, Rev. Mod. Phys. 93 041002 (2021).
  * [6] T. Morimoto, S. Kitamura, and N. Nagaosa, J. Phys. Soc. Jpn. 92, 072001 (2023).
  * [7] R. W. Boyd, Nonlinearoptics (Academic press, London, 2003).
  * [8] N. Bloembergen: Nonlinear optics (World Scientific, Singapore, 1996).
  * [9] P. J. Sturman and V. M. Fridkin: Photovoltaic and Photo-refractive Effects in Noncentrosymmetric Materials (CRC Press, Philadelphia, 1992), Vol. 8.
  * [10] W. Nie, H. Tsai, R. Asadpour, J. -C. Blancon, A. J. Neukirch, G. Gupta, J. J. Crochet, M. Chhowalla, S. Tretiak, M. A. Alam, H. -L. Wang, and A. D. Mohite: Science 347 522 (2015).
  * [11] D. Shi, V. Adinolfi, R. Comin, M. Yuan, E. Alarousu, A. Buin, Y. Chen, S. Hoogland, A. Rothenberger, K. Katsiev, Y. Losovyj, X. Zhang, P. A. Dowben, O. F. Mohammed, E. H. Sargent, and O. M. Bakr, Science 347 519 (2015).
  * [12] D. W. de Quilettes, S. M. Vorpahl, S. D. Stranks, H. Nagaoka, G. E. Eperon, M. E. Ziffer, H. J. Snaith, and D. S. Ginger: Science 348 683 (2015).
  * [13] S. Kitamura, N. Nagaosa, and T. Morimoto: Commun. Phys. 3 63 (2020).
  * [14] S. Kitamura, N. Nagaosa, and T. Morimoto: Phys. Rev. B 102 245141 (2020).
  * [15] S. Takayoshi, J. Wu, T. Oka, SciPost Phys. 11, 075 (2021)
  * [16] Y. Suzuki, Phys. Rev. B 105, 075201 (2022).
  * [17] Inti Sodemann and Liang Fu, Phys. Rev. Lett. 115, 216806 (2015).
  * [18] K. Nasu, Photoinduced phase transitions (World Scientific, 2004).
  * [19] D. N. Basov, R. D. Averitt, D. van der Marel, M. Dressel, and K. Haule, Rev. Mod. Phys. 83, 471 (2011).
  * [20] T. Oka and H. Aoki: Phys. Rev. B 79 081406 (2009).
  * [21] T. Kitagawa, T. Oka, A. Brataas, L. Fu, and E. Demler Phys. Rev. B 84, 235108 (2011).
  * [22] J. W. Mclver, B. Schulte, F.-U. Stein, T. Matsuyama, G. Jotzu, G. Meier and A. Cavalleri, Nat. Phys 16 38–41 (2020).
  * [23] A. Kirilyuk, A. V. Kimel, and T. Rasing, Rev. Mod. Phys. 82, 2731 (2010).
  * [24] S. Takayoshi, H. Aoki, and T. Oka, Phys. Rev. B 90, 085150 (2014).
  * [25] S. Takayoshi, M. Sato, and T. Oka, Phys. Rev. B 90, 214413 (2014).
  * [26] M. Sato, S. Takayoshi, and T. Oka, Phys. Rev. Lett. 117, 147202 (2016).
  * [27] U. Weiss, Quantum Dissipative Systems (World Scientific, 1999).
  * [28] A. Barreiro, M. Lazzeri, J. Moser, F. Mauri, and A. Bachtold, Phys. Rev. Lett. 103, 076601 (2009).
  * [29] N. Vandecasteele, A. Barreiro, M. Lazzeri, A. Bachtold, and F. Mauri, Phys. Rev. B 82, 045416 (2010).
  * [30] T. Fang, A. Konar, H. Xing, and D. Jena, Phys. Rev. B 84, 125450 (2011).
  * [31] J. Li and J. E. Han, Phys. Rev. B 97, 205412 (2018).
  * [32] S. Okamoto, Phys. Rev. B 76, 035105 (2007).
  * [33] N. Sugimoto, S. Onoda, and N. Nagaosa, Phys. Rev. B 78, 155104 (2008).
  * [34] F. Heidrich-Meisner, I. González, K. A. Al-Hassanieh, A. E. Feiguin, M. J. Rozenberg, and E. Dagotto, Phys. Rev. B 82, 205110 (2010).
  * [35] T. N. Ikeda, K. Chinzei, M. Sato, SciPost Phys. Core 4, 033 (2021).
  * [36] A. Blandin, A. Nourtier, and D. W. Hone, J. Phys. (Paris) 37, 369 (1976).
  * [37] A. P. Jauho, N. S. Wingreen, and Y. Meir, Phys. Rev. B 50, 5528 (1994).
  * [38] S. Onoda, N. Sugimoto and N. Nagaosa, Progress of Theoretical Physics 116, 61 (2006).
  * [39] J. E. Han,Phys. Rev. B 87, 085119 (2013)
  * [40] J. E. Sipe and Ed Ghahramani, Phys. Rev. B 48, 11705 (1993).
  * [41] J. E. Sipe and A. I. Shkrebtii, Phys. Rev. B 61, 5337 (2000).
  * [42] G. B. Ventura, D. J. Passos, J. M. B. Lopes dos Santos, J. M. Viana Parente Lopes, and N. M. R. Peres, Phys. Rev. B 96, 035431 (2017).
  * [43] D. J. Passos, G. B. Ventura, J. M. Viana Parente Lopes, J. M. B. Lopes dos Santos, and N. M. R. Peres, Phys. Rev. B 97, 235446 (2018).
  * [44] Y. Michishita and R. Peters: Phys. Rev. B 103 195133 (2021).
  * [45] H. -P. Breuer and F. Petruccione, The Theory of Open Quantum Systems (Oxford University Press, Oxford, 2007).
  * [46] R. Alicki, and K. Lendi, Quantum Dynamical Semigroups and Applications (Springer, 2007).
  * [47] C. Zener: Proc. R. Soc. Lond. A 137 696 (1932).
  * [48] J. P. Daivis and P. Pechukas, J. Chem. Phys. 64 3129 (1976).
  * [49] A. Dykhne, J. Exptl. Theoret. Phys. (USSR) 41, 1324 (1961) [Sov. Phys.JETP 14, 1 (1962)].
  * [50] S. A. Sato, P. Tang, M. A. Sentef, U. De. Giovannini, H. Hübener, and A. Rubio, New J. Phys. 21 093005 (2019).
  * [51] M. Nuske, L. Broers, B. Schulte, G. Jotzu, S. A. Sato, A. Cavalleri, A. Rubio, J. W. McIver, and L. Mathey Phys. Rev. Research 2, 043408 (2020).
  * [52] S. A Sato and A. Rubio, New J. Phys. 23 063047 (2021).
  * [53] The explicit form reads $\Theta_{\alpha}(t)=\int_{t_{0}}^{t}dt^{\prime}\,[\varepsilon_{k\alpha}(t^{\prime})+E\matrixelement{u_{\alpha,k-Et^{\prime}}}{i\partial_{k}}{u_{\alpha,k-Et^{\prime}}}]$.
  * [54] C. D. Grandi and A. Polkovnikov, Adiabatic perturbation theory: From Landau-Zener problem to quenching through a quantum critical point, Quantum Quenching, Annealing and Computation (Springer, Berlin, 2010), pp. 75–114.
  * [55] M. Büttiker,Phys. Rev. Lett. 57, 1761 (1986).
  * [56] L. Del Re, B. Rost, A. F. Kemper, and J. K. Freer-icks, Phys. Rev. B 102, 125112 (2020).
  * [57] Y. Michishita and N. Nagaosa: Phys. Rev. B 106 125114 (2022).
  * [58] M. Yamaguchi, T. Yuge, and T. Ogawa Phys. Rev. E 95, 012136 (2017).
  * [59] See Supplemental Material for the derivation of the dissipation term $\mathcal{D}_{n}(t)$.

Supplemental Material: Unexpected linear conductivity in Landau-Zener model:
limitations and improvements of the relaxation time approximation in the
quantum master equation

Ibuki Terada Sota Kitamura Hiroshi Watanabe Hiroaki Ikeda

##  Appendix A Supplemental Material

##  Appendix B Time evolution operator of isolated two-band systems

We here introduce the snapshot basis [13, 14] as a solution of the time-
dependent Schrödinger equation in the adiabatic limit (without the fermionic
reservoir) [48, 49], in order to systematically expand the quantum master
equation (QME) with respect to the external field (known as adiabatic
perturbation theory).

According to the adiabatic theorem, the solution of the time-dependent
Schrödinger equation in the adiabatic limit can be written as the
instantaneous eigenstate of the system’s Hamiltonian $H(k-Et)$ multiplied by a
phase factor,

$\ket{\Phi_{k\alpha}(t)}=\ket{u_{k\alpha}(t)}e^{-i\Theta_{\alpha}(t)}.$ (S1)

Here we define the instantaneous eigenstate
$\ket{u_{k\alpha}(t)}=\ket{u_{\alpha,k-Et}}$ with
$H(k)|u_{\alpha,k}\rangle=\varepsilon_{k\alpha}|u_{\alpha,k}\rangle$, and the
phase factor $\Theta_{\alpha}(t)$ given by

$\displaystyle\Theta_{\alpha}(t)$
$\displaystyle=\int_{t_{0}}^{t}dt^{\prime}\,[\varepsilon_{k\alpha}(t^{\prime})+EA_{\alpha\alpha}(t^{\prime})]$
(S2)
$\displaystyle=\int_{t_{0}}^{t}dt^{\prime}\,[\varepsilon_{k\alpha}(t^{\prime})+E\matrixelement{u_{k\alpha}(t^{\prime})}{i\partial_{k}}{u_{k\alpha}(t^{\prime})}],$
(S3)

which is introduced such that the diagonal element of the time-dependent
Schrödinger equation vanishes. The label $\alpha=\pm$ specifies the band
index.

We expand the one-particle wave function $\ket{\Psi_{k}(t)}$ by the snapshot
basis introduced above, as

$\displaystyle\ket{\Psi_{k}(t)}=\sum_{\alpha=\pm}a_{k\alpha}(t)\ket{\Phi_{k\alpha}(t)}.$
(S4)

Then the time-dependent Schrödinger equation
$i\partial_{t}\ket{\Psi_{k}(t)}=H(k-Et)\ket{\Psi_{k}(t)}$ can be rewritten as

$\displaystyle i\frac{d}{dt}\left(\begin{matrix}a_{k+}(t)\\\
a_{k-}(t)\end{matrix}\right)$
$\displaystyle=\mathcal{W}_{k}(t)\left(\begin{matrix}a_{k+}(t)\\\
a_{k-}(t)\end{matrix}\right)$ (S5)
$\displaystyle=\left(\begin{matrix}0&W_{k}(t)\\\
W_{k}^{\ast}(t)&0\end{matrix}\right)\left(\begin{matrix}a_{k+}(t)\\\
a_{k-}(t)\end{matrix}\right)$ (S6)
$\displaystyle=\left(\begin{matrix}0&E\matrixelement{u_{k+}(t)}{i\partial_{k}}{u_{k-}(t)}e^{i\Theta_{+}(t)-i\Theta_{-}(t)}\\\
E\matrixelement{u_{k-}(t)}{i\partial_{k}}{u_{k+}(t)}e^{i\Theta_{-}(t)-i\Theta_{+}(t)}&0\end{matrix}\right)\left(\begin{matrix}a_{k+}(t)\\\
a_{k-}(t)\end{matrix}\right).$ (S7)

The formal solution of Eq. (S5) is
$\left(a_{k+}(t),a_{k-}(t)\right)^{T}=U(t,t_{0})\left(a_{k+}(t_{0}),a_{k-}(t_{0})\right)^{T}$,
where we introduce the time evolution operator (of the isolated system) in the
snapshot basis,

$U(t,t_{0})=\mathcal{T}\exp\left[-i\int_{t_{0}}^{t}\mathcal{W}_{k}(t^{\prime})dt^{\prime}\right],$
(S8)

with the time-ordered product operator $\mathcal{T}$. One can show that
$U(t,t_{0})\to 1$ for arbitrary $t$, when the system is gapped and the
external field $E$ is infinitesimally small.

Here we also introduce the snapshot basis in the second-quantized form. We
define the second-quantized Hamiltonian as
$\tilde{H}(k,t)=\sum_{\sigma\sigma^{\prime}}\tilde{c}_{k\sigma}^{\dagger}(t)[H(k-Et)]_{\sigma\sigma^{\prime}}\tilde{c}_{k\sigma^{\prime}}(t)$
with the field operator $\tilde{c}_{k\sigma}(t)$ in the Heisenberg
representation. Then the snapshot basis for the field operator is defined by
[14]

$\tilde{\psi}_{k\alpha}(t)=\sum_{\sigma}\innerproduct{\Phi_{k\alpha}(t)}{\sigma}\tilde{c}_{k\sigma}(t).$
(S9)

The equation of motion for the snapshot field operator is obtained using
$i\partial_{t}\tilde{c}_{k\sigma}(t)=[\tilde{c}_{k\sigma}(t),\tilde{H}(k,t)]$
as

$i\partial_{t}\tilde{\psi}_{k\alpha}(t)=\sum_{\beta}[\mathcal{W}_{k}(t)]_{\alpha\beta}\tilde{\psi}_{k\beta}(t),$
(S10)

i.e., the first-quantized coefficients $a_{k\alpha}(t)$ and the second-
quantized operators $\tilde{\psi}_{k\alpha}(t)$ obey the same equation of
motion. It immediately follows that the snapshot field operator satisfies
$\left(\tilde{\psi}_{k+}(t),\tilde{\psi}_{k-}(t)\right)^{T}=U(t,t_{0})\left(\tilde{\psi}_{k+}(t_{0}),\tilde{\psi}_{k-}(t_{0})\right)^{T}$.

##  Appendix C Density matrix within the RTA in the low $E$ regime

Here we calculate the density matrix within the RTA at the low electric field.
Diagonal and offdiagonal elements of Eq. (1) with Eq. (2) read

$\displaystyle\frac{d[\rho_{k}(t)]_{\alpha\alpha}}{dt}$ $\displaystyle=-\alpha
2\text{Im}\left[W^{\ast}_{k}(t)[\rho_{k}(t)]_{+-}\right]-\frac{[\rho_{k}(t)]_{\alpha\alpha}-f_{D}(\varepsilon_{k\alpha}(t))}{\tau_{1}},$
(S11) $\displaystyle\frac{d[\rho_{k}(t)]_{+-}}{dt}$
$\displaystyle=-iW_{k}(t)\left[[\rho_{k}(t)]_{--}-[\rho_{k}(t)]_{++}\right]-\frac{[\rho_{k}(t)]_{+-}}{\tau_{2}}.$
(S12)

First (second) equation can be formally integrated when the offdiagonal
(diagonal) matrix element is regarded as a source term. We obtain

$\displaystyle[\rho_{k}(t)]_{\alpha\alpha}$
$\displaystyle=\int_{-\infty}^{t}\left\\{\frac{f_{D}(\varepsilon_{k\alpha}(s))}{\tau_{1}}-\alpha
2\text{Im}\left[W^{\ast}_{k}(s)[\rho_{k}(s)]_{+-}\right]\right\\}e^{-(t-s)/\tau_{1}}ds,$
(S13) $\displaystyle[\rho_{k}(t)]_{+-}$
$\displaystyle=-i\int_{-\infty}^{t}W_{k}(s)\left[[\rho_{k}(s)]_{--}-[\rho_{k}(s)]_{++}\right]e^{-(t-s)/\tau_{2}}ds,$
(S14)

where the dependence on the initial condition vanishes as we take the initial
time as $t_{0}\rightarrow-\infty$. We can expand the first term in Eq. (S13)
as

$\int_{-\infty}^{t}\frac{f_{D}(\varepsilon_{k\alpha}(s))}{\tau_{1}}e^{-(t-s)/\tau_{1}}ds=f_{D}(\varepsilon_{k\alpha}(t))+E\tau_{1}\frac{\partial
f_{D}(\varepsilon_{k\alpha}(t))}{\partial k}+\mathcal{O}(E^{2}),$ (S15)

by performing the integral by parts. On the other hand, the second term in Eq.
(S13) has an $\mathcal{O}(E^{2})$ contribution because
$[\rho_{k}(s)]_{+-}=O(E)$ as we see below. Therefore, the diagonal part of the
density matrix is written as

$[\rho_{k}(t)]_{\alpha\alpha}\simeq
f_{D}(\varepsilon_{k\alpha}(t))+E\tau_{1}\frac{\partial
f_{D}(\varepsilon_{k\alpha}(t))}{\partial k}$ (S16)

up to the first order of $E$, which corresponds to the Boltzmann transport
theory. In a similar manner, we can expand Eq. (S14) with respect to $E$,
using the relation

$\displaystyle e^{i\Theta_{+}(s)-i\Theta_{-}(s)+s/\tau_{2}}$
$\displaystyle=e^{i\int_{t_{0}}^{s}dt^{\prime}\,[\Delta_{k}(t^{\prime})+E(A_{++}(t^{\prime})-A_{--}(t^{\prime}))]+s/\tau_{2}}$
(S17)
$\displaystyle=e^{i\int_{t_{0}}^{s}dt^{\prime}\,E(A_{++}(t^{\prime})-A_{--}(t^{\prime}))}\frac{\frac{d}{ds}e^{i\int_{t_{0}}^{s}dt^{\prime}\,\Delta_{k}(t^{\prime})+s/\tau_{2}}}{i\Delta_{k}(s)+\tau_{2}^{-1}}$
(S18)

and performing the integral by parts. The result reads

$[\rho_{k}(t)]_{+-}\simeq-\dfrac{W_{k}(t)}{\Delta_{k}(t)-i\tau_{2}^{-1}}\left[f_{D}(\varepsilon_{k-}(t))-f_{D}(\varepsilon_{k+}(t))\right]$
(S19)

at the first order of $E$ [Note that $W_{k}(t)=\mathcal{O}(E)$].

On the basis of above results (S16) and (S19), we obtain the linear
conductivity $\sigma^{\rm{RTA}}$,

$\displaystyle\sigma^{\rm{RTA}}$
$\displaystyle=\sigma^{\rm{RTA}}_{\rm{intra}}+\sigma^{\rm{RTA}}_{\rm{inter}},$
(S20) $\displaystyle\sigma^{\rm{RTA}}_{\rm{intra}}$
$\displaystyle=\tau_{1}\sum_{\alpha=\pm}\int\frac{dk}{2\pi}\left(\frac{\partial\varepsilon_{k\alpha}}{\partial
k}\right)^{2}\left(-\frac{\partial
f_{D}}{\partial\varepsilon_{k\alpha}}\right),$ (S21)
$\displaystyle\sigma^{\rm{RTA}}_{\rm{inter}}$
$\displaystyle=2\tau^{-1}_{2}\int\frac{dk}{2\pi}\frac{|\langle
u_{+,k}|i\partial_{k}|u_{-,k}\rangle|^{2}\Delta_{k}}{\Delta_{k}^{2}+\tau_{2}^{-2}}\delta
f_{k}.$ (S22)

In insulating systems, the intraband contribution
$\sigma^{\rm{RTA}}_{\rm{intra}}$ vanishes because they have no Fermi surface.
On the other hand, the interband contribution $\sigma^{\rm{RTA}}_{\rm{iner}}$
has finite value even in insulating systems and grows roughly in proportion to
the interband damping $\tau_{2}^{-1}$.

The damping $i\tau_{2}$ in the denominator of Eq. (S19) also affects the Hall
conductivity. The above formulation in 1D can be easily extended to 2D, with
which the Hall conductivity can be calculated as

$\displaystyle\frac{\sigma_{xy}-\sigma_{yx}}{2}$
$\displaystyle=\frac{e^{2}}{\hbar}\int\frac{d^{2}\bm{k}}{(2\pi)^{2}}\mathrm{Im}\left[\frac{\langle
u_{+,\bm{k}}|i\partial_{k_{x}}|u_{-,\bm{k}}\rangle\langle
u_{-,\bm{k}}|i\partial_{k_{y}}|u_{+,\bm{k}}\rangle\Delta_{\bm{k}}}{\Delta_{\bm{k}}-i\tau_{2}^{-1}}-\frac{\langle
u_{+,\bm{k}}|i\partial_{k_{y}}|u_{-,\bm{k}}\rangle\langle
u_{-,\bm{k}}|i\partial_{k_{x}}|u_{+,\bm{k}}\rangle\Delta_{\bm{k}}}{\Delta_{\bm{k}}-i\tau_{2}^{-1}}\right]\delta
f_{\bm{k}}$ (S23)
$\displaystyle=\frac{e^{2}}{\hbar}\int\frac{d^{2}\bm{k}}{(2\pi)^{2}}\frac{\Delta_{\bm{k}}^{2}}{\Delta_{\bm{k}}^{2}+\tau_{2}^{-2}}2\mathrm{Im}\left[\langle
u_{+,\bm{k}}|i\partial_{k_{x}}|u_{-,\bm{k}}\rangle\langle
u_{-,\bm{k}}|i\partial_{k_{y}}|u_{+,\bm{k}}\rangle\right]\delta f_{\bm{k}}$
(S24)
$\displaystyle=-\frac{e^{2}}{h}\sum_{\alpha}\int\frac{d^{2}\bm{k}}{2\pi}\frac{\Delta_{\bm{k}}^{2}}{\Delta_{\bm{k}}^{2}+\tau_{2}^{-2}}\Omega_{\alpha}(\bm{k})f_{D}(\varepsilon_{\bm{k}\alpha}),$
(S25)

where we introduce the Berry curvature
$\Omega_{\alpha}(\bm{k})=2\mathrm{Im}[\langle\partial_{k_{x}}u_{\alpha,\bm{k}}|\partial_{k_{y}}u_{\alpha,\bm{k}}\rangle]$.

We recover the well-known TKNN formula for the anomalous Hall effect by taking
$\tau_{2}^{-1}\to 0$, which is known to be quantized in insulating systems at
zero temperature. While the nonzero damping $\tau_{2}^{-1}\neq 0$ leads to the
deviation from the quantized value in the RTA, this is also an artifact of the
approximation as in the longitudinal conductivity derived above.

##  Appendix D The derivation of the quantum master equation

To reveal the dynamics of the two-band systems coupled to the fermionic
reservoir, we here derive the QME in the snapshot basis microscopically. The
QME determines the reduced density operator
$\hat{\rho}(t)=\mathrm{Tr_{B}}\hat{\rho}_{\rm{tot}}(t)$, where
$\hat{\rho}_{\rm{tot}}(t)$ is the density matrix of the total system. If
Hamiltonian is written as the sum of the momentum subspace
$\hat{H}=\sum_{k}\hat{H}(k)$, we can construct the QME for the single electron
with $k$. We consider

$\begin{split}\frac{d\tilde{\rho}_{k}(t)}{dt}=-\int_{-\infty}^{t}\mathrm{Tr_{B}}\left[\tilde{H}_{I,k}(t),\left[\tilde{H}_{I,k}(s),\tilde{\rho}_{k}(t)\otimes\tilde{\rho}_{B}(t_{0})\right]\right]ds,\end{split}$
(S26)

where $\tilde{H}_{I,k}$ is the interaction Hamiltonian between the system and
the bath. Here the tilde symbol denotes the interaction picture based on the
system decoupled from the reservoir. Equation (S26) is known as the Redfield
equation [45]. We note that the following two approximations are employed in
deriving the Redfield equation; (1) The Born approximation, i.e., the
assumption that the coupling between the system and the bath is weak. (2) The
Markov approximation, i.e., the assumption that the dynamics has no memory
effect.

For the fermionic reservoir employed in Ref. [14], the interaction Hamiltonian
is given by

$\tilde{H}_{I,k}(t)=\sum_{\sigma,p}V_{p}\tilde{b}^{{\dagger}}_{k\sigma
p}(t)\tilde{c}_{k\sigma}(t)+\text{H.c.},$ (S27)

where $\tilde{c}_{k\sigma}(t)$ is a fermion annihilation operator of an
electron with $k$ and pseudospin $\sigma$, while $\tilde{b}_{k\sigma
p}(t)=\tilde{b}_{k\sigma p}(0)e^{-i\omega_{p}t}$ is an annihilation operator
of an electron in the fermionic reservoir whose mode energy is $\omega_{p}$.

Let us compute the partial trace in Eq. (S26) using the explicit expression
Eq. (S27). Before that, we introduce the snapshot basis for the field operator
$\tilde{\psi}_{k\alpha}(t)$ defined as Eq. (S9), in order to express the $s$
integral in a form convenient for the adiabatic perturbation theory. We can
rewrite the QME as

$\begin{split}\frac{d\tilde{\rho}_{k}(t)}{dt}&=\sum_{\alpha,\beta}\int_{0}^{\infty}ds\
\mathcal{C}_{\rm
p}(s)\innerproduct{\Phi_{k\beta}(t-s)}{\Phi_{k\alpha}(t)}\left[\tilde{\psi}^{{\dagger}}_{k\beta}(t-s)\tilde{\rho}_{k}(t),\tilde{\psi}_{k\alpha}(t)\right]+\text{H.c.}\\\
&+\sum_{\alpha,\beta}\int_{0}^{\infty}ds\ \mathcal{C}_{\rm
h}(s)\innerproduct{\Phi_{k\alpha}(t)}{\Phi_{k\beta}(t-s)}\left[\tilde{\psi}_{k\beta}(t-s)\tilde{\rho}_{k}(t),\tilde{\psi}^{{\dagger}}_{k\alpha}(t)\right]+\text{H.c.},\end{split}$
(S28)

where we introduced the particle/hole time correlation function of the
fermionic reservoir [56, 14]

$\displaystyle\mathcal{C}_{\rm p}(s)$
$\displaystyle=\sum_{p}|V_{p}|^{2}\text{Tr}_{\text{B}}\left[\tilde{b}_{k\sigma
p}^{\dagger}(t)\tilde{b}_{k\sigma p}(t-s)\tilde{\rho}_{B}(t_{0})\right],$
(S29) $\displaystyle\mathcal{C}_{\rm h}(s)$
$\displaystyle=\sum_{p}|V_{p}|^{2}\text{Tr}_{\text{B}}\left[\tilde{b}_{k\sigma
p}(t)\tilde{b}_{k\sigma p}^{\dagger}(t-s)\tilde{\rho}_{B}(t_{0})\right].$
(S30)

Using the broadband condition Eq. (12), we can explicitly calculate the
correlation function as

$\displaystyle\mathcal{C}_{\rm h}(s)=\mathcal{C}_{\rm p}(s)$
$\displaystyle=2\Gamma\int\frac{d\omega}{2\pi}e^{i\omega
s}f_{D}(\omega)=\Gamma\left(\delta(s)-i\frac{T}{\sinh(\pi Ts)}\right),$ (S31)

which decays with a decay time $\tau_{B}=1/\pi T$ ($k_{B}=1$) [56].
Substituting Eq. (S31) into Eq. (S28), we obtain

$\begin{split}\frac{d\tilde{\rho}_{k}(t)}{dt}&=\Gamma\sum_{\alpha}\left(\tilde{\psi}^{{\dagger}}_{k\alpha}(t)\tilde{\rho}_{k}(t)\tilde{\psi}_{k\alpha}(t)+\tilde{\psi}_{k\alpha}(t)\tilde{\rho}_{k}(t)\tilde{\psi}^{{\dagger}}_{k\alpha}(t)-\tilde{\rho}_{k}(t)\right)\\\
&-\Gamma\sum_{\alpha,\beta}\int_{0}^{\infty}ds\ \frac{iT}{\sinh(\pi
Ts)}\innerproduct{\Phi_{k\alpha}(t)}{\Phi_{k\beta}(t-s)}\left[\left\\{\tilde{\rho}_{k}(t),\tilde{\psi}_{k\beta}(t-s)\right\\},\tilde{\psi}^{{\dagger}}_{k\alpha}(t)\right]+\text{H.c.}\end{split}$
(S32)

Note that we used $\int_{0}^{\infty}ds\delta(s)=1/2$ for the first term.

Due to the decaying nature of the correlation function $\mathcal{C}_{\rm
p/h}(s)$, the integrand of Eq. (S32) only need to be accurate in
$s\lesssim\tau_{B}$. Here we assume that the change in the momentum $\Delta
k=Es$ is small enough in $s\lesssim\tau_{B}$, with which one can approximate
$\ket{\Phi_{k\alpha}(t-s)}$ and $\tilde{\psi}_{k\alpha}(t-s)$ using the Taylor
expansion.

Let us start with the expansion of the phase factor
$e^{-i\Theta_{\alpha}(t)}=e^{-i\int_{t_{0}}^{t}dt^{\prime}\left[\varepsilon_{k\alpha}(t^{\prime})+EA_{\alpha\alpha}(t^{\prime})\right]dt^{\prime}}$.
The exponent at $t-s$ is expanded as

$-i\int_{t_{0}}^{t-s}dt^{\prime}\left[\varepsilon_{k\alpha}(t^{\prime})+EA_{\alpha\alpha}(t^{\prime})\right]dt^{\prime}=-i\int_{t_{0}}^{t}dt^{\prime}\left[\varepsilon_{k\alpha}(t^{\prime})+EA_{\alpha\alpha}(t^{\prime})\right]dt^{\prime}+i\varepsilon_{k\alpha}(t)s+iEsA_{\alpha\alpha}(t)+\frac{iEs^{2}}{2}\frac{\partial\varepsilon_{k\alpha}(t)}{\partial
k}+\cdots.$ (S33)

The third and subsequent terms vanish in the adiabatic limit, so that here we
only take account of them in the leading order in $\Delta k=Es$. Then the
phase factor is approximated as

$e^{-i\Theta_{\alpha}(t-s)}\simeq\left(1+iEsA_{\alpha\alpha}(t)\right)e^{-i\Theta_{\alpha}(t)+i\varepsilon_{k\alpha}(t)s}$
(S34)

where we use $e^{iEsA_{\alpha\alpha}(t)}\simeq 1+iEsA_{\alpha\alpha}(t)$.
Because the factor $e^{+i\varepsilon_{k\alpha}(t)s}$ plays an important role
in the later discussion, we here call this approximation the dynamical phase
approximation (DPA).

In a similar manner, we can expand the instantaneous eigenstate
$\ket{u_{k\alpha}(t-s)}$ as

$\displaystyle\ket{u_{k\alpha}(t-s)}$
$\displaystyle\simeq\ket{u_{k\alpha}(t)}-s\frac{d}{dt}\ket{u_{k\alpha}(t)}$
(S35)
$\displaystyle=\ket{u_{k\alpha}(t)}-iEs\sum_{\beta}A_{\beta\alpha}(t)\ket{u_{k\beta}(t)},$
(S36)

where the terms with the second or higher order of $s$ were dropped. With
these, we can evaluate the snapshot basis
$\ket{\Phi_{k\alpha}(t)}=\ket{u_{k\alpha}(t)}e^{-i\Theta_{\alpha}(t)}$ at
$t-s$ with the DPA as

$\begin{split}\ket{\Phi_{k\alpha}(t-s)}&\simeq\left(1+iEsA_{\alpha\alpha}\right)e^{-i\Theta_{\alpha}(t)+i\varepsilon_{k\alpha}(t)s}\left(\ket{u_{k\alpha}(t)}-iEs\sum_{\beta}A_{\beta\alpha}(t)\ket{u_{k\beta}(t)}\right)\\\
&\simeq
e^{+i\varepsilon_{k\alpha}(t)s}\left(\ket{\Phi_{k\alpha}(t)}-is\sum_{\beta}\left[\mathcal{W}_{k}(t)\right]_{\beta\alpha}\ket{\Phi_{k\beta}(t)}\right).\end{split}$
(S37)

Next, let us consider the snapshot basis operator
$\tilde{\psi}_{k\alpha}(t-s)$ with the DPA. Since the snapshot basis operator
evolves as
$\tilde{\psi}_{k\alpha}(t)=\sum_{\beta}[U(t,t^{\prime})]_{\alpha\beta}\tilde{\psi}_{k\beta}(t^{\prime})$
with Eq. (S8) [14], the quantity we need to approximate here is

$\left[U(t-s,t)\right]_{\alpha\beta}=\delta_{\alpha\beta}-i\int_{0}^{-s}dt_{1}\left[\mathcal{W}_{k}(t+t_{1})\right]_{\alpha\beta}+(-i)^{2}\int_{0}^{-s}dt_{1}\int_{0}^{t_{1}}dt_{2}\sum_{\gamma}\left[\mathcal{W}_{k}(t+t_{1})\right]_{\alpha\gamma}\left[\mathcal{W}_{k}(t+t_{2})\right]_{\gamma\beta}+\cdots.$
(S38)

In order to reproduce the first order of $Es$, it is sufficient to truncate at
the second term and approximate the integrand as

$\left[\mathcal{W}_{k}(t+t_{1})\right]_{\alpha\beta}\simeq\left[\mathcal{W}_{k}(t)\right]_{\alpha\beta}e^{\alpha
i\Delta_{k}(t)t_{1}},$ (S39)

according to Eq. (S34). However, this approximation turns out to show
unphysical behavior due to the broken unitarity of $U(t,t^{\prime})$ upon the
truncation. Also, we note that the norm of $\ket{\Phi_{k\alpha}(t-s)}$ is not
conserved in this treatment.

To remedy this problem, we here include the third term in Eq. (S38) while
retaining the use of Eq. (S39), which results in

$\left[U(t-s,t)\right]_{\alpha\beta}\simeq\delta_{\alpha\beta}-\alpha\left[\mathcal{W}_{k}(t)\right]_{\alpha\beta}\frac{e^{-\alpha
i\Delta_{k}(t)s}-1}{\Delta_{k}(t)}+\frac{|W_{k}(t)|^{2}}{\Delta_{k}(t)}\left(i\alpha
s+\frac{e^{-\alpha
i\Delta_{k}(t)s}-1}{\Delta_{k}(t)}\right)\delta_{\alpha\beta}.$ (S40)

With this approximation, the transformed operator

$\displaystyle\tilde{\psi}_{k\alpha}^{\prime}(t,s)$
$\displaystyle=\sum_{\beta}\innerproduct{\Phi_{k\alpha}(t)}{\Phi_{k\beta}(t-s)}\tilde{\psi}_{k\beta}(t-s)$
(S41)

appearing in Eq. (S32) is approximated as

$\displaystyle\tilde{\psi}_{k\alpha}^{\prime}(t,s)$
$\displaystyle\simeq\sum_{\beta}e^{i\varepsilon_{k\alpha}(t)s}\left[\delta_{\alpha\beta}+\left(-\alpha\left[\mathcal{W}_{k}(t)\right]_{\alpha\beta}+\frac{|W_{k}(t)|^{2}}{\Delta_{k}(t)}\delta_{\alpha\beta}\right)\left(i\alpha
se^{-\alpha i\Delta_{k}(t)s}+\frac{e^{-\alpha
i\Delta_{k}(t)s}-1}{\Delta_{k}(t)}\right)\right]\tilde{\psi}_{k\beta}(t),$
(S42)

which satisfies

$\displaystyle\left\\{\tilde{\psi}_{k\alpha}^{\prime}(t,s),\tilde{\psi}_{k\beta}^{\prime{\dagger}}(t,s)\right\\}$
$\displaystyle=\left(1+|W_{k}(t)|^{2}s^{2}\right)\delta_{\alpha\beta}+\mathcal{O}(E^{3}).$
(S43)

This implies that the unitary nature of the time evolution of the isolated
system is maintained up to the first order of $s$ and the second order of $E$.

With the use of Eq. (S42) and

$\displaystyle\int_{0}^{\infty}\frac{T\sin(\varepsilon s)}{\sinh(\pi Ts)}ds$
$\displaystyle=-f_{D}(\varepsilon)+\frac{1}{2},$ (S44)
$\displaystyle-i\int_{0}^{\infty}\frac{Ts^{l}e^{-i\varepsilon s}}{\sinh(\pi
Ts)}ds$
$\displaystyle=i^{l}\frac{\partial^{l}}{\partial\varepsilon^{l}}\left[f_{D}(\varepsilon)+\frac{i}{\pi}\mathrm{Re}\Psi\left(\frac{1}{2}+\frac{i\varepsilon}{2\pi
T}\right)\right]\quad(l\geq 1),$ (S45)

we can perform the integral in Eq. (S32). Here $\Psi$ is the digamma function.
The resultant equation is given by

$\begin{split}\frac{d\tilde{\rho}_{k}(t)}{dt}&\simeq\tilde{\mathcal{D}}_{0}(t)+\tilde{\mathcal{D}}_{1}(t)+\tilde{\mathcal{D}}_{2}(t),\\\
\tilde{\mathcal{D}}_{0}(t)&=2\Gamma\sum_{\alpha}f_{D}(\varepsilon_{k\alpha}(t))\left(\tilde{\psi}^{{\dagger}}_{k\alpha}(t)\tilde{\rho}_{k}(t)\tilde{\psi}_{k\alpha}(t)-\frac{1}{2}\left\\{\tilde{\rho}_{k}(t),\tilde{\psi}_{k\alpha}(t)\tilde{\psi}^{{\dagger}}_{k\alpha}(t)\right\\}\right)\\\
&+2\Gamma\sum_{\alpha}(1-f_{D}(\varepsilon_{k\alpha}(t)))\left(\tilde{\psi}_{k\alpha}(t)\tilde{\rho}_{k}(t)\tilde{\psi}^{{\dagger}}_{k\alpha}(t)-\frac{1}{2}\left\\{\tilde{\rho}_{k}(t),\tilde{\psi}^{{\dagger}}_{k\alpha}(t)\tilde{\psi}_{k\alpha}(t)\right\\}\right),\\\
\tilde{\mathcal{D}}_{1}(t)&=2\Gamma\sum_{\alpha\beta}\frac{\left[\mathcal{W}_{k}(t)\right]_{\alpha\beta}}{\Delta_{k}(t)}\left(\delta
f_{k}(t)+\frac{\Delta_{k}(t)}{2}\frac{\partial}{\partial\varepsilon}\left[f_{D}(\varepsilon_{k,-}(t))+f_{D}(\varepsilon_{k,+}(t))-\frac{i\alpha}{\pi}\delta\Psi_{k}(t)\right]\right)\\\
&\times\left(\tilde{\psi}_{k\beta}(t)\tilde{\rho}_{k}(t)\tilde{\psi}^{{\dagger}}_{k\alpha}(t)-\tilde{\psi}^{{\dagger}}_{k\alpha}(t)\tilde{\rho}_{k}(t)\tilde{\psi}_{k\beta}(t)-\left\\{\tilde{\rho}_{k}(t),\tilde{\psi}^{{\dagger}}_{k\alpha}(t)\tilde{\psi}_{k\beta}(t)\right\\}\right),\\\
\tilde{\mathcal{D}}_{2}(t)&=-2\Gamma\sum_{\alpha}\alpha\frac{|W_{k}(t)|^{2}}{\Delta_{k}^{2}(t)}\left(\delta
f_{k}(t)+\Delta_{k}(t)\frac{\partial
f_{D}(\varepsilon_{k,-\alpha}(t))}{\partial\varepsilon}\right)\\\
&\times\left(\tilde{\rho}_{k}(t)+\tilde{\psi}_{k\alpha}(t)\tilde{\rho}_{k}(t)\tilde{\psi}^{{\dagger}}_{k\alpha}(t)-\tilde{\psi}^{{\dagger}}_{k\alpha}(t)\tilde{\rho}_{k}(t)\tilde{\psi}_{k\alpha}(t)-\left\\{\tilde{\rho}_{k}(t),\tilde{\psi}^{{\dagger}}_{k\alpha}(t)\tilde{\psi}_{k\alpha}(t)\right\\}\right),\end{split}$
(S46)

where
$\delta\Psi_{k}(t)=\mathrm{Re}\Psi(\frac{1}{2}-\frac{i\beta\varepsilon_{k-}(t)}{2\pi})-\mathrm{Re}\Psi(\frac{1}{2}-\frac{i\beta\varepsilon_{k+}(t)}{2\pi})$.

We redefine the density matrix in the snapshot basis as
$[\rho_{k}(t)]_{\alpha\beta}=\mathrm{Tr}[\tilde{\psi}^{{\dagger}}_{k\beta}(t)\tilde{\psi}_{k\alpha}(t)\tilde{\rho}_{k}(t)]$.
We obtain the equation of motion for the matrix element of $\rho_{k}(t)$ as

$\frac{d[\rho_{k}(t)]_{\alpha\beta}}{dt}=-i\left[[\mathcal{W}_{k}(t),\rho_{k}(t)]\right]_{\alpha\beta}+\mathrm{Tr}\left[\tilde{\psi}^{{\dagger}}_{k\beta}(t)\tilde{\psi}_{k\alpha}(t)\frac{d\tilde{\rho}_{k}(t)}{dt}\right].$
(S47)

The second term of the right hand side in Eq. (S47) corresponds to the
dissipation term. Substituting Eq. (S46) into Eq. (S47), we finally arrive at

$\displaystyle[{\mathcal{D}}_{0}(t)]_{\alpha\beta}$
$\displaystyle=-2\Gamma([\rho_{k}(t)]_{\alpha\beta}-f_{D}(\varepsilon_{k\alpha}(t))\delta_{\alpha\beta}),$
(S48) $\displaystyle{}[{\mathcal{D}}_{1}(t)]_{\alpha\beta}$
$\displaystyle=-2\Gamma\frac{\left[\mathcal{W}_{k}(t)\right]_{\alpha\beta}}{\Delta_{k}(t)}\left(\delta
f_{k}(t)+\frac{\Delta_{k}(t)}{2}\frac{\partial}{\partial\varepsilon}\left[f_{D}(\varepsilon_{k,-}(t))+f_{D}(\varepsilon_{k,+}(t))-\frac{i\alpha}{\pi}\delta\Psi_{k}(t)\right]\right),$
(S49) $\displaystyle{}[{\mathcal{D}}_{2}(t)]_{\alpha\beta}$
$\displaystyle=2\Gamma\alpha\frac{|W_{k}(t)|^{2}}{\Delta_{k}^{2}(t)}\left(\delta
f_{k}(t)+\Delta_{k}(t)\frac{\partial
f_{D}(\varepsilon_{k,-\alpha}(t))}{\partial\varepsilon}\right)\delta_{\alpha\beta}$
(S50)

where $\mathcal{D}_{n}$ is the $n$-th order of $E$ term in the dissipation
term. In the presence of the particle-hole symmetry, the derivative term in
$\mathcal{D}_{1}$ vanishes since
$f_{D}(\varepsilon_{k+}(t))+f_{D}(\varepsilon_{k-}(t))=1$ and
$\delta\Psi_{k}(t)=0$.