The driven-Markovian master equation based on the Lewis-Riesenfeld invariants theory

S. L. Wu [email protected] School of Physics and Materials Engineering, Dalian Nationalities University, Dalian 116600 China X. L. Huang School of Physics and Electronic Technology, Liaoning Normal University, Dalian 116029, China X. X. Yi [email protected] Center for Quantum Sciences and School of Physics, Northeast Normal University, Changchun 130024, China

Abstract

We derive a Markovian master equation for driven open quantum systems based on the Lewis-Riesenfeld invariants theory, which is available for arbitrary driving protocols.The role of the Lewis-Riesenfeld invariants is to help us bypass the time-ordering obstacle in expanding the propagator of the free dynamics, such that the Lindblad operators in our driven-Markovian master equation can be determined easily. We also illustrate that, for the driven open quantum systems, the spontaneous emission and the thermal excitation induce the transitions between eigenstates of the Lewis-Riesenfeld invariant, but not the system Hamiltonian’s. As an example, we present the driven-Markovian master equation for a driven two-level system coupled to a heat reservoir. By comparing to the exactly solvable models, the availability of the driven-Markovian master equation is verified. Meanwhile, the adiabatic limit and inertial limit of the driven-Markovian master equation are also discussed, which result in the same Markovian master equations as those presented before in the corresponding limits.

pacs:

03.67.-a, 03.65.Yz, 05.70.Ln, 05.40.Ca

I Introduction

Any physical system in nature, no matter classical or quantum, couples to its surroundings exchanging energy and matter to make it as an open system. The theory of open quantum systems aims at providing a concise manner to describe the dynamics of the primary systemBreuer2007 . For open quantum systems satisfying the Born-Markov approximationDavies1974 , the Gorini-Kossakowski- Lindblad-Sudarshan (GKLS) Markovian master equation gives a general completely positive and trace preserving map of the reduced dynamics Gorini1976 ; Lindblad1976 . In the original derivation, it is assumed that the system Hamiltonian is time-independent. The coupling to the environment induces transitions between the static eigenstates of the system Hamiltonian. For the open systems with time-dependent external drives, the GKLS master equations have been derived in and beyond the adiabatic limit, which leads to the adiabatic Davies1978 ; Albash2012 ; Childs2001 ; Kamleitner2013 ; Sarandy2004 and non-adiabatic Yamaguchi2017 ; Dann2018 ; Potts2021 Markovian master equation.

For the driven open quantum systems without memory of the driving protocol, a non-adiabatic Markovian master equation (NAME) has been derived Dann2018 . The Lindblad operators in this non-adiabatic master equation are eigenoperators of the propagator of the free dynamics which associates with the system Hamiltonian. Generally speaking, these eigenoperators can be determined by representing the dynamics in the operator space which is also known as the Liouvillian space Albert2016 ; Scopa2019 . However, it is difficult to give these eigenoperators explicitly, and these eigenoperators absent clear physical meanings . For second best thing, a method based on the inertial theorem is proposedDann2020 . The inertial theorem relies on a priori decomposition of the Liouvillian superoperator into a rapid-changing scalar parameter and a slow-changing superoperator, which is equivalent to additional restrictions on the driving protocol Dann2021 . In this way, the Lindblad operators can be obtained explicitly, if the inertial limit is reachedDann2020 ; Dann2021 ; Dann2019 . At the same time, effective numerical methods to simulate the driven open quantum system dynamics are also proposed Hwang2012 ; Stockburger2016 , but may provide less structural insight into the dynamics.

For the open quantum systems with a static Hamiltonian, the population transitions caused by decoherence occur between the eigenstates of the static Hamiltonian Breuer2007 . Hence, even if the Markovian equation can not be given explicitly, we may formulate a phenomenological master equation due to the clearly physical meanings of the transitions Ozaki2021 ; Taran2019 . However, it is totally different for the driven open quantum systems with a non-adiabatic driving protocol, since there is no such a physical meaning for the decoherence-induced transitions. Thus, it is natural to ask if there is a simple method to formulate a Markovian master equation for a driven open quantum system with arbitrary driving protocols as we used in formulating the Markovian master equation with the static Hamiltonian?

In this paper, we derive a driven-Markovian master equation (DMME) for arbitrary driving protocols by using the Lewis-Riesenfeld invariants (LRIs) Lewis1968 , which is easy to formulate and has a clear meaning of the decoherence-induced transitions. Since the solution of the Shrödinger equation with a time-dependent Hamiltonian can be expressed as a superposition of the eigenstates of the Lewis-Riesenfeld invariant with constant amplitudes Lewis1969 , the unitary operator corresponding to the free propagator can be written down explicitly. On the other hand, if the the timescale for the driving protocol, also known as ”the non-adiabatic timescale”, approaches to or is larger than the reservoir correlation time, the memory effect of the driving protocol can not be neglected. By using the Fourier transformation and its inverse transformation, we collect the frequency contribution caused by the driving memory effect into the one-side Fourier transformation of the reservoir correlation function. Thus, the driving memory effect can be included in the DMME with a concise manner. As we shall see, the transitions induced by the coupling to the environment only occur between the invariant’s eigenstates. Therefore, it is quiet straightforward to formulate a DMME for a general driven open system, if its LRIs are known.

This paper is organized as follows. In Sec. II, we present the general formula of the DMME based on the Lewis-Riesenfeld invariants theory. The memory of the driven protocol is encoded in the decoherence rates and the Lamb shifts of the DMME. Then, we apply this DMME to the two-level system with time-dependent driving fields in Sec. III, and the corresponding adiabatic and inertial limits are discussed. We also derived exact dynamics for the driven two-level system, which will help to illustrate validity and availability of the DMME. Finally, the conclusions are given in Sec. IV.

II General Formalism

In this section, we apply the LRIs theory to present the DMME with explicit mathematical and physical meaning. Consider the dynamics of the total system which is governed by the Hamiltonian

H(t)=H_{\text{s}}(t)+H_{\text{B}}+H_{\text{I}}.

$H_{\text{s}}(t)$ stands for the system Hamiltonian; the reservoir is represented by the Hamiltonian

H_{\text{B}}=\sum_{k}\omega_{k}b_{k}^{\dagger}b_{k},

where $b_{k}$ and $\omega_{k}$ are, respectively, the annihilation operator and the eigenfrequency of the $k$ -th mode of the reservoir. In the following, the natural units $\hbar=c=1$ are used throughout. We assume that the system-reservoir interaction Hamiltonian is given by

H_{\text{I}}=\sum_{k}\textsl{g}_{k}A_{k}\otimes B_{k}.

$A_{k}$ and $B_{k}$ are the Hermitian operators of the system and reservoir, respectively. $g_{k}$ stands for the coupling strength.

The von Neumann equation for the density operator of the total system in the interaction picture reads

\partial_{t}\tilde{\rho}(t)=-i\left[\tilde{H_{\text{I}}}(t),\tilde{\rho}(t)\right],

where $\tilde{\rho}(t)$ denotes the density operator of the total system in the interaction picture, and a similar notations is applied for the other system and reservoir operators. By assuming the weak system-reservoir coupling (the Born approximation), we obtain the Born equation for the system density operator $\tilde{\rho}_{\text{s}}(t)$ ,

	$\displaystyle\partial_{t}$		$\displaystyle\tilde{\rho}_{\text{s}}(t)=$
			$\displaystyle-\int_{0}^{t}d\tau\text{Tr}_{\text{B}}\left\{\left[\tilde{H_{\text{I}}}(t),\left[\tilde{H_{\text{I}}}(t-\tau),\tilde{\rho}_{\text{s}}(t-\tau)\otimes\tilde{\rho}_{\text{B}}\right]\right]\right\}.$

Here, we have assumed that $\text{Tr}_{\text{B}}\{\tilde{H_{\text{I}}}(t)\tilde{\rho}(0)\}=0$ , and the initial state of the total system can be written as $\tilde{\rho}(0)=\tilde{\rho}_{\text{s}}(0)\otimes\tilde{\rho}_{\text{B}}$ . The Born approximation assumes that the coupling between the system and the reservoir is weak, such that the influence of the system on the reservoir is small. If the system evolution time is much larger than the reservoir correlation time $\tau_{\text{B}}$ , we can replace $\tilde{\rho}_{\text{s}}(t-\tau)$ by $\tilde{\rho}_{\text{s}}(t)$ and the integral limits can be extended to $\infty$ , which is known as the Markovian approximation. In such a case, the dynamics governed by the following Redfield master equation within the Born-Markovian approximation Redfield1965 ,

	$\displaystyle\partial_{t}$		$\displaystyle\tilde{\rho}_{\text{s}}(t)=$		(1)
			$\displaystyle-\int_{0}^{\infty}d\tau\text{Tr}_{\text{B}}\left\{\left[\tilde{H_{\text{I}}}(t),\left[\tilde{H_{\text{I}}}(t-\tau),\tilde{\rho}_{\text{s}}(t)\otimes\tilde{\rho}_{\text{B}}\right]\right]\right\}.$		(1)

For an operator $A$ of the system, the corresponding operator in the interaction picture can be connected by an unitary transformation, i.e.,

\displaystyle\tilde{A}_{k}(t)=\hat{\mathcal{U}}_{\text{s}}(t)A_{k}=U_{\text{s}}^{\dagger}(t)A_{k}U_{\text{s}}(t).

(2)

$U_{\text{s}}(t)$ describing the free propagator of the system satisfies a Schrödinger equation for the time-dependent system Hamiltonian

i\partial_{t}U_{\text{s}}(t)=H_{\text{s}}(t)U_{\text{s}}(t),\>U_{\text{s}}(0)=I,

(3)

which results in $U_{\text{s}}(t)=\mathcal{T}\exp\left(-i\int_{0}^{t}\text{d}\tau\,H_{\text{s}}(\tau)\right)$ with the time-ordering operator $\mathcal{T}$ .

To reduce Eq.(1), a set of eigenoperator of the superoperator $\mathcal{\hat{U}}_{\text{s}}(t)$ are needed, where the eigenoperators satisfy $\tilde{F}_{j}(t)=\hat{\mathcal{U}}_{\text{s}}(t)\tilde{F}_{j}(0)=\lambda_{j}(t)\tilde{F}_{j}(0)$ Dann2018 . However, it is difficult to solve the eigenequation of the superoperator $\mathcal{\hat{U}}_{\text{s}}(t)$ directly. To overcome this difficulty, the inertial theorem has been used to obtain an approximative solution of $\hat{\mathcal{U}}_{\text{s}}(t)$ Dann2021 ; Dann2020 . However, the solution is accurate only under that the inertial parameter is small, which requires a slow acceleration of the drive.

In fact, the free propagator of the system can be obtained directly, if the Lewis-Riesenfeld invariants for the system Hamiltonian $H_{\text{s}}(t)$ are known. A LRI $I_{\text{s}}(t)$ for the systems with the Hamiltonian $H_{\text{s}}(t)$ is a Hermitian operator which obeys an equation in the Schrödinger picture Lewis1969

i\partial_{t}I_{\text{s}}(t)-\left[H_{\text{s}}(t),I_{\text{s}}(t)\right]=0.

(4)

For the closed dynamics of the systems, a general solution of the Schrödinger equation can be written as

|\Psi(t)\rangle=\sum_{n=1}^{N}c_{n}\text{e}^{i\alpha_{n}(t)}|\psi_{n}(t)\rangle.

(5)

Here, $|\psi_{n}(t)\rangle$ is the $n$ -th eigenstate of the LRI with a real constant eigenvalue $\lambda_{n}$ , i.e., $I_{\text{s}}(t)|\psi_{n}(t)\rangle=\lambda_{n}|\psi_{n}(t)\rangle$ , $\{c_{n}\}$ are time-independent amplitudes, and the Lewis-Riesenfeld phases are defined as Lewis1969

\alpha_{n}(t)=\int_{0}^{t}\langle\psi_{n}(\tau)|\left(i\partial_{\tau}-H_{\text{s}}(\tau)\right)|\psi_{n}(\tau)\rangle\,\text{d}\tau.

(6)

Therefore, the solution of Eq.(3) can be expressed by means of the eigenstates of the LRIs,

U_{\text{s}}(t)=\sum_{n}\text{e}^{i\alpha_{n}(t)}|\psi_{n}(t)\rangle\langle\psi_{n}(0)|.

(7)

The LRIs theory was designed to investigate the time evolution of dynamical systems with an explicitly time-dependent Hamiltonian Lewis1969 . The invariants comply with the following properties: (i) The expectation values of the LRIs are constant. (ii) The eigenvalues of a LRI are constant, while the eigenstates are time dependent. (iii) Any time-dependent Hermitian operator which satisfies Eq. (4) is a LRI for closed quantum systems. Each LRI corresponds to a symmetry of the closed quantum system. Thereafter, the LRIs are successfully applied to investigate time-dependent problems in quantum mechanics Pedrosa2009 such as the Berry phase Sarandy2007 , the connection between quantum theory and classical theory Schuch2006 , and the quantum control Chen2011 . At the same time, the method to construct the LRIs for various quantum systems has been explored, for instance, the harmonic oscillator system Lewis1969 , the few-level systems Monteoliva1994 ; Nakahara2012 , the pseudo-Hermitian systemSimeonov2016 , and the open fermionic systemsKim2000 . Also a general method for constructing LRIs has been proposedPonte2018 .

By means of the explicit formula of the free evolution operator $U_{\text{s}}(t)$ (Eq.(7)), the system operator Eq.(2) in the interaction picture can be rewritten as

	$\displaystyle\tilde{A}_{k}(t)$	$\displaystyle=$	$\displaystyle U_{\text{s}}^{\dagger}(t)A_{k}U_{\text{s}}(t)$		(8)
		$\displaystyle=$	$\displaystyle\sum_{n,m}\text{e}^{i\theta_{mn}^{k}(t)}\xi_{mn}^{k}(t)\tilde{F}_{mn}$		(8)

with

\displaystyle\theta_{mn}^{k}(t)=\alpha_{n}(t)-\alpha_{m}(t)+\text{Arg}\left(\langle\psi_{m}(t)|A_{k}|\psi_{n}(t)\rangle\right)

(9)

and $\xi_{mn}^{k}(t)=|\langle\psi_{m}(t)|A_{k}|\psi_{n}(t)\rangle|$ which satisfy $\theta_{mn}^{k}(t)\in\mathbb{R}$ and $\xi_{mn}^{k}(t)>0.$ The time-independent operators $\tilde{F}_{mn}=|\psi_{m}(0)\rangle\langle\psi_{n}(0)|$ denotes one of Lindblad operators in the interaction picture. Since $\tilde{A}_{k}(t)$ are Hermitian operators, it yields

\tilde{A}_{k}(t)=\sum_{n^{\prime},m^{\prime}}\text{e}^{-i\theta_{m^{\prime}n^{\prime}}^{k}(t)}\xi_{m^{\prime}n^{\prime}}^{k}(t)\tilde{F}_{m^{\prime}n^{\prime}}^{\dagger}.

(10)

Any $\tilde{F}_{m^{\prime}n^{\prime}}^{\dagger}$ contains in the operator set $\left\{\tilde{F}_{mn}\right\}$ which can be used to expand the corresponding Liouvillian spacePetrosky1997 . Substituting Eqs.(8) and (10) into Eq.(1), we can express the Markovian master equation as

	$\displaystyle\partial_{t}\tilde{\rho}_{\text{s}}(t)$	$\displaystyle=$	$\displaystyle\sum_{m,m^{\prime},n,,n^{\prime}}\Gamma_{mn,m^{\prime}n^{\prime}}(t)\left(\tilde{F}_{mn}\tilde{\rho}_{\text{s}}(t)\tilde{F}_{m^{\prime}n^{\prime}}^{\dagger}\right.$		(11)
			$\displaystyle\left.-\tilde{F}_{m^{\prime}n^{\prime}}^{\dagger}\tilde{F}_{mn}\tilde{\rho}_{\text{s}}(t)\right)+\text{H.c.}$		(11)

with

$\displaystyle\Gamma_{mn,m^{\prime}n^{\prime}}$	$\displaystyle(t)=\sum_{k,k^{\prime}}\textsl{g}_{k}\textsl{g}_{k^{\prime}}$
$\displaystyle\times$	$\displaystyle\int_{0}^{\infty}\text{d}s\xi_{m^{\prime}n^{\prime}}^{k^{\prime}}(t)\xi_{mn}^{k}(t-s)\text{e}^{i\left(\theta_{mn}^{k}(t-s)-\theta_{m^{\prime}n^{\prime}}^{k^{\prime}}(t)\right)}$
$\displaystyle\times$	$\displaystyle\text{Tr}_{\text{B}}\left\{\tilde{B}_{k^{\prime}}(t)\tilde{B}_{k}(t-s)\rho_{B}\right\},$	(12)

where H.c. denotes the Hermitian conjugated expression and $\tilde{B}_{k^{\prime}}(t)$ is the bath operator in the interaction picture.

As shown in Eq.(11), there is memory effect of the driving protocol, which contains in $\xi_{mn}^{k}(t-s)$ and $\theta_{mn}^{k}(t-s)$ . At first, by means of the Taylor expansion, the phase $\theta_{mn}^{k}(t-s)$ can be written as

$\displaystyle\theta_{mn}^{k}(t-s)$	$\displaystyle=$	$\displaystyle\theta_{mn}^{k}(t)+\partial_{s}\theta_{mn}^{k}(t-s)\|_{s=0}s$	(13)
		$\displaystyle+\sum_{l=2}^{\infty}\frac{1}{l!}\partial_{s}^{l}\theta_{mn}^{k}(t-s)s^{l}$
	$\displaystyle\equiv$	$\displaystyle\theta_{mn}^{k}(t)+\alpha_{mn}^{k}(t)s+\varTheta_{mn}^{k}(t,t-s),$

where

	$\displaystyle\varTheta_{mn}^{k}(t,t-s)$	$\displaystyle=$	$\displaystyle\theta_{mn}^{k}(t-s)-\theta_{mn}^{k}(t)-\alpha_{mn}^{k}(t)s$		(14)
		$\displaystyle=$	$\displaystyle\int_{t-s}^{t}(\alpha_{mn}^{k}(\tau)-\alpha_{mn}^{k}(t))d\tau$		(14)

is a function of $t$ and $t-s$ with $\alpha_{mn}^{k}(t)=-\partial_{t}\theta_{mn}^{k}(t)$ . With the consideration of $\text{e}^{i\varTheta_{mn}^{k}}=\cos\varTheta_{mn}^{k}+i\sin\varTheta_{mn}^{k}$ , Eq.(12) becomes

$\displaystyle\Gamma$	${}_{mn,m^{\prime}n^{\prime}}(t)$	(15)
	$\displaystyle=\sum_{k,k^{\prime}}\textsl{g}_{k}\textsl{g}_{k^{\prime}}\xi_{m^{\prime}n^{\prime}}^{k^{\prime}}(t)\text{e}^{i\left(\theta_{mn}^{k}(t)-\theta_{m^{\prime}n^{\prime}}^{k^{\prime}}(t)\right)}$
	$\displaystyle\times\int_{0}^{\infty}\text{d}s\left(\Xi_{mn}^{\text{c},k}(t,t-s)+i\Xi_{mn}^{\text{s},k}(t,t-s)\right)$
	$\displaystyle\times\text{Tr}_{\text{B}}\left\{\tilde{B}_{k^{\prime}}(t)\tilde{B}_{k}(t-s)\rho_{B}\right\}\text{e}^{i\alpha_{mn}^{k}(t)s},$

where $\Xi_{mn}^{\text{c},k}(t,t-s)=\xi_{mn}^{k}(t-s)\cos\varTheta_{mn}^{k}(t,t-s)$ and $\Xi_{mn}^{\text{s},k}(t,t-s)=\xi_{mn}^{k}(t-s)\sin\varTheta_{mn}^{k}(t,t-s)$ . For $\Xi_{mn}^{\text{c(s)},k}(t,t-s)$ , we take the Fourier expansion with respect to $t-s$ ,

\Xi_{mn}^{\text{c(s)},k}(t,t-s)=\frac{1}{\sqrt{2\pi}}\int_{-\infty}^{+\infty}d\omega_{\xi}\bar{\Xi}_{mn}^{\text{c(s)},k}(t,\omega_{\xi})\text{e}^{i\omega_{\xi}(t-s)}

(16)

with $\bar{\Xi}_{mn}^{\text{c(s)},k}(t,\omega_{\xi})=\frac{1}{\sqrt{2\pi}}\int_{-\infty}^{+\infty}\Xi_{mn}^{c(s),k}(t,\tau)\text{e}^{-i\omega_{\xi}\tau}d\tau$ . By substituting Eq. (16) into Eq.(15), it yields

	$\displaystyle\Gamma$		${}_{mn,m^{\prime}n^{\prime}}(t)$
			$\displaystyle=\frac{1}{\sqrt{2\pi}}\sum_{k,k^{\prime}}\textsl{g}_{k}\textsl{g}_{k^{\prime}}\xi_{m^{\prime}n^{\prime}}^{k^{\prime}}(t)\text{e}^{i\left(\theta_{mn}^{k}(t)-\theta_{m^{\prime}n^{\prime}}^{k^{\prime}}(t)\right)}$
			$\displaystyle\times\int_{-\infty}^{+\infty}d\omega_{\xi}\left(\bar{\Xi}_{mn}^{\text{c},k}(t,\omega_{\xi})+i\bar{\Xi}_{mn}^{\text{s},k}(t,\omega_{\xi})\right)\text{e}^{i\omega_{\xi}t}$
			$\displaystyle\times\bar{\varLambda}_{kk^{\prime}}(\alpha_{mn}^{k}-\omega_{\xi}),$

where $\bar{\varLambda}_{kk^{\prime}}$ is the one-sided Fourier transform of the instantaneous reservoir correlation function

\displaystyle\bar{\varLambda}_{kk^{\prime}}(\alpha)

\displaystyle=

\displaystyle\int_{0}^{\infty}\text{d}s\text{e}^{i\alpha s}\text{Tr}_{\text{B}}\left\{\tilde{B}_{k^{\prime}}(t)\tilde{B}_{k}(t-s)\rho_{B}\right\}.

(17)

with $\alpha=\alpha_{mn}^{k}-\omega_{\xi}$ . It is convenient to decompose $\bar{\varLambda}_{kk^{\prime}}$ into a real and imaginary part, i.e.,

\bar{\varLambda}_{kk^{\prime}}(\alpha)=\bar{\varLambda}_{kk^{\prime}}^{\text{R}}(\alpha)+i\,\bar{\varLambda}_{kk^{\prime}}^{\text{I}}(\alpha),

where $\bar{\Lambda}_{kk^{\prime}}^{\text{I}}(\alpha)=-\frac{i}{2}\left(\bar{\Lambda}_{kk^{\prime}}(\alpha)-\bar{\Lambda}_{kk^{\prime}}^{*}(\alpha)\right)$ is a Hermitian matrix and $\bar{\Lambda}_{kk^{\prime}}^{\text{R}}(\alpha)$ can be written as

\bar{\varLambda}_{kk^{\prime}}^{\text{R}}(\alpha)=\frac{1}{2}\int_{-\infty}^{\infty}\text{d}s\text{e}^{i\alpha s}\text{Tr}_{\text{B}}\left\{\tilde{B}_{k}(s)\tilde{B}_{k^{\prime}}(0)\rho_{B}\right\}.

We divide $\Gamma_{mn,m^{\prime}n^{\prime}}$ into the real and imaginary parts, i.e.,

$\displaystyle\Gamma_{mn,m^{\prime}n^{\prime}}(t)$	$\displaystyle=$	$\displaystyle\frac{1}{\sqrt{2\pi}}\sum_{k,k^{\prime}}\textsl{g}_{k}\textsl{g}_{k^{\prime}}\xi_{m^{\prime}n^{\prime}}^{k^{\prime}}(t)\text{e}^{i\left(\theta_{mn}^{k}(t)-\theta_{m^{\prime}n^{\prime}}^{k^{\prime}}(t)\right)}\int_{-\infty}^{+\infty}d\omega_{\xi}\text{e}^{i\omega_{\xi}t}$	(18)
	$\displaystyle\times$	$\displaystyle\left(\left(\bar{\Xi}_{mn}^{\text{c},k}(t,\omega_{\xi})\bar{\varLambda}_{kk^{\prime}}^{\text{R}}(\alpha_{mn}^{k}-\omega_{\xi})-\bar{\Xi}_{mn}^{\text{s},k}(t,\omega_{\xi})\bar{\varLambda}_{kk^{\prime}}^{\text{I}}(\alpha_{mn}^{k}-\omega_{\xi})\right)\right.$
	$\displaystyle+$	$\displaystyle\left.i\!\left(\bar{\Xi}_{mn}^{\text{c},k}(t,\omega_{\xi})\bar{\varLambda}_{kk^{\prime}}^{\text{I}}(\alpha_{mn}^{k}-\omega_{\xi})+\bar{\Xi}_{mn}^{\text{s},k}(t,\omega_{\xi})\bar{\varLambda}_{kk^{\prime}}^{\text{R}}(\alpha_{mn}^{k}-\omega_{\xi})\right)\right).$

According to the convolution theorem of the Fourier transformation, we can transform the integral in $\Gamma_{mn,m^{\prime}n^{\prime}}$ with respect to $\omega_{\xi}$ into a convolution of the time-domain integral, which leads to

$\displaystyle\Gamma_{mn,m^{\prime}n^{\prime}}(t)$	$\displaystyle=$	$\displaystyle\frac{1}{2\pi}\sum_{k,k^{\prime}}\textsl{g}_{k}\textsl{g}_{k^{\prime}}\xi_{m^{\prime}n^{\prime}}^{k^{\prime}}(t)\text{e}^{i\left(\theta_{mn}^{k}(t)-\theta_{m^{\prime}n^{\prime}}^{k^{\prime}}(t)\right)}\int_{-\infty}^{+\infty}ds^{\prime}$
	$\displaystyle\times$	$\displaystyle\left(\left(\Xi_{mn}^{\text{c},k}(t,t-s^{\prime})\Lambda_{kk^{\prime}}^{\text{R}}(\alpha_{mn}^{k},s^{\prime})-\Xi_{mn}^{\text{s},k}(t,t-s^{\prime})\Lambda_{kk^{\prime}}^{\text{I}}(\alpha_{mn}^{k},s^{\prime})\right)\right.$
	$\displaystyle+$	$\displaystyle\left.i\!\left(\Xi_{mn}^{\text{c},k}(t,t-s^{\prime})\Lambda_{kk^{\prime}}^{\text{I}}(\alpha_{mn}^{k},s^{\prime})+\Xi_{mn}^{\text{s},k}(t,t-s^{\prime})\Lambda_{kk^{\prime}}^{\text{R}}(\alpha_{mn}^{k},s^{\prime})\right)\right)$

with $\Lambda_{kk^{\prime}}^{\text{R(I)}}(\alpha_{mn}^{k},s^{\prime})=\frac{1}{\sqrt{2\pi}}\int_{-\infty}^{+\infty}d\omega_{\xi}\text{e}^{i\omega_{\xi}s^{\prime}}\bar{\varLambda}_{kk^{\prime}}^{\text{R(I)}}(\alpha_{mn}^{k}-\omega_{\xi})$ .

The system evolution time $\tau_{\text{s}}$ is the typical timescale of the intrinsic evolution of the system, which is defined by a typical value for the inverse of the instantaneous frequency differences involved, i.e. $\tau_{\text{s}}\propto|\alpha_{mn}^{k}(t)-\alpha_{m^{\prime}n^{\prime}}^{k^{\prime}}(t)|^{-1}$ for $\alpha_{mn}^{k}(t)\neq\alpha_{m^{\prime}n^{\prime}}^{k^{\prime}}(t)$ . If $\tau_{\text{s}}$ is small compared to the relaxation time $\tau_{\text{R}}$ , the non-secular terms in the DMME with $\alpha_{mn}^{k}(t)\neq\alpha_{m^{\prime}n^{\prime}}^{k^{\prime}}(t)$ may be neglected, which is known as the secular approximation. Thus we have the DMME within the secular approximation in the interaction picture,

	$\displaystyle\partial_{t}\tilde{\rho}_{\text{s}}$		$\displaystyle(t)=-i\left[\tilde{H}_{\text{LS}}(t),\tilde{\rho}_{\text{s}}(t)\right]+\sum_{m,m^{\prime},n,,n^{\prime}}\Gamma_{mn,m^{\prime}n^{\prime}}^{\text{R}}(t)$		(19)
			$\displaystyle\times\left(\tilde{F}_{mn}\tilde{\rho}_{\text{s}}(t)\tilde{F}_{m^{\prime}n^{\prime}}^{\dagger}-\frac{1}{2}\left\{\tilde{F}_{m^{\prime}n^{\prime}}^{\dagger}\tilde{F}_{mn},\tilde{\rho}_{\text{s}}(t)\right\}\right),$		(19)

in which $\tilde{H}_{\text{LS}}(t)=\sum_{m,n,m^{\prime},n^{\prime}}\Gamma_{mn,m^{\prime}n^{\prime}}^{\text{I}}(t)\tilde{F}_{m^{\prime}n^{\prime}}^{\dagger}\tilde{F}_{mn}$ is the Lamb shift Hamiltonian and $\Gamma_{mn,m^{\prime}n^{\prime}}^{\text{R(I)}}(t)$ denotes the real (imaginary) part of $\Gamma_{mn,m^{\prime}n^{\prime}}(t)$ ,

		$\displaystyle\Gamma$	${}_{mn,m^{\prime}n^{\prime}}^{\text{R}}(t)=\frac{1}{\pi}\sum_{k}\textsl{g}_{k}^{2}\xi_{m^{\prime}n^{\prime}}^{k}(t)\int_{-\infty}^{+\infty}ds^{\prime}$
			$\displaystyle\times\left(\Xi_{mn}^{\text{c},k}(t-s^{\prime})\Lambda_{kk}^{\text{R}}(\alpha_{mn}^{k},s^{\prime})-\Xi_{mn}^{\text{s},k}(t-s^{\prime})\Lambda_{kk}^{\text{I}}(\alpha_{mn}^{k},s^{\prime})\right),$

		$\displaystyle\Gamma$	${}_{mn,m^{\prime}n^{\prime}}^{\text{I}}(t)=\frac{1}{2\pi}\sum_{k}\textsl{g}_{k}^{2}\xi_{m^{\prime}n^{\prime}}^{k}(t)\int_{-\infty}^{+\infty}ds^{\prime}$
			$\displaystyle\times\left(\Xi_{mn}^{\text{c},k}(t-s^{\prime})\Lambda_{kk}^{\text{I}}(\alpha_{mn}^{k},s^{\prime})+\Xi_{mn}^{\text{s},k}(t-s^{\prime})\Lambda_{kk}^{\text{R}}(\alpha_{mn}^{k},s^{\prime})\right).$

In fact, although we admit non-adiabatic change of the driving protocol, the DMME presented in Eq.(19) describes a Markovian dynamics if $\Gamma_{mn,m^{\prime}n^{\prime}}^{\text{R}}(t)\geq 0$ for $\forall t$ . The memory effects of the driving protocol is explicitly encoded into a convolution with the reservoir correlation function. In order to connect to the previous resultsDann2018 ; Dann2019 , we may assume that the change of the phase $\theta_{mn}^{k}(t)$ is slow comparing to the reservoir correlation decay rate. Thus, there is a typical timescale $\tau_{\text{d}}$ , called “the non-adiabatic phase timescale”, defined asDann2018

\tau_{\text{d}}\equiv\text{Min}_{m,n,k,t}\left\{\frac{\partial_{t}\theta_{mn}^{k}(t)}{\partial_{t}^{2}\theta_{mn}^{k}(t)}\right\},

which is related to the change of the phase in the driving protocol. Thus, the assumption of the slow changing phase is equivalent to require that the reservoir correlation time $\tau_{\text{B}}$ has to be much smaller than the non-adiabatic phase timescale $\tau_{\text{d}}$ , i.e., $\tau_{\text{B}}\ll\tau_{\text{d}}$ . For $s\in\left[0,\tau_{\text{B}}\right]$ and $s\ll t$ , the terms up to second order in Eq.(13) can be ignored, i.e., $\varTheta_{mn}^{k}=0,$ so that $\Xi_{mn}^{\text{c},k}(t,t-s)=\xi_{mn}^{k}(t-s)$ and $\Xi_{mn}^{\text{s},k}(t,t-s)=0$ . Therefore, the real (imaginary) part of $\Gamma_{mn,m^{\prime}n^{\prime}}(t)$ becomes

$\displaystyle\Gamma_{mn,m^{\prime}n^{\prime}}^{\text{R}}(t)$	$\displaystyle=$	$\displaystyle\frac{1}{\pi}\sum_{k}\textsl{g}_{k}^{2}\xi_{m^{\prime}n^{\prime}}^{k}(t)\int_{-\infty}^{+\infty}d\omega_{\xi}\bar{\xi}_{mn}^{k}(\omega_{\xi})$
		$\displaystyle\times\bar{\varLambda}_{kk}^{\text{R}}(\alpha_{mn}^{k}-\omega_{\xi})\text{e}^{i\omega_{\xi}t}.$
$\displaystyle\Gamma_{mn,m^{\prime}n^{\prime}}^{\text{I}}(t)$	$\displaystyle=$	$\displaystyle\frac{1}{2\pi}\sum_{k}\textsl{g}_{k}^{2}\xi_{m^{\prime}n^{\prime}}^{k}(t)\int_{-\infty}^{+\infty}d\omega_{\xi}\bar{\xi}_{mn}^{k}(\omega_{\xi})$
		$\displaystyle\bar{\varLambda}_{kk}^{\text{I}}(\alpha_{mn}^{k}-\omega_{\xi})\text{e}^{i\omega_{\xi}t}.$

If the change of $\xi_{mn}^{k}(t)$ is much smaller than the instantaneous frequency, i.e., $\omega_{\xi}\ll\alpha_{mn}^{k}$ , we immediately obtain

	$\displaystyle\Gamma_{mn,m^{\prime}n^{\prime}}^{\text{R}}(t)$	$\displaystyle=$	$\displaystyle 2\sum_{k}\textsl{g}_{k}^{2}\xi_{m^{\prime}n^{\prime}}^{k}(t)\xi_{mn}^{k}(t)\bar{\varLambda}_{kk}^{\text{R}}(\alpha_{mn}^{k}),$
	$\displaystyle\Gamma_{mn,m^{\prime}n^{\prime}}^{\text{I}}(t)$	$\displaystyle=$	$\displaystyle\sum_{k}\textsl{g}_{k}^{2}\xi_{m^{\prime}n^{\prime}}^{k}(t)\xi_{mn}^{k}(t)\bar{\varLambda}_{kk}^{\text{I}}(\alpha_{mn}^{k}),$		(22)

which lead to the non-adiabatic Markovian master equation given in Ref.Dann2018 .

The DMME presented in Eq.(19) does not contain any approximation on the driving protocol. It is interesting that both the real and imaginary parts of the one-sided Fourier transform of the instantaneous reservoir correlation function $\bar{\varLambda}_{kk}(\alpha)$ are involved in both the Lamb shift and the decoherence (see Eqs. (II) and (II)). For the Markovian master equation with the static Hamiltonian and the time-dependent Hamiltonian satisfying $\tau_{\text{B}}\ll\tau_{\text{d}}$ , $\bar{\varLambda}_{kk}^{\text{R}}(\alpha)$ only contributes to the decoherence process, while $\bar{\varLambda}_{kk}^{\text{I}}(\alpha)$ just appears in the Lamb shift. As a result, the positive decoherence rates may not be ensured and additional energy level shifts can be observed in the driven open quantum systems for the timescale $\tau_{\text{B}}\sim\tau_{\text{d}}$ .

In the DMME (Eq.(19)), the jump operator $\tilde{F}_{mn}$ denotes a transition from the state $|\psi_{n}(0)\rangle$ to another one $|\psi_{m}(0)\rangle$ . In other words, the transitions caused by the decoherence occur between eigenstates of the LRI. Based on the Lewis-Riesenfeld phase (Eq.(6)), the instantaneous frequency $\alpha_{mn}^{k}$ can be divided into three parts, i.e.,

$\displaystyle\alpha_{mn}^{k}$	$\displaystyle=$	$\displaystyle-\left(\langle\psi_{m}(t)\|H_{\text{s}}(t)\|\psi_{m}(t)\rangle-\langle\psi_{n}(t)\|H_{\text{s}}(t)\|\psi_{n}(t)\rangle\right)$	(23)
		$\displaystyle+i\left(\langle\psi_{m}(t)\|\partial_{t}\|\psi_{m}(t)\rangle-\langle\psi_{n}(t)\|\partial_{t}\|\psi_{n}(t)\rangle\right)$
		$\displaystyle-\partial_{t}\text{Arg}\left(\langle\psi_{m}(t)\|A_{k}\|\psi_{n}(t)\rangle\right).$

The first term in Eq.(23) attributes to a difference between the energy average values of the eigenstates $|\psi_{n}(t)\rangle$ and $|\psi_{m}(t)\rangle$ . The second term is a geometric contribution from the time-dependent eigenstates, while the third term comes from the phase changing rate in the transitions caused by the interaction Hamiltonian. In the adiabatic limit, the eigenstates of the LRI are the eigenstates of the system Hamiltonian, and the adiabatic condition must be satisfied. Thus, the last two terms are no contributions to the instantaneous frequency, while the first term becomes the energy gap between the $n$ -th and the $m$ -th energy levels, which leads to the adiabatic master equation given in Ref. Albash2012 ; Kamleitner2013 .

III The Driven Open Two-Level System

In this section, we apply the general formulism to a driven two-level system which couples with a heat reservoir. Here, we consider that the driven two-level system Hamiltonian in a laser adapted interaction picture takes the form Chen2010

H_{\text{s}}(t)=\Delta(t)\sigma_{z}+\Omega(t)\sigma_{x},

(24)

where $\Delta(t)=\omega_{0}(t)-\omega_{L}$ is the time-dependent detuning with the time-dependent Rabi frequency $\omega_{0}(t)$ and a constant laser frequency $\omega_{L}$ ; $\Omega(t)$ is time-dependent driven field. The heat reservoir can be represented by the reservoir Hamiltonian

H_{\text{B}}=\sum_{k}\Omega_{k}b_{k}^{\dagger}b_{k}

with $\Omega_{k}=\omega_{k}-\omega_{L}$ , where $b_{k}$ and $\omega_{k}$ are the annihilation operator and the eigen-frequency of the $k$ -th mode of the reservoir Shen2014 . Without loss of generality, the interaction Hamiltonian is selected as

H_{\text{I}}=\sum_{j=x,y}A^{j}\otimes B^{j},

(25)

where the system and bath operators are

	$\displaystyle A^{x}=\sigma_{x},\>B^{x}=\sum_{k}g_{k}^{x}(b_{k}^{\dagger}+b_{k}),$
	$\displaystyle A^{y}=\sigma_{y},\>B^{y}=\sum_{k}ig_{k}^{y}(b_{k}-b_{k}^{\dagger}).$		(26)

For the two-level system governed by the Hamiltonian Eq.(24), the LRIs have been explored before Chen2011 ; Lai1996 . Here, we write the LRIs of the two-level system in form of the spectrum decomposition

\displaystyle I_{\text{s}}(t)

\displaystyle=

\displaystyle\sum_{k=1,2}\pm\Omega_{\text{I}}|\psi_{k}(t)\rangle\langle\psi_{k}(t)|,

(27)

where $\pm\Omega_{\text{I}}$ are constant eigenvalues and

	$\displaystyle\|\psi_{1}(t)\rangle$	$\displaystyle=$	$\displaystyle\left(\cos\eta(t)\text{e}^{i\zeta(t)},\sin\eta(t)\right)^{\text{T}},$
	$\displaystyle\|\psi_{2}(t)\rangle$	$\displaystyle=$	$\displaystyle\left(\sin\eta(t)\text{e}^{i\zeta(t)},-\cos\eta(t)\right)^{\text{T}},$		(28)

are the eigenstates of the LRI (Eq.(27)), correspondingly. Inserting Eqs. (24) and (27) into Eq.(4), the parameters $\eta(t)$ and $\zeta(t)$ needs to satisfy the following differential equation

	$\displaystyle\partial_{t}\eta=\Omega\,\sin\zeta,$
	$\displaystyle\sin 2\eta\,\left(2\Delta+\partial_{t}\zeta\right)=2\,\Omega\,\cos 2\eta\,\cos\zeta.$		(29)

In what follows, we identify the system operator $\tilde{A}^{x(y)}(t)$ based on the Eq.(8). On the one hand, $\langle\psi_{m}(t)|A^{x}|\psi_{n}(t)\rangle$ can be obtained by means of Eqs.(26) and (28)

$\displaystyle A_{11}^{x}$	$\displaystyle=$	$\displaystyle\sin 2\eta\,\cos\zeta\,\text{e}^{i\varphi_{11}^{x}},$
$\displaystyle A_{12}^{x}$	$\displaystyle=$	$\displaystyle\sqrt{1-\sin^{2}2\eta\,\cos^{2}\zeta}\,\text{e}^{i\varphi_{12}^{x}},$
$\displaystyle A_{21}^{x}$	$\displaystyle=$	$\displaystyle\sqrt{1-\sin^{2}2\eta\,\cos^{2}\zeta}\,\text{e}^{i\varphi_{21}^{x}},$
$\displaystyle A_{22}^{x}$	$\displaystyle=$	$\displaystyle\sin 2\eta\,\cos\zeta\,\text{e}^{i\varphi_{22}^{x}},$	(30)

so do $\langle\psi_{m}(t)|A^{y}|\psi_{n}(t)\rangle$ , i.e.,

$\displaystyle A_{11}^{y}$	$\displaystyle=$	$\displaystyle\sin 2\eta\,\sin\zeta\,\text{e}^{i\varphi_{11}^{y}},$
$\displaystyle A_{12}^{y}$	$\displaystyle=$	$\displaystyle\sqrt{1-\sin^{2}2\eta\,\sin^{2}\zeta}\,\text{e}^{i\varphi_{12}^{y}},$
$\displaystyle A_{21}^{y}$	$\displaystyle=$	$\displaystyle\sqrt{1-\sin^{2}2\eta\,\sin^{2}\zeta}\,\text{e}^{i\varphi_{21}^{y}},$
$\displaystyle A_{22}^{y}$	$\displaystyle=$	$\displaystyle\sin 2\eta\,\sin\zeta\,\text{e}^{i\varphi_{22}^{y}},$	(31)

in which the phases are $\varphi_{11}^{x}=0$ , $\varphi_{22}^{x}=\pi$ , $\varphi_{11}^{y}=\pi$ , $\varphi_{22}^{y}=0$ ,

	$\displaystyle\tan\varphi_{12}^{x}$	$\displaystyle=$	$\displaystyle-\frac{\sin\zeta}{\cos 2\eta\,\cos\zeta},\,\tan\varphi_{21}^{x}=\frac{\sin\zeta}{\cos 2\eta\,\cos\zeta},$
	$\displaystyle\tan\varphi_{12}^{y}$	$\displaystyle=$	$\displaystyle\frac{\cos\zeta}{\cos 2\eta\,\sin\zeta},\,\tan\varphi_{21}^{y}=-\frac{\cos\zeta}{\cos 2\eta\,\sin\zeta},$

correspondingly. After substituting Eq.(28) into Eq.(6), we can obtain the Lewis-Riesenfeld phases,

	$\displaystyle\alpha_{1}$	$\displaystyle=$	$\displaystyle\int_{0}^{t}\text{d}\tau\left(-\partial_{\tau}\zeta\cos^{2}\eta-{\Delta}\cos 2\eta-\Omega\cos\zeta\sin 2\eta\right),$
	$\displaystyle\alpha_{2}$	$\displaystyle=$	$\displaystyle\int_{0}^{t}\text{d}\tau\left(-\partial_{\tau}\zeta\sin^{2}\eta+{\Delta}\cos 2\eta+\Omega\cos\zeta\sin 2\eta\right).$

Thus, the propagator of the free dynamics for the driven two-level system with the system Hamiltonian Eq.(24) can be written down explicitly according to Eq.(7). From Eqs. (30) and (31), we have

$\displaystyle\xi_{11}^{x}$	$\displaystyle=$	$\displaystyle\xi_{22}^{x}=\sin 2\eta\,\cos\zeta,$
$\displaystyle\xi_{12}^{x}$	$\displaystyle=$	$\displaystyle\xi_{21}^{x}=\sqrt{1-\sin^{2}2\eta\,\cos^{2}\zeta},$
$\displaystyle\xi_{11}^{y}$	$\displaystyle=$	$\displaystyle\xi_{22}^{y}=\sin 2\eta\,\sin\zeta,$
$\displaystyle\xi_{12}^{y}$	$\displaystyle=$	$\displaystyle\xi_{21}^{y}=\sqrt{1-\sin^{2}2\eta\,\sin^{2}\zeta}.$	(32)

and

$\displaystyle\theta_{12}^{j}$	$\displaystyle=$	$\displaystyle\alpha_{2}-\alpha_{1}+{\varphi}_{12}^{j},$
$\displaystyle\theta_{21}^{j}$	$\displaystyle=$	$\displaystyle\alpha_{1}-\alpha_{2}+{\varphi}_{21}^{j},$
$\displaystyle\theta_{11}^{j}$	$\displaystyle=$	$\displaystyle\varphi_{11}^{j},\,\theta_{22}^{j}=\varphi_{22}^{j},$

for $j=x,y$ , which result in the instantaneous frequencies as

$\displaystyle\alpha_{12}^{x}$	$\displaystyle=$	$\displaystyle-\alpha_{21}^{x}$
	$\displaystyle=$	$\displaystyle-\partial_{\tau}\zeta\,\cos 2\eta-2\Delta\cos 2\eta-2\Omega\cos\zeta\sin 2\eta$
		$\displaystyle+\frac{\partial_{t}\eta\,\sin 2\eta\sin 2\zeta+\partial_{t}\zeta\,\cos 2\eta}{1-{\sin^{2}2\,\eta}\,{\cos^{2}\zeta}},$
$\displaystyle\alpha_{12}^{y}$	$\displaystyle=$	$\displaystyle-\alpha_{21}^{y}$
	$\displaystyle=$	$\displaystyle-\partial_{\tau}\zeta\,\cos 2\eta-2\Delta\cos 2\eta-2\Omega\cos\zeta\sin 2\eta$
		$\displaystyle-\frac{\partial_{t}\eta\,\sin 2\eta\,\sin 2\zeta-\partial_{t}\zeta\,\cos 2\eta}{1-\sin^{2}2\eta\,\sin^{2}\zeta},$
$\displaystyle\alpha_{11}^{x(y)}$	$\displaystyle=$	$\displaystyle\alpha_{22}^{x(y)}=0.$	(33)

Therefore, the system operators $\tilde{A}^{x(y)}(t)$ are determined by taking $\xi_{mn}^{j}$ , $\theta_{mn}^{j}$ and $\alpha_{mn}^{j}$ into Eq.(8).

Based on the parameters provided above, we can obtain the Lamb shifts and the decoherence rates via Eqs. (II) and (II). Firstly, according to Eq.(13), we have

\displaystyle\varTheta_{mn}^{x(y)}(t,t-s)=\theta_{mn}^{x(y)}(t-s)-\theta_{mn}^{x(y)}(t)-\alpha_{mn}^{x(y)}(t)s,

which yields

\displaystyle\varTheta_{12}^{x(y)}(t,t-s)

\displaystyle=\int_{t-s}^{t}\text{d}\tau\,\left(\alpha_{12}^{x(y)}(\tau)-\alpha_{12}^{x(y)}(t)\right)

and $\varTheta_{11}^{x(y)}(t,t-s)=\varTheta_{22}^{x(y)}(t,t-s)=0$ . Secondly, let us take the reservoir to be in an equilibrium state at temperature $T_{R}$ . The correlation functions of the heat reservoir operators satisfy

$\displaystyle\text{Tr}_{\text{B}}\left\{b_{k^{\prime}}b_{k}^{\dagger}\rho_{B}\right\}$	$\displaystyle=$	$\displaystyle\delta_{k^{\prime}k}(1+N_{k}),$
$\displaystyle\text{Tr}_{\text{B}}\left\{b_{k^{\prime}}^{\dagger}b_{k}\rho_{B}\right\}$	$\displaystyle=$	$\displaystyle\delta_{k^{\prime}k}N_{k},$
$\displaystyle\text{Tr}_{\text{B}}\left\{b_{k^{\prime}}b_{k}\rho_{B}\right\}$	$\displaystyle=$	$\displaystyle 0,$
$\displaystyle\text{Tr}_{\text{B}}\left\{b_{k^{\prime}}^{\dagger}b_{k}^{\dagger}\rho_{B}\right\}$	$\displaystyle=$	$\displaystyle 0,$	(34)

where $N_{k}=\left(\exp(\omega_{k}/T_{R})-1\right)^{-1}$ denotes the Planck distribution with the reservoir temperature $T_{R}$ . In continuum limit, the sum over $(\textsl{g}_{k}^{x(y)})^{2}$ can be replaced by an integral

\displaystyle\sum_{k}\left(\textsl{g}_{k}^{x(y)}\right)^{2}\rightarrow\int_{0}^{\infty}\text{d}\omega_{k}J^{x(y)}(\omega_{k})

(35)

with the spectral density function $J^{x(y)}(\omega_{k})$ . Inserting Eq.(26) into Eq.(17), it yields

$\displaystyle\bar{\Lambda}^{x(y)}(\alpha)$	$\displaystyle\equiv$	$\displaystyle\sum_{k,k^{\prime}}g_{k}^{x(y)}g_{k^{\prime}}^{x(y)}\bar{\Lambda}_{kk^{\prime}}^{x(y)}(\alpha)$	(36)
	$\displaystyle=$	$\displaystyle\int_{0}^{\infty}\text{d}\Omega_{k}J^{x(y)}(\Omega_{k}+\omega_{L})\left(N_{k}\int_{0}^{\infty}\text{d}s\text{e}^{i\left(\alpha+\Omega_{k}\right)s}\right.$
		$\displaystyle\left.+(N_{k}+1)\int_{0}^{\infty}\text{d}s\text{e}^{i\left(\alpha-\Omega_{k}\right)s}\right).$

with $\alpha=\alpha_{mn}^{x(y)}-\omega_{\xi}$ . On making use of the formula

\displaystyle\int_{0}^{\infty}\text{d}s\text{e}^{-i\varepsilon s}=\pi\delta(\varepsilon)-i\text{P}\frac{1}{\varepsilon}

(37)

with the Cauchy principal value P, we finally arrive at

\bar{\Lambda}^{x(y)}(\alpha)=\bar{\Lambda}^{\text{R},x(y)}(\alpha)+i\bar{\Lambda}^{\text{I},x(y)}(\alpha),

where

\bar{\Lambda}^{\text{R},x(y)}(\alpha)=\gamma_{0}(\alpha)\left(N(\alpha+\omega_{L})+1\right)

and

			$\displaystyle\bar{\Lambda}^{\text{I},x(y)}(\alpha)$
	$\displaystyle=$		$\displaystyle\text{P}\left[\int_{0}^{\infty}\text{d}\omega_{k}J^{x(y)}(\omega_{k})\left[\frac{N(\omega_{k})+1}{\alpha+\omega_{L}-\omega_{k}}+\frac{N(\omega_{k})}{\alpha-\omega_{L}+\omega_{k}}\right]\right].$

with $\gamma_{0}(\alpha)=\pi J(\alpha)$ . After inserting $\bar{\Xi}_{mn}^{\text{c(s)},x(y)}(t,\omega_{\xi})$ and $\bar{\Lambda}^{\text{R(I)},x(y)}(\alpha)$ into Eq.(18) and taking the inverse Fourier transformation respect to $\omega_{\xi}$ , the Lamb shifts and the decoherence rates can be obtained.

Without any restriction on the driving protocol, the dynamics of the driven two-level system is governed by the following DMME in the interaction picture,

\displaystyle\mathcal{\tilde{L}}\tilde{\rho}_{\text{s}}(t)=-i\left[\tilde{H}_{\text{LS}}(t),\tilde{\rho}_{\text{s}}(t)\right]+\mathcal{D}^{\text{R}}\tilde{\rho}_{\text{s}}(t)+\mathcal{D}^{\text{D}}\tilde{\rho}_{\text{s}}(t),

(38)

with the Lamb shifts $\tilde{H}_{\text{LS}}(t)=\sum_{j,mn}\Gamma^{\text{I},j}_{mn}(t)\tilde{F}_{mn}^{\dagger}\tilde{F}_{mn}$ . According to Eq.(33), the instantaneous frequency is degenerate for $mn=\{11,\,22\}$ , which indicates a dephasing process on $|\varphi_{1}\rangle$ and $|\varphi_{2}\rangle$ . Therefore, we divide the Lindblandian into two parts with the dissipators

$\displaystyle\mathcal{D}^{\text{R }}\tilde{\rho}_{\text{s}}=$	$\displaystyle\sum_{mn=12,21}\sum_{j=x,y}\Gamma^{\text{R},j}_{mn,mn}(t)$	(39)
	$\displaystyle\times\left(\tilde{F}_{mn}\tilde{\rho}_{\text{s}}\tilde{F}_{mn}^{\dagger}-\frac{1}{2}\left\{\tilde{F}_{mn}^{\dagger}\tilde{F}_{mn},\tilde{\rho}_{\text{s}}\right\}\right),$	(39)
$\displaystyle\mathcal{D}^{\text{D}}\tilde{\rho}_{\text{s}}=$	$\displaystyle\sum_{mn=11,22}^{m^{\prime}n^{\prime}=11,22}\sum_{j=x,y}\Gamma^{\text{R},j}_{mn,m^{\prime}n^{\prime}}(t)\text{e}^{i\left(\theta_{mn}^{j}(t)-\theta_{m^{\prime}n^{\prime}}^{j}(t)\right)}$	(40)
	$\displaystyle\times\left(\tilde{F}_{m^{\prime}n^{\prime}}\tilde{\rho}_{\text{s}}\tilde{F}_{mn}^{\dagger}-\frac{1}{2}\left\{\tilde{F}_{m^{\prime}n^{\prime}}^{\dagger}\tilde{F}_{mn},\tilde{\rho}_{\text{s}}\right\}\right),$	(40)

which correspond to the energy dissipation and the dephasing processes, respectively. Here, we have used the fact $\alpha_{12}(t)=-\alpha_{21}(t)$ , so that the terms with $mn\neq m^{\prime}n^{\prime}$ in Eq.(39) vanish because of the secular approximation. It is noteworthy that the dephasing rates in Eq.(40) satisfy $\Gamma^{\text{R},x(y)}_{11,11}=\Gamma^{\text{R},x(y)}_{22,22}=-\Gamma^{\text{R},x(y)}_{11,22}=-\Gamma^{\text{R},x(y)}_{22,11}$ , due to $\theta_{11}^{x}=\theta_{22}^{y}=\pi$ , $\theta_{22}^{x}=\theta_{11}^{y}=0$ . By introducing a Hermitian operator of the interaction picture $\tilde{\Sigma}_{z}=\tilde{F}_{22}-\tilde{F}_{11}$ , the dephasing term in Eq.(38) can be rewritten as

\displaystyle\mathcal{D}^{\text{D}}\tilde{\rho}_{\text{s}}=\Gamma^{\text{R}}_{d}(t)\left(\tilde{\Sigma}_{z}\tilde{\rho}_{\text{s}}\tilde{\Sigma}_{z}^{\dagger}-\frac{1}{2}\left\{\tilde{\Sigma}_{z}^{\dagger}\tilde{\Sigma}_{z},\tilde{\rho}_{\text{s}}\right\}\right),

with $\Gamma^{\text{R}}_{d}(t)=\Gamma^{\text{R},x}_{11,11}+\Gamma^{\text{R},y}_{11,11}$ . We further define two operators $\tilde{\Sigma}_{+}\equiv\tilde{F}_{21}$ and $\tilde{\Sigma}_{-}\equiv\tilde{F}_{12}$ , which fulfills

\tilde{\Sigma}_{+}=\tilde{\Sigma}_{-}^{\dagger},\,\left[\tilde{\Sigma}_{z},\tilde{\Sigma}_{+}\right]=\frac{1}{2}\tilde{\Sigma}_{+},\,\left[\tilde{\Sigma}_{z},\tilde{\Sigma}_{-}\right]=-\frac{1}{2}\tilde{\Sigma}_{-}.

Thus, the dissipative term as shown Eq.(39) can be reproduced as

	$\displaystyle\mathcal{D}^{\text{R }}\tilde{\rho}_{\text{s}}$	$\displaystyle=$	$\displaystyle\Gamma^{\text{R}}_{+}(t)\left(\tilde{\Sigma}_{+}\tilde{\rho}_{\text{s}}\tilde{\Sigma}_{-}-\frac{1}{2}\left\{\tilde{\Sigma}_{-}\tilde{\Sigma}_{+},\tilde{\rho}_{\text{s}}\right\}\right)$
		$\displaystyle+$	$\displaystyle\Gamma^{\text{R}}_{-}(t)\left(\tilde{\Sigma}_{-}\tilde{\rho}_{\text{s}}\tilde{\Sigma}_{+}-\frac{1}{2}\left\{\tilde{\Sigma}_{+}\tilde{\Sigma}_{-},\tilde{\rho}_{\text{s}}\right\}\right),$

with $\Gamma^{\text{R}}_{+}\equiv\Gamma^{\text{R},x}_{21}+\Gamma^{\text{R},y}_{21}$ and $\Gamma^{\text{R}}_{-}\equiv\Gamma^{\text{R},x}_{12}+\Gamma^{\text{R},y}_{12}$ . Transforming back to the Schrödinger picture, we finally arrive at the DMME,

$\displaystyle\partial_{t}\rho_{\text{s}}$	$\displaystyle=$	$\displaystyle\mathcal{L}(t)\rho_{\text{s}}$	(41)
	$\displaystyle=$	$\displaystyle-{i}\left[H_{\text{s}}(t)+H_{\text{LS}}(t),\rho_{\text{s}}(t)\right]$
	$\displaystyle+$	$\displaystyle\Gamma^{\text{R}}_{+}(t)\left(\Sigma_{+}{\rho}_{\text{s}}(t)\Sigma_{-}-\frac{1}{2}\left\{\Sigma_{-}\Sigma_{+},{\rho}_{\text{s}}(t)\right\}\right)$
	$\displaystyle+$	$\displaystyle\Gamma^{\text{R}}_{-}(t)\left(\Sigma_{-}{\rho}_{\text{s}}(t)\Sigma_{+}-\frac{1}{2}\left\{\Sigma_{+}\Sigma_{-},{\rho}_{\text{s}}(t)\right\}\right)$
	$\displaystyle+$	$\displaystyle\Gamma^{\text{R}}_{d}(t)\left[\Sigma_{z},\left[{\rho}_{\text{s}}(t),\Sigma_{z}\right]\right].$

with the time-dependent Lindblad operators $\Sigma_{k}=U_{\text{s}}(t)\tilde{\Sigma}_{k}U_{\text{s}}^{\dagger}(t)$ for $k=+,-,z$ , and the Lamb shift $H_{\text{LS}}(t)=U_{\text{s}}(t)\tilde{H}_{\text{LS}}(t)U_{\text{s}}^{\dagger}(t)$ .

III.1 The Adiabatic Limit

In the adiabatic limit, the corresponding LRIs satisfy $\left[H_{\text{s}}(t),I_{\text{s}}(t)\right]=0$ , and share the same eigenstates to the system Hamiltonian. According to Eq.(29), if $\partial_{t}\eta(t)=\partial_{t}\zeta(t)=0$ , it yields $\sin\zeta=0$ and $\tan 2\eta=\Omega/\Delta$ . Thus, we can write down the eigenstates of the system Hamiltonian (Eq.(24)) in form of Eq.(28) with

\displaystyle\zeta=0,\,\eta=\arccos\left(-\frac{\sqrt{2}}{2}\sqrt{\frac{\sqrt{\Delta^{2}+\Omega^{2}}-\Delta}{\sqrt{\Delta^{2}+\Omega^{2}}}}\right).

(42)

It can be verified that

H_{\text{s}}(t)|\psi_{i}(t)\rangle=\epsilon_{i}(t)|\psi_{i}(t)\rangle

with the eigenvalues of the system Hamiltonian $\epsilon_{1,2}(t)=\mp\sqrt{\Delta^{2}+\Omega^{2}}/2$ . In such a case, the propagator can be represented in terms of the instantaneous eigenstates of the system Hamiltonian as shown in Eq.(7). The phases in the propagator come back to a sum of the geometric phases and the dynamical phases.

Here, we consider the situation where $\textsl{g}_{k}^{x}=0$ in the interaction Hamiltonian Eq.(25) for all $k$ . Thus, the expansion coefficients in Eq.(8) are

	$\displaystyle\xi_{11}^{y}$	$\displaystyle=$	$\displaystyle\xi_{22}^{y}=\|\sin 2\eta\,\sin\zeta\|=0,$
	$\displaystyle\xi_{12}^{y}$	$\displaystyle=$	$\displaystyle\xi_{21}^{y}=\sqrt{1-\sin^{2}2\eta\,\sin^{2}\zeta}=1,$

with the phase $\varphi_{12}^{y}=-\varphi_{21}^{y}=\pi/2$ . Due to $\partial_{t}\zeta\rightarrow 0$ , the geometric phases in $\alpha_{1}$ and $\alpha_{2}$ are much smaller than the coresponding dynamical phases, so that the phase in Eq.(8) reads

	$\displaystyle\theta_{12}^{y}$	$\displaystyle=$	$\displaystyle\alpha_{2}-\alpha_{1}$
		$\displaystyle=$	$\displaystyle-\int_{0}^{t}\text{d}\tau\sqrt{\Delta(\tau)^{2}+\Omega(\tau)^{2}}+\pi,$

and $\theta_{21}^{y}=-\theta_{12}^{y}$ , whose derivatives are

	$\displaystyle\alpha_{12}^{y}$	$\displaystyle=$	$\displaystyle\sqrt{\Delta(t)^{2}+\Omega(t)^{2}},$
	$\displaystyle\alpha_{21}^{y}$	$\displaystyle=-$	$\displaystyle\sqrt{\Delta(t)^{2}+\Omega(t)^{2}},$

respectively.

In the adiabatic limits, the reservoir correlation time $\tau_{\text{B}}$ is much smaller than the non-adiabatic timescale of the driving protocol $\tau_{\text{d}}$ , i.e., $\tau_{\text{B}}\ll\tau_{\text{d}}$ Dann2021 . Thus the Lamb shifts and the decoherence rates can be obtained from Eq. (22). By considering Eq.(36), it yields

\displaystyle\Gamma_{mn}^{\text{R},y}(t)

\displaystyle=

\displaystyle 2\gamma_{0}(\alpha_{mn})\left(N(\alpha_{mn})+1\right),

(43)

with $mn=12,\,21$ . Note that No matter $\Delta(t)$ and $\Omega(t)$ are either positive or negative, $\alpha_{12}^{y}$ ( $\alpha_{21}^{y}$ ) is always positive (negative), and the Planck distribution satisfies $N(-\alpha_{mn})=-\left(N(\alpha_{mn})+1\right)$ . Therefore, the adiabatic Markovian master equation (AME) for the driven open two-level system can be written as Albash2012 ; Kamleitner2013

$\displaystyle\partial_{t}\rho_{\text{s}}$	$\displaystyle=$	$\displaystyle\mathcal{L}(t)\rho_{\text{s}}$	(44)
	$\displaystyle=$	$\displaystyle-{i}\left[H_{\text{s}}(t)+H_{\text{LS}}(t),\rho_{\text{s}}(t)\right]$
	$\displaystyle+$	$\displaystyle 2\gamma_{0}(N+1)\left(\Sigma_{-}{\rho}_{\text{s}}(t)\Sigma_{+}-\frac{1}{2}\left\{\Sigma_{+}\Sigma_{-},{\rho}_{\text{s}}(t)\right\}\right)$
	$\displaystyle+$	$\displaystyle 2\gamma_{0}N\left(\Sigma_{+}{\rho}_{\text{s}}(t)\Sigma_{-}-\frac{1}{2}\left\{\Sigma_{-}\Sigma_{+},{\rho}_{\text{s}}(t)\right\}\right).$

III.2 The Inertial Limit

A NAME based on the inertial theorem has been proposedDann2018 ; Dann2020 ; Dann2021 , in which the free propagator of the closed quantum system is determined by decomposing the dynamical generator in the Hilbert-Schmidt space into a rapidly changed scalar function and an adiabatically changed matrix. In this subsection, we illustrate that the NAME based on the inertial theorem is the DMME as shown in Eq.(41) in the inertial limits.

Besides the system Hamiltonian given by Eq.(24), two following additional operators are needed to determine the free propagator, which areDann2021

\displaystyle L(t)=\Omega(t)\sigma_{z}-\Delta(t)\sigma_{x},\>C(t)=\bar{\Omega}(t)\sigma_{z}

(45)

with $\bar{\Omega}(t)=\sqrt{\Omega^{2}(t)+\Delta^{2}(t)}$ . We may construct the Liouvillian vector as $\vec{v}=\{H_{\text{s}}(t),L(t),C(t)\}$ . The inertial theorem requires that the adiabatic parameters for $H_{\text{s}}(t)$ is constant Dann2021 , i.e.,

\displaystyle\mu=\frac{\Omega(t)\partial_{t}\Delta(t)-\Delta(t)\partial_{t}\Omega(t)}{2\,\bar{\Omega}^{3}(t)}\equiv\text{cons.}.

Here, we call the dynamics which satisfies the requirement of the inertial theorem as the dynamics in the inertial limit. Under the inertial limit mentioned above, we have $\partial_{t}\vec{w}(t)=-i\bar{\Omega}(t)\mathcal{B}(\mu)\vec{w}(t)$ with $\vec{w}(t)=\frac{\bar{\Omega}(0)}{\bar{\Omega}(t)}\vec{v}(t)$ and

\mathcal{B}(\mu)=i\left(\begin{array}[]{ccc}0&\mu&0\\ -\mu&0&1\\ 0&-1&0\\ \end{array}\right).

As a result, by calculating the eigenstates of $\mathcal{B}(\mu)$ , we have the eigenoperators of the free propagator

$\displaystyle\Sigma_{x}$	$\displaystyle=$	$\displaystyle\frac{1}{2\kappa^{2}\bar{\Omega}(t)}[-\mu H_{\text{s}}(t)-i\kappa L(t)+C(t)],$
$\displaystyle\Sigma_{y}$	$\displaystyle=$	$\displaystyle\frac{1}{2\kappa^{2}\bar{\Omega}(t)}[-\mu H_{\text{s}}(t)+i\kappa L(t)+C(t)],$
$\displaystyle\Sigma_{z}$	$\displaystyle=$	$\displaystyle\frac{1}{\kappa\bar{\Omega}(t)}[H_{\text{s}}(t)+\mu C(t)],$	(46)

with $\kappa=\sqrt{1+\mu^{2}}$ , which are the Lindblad operators in the NAME based on the inertial theoremDann2020 .

As discussed in Sec.II, the Lindblad operators can be determined by the eigenstates of the LRIs according to Eq.(8). For the driven open quantum systems with the system Hamiltonian Eq.(24), the eigenstates of $\Sigma_{z}$ (Eq.(46)) must be the eigenstates of the LRIs defined in Eq.(27), since $\Sigma_{z}$ is Hermitian. In order to verify this correspondence, we check whether the eigenstates of $\Sigma_{z}$ fulfill the differential equation for the parameters of the LRIs as shown in Eq.(29). Substituting Eqs.(24) and (45) into Eq.(46), the eigenstates of $\Sigma_{z}$ are obtained straightforwardly,

	$\displaystyle\|\varphi_{1}\rangle=\left(\begin{array}[]{c}\frac{\left(i\,\mu\bar{\Omega}-\Omega\right)}{\sqrt{2\,\kappa\bar{\Omega}\left(\Delta+\kappa\bar{\Omega}\right)}}\\ \frac{\sqrt{\Delta+\kappa\bar{\Omega}}}{\sqrt{2\,\kappa\bar{\Omega}}}\end{array}\right),$		(49)
	$\displaystyle\|\varphi_{2}\rangle=\left(\begin{array}[]{c}\frac{\left(\Omega-i\,\mu\bar{\Omega}\right)}{\,\sqrt{2\,\kappa\bar{\Omega}\left(\kappa\bar{\Omega}-\Delta\right)}}\\ \frac{\sqrt{\kappa\bar{\Omega}-\Delta}}{\sqrt{2\,\kappa\bar{\Omega}}}\end{array}\right),$		(52)

which can be parameterized as Eq.(28) with

\displaystyle\zeta=-\arctan\left(\frac{\mu\bar{\Omega}}{\Omega}\right),\,\eta=\arccos\left(-\frac{\sqrt{2}}{2}\sqrt{\frac{\kappa\bar{\Omega}-\Delta}{\kappa\bar{\Omega}}}\right).

The time derivatives of $\eta(t)$ and $\zeta(t)$ read

	$\displaystyle\partial_{t}\zeta$	$\displaystyle=$	$\displaystyle-\frac{2\mu^{2}\bar{\Omega}^{2}\Delta}{\mu^{2}\bar{\Omega}^{2}+\Omega^{2}}-\frac{\Omega\bar{\Omega}}{\mu^{2}\bar{\Omega}^{2}+\Omega^{2}}\partial_{t}\mu,$
	$\displaystyle\partial_{t}\eta$	$\displaystyle=$	$\displaystyle-\frac{\mu\Omega\bar{\Omega}}{\sqrt{\mu^{2}\bar{\Omega}^{2}+\Omega^{2}}}+\frac{\mu\Delta}{2\kappa^{2}\sqrt{\mu^{2}\bar{\Omega}^{2}+\Omega^{2}}}\,\partial_{t}\mu.$

By taking $\zeta$ , $\eta$ and their time-derivatives into Eq.(29), it can be verify that, the differential equations Eq.(29) hold in the inertial limits, i.e., $\partial_{t}\mu/\sqrt{\mu^{2}\bar{\Omega}^{2}+\Omega^{2}}\rightarrow 0$ . In other words, $|\varphi_{1}\rangle$ and $|\varphi_{2}\rangle$ are the eigenstates of a inertial LRI which requires $\partial_{t}\mu=0$ .

In the following, we derive the inertial Markovian master equation according to the inertial LRI. Here, we still consider that $\textsl{g}_{k}^{x}=0$ for all $k$ in the interaction Hamiltonian $H_{\text{I}}$ . By inserting Eq.(52) into Eq.(10), it yields the expanding coefficients

\displaystyle\xi_{11}=\xi_{22}=\frac{\mu}{\kappa},\,\xi_{12}=\xi_{21}=\frac{1}{\kappa},

(53)

and the phases

	$\displaystyle\theta_{12}$	$\displaystyle=$	$\displaystyle-\theta_{21}$
		$\displaystyle=$	$\displaystyle-\int_{0}^{t}\text{d}\tau\,\frac{2\kappa\bar{\Omega}(\tau)\Omega^{2}(\tau)}{\mu^{2}\bar{\Omega}^{2}(\tau)+\Omega^{2}(\tau)}+\varphi_{12}(t),$

with $\varphi_{12}=\arctan(\kappa\Omega/\mu\Delta)$ . After some simple algebra, we obtain the instantaneous frequency with a concise representation

\displaystyle\alpha_{12}=-\alpha_{21}=2\kappa\bar{\Omega}(t).

If we assume that the non-adiabatic timescale $\tau_{\text{d}}$ is great shorter than the reservoir correlation time $\tau_{\text{B}}$ , we can determine the Lamb shifts and the decoherence rates according to Eq.(22). By the same procedure in Sec.III.1, we obtain the inertial master equation , which has precisely same representation as shown in Ref. Dann2020 .

III.3 Comparison to the Exactly Solvable Models: the Dephasing Model

Here, we further compare to another toy model, which is exactly solved in the interaction picture. With the same system Hamiltonian as shown in Eq.(24), i.e.,

H_{\text{s}}(t)=\Delta(t)\sigma_{z}+\Omega(t)\sigma_{x},

we consider a time-dependent interaction Hamiltonian

\displaystyle H_{\text{I}}=A(t)\otimes B,

where the system and reservoir operators are

	$\displaystyle A(t)$	$\displaystyle=$	$\displaystyle\sin 2\eta\cos\zeta\sigma_{x}+\sin 2\eta\sin\zeta\sigma_{y}+\cos 2\eta\sigma_{z},$		(54)
	$\displaystyle B$	$\displaystyle=$	$\displaystyle\sum_{k}g_{k}\left(b_{k}^{\dagger}+b_{k}\right).$

$\eta$ and $\zeta$ are time-dependent parameters in the eigenstates of the LRI (Eq.(28)), which are governed by Eq.(29). In this way, the Hamiltonian in the interaction picture can be written as

\tilde{H}_{\text{I}}=\sigma_{z}\otimes\sum_{k}g_{k}\left(b_{k}^{\dagger}\text{e}^{i\Omega_{k}t}+b_{k}\text{e}^{-i\Omega_{k}t}\right).

By using an unitary transformation

\tilde{V}=\exp\left[\frac{1}{2}\sigma_{z}\sum_{k}\left(\gamma_{k}b_{k}^{\dagger}-\gamma_{k}^{*}b_{k}\right)\right]

(55)

with $\gamma_{k}=2g_{k}\left(1-\text{e}^{i\Omega_{k}t}\right)/\Omega_{k}$ , the reservoir and two-level system decouple to each other, which leads to an exactly solvable dynamics of the open two-level systemBreuer2007 ; Uhrig2007 ; Yi2000 . The detailed derivation can be found in Appendix A.

Putting aside the exact dynamics of this toy model, we derive the DMME for the open two-level system. Taking Eqs. (28) and (54) into Eq.(8), it follows that the amplitudes are $\xi_{11}=\xi_{22}=1$ and $\xi_{12}=\xi_{21}=0$ while the phases and instantaneous frequency read $\theta_{11}=\pi$ , $\theta_{22}=0$ , and $\alpha_{11}=\alpha_{22}=0$ . Thus we have $\Gamma_{12,12}=\Gamma_{21,21}=0$ and $\Gamma_{11,11}=\Gamma_{22,22}=-\Gamma_{11,22}=-\Gamma_{22,11}\equiv\Gamma_{\text{D}}$ according to Eq.(15). In view of the reservoir correlation functions Eq.(34), $\Gamma_{\text{D}}$ may be taken the following form

	$\displaystyle\Gamma_{\text{D}}$	$\displaystyle=$	$\displaystyle\lim_{t\rightarrow\infty}\int_{0}^{\infty}d\Omega_{k}J\left(\Omega_{k}\right)$
			$\displaystyle\times\int_{0}^{t}ds\left[\left(2N_{k}+1\right)\cos\Omega_{k}s-i\sin\Omega_{k}s\right].$

In case of the zero reservoir temperature, i.e., $N_{k}=0$ , it yields

	$\displaystyle\Gamma_{\text{D}}$	$\displaystyle=$	$\displaystyle\lim_{t\rightarrow\infty}\left\{\frac{\kappa\,\Omega_{c}^{2}t}{\Omega_{c}^{2}t^{2}+1}+i\!\frac{\kappa\,\Omega_{c}^{3}t^{2}}{\Omega_{c}^{2}t^{2}+1}\right\}$
		$\displaystyle=$	$\displaystyle i\,\kappa\Omega_{c},$

in which the following spectral density has been used Breuer2007

J(\Omega_{k})=\kappa\Omega_{k}\exp\left(-\frac{\Omega_{k}}{\Omega_{c}}\right)

(56)

with the cut-off frequency $\Omega_{c}$ and a dimensionless coupling rate $\kappa$ . As we see, under the Markovian approximation, the real part of $\Gamma_{\text{D}}$ vanishes, which leads to a meaningless DMME. Hence, we restrict the upper limit of the integration over $s$ to be $t$ , but not $\infty$ . The DMME in the interaction picture can be written as

\partial_{t}\tilde{\rho}_{\text{s}}(t)=-i\left[\tilde{H}_{\text{LS}}(t),\tilde{\rho}_{\text{s}}(t)\right]+\mathcal{D}^{\text{D}}\tilde{\rho}_{\text{s}}(t)

with the Lamb shift Hamiltonian $\tilde{H}_{\text{LS}}(t)=\Gamma_{\text{D}}^{\text{I}}(t)\tilde{\Sigma}_{z}^{\dagger}\tilde{\Sigma}_{z}$ and the dissipator

\mathcal{D}^{\text{D}}\tilde{\rho}_{\text{s}}=\Gamma_{\text{D}}^{\text{R}}(t)\left(\tilde{\Sigma}_{z}\tilde{\rho}_{\text{s}}\tilde{\Sigma}_{z}^{\dagger}-\frac{1}{2}\left\{\tilde{\Sigma}_{z}^{\dagger}\tilde{\Sigma}_{z},\tilde{\rho}_{\text{s}}\right\}\right).

The Lamb shift strength and the dephasing rate are

\Gamma_{\text{D}}^{\text{I}}(t)=\frac{\kappa\,\Omega_{c}^{3}t^{2}}{\Omega_{c}^{2}t^{2}+1},\>\Gamma_{\text{D}}^{\text{R}}(t)=\frac{\kappa\,\Omega_{c}^{2}t}{\Omega_{c}^{2}t^{2}+1},

while the Lindblad operator is $\tilde{\Sigma}_{z}=\tilde{F}_{22}-\tilde{F}_{11}$ .

Refer to caption — Figure 1: The main values of the Pauli operators as a function of the dimensionless time $\Omega_{c}t$ for the dynamics given by the exact solution (red dashed lines), the DMME without the Lamb shift (blue solid lines), and the DMME with the Lamb shift (green dotted lines). The parameters are chosen as $\Omega_{0}=\Delta=\omega_{c}$ , $\Omega_{c}=20\omega_{c}$ , and $\kappa=1$ . We set $\omega_{c}=1$ as an unity of the other parameters.

The numerical results of the main values of the Pauli operators are plotted in FIG.(1). The initial state of the system is prepared on $\rho_{\text{s}}(0)=(|1\rangle+|0\rangle)(\langle 1|+\langle 0|)/2$ . The red dashed lines are the results given by the exact solution, while the blue solid lines are the results associated with the DMME.We consider a driving protocol with a constant detuning $\Delta$ and a driving field with a time-dependent strength

\displaystyle\Omega=\Omega_{0}\sin^{2}(\omega_{c}t)

with constant $\Omega_{0}$ and $\omega_{c}$ . In FIG.(1), we observe a strikingly good agreement between the DMME (blue solid lines) and the exact solution ( red dashed lines). In fact, the high uniformity does not depend on the choice of the driving rate and the driving protocol. It is easy to verify

\Gamma_{\text{e}}(t)=\int_{0}^{t}d\tau\Gamma_{\text{D}}^{\text{R}}(\tau),

which results in the same deocherence process for both the DMME and the exact solution. On the other hand, due to $\tilde{\Sigma}_{z}^{\dagger}\tilde{\Sigma}_{z}=I$ ( $I$ is an identity operator), the Lamb shift does not affect the evolution of the two level system.

III.4 Dissipative Landau-Zener transition

In this subsection, we consider the dissipative Landau-Zener problem, in which the two-level system couples to a reservoir at zero temperature. The exact transition probabilities for such a model have been presentedWubs2006 ; Saito2007 . Here, we simulate the open system dynamics of the dissipative Landau-Zener problem by the DMME, and show that the transition probability given by the DMME almost coincides with the exact one. Meanwhile, a clear physical explanation is also presented.

The dissipative Landau-Zener problem is a scattering problem in the restricted sense that changes in the two-level system’s state will occur only during a finite time interval around $t=0$ . The two-level system’s Hamiltonian has the same form as Eq.(24), i.e.,

H_{\text{s}}(t)=\Delta(t)\sigma_{z}+\Omega(t)\sigma_{x},

(57)

with $\Delta(t)=vt/2$ and $\Omega=\Omega_{0}/2$ , where $v$ is the constant sweep velocity and $\Omega_{0}$ denotes the real intrinsic interaction amplitude between the diabatic states $|1\rangle$ and $|0\rangle$ . $\sigma_{x}=|1\rangle\langle 0|+|0\rangle\langle 1|$ and $\sigma_{z}=|1\rangle\langle 1|-|0\rangle\langle 0|$ stand for the $x$ and $z$ components of the Pauli operators. The eigenstates of the LRIs are still given by Eq.(28), in which $\eta(t)$ and $\zeta(t)$ can be obtained via solving the differential equations Eq.(29) by the help of the system Hamiltonian Eq.(57). Further, we assume that the interaction Hamiltonian takes the same form as Eq.(25) with $g_{k}^{y}=0$ for all $k$ , so that we obtain

\displaystyle H_{\text{I}}=\sigma_{x}\otimes\sum_{k}g_{k}^{x}(b_{k}^{\dagger}+b_{k}),

(58)

where $b_{k}$ and $g_{k}$ are, respectively, the annihilation operator and the coupling strength of the $k$ -th mode of the reservoir.

With the setting above, the DMME for the dissipative Landau-Zener problem can be written as the same form as Eq.(41) with $\Gamma_{mn}^{\text{R(I)},y}=0$ . For simplifying our discussion, we neglect the Lamb shifts ( $H_{\text{LS}}(t)=0$ ), and consider that the reservoir correlation time $\tau_{\text{B}}$ have to much smaller than the non-adiabatic phase timescale $\tau_{\text{d}}$ , i.e., $\tau_{\text{B}}\ll\tau_{\text{d}}$ . According to Eq.(22) and by using the relations Eqs. (35) and (37), the decoherence rates become

	$\displaystyle\Gamma^{\text{R}}_{-}(t,\alpha_{12}^{y})$	$\displaystyle=$	$\displaystyle\pi(\xi_{12}^{y}(t))^{2}J(\alpha_{12}^{y})N(\alpha_{12}^{y}),$
	$\displaystyle\Gamma^{\text{R}}_{+}(t,\alpha_{21}^{y})$	$\displaystyle=$	$\displaystyle\pi(\xi_{21}^{y}(t))^{2}J(\alpha_{21}^{y})N(\alpha_{21}^{y}),$		(59)

where $\xi_{12}^{y}(t)$ and $\xi_{21}^{y}(t)$ are given by Eq.(32). $J(\alpha)$ stands for the spectral density associated with the instantaneous frequency $\alpha_{mn}^{y}$ , which can be selected as

\displaystyle J(\alpha_{mn}^{y})=\kappa\alpha_{mn}^{y}\exp\left(-\frac{|\alpha_{mn}^{y}|}{\Omega_{c}}\right)

(60)

with the cut-off frequency $\Omega_{c}$ and a dimensionless coupling rate $\kappa$ . $N(\alpha_{mn}^{y})=\left(\exp(\alpha_{mn}^{y}/T_{R})-1\right)^{-1}$ denotes the Planck distribution with the reservoir temperature $T_{R}$ . Since the Planck distribution satisfies $N(-\alpha_{mn}^{y})=-\left(N(\alpha_{mn}^{y})+1\right)$ and $\alpha_{12}^{y}=-\alpha_{21}^{y}$ is always fulfilled, the instantaneous frequency $\alpha_{12}^{y}$ determines the transition direction caused by decoherence. When the reservoir is at zero temperature, it is easy to illustrate that, if $\alpha_{12}^{y}>0$ , we have $N(\alpha_{12}^{y})=1$ and $N(\alpha_{21}^{y})=0$ , which implies a decay from $|\psi_{2}(t)\rangle$ to $|\psi_{1}(t)\rangle$ ; if $\alpha_{12}^{y}<0$ , i.e., $\alpha_{21}^{y}>0$ , we have $N(\alpha_{12}^{y})=0$ and $N(\alpha_{21}^{y})=1$ , which leads to a decay from $|\psi_{1}(t)\rangle$ to $|\psi_{2}(t)\rangle$ . Since the instantaneous frequency $\alpha_{12}^{y}$ is time-dependent, the DMME may take the following form,

\displaystyle\partial_{t}\rho_{\text{s}}(t)

\displaystyle=

\displaystyle-{i}\left[H_{\text{s}}(t),\rho_{\text{s}}(t)\right]+\mathcal{D}\rho_{\text{s}}(t)

(61)

with

\displaystyle\mathcal{D}\rho_{\text{s}}=\begin{cases}\Gamma\left(\Sigma_{-}\rho_{\text{s}}\Sigma_{+}-\frac{1}{2}\left\{\Sigma_{+}\Sigma_{-},\rho_{\text{s}}\right\}\right)&\text{for}\,\alpha_{12}^{y}>0\\ \Gamma\left(\Sigma_{+}\rho_{\text{s}}\Sigma_{-}-\frac{1}{2}\left\{\Sigma_{-}\Sigma_{+},\rho_{+}\right\}\right)&\text{for}\,\alpha_{12}^{y}<0,\end{cases}

where $\Gamma=\kappa\pi(\xi_{12}^{y}(t))^{2}|\alpha_{12}^{y}|\exp\left(-{|\alpha_{12}^{y}|}/{\Omega_{c}}\right)$ is the decoherence rate depending on $t$ and $|\alpha_{12}^{y}|$ .

For such a dissipative Landau-Zener problem, the exact transition probability $P_{11}$ reads Wubs2006 ; Saito2007

\displaystyle P_{11}=\exp\left(-\frac{\pi W^{2}}{2v}\right)

(62)

with

\displaystyle W^{2}=\Omega_{0}^{2}+\sum_{k}\left(\frac{g_{k}^{x}}{2}\right)^{2}.

(63)

$P_{11}$ denotes the probability for a transition to the diabatic state $|1\rangle$ , if the initial state of the two level system is prepared in $|\psi_{1}(-\infty)\rangle=|1\rangle$ . More directly speaking, $P_{11}=\langle 1|\rho_{\text{s}}(\infty)|1\rangle$ is the population on the diabatic state $|1\rangle$ at $t=\infty$ . Considering the relation Eq.(35), i.e.,

\sum_{k}(g_{k}^{x})^{2}\rightarrow\int_{0}^{\infty}d\omega_{k}J(\omega_{k})

with the spectral density given by Eq.(60), we have

\displaystyle W^{2}=\Omega_{0}^{2}+\frac{\kappa}{4}\Omega_{c}^{2}.

(64)

In case of the closed system, the exact transition probability can be evaluated by means of $|{\langle 1|\psi_{1}(\infty)\rangle}|^{2}$ , since $|\psi_{1}(t)\rangle$ denotes the quantum state evolution for the closed system dynamics with the initial state $|1\rangle$ .

In FIG. 2 (a), we plot the evolution of the population on the diabatic state $|1\rangle$ for both the dissipative case $\rho_{11}\equiv\langle 1|\rho_{\text{s}}(t)|1\rangle$ (the blue solid line) and the closed case $|\langle 1|\psi_{1}(t)\rangle|^{2}$ (the yellow dashed line). Since we have selected $\Omega_{0}^{2}/v=4$ , the adiabatic condition for the Landau-Zener transition is close to be satisfied. Therefore, both the dissipative case and the closed case give an almost similar dynamical evolution. On the one hand, the instantaneous frequency $\alpha_{12}^{y}$ is always positive at the most of time as shown FIG. 2 (b). Hence, the instantaneous eigenstate $|\psi_{1}(t)\rangle$ must be the instantaneous steady state of the DMME (Eq.(61)). Therefore, in the adiabatic limits, the quantum state $\rho_{\text{s}}(t)$ must follow with $|\psi_{1}(t)\rangle$ . Also we find that $\alpha_{12}^{y}<0$ at a small interval around $t=0$ . At this time, the instantaneous steady state becomes $|\psi_{2}(t)\rangle$ , so that the dissipative dynamical evolution divides from the closed case. On the other hand, it can be observed that the transition probability given by the DMME is very close to the exact one $P_{11}$ (Eq.(62)). This can be understood as follows: In the adiabatic limits, i.e., $\Omega_{0}^{2}/v\gg 1$ , we have $\Omega_{0}^{2}\gg\sum_{k}(g_{k}^{x})^{2}$ , when the Born approximation (the weak coupling approximation) is satisfied. Thus, the influence of the dissipation on the transition probability $P_{11}$ is inapparent.

The non-adiabatic case is presented in FIG. 3, in which we have selected $\Omega_{0}^{2}/v=0.04$ . For the closed system case, since the adiabatic condition can not be fulfilled, the quantum state can not follow the eigenstates of the Hamiltonian into the diabatic state $|0\rangle$ , which are illustrated by the yellow dashed line in FIG. 3 (a). When the coupling to the reservoir is considered, the result is entirely different. We can observe that the numerical result given by the DMME is in good agreement with the exact transition probability (the red dotted line). For the exact transition probability (Eq.(62)), due to $\Omega^{2}/v\ll 1$ , the effect of the dissipation on the Landau-Zener transition is dominate, which reduces $P_{11}$ evidently.The decrease of $P_{11}$ can be understood by means of the DMME (Eq.(61)). As shown in FIG. 3 (b), the instantaneous frequency $\alpha_{12}^{y}(t)$ is positive at $t<0$ , so that the instantaneous steady state of the DMME is the eigenstate of the LRI $|\psi_{1}(t)\rangle$ . When $t>0$ , the instantaneous frequency $\alpha_{12}^{y}(t)$ becomes negative. Thus, the instantaneous steady state of the DMME is turned over to the other eigenstate of the LRI, i.e., $|\psi_{2}(t)\rangle$ , for $t>0$ . Therefore, the population will decays into $|\psi_{2}(t)\rangle$ gradually, and the final transition probability becomes $\rho_{11}(\infty)=|\langle 1|\psi_{2}(\infty)\rangle|$ .

IV Conclusion

The driven-Markovian master equation is derived by using the LRI theory within the Born-Markovian approximation in this paper. Since the unitary operator associated with the free propagator of the quantum system can be decomposed by the eigenstates of the LRI, our derivation overcomes the time-ordering obstacle in writing down a exact formula of the propagator for the free dynamics. Due to the rapid changing of the driving protocols, the non-adiabatic timescale may approach to, even is larger than, the reservoir correlation time, which leads to the memory effect of the driving protocols. The DMME presented here includes this memory effect, which leads to additional Lamb shifts and decoherence terms. Therefore, the DMME does not contains any constraint on the driving protocols, such as the adiabatic or inertial approximation Albash2012 ; Dann2021 .

According to the DMME, the transitions of the driven open quantum system occur on the eigenstates of the LRI, but not on the instantaneous eigenstates of the system Hamiltonian. This is very practical in determining the Lindblad operators in the DMME, if the LRIs are knownDann2018 , which is illustrated by the example of the driven two-level system. Similar to the Makovian master equation with a static Hamiltonian, both the energy relaxation and dephasing processes emerge in the dynamics of the driven two-level system. But the decoherence rates and the Lindblad operators are time-dependent, which implies a time-dependent steady state. Such a time-dependent steady state is an important candidate in the quantum state engineering of open quantum systems. What is more, if the reservoir is at low, or ultra-low, temperature, the steady state is close to a pure state which is one of the eigenstates of the LRI. Therefore, the inverse engineering method based on the LRIsChen2010 ; Chen2011 ; Wu2021 ; Tobalina2020 will be a more promising controlling method than the others in the field of the shortcuts to adiabaticity Campo2013 ; Odelin2019 .

Here, we would like to emphasise that the dynamics of the driven open quantum systems is closely connected to the symmetry of the system, which is contained in the LRIs of the corresponding system Hamiltonian. For instance, the DMME given in Eq.(41) can describe the dynamics only for open two-level systems with the system Hamiltonian like Eq.(24). For different types of the system Hamiltonian, the open system dynamics must be different due to distinct symmetries of the driven system. Therefore, to present the LRIs of a particular system with a certain driving protocol is a very essential task in applying the general formula of the DMME Ponte2018 . For driven two qubits system, a Lie-algebraic classification and detailed construction of the LRIs has been exploitedNakahara2012 . However, for the many body and multi-level system, this is a hard taskJezek1988 . Fortunately, there are many useful (semi-) simple subalgebras in a complicated Lie algebraSlansky1981 , which corresponds to many available driving potocols. Therefore, the DMME based on the LRI theory have broad potential applications in the quantum information processWild2016 ; Albash2018 and the quantum thermodynamicsGuarnieri2019 ; Vu2022 .

This work is supported by the National Natural Science Foundation of China (NSFC) under Grants Nos. 12075050, 12175033, 12205037, and National Key R&D Program of China (Grant No. 2021YFE0193500). We thank Prof. J. H. An for pointing out that the derivation of exact solution in the example ”Comparison to the Exactly Solvable Models I: the Dissipative Model” is wrong, which does not affect our main results.

Appendix A The Exact Evolution of the Dephasing Model

We start our derivation form taking an unitary transformation $U_{\text{s}}$ (Eq.(7)) on the total Hamiltonian, which leads to

	$\displaystyle h$	$\displaystyle=$	$\displaystyle U_{\text{s}}^{\dagger}HU_{\text{s}}$
		$\displaystyle=$	$\displaystyle\sigma_{z}\otimes\sum_{k}g_{k}\left(b_{k}^{\dagger}+b_{k}\right)+\sum_{k}\Omega_{k}b_{k}^{\dagger}b_{k}.$

Here an interaction Hamiltonian as Eq.(54) has been considered. It follows that the above Hamiltonian can be exactly diagonalizated by means of the unitary operator defined as $V=\exp\left(\sigma_{z}\otimes K\right)$ with $K=\sum_{k}\frac{g_{k}}{\Omega_{k}}\left(b_{k}^{\dagger}-b_{k}\right)$ Breuer2007 ; Uhrig2007 ,

	$\displaystyle\tilde{h}$	$\displaystyle=$	$\displaystyle VhV^{\dagger}$
		$\displaystyle=$	$\displaystyle\sum_{k}\Omega_{k}b_{k}^{\dagger}b_{k}-\sum_{k}\frac{g_{k}^{2}}{\Omega_{k}},$

in which the two-level system decouples to the heat reservoir. Mean while, it is easy to verify that $V\sigma_{x,y}V=\sigma_{x,y}$ . In the rotating frame given by $U_{\text{s}}$ , the evolution of the system quantum state can be written down formally

\tilde{\rho}(t)=U_{\text{eff}}(t)\rho(0)U_{\text{eff}}^{\dagger}(t),

where $U_{\text{eff}}(t)=\exp(-i\int_{0}^{t}h(s)ds)$ is the evolution operator for the total system, and $\rho(t)=U_{\text{s}}\tilde{\rho}(t)U_{\text{s}}^{\dagger}$ . Thus, an useful relation for $V$ and $U_{\text{eff}}(t)$ reads

	$\displaystyle U_{\text{eff}}^{\dagger}(t)V^{n}U_{\text{eff}}(t)$	$\displaystyle=$	$\displaystyle V^{\dagger}\left(VU_{\text{eff}}^{\dagger}(t)V^{\dagger}\right)V^{n}\left(VU_{\text{eff}}(t)V^{\dagger}\right)V$		(65)
		$\displaystyle=$	$\displaystyle V^{\dagger}\exp\left(n\sigma_{z}\otimes K(t)\right)V,$		(65)

with $K(t)=\exp\left(iH_{\text{B}}t\right)K\exp\left(-iH_{\text{B}}t\right)$ .

Let us consider an initial product state $\rho(0)=\rho_{\text{s}}(0)\otimes\rho_{\text{B}}$ with a heat equilibrium state at temperature $T_{\text{R}}$ and a system strate as

\rho_{s}(0)=\frac{1}{2}(I+\sum_{n}r_{n}\sigma_{n}).

Here, $\sigma_{n}$ are Pauli operators and $r_{n}$ is the corresponding component of Bloch vector satisfing $r_{n}=\text{Tr}_{\text{s}}\left\{\rho_{s}\sigma_{n}\right\}$ . Therefore, for obtain exact evolution of the quantum state, we needs to calculate the main values of the Pauli operators. For the x-component in the rotating frame given by $U_{\text{s}}$ , we find

$\displaystyle\tilde{r}_{x}(t)$	$\displaystyle=$	$\displaystyle\text{Tr}\left\{\tilde{\rho}(t)\sigma_{x}\right\}$
	$\displaystyle=$	$\displaystyle\text{Tr}\left\{\rho(0)U_{\text{eff}}^{\dagger}(t)\sigma_{x}U_{\text{eff}}(t)\right\}$
	$\displaystyle=$	$\displaystyle\text{Tr}\left\{\rho(0)U_{\text{eff}}^{\dagger}(t)V\sigma_{x}VU_{\text{eff}}(t)V^{\dagger}V\right\},$

where $V\sigma_{x}V=\sigma_{x}$ has been used. Since $VU_{\text{eff}}(t)V^{\dagger}=\exp(-i\int_{0}^{t}\tilde{h}(s)ds),$ it yields

	$\displaystyle\tilde{r}_{x}(t)$	$\displaystyle=$	$\displaystyle\text{Tr}\left\{\rho(0)U_{\text{eff}}^{\dagger}(t)V^{2}U_{\text{eff}}(t)V^{\dagger}\sigma_{x}V\right\}$
		$\displaystyle=$	$\displaystyle\text{Tr}\left\{\rho(0)U_{\text{eff}}^{\dagger}(t)V^{2}U_{\text{eff}}(t)V^{\dagger 2}\sigma_{x}\right\}.$

By considering Eq.(65), we arrive at

$\displaystyle\tilde{r}_{x}(t)$	$\displaystyle=$	$\displaystyle\text{Tr}\left\{\rho(0)V^{\dagger}\exp\left(2\sigma_{z}\otimes K(t)\right)V^{\dagger}\sigma_{x}\right\}$
	$\displaystyle=$	$\displaystyle\text{Tr}\left\{\left(\sigma_{x}\rho_{\text{s}}(0)\otimes\rho_{\text{B}}\right)\exp\left(2\sigma_{z}\otimes\left(K(t)-K(0)\right)\right)\right\}$
	$\displaystyle=$	$\displaystyle\langle 1\|\sigma_{x}\rho_{\text{s}}(0)\|1\rangle\left\langle\exp\left(2\left(K(t)-K(0)\right)\right)\right\rangle$
		$\displaystyle+\langle 0\|\sigma_{x}\rho_{\text{s}}(0)\|0\rangle\left\langle\exp\left(-2\left(K(t)-K(0)\right)\right)\right\rangle,$

where $\sigma_{z}|1\rangle=|1\rangle$ and $\sigma_{x}|0\rangle=-|0\rangle$ have been used, and $\left\langle\exp\left(\pm 2\left[K(t)-K(0)\right]\right)\right\rangle=\text{Tr}_{\text{B}}\left\{\rho_{\text{B}}\exp\left(\pm 2\left[K(t)-K(0)\right]\right)\right\}$ are the Wigner characteristic function of the reservoir mode $k$ . It can be easily determined by noting that it represents a Gaussian function, which immediately leads to

			$\displaystyle\left\langle\exp\left(\pm 2\left(K(t)-K(0)\right)\right)\right\rangle$
		$\displaystyle=$	$\displaystyle\exp\left(2\left\langle\left(K(t)-K(0)\right)^{2}\right\rangle\right)$
		$\displaystyle=$	$\displaystyle\exp\left(-4\sum_{k}\frac{g_{k}^{2}}{\Omega_{k}^{2}}\left\langle 2b_{k}^{\dagger}b_{k}+1\right\rangle\left(1-\cos\Omega_{k}t\right)\right).$

We now perform the continuum limit of the bath modes. Introducing the density $f(\Omega_{k})$ of the modes of frequency $\Omega_{k}$ and defining the spectral density asBreuer2007

J(\Omega_{k})=4f(\Omega_{k})g_{k}^{2},

we can write down the decoherence function

	$\displaystyle\Gamma_{\text{e}}(t)$	$\displaystyle=$	$\displaystyle\int_{0}^{\infty}d\Omega_{k}\frac{J(\Omega_{k})}{\Omega_{k}^{2}}\left(2N_{k}+1\right)\left(1-\cos\Omega_{k}t\right)$
		$\displaystyle=$	$\displaystyle\frac{1}{2}\ln\left(1+\Omega_{c}^{2}t^{2}\right),$

where the spectral density Eq.(56) has been used. Therefore, we obtain

\tilde{r}_{x}(t)=r_{x}(0)\exp\left(-\Gamma_{\text{e}}(t)\right),

with $N_{k}=\left\langle b_{k}^{\dagger}b_{k}\right\rangle$ . With the same procedure, we can determine the y-component of the Bloch vector, which present a similar expression as $\tilde{r}_{x}$ ,

\tilde{r}_{y}(t)=r_{y}(0)\exp\left(-\Gamma_{\text{e}}(t)\right).

Due to $\left[\sigma_{z},U_{\text{eff}}(t)\right]=0$ , we have $\tilde{r}_{z}(t)=r_{z}(0)$ . Thus, the exact quantum state evolution of the driven two-level system in the Schrödinger picture can be obtained by using the unitary transformation $U_{\text{s}}(t),$

	$\displaystyle\rho(t)$	$\displaystyle=$	$\displaystyle U_{\text{s}}(t)\tilde{\rho}(t)U_{\text{s}}^{\dagger}(t)$
		$\displaystyle=$	$\displaystyle\frac{1}{2}(I+\sum_{n}\tilde{r}_{n}(t)U_{\text{s}}(t)\sigma_{n}U_{\text{s}}^{\dagger}(t)).$

References

(1) H. P. Breuer and F. Petruccione, The Theory of Open Quantum Systems (Oxford University Press, New York, 2007).
(2) E. B. Davies, Commun. Math. Phys. 39, 91 (1974).
(3) V. Gorini, A. Kossakowski, and E. C. G. Sudarshan, J. Math. Phys. 17, 821 (1976).
(4) G. Lindblad, Commun. Math. Phys. 48, 119 (1976).
(5) E. Davies and H. Spohn, J. Stat. Phys. 19, 511 (1978).
(6) T. Albash, S. Boixo, D. A. Lidar, and P. Zanardi, New J. Phys. 14, 123016 (2012).
(7) A. M. Childs, E. Farhi, and J. Preskill, Phys. Rev. A 65, 012322 (2001).
(8) I. Kamleitner, Phys.Rev.A 87, 042111 (2013).
(9) M. S. Sarandy, L. A. Wu, and D. A. Lidar, Quant. Info. Proc. 3, 331 (2004).
(10) M. Yamaguchi, T. Yuge, and T. Ogawa, Phys. Rev. E 95, 012136 (2017).
(11) R. Dann, A. Levy, and R. Kosloff, Phys. Rev. A 98, 052129 (2018).
(12) P. P. Potts, A. A. Sand Kalaee, and A. Wacker, New J. Phys. 23, 123013 (2021).
(13) V. V. Albert, B. Bradlyn, M. Fraas, and L. Jiang, Phys. Rev. X 6, 041031 (2016).
(14) S. Scopa, G. T. Landi, A. Hammoumi, and D. Karevski, Phys. Rev. A 99, 022105 (2019).
(15) R. Dann, A. Tobalina, and R. Kosloff, Phys. Rev. A 101, 052102 (2020).
(16) R. Dann, and R. Kosloff, Phys. Rev. Research 3, 013064 (2021).
(17) R. Dann, A. Tobalina, and R. Kosloff, Phys. Rev. Lett. 122, 250402 (2019).
(18) B. Hwang, and H. S. Goan, Phys. Rev. A 85, 032321 (2012).
(19) J. T. Stockburger, Euro. Phys. Lett. 115, 40010 (2016).
(20) S. Ozaki and H. Nakazato, Phys. Rev. A 103, 053713 (2021).
(21) G. Taran, E. Bonet, and W. Wernsdorfer, Phys. Rev. B 99, 180408(R) (2019).
(22) H. R. Lewis Jr., J. Math. Phys. 9, 1976 (1968).
(23) H. R. Lewis Jr. and W. B. Riesenfeld, J. Math. Phys. 10, 1458 (1969).
(24) A. G. Redfield, Advances in Magnetic and Optical Resonance (Elsevier, 1965).
(25) I. A. Pedrosa and A. Rosas, Phys. Rev. Lett. 103, 010402 (2009).
(26) M. S. Sarandy, E. I. Duzzioni, and M. H. Y. Moussa, Phys. Rev. A 76, 052112 (2007).
(27) D. Schuch and M. Moshinsky, Phys. Rev. A 73, 062111 (2006).
(28) X. Chen, E. Torrontegui, and J. G. Muga, Phys. Rev. A 83, 062116 (2011).
(29) D. B. Monteoliva, H. J. Korsch, and J, A, Nunez, J. Phys. A 27, 6897 (1994).
(30) U. Güngördü, Y. Wan, M. A. Fasihi, and M. Nakahara, Phys. Rev. A 86, 062312 (2012).
(31) L. S. Simeonov and N. V. Vitanov, Phys. Rev. A 93, 012123 (2016).
(32) S.P. Kim, A. E. Santana, F.C, Khanna, Phys. Lett. A 272, 46 (2000).
(33) M. A. de Ponte, P. M. Consoli, and M. H. Y. Moussa, Phys. Rev. A 98, 032102 (2018).
(34) T. Petrosky and I. Prigogine, The Liouville Space Extension of Quantum Mechanics, edited by I. Prigogine and S. A. Rice (John Wiley & Sons, New York, 1997).
(35) X. Chen, I. Lizuain, A. Ruschhaupt, D. Guéry-Odelin, and J. G. Muga, Phys.Rev.Lett. 105, 123003 (2010).
(36) H. Z. Shen, M. Qin, X. M. Xiu, and X. X. Yi, Phys. Rev. A 89, 062113 (2014).
(37) Y. Z. Lai, J. Q. Liang, H. J. W. Müller-Kirsten, and J.-G. Zhou, Phys. Rev. A 53, 3691 (1996).
(38) G. Uhrig, Phys. Rev. Lett. 98, 100504 (2007).
(39) X. X. Yi and G. C. Guo, Phys. Rev. A 62, 062312 (2000).
(40) Martijn Wubs, Keiji Saito, Sigmund Kohler, Peter Hänggi, and Yosuke Kayanuma, Phys. Rev. Lett. 97, 200404 (2006).
(41) Keiji Saito, Martijn Wubs, Sigmund Kohler, Yosuke Kayanuma, and Peter Hänggi, Phys. Rev. B 75, 214308 (2007).
(42) S. L. Wu, W. Ma, X. L. Huang, and X. X. Yi, Phys. Rev. Applied 16, 044028 (2021).
(43) A. Tobalina, E. Torrontegui, I. Lizuain, M. Palmero, and J. G. Muga, Phys. Rev. A 102, 063112 (2020).
(44) A. del Campo, Phys. Rev. Lett. 111, 100502 (2013).
(45) D. Guéry-Odelin, A. Ruschhaupt, A. Kiely, E. Torrontegui, S. Martínez-Garaot, and J. G. Muga, Rev. Mod. Phys. 91, 045001 (2019).
(46) E. S. Hernández and D. M. Jezek, Phys. Rev. A 38, 4455 (1988).
(47) R. Slansky, Phys. Rep. 79, 1 (1981).
(48) D. S. Wild, S. Gopalakrishnan, M. Knap, N. Y. Yao, and M. D. Lukin, Phys. Rev. Lett. 117, 150501 (2016).
(49) T. Albash and D. A. Lidar, Rev. Mod. Phys. 90, 015002 (2018).
(50) G. Guarnieri, G. T. Landi, S. R. Clark, and J. Goold, Phys. Rev. Research 1, 033021 (2019).
(51) T. VanVu, and K. Saito, Phys. Rev. Lett. 128, 140602 (2022).