Markovian and non-Markovian master equations versus an exactly solvable model of a qubit in a cavity

Zihan Xia Juan Garcia-Nila Center for Quantum Information Science & Technology Department of Electrical & Computer Engineering Daniel A. Lidar Center for Quantum Information Science & Technology Department of Electrical & Computer Engineering Department of Physics & Astronomy Department of Chemistry
University of Southern California, Los Angeles, CA 90089, USA

Abstract

Quantum master equations are commonly used to model the dynamics of open quantum systems, but their accuracy is rarely compared with the analytical solution of exactly solvable models. In this work, we perform such a comparison for the damped Jaynes-Cummings model of a qubit in a leaky cavity, for which an analytical solution is available in the one-excitation subspace. We consider the non-Markovian time-convolutionless master equation up to the second (Redfield) and fourth orders as well as three types of Markovian master equations: the coarse-grained, cumulant, and standard rotating-wave approximation (RWA) Lindblad equations. We compare the exact solution to these master equations for three different spectral densities: impulse, Ohmic, and triangular. We demonstrate that the coarse-grained master equation outperforms the standard RWA-based Lindblad master equation for weak coupling or high qubit frequency (relative to the spectral density high-frequency cutoff $\omega_{c}$ ), where the Markovian approximation is valid. In the presence of non-Markovian effects characterized by oscillatory, non-decaying behavior, the TCL approximation closely matches the exact solution for short evolution times (in units of $\omega_{c}^{-1}$ ) even outside the regime of validity of the Markovian approximations. For long evolution times, all master equations perform poorly, as quantified in terms of the trace-norm distance from the exact solution. The fourth-order time-convolutionless master equation achieves the top performance in all cases. Our results highlight the need for reliable approximation methods to describe open-system quantum dynamics beyond the short-time limit.

I Introduction

The study of open quantum systems presents both conceptual and technical challenges due to the complexity and high dimensionality of the environment, or bath. Exact analytical solutions describing the joint system-bath evolution are rarely attainable, necessitating the development of approximation methods to capture the reduced system dynamics [1, 2, 3]. To address this challenge, various approaches have been developed to derive master equations that describe the system’s evolution using the reduced density matrix, the best known of which is the Markovian Lindblad [or Gorini-Kossakowski-Lindblad-Sudharshan (GKLS)] equation [4, 5]. Numerous other master equations have been derived, some of which include non-Markovian effects. In some cases, rigorous error bounds have been derived that quantify the deviation between the solutions of known master equations and the exact solution [6]. However, these tend to be rather loose. Hence, it is desirable to compare the predictions of various master equations to non-trivial examples of exactly solvable open-system problems. This has been done, e.g., for the central spin model [7, 8].

This is the goal of the present work, where we study the damped Jaynes-Cummings model of a qubit inside a leaky cavity [9]. In this model, the qubit system interacts with the cavity electromagnetic field through the dipole approximation [10]. The qubit decay rate can be associated with experimentally measurable parameters such as the dipole moment and the energy gap [11]. Experimental proposals for simulating the spin-boson model with an Ohmic spectral density using superconducting circuits have been previously discussed [12, 13, 14]. We solve this model analytically assuming the zero temperature limit and the $1$ -excitation subspace of the joint qubit-cavity system, where the cavity is populated by at most a single photon. This is similar to previous studies assuming a Lorentzian spectral density [15, 16], but we do this for three different spectral densities: impulse, Ohmic, and triangular (formally defined below). These choices are motivated by experiments involving condensed matter systems such as superconducting qubits interacting with bosonic modes [17], rather than the original quantum-optical setting of an atom in a cavity that inspired the damped Jaynes-Cummings model. We then compare the exact solution to a number of different Markovian and non-Markovian master equations. We find that in the weak-coupling limit and for large qubit frequencies relative to the spectral density cutoff, where the Markovian approximation holds, the Lindblad equation derived using the coarse-graining approach [18] is more accurate than the standard rotating wave approximation based Lindblad equation. In the non-Markovian regime, the time-convolutionless master equation proves to be accurate in approximating the exact solution for relatively short evolution times. All the approximation methods we consider struggle at long evolution times.

This paper is structured as follows. In Section II, we introduce the damped Jaynes-Cummings model and derive the exact solutions for all three spectral densities. In Section III, we explore various Markovian approximation methods in the context of the damped Jaynes-Cummings model. We begin with the coarse-grained Lindblad equation (CG-LE) in Section III.1, then the cumulant LE (C-LE) in Section III.2, and finally the commonly used rotating-wave approximation (RWA)-based LE (RWA-LE, Section III.3). We apply the time-convolutionless (TCL) approach to second (TCL2, also known as the Redfield equation) and fourth orders (TCL4) in Section III.4. Section IV is the heart of this work, where we present comparisons between the exact solutions and the various approximation methods. This includes a comparison of the exact solution with Markovian and TCL approximations for the Ohmic spectral density, a comparison of the exact solution with the Markovian CG-LE and RWA-LE, and finally a comparison with the non-Markovian TCL models for the impulse and triangular spectral densities. We summarize our findings and present our conclusions in Section V. A variety of technical details that complement the main text are presented in the Appendices.

Readers who are already familiar with the different types of master equations may choose to skip Section III. All our key analytical results are conveniently accessible via Table 1, which provides the corresponding equation numbers. Readers who are interested primarily in the results of the comparison between the exact model results and the various master equations may choose to skip ahead to Section IV and focus on the graphs presented there.

II Exact Dynamics

II.1 General open system setup

The total Hamiltonian of the system and the bath is given by

\displaystyle H=H_{0}+H_{SB}\ ,

(1)

where $H_{0}=H_{S}\otimes I_{B}+I_{S}\otimes H_{B}$ with $H_{S}$ and $H_{B}$ being the pure system and bath Hamiltonians, respectively, $I$ the identity operator.

We move to the interaction picture, where all operators transform according to:

\displaystyle X\mapsto\tilde{X}(t)=e^{iH_{0}t}Xe^{-iH_{0}t}\ .

(2)

The dynamics of the total system in the interaction picture are governed by the Liouville-von Neumann equation:

\displaystyle\frac{d\tilde{\rho}_{SB}}{dt}=-i[\tilde{H}_{SB},\tilde{\rho}_{SB}]\ ,

(3)

where $\tilde{\rho}_{SB}$ is the density matrix of the total system acting on the Hilbert space $\mathcal{H}_{SB}=\mathcal{H}_{S}\otimes\mathcal{H}_{B}$ . The joint system-bath state is thus given by:

\displaystyle\tilde{\rho}_{SB}(t)=\tilde{U}(t)\tilde{\rho}_{SB}(0)\tilde{U}^{\dagger}(t)\ ,

(4)

where the unitary evolution operator is:

\displaystyle\tilde{U}(t)=T_{+}\exp(-i\int_{0}^{t}\tilde{H}_{SB}(t^{\prime})dt^{\prime})\ .

(5)

and $T_{+}$ denotes Dyson time-ordering.

The solution can equivalently be expressed as a Dyson series by integrating and iterating Eq. 3:

		$\displaystyle\tilde{\rho}_{SB}(t)=\rho_{SB}(0)+\sum_{n=1}^{\infty}(-i)^{n}\int_{0}^{t}dt_{1}\int_{0}^{t_{1}}dt_{2}\cdots$		(6)
		$\displaystyle\ \ \int_{0}^{t_{n-1}}dt_{n}[\tilde{H}_{SB}(t_{1}),[\tilde{H}_{SB}(t_{2}),...[\tilde{H}_{SB}(t_{n}),\rho_{SB}(0)]]....]\ .$

The state of the system is given by the reduced density operator:

\displaystyle\tilde{\rho}(t)=\mathrm{Tr}_{B}[\tilde{\rho}_{SB}(t)]\ ,

(7)

where $\mathrm{Tr}_{B}$ denotes the partial trace over the bath state.

II.2 Model of a qubit in a leaky cavity

We analyze the dynamics of a single qubit inside a leaky cavity, coupled to a bosonic bath at zero temperature. By working in the single-excitation subspace this model becomes analytically solvable, and we closely follow the solution method of Refs. [10, 15] (see also Ref. [2]). However, we consider different bath spectral densities.

The system Hamiltonian, $H_{S}$ , can be expressed as

\displaystyle H_{S}=\Omega_{0}|{1}\rangle\!\langle 1|=\Omega_{0}\sigma_{+}\sigma_{-}\ .

(8)

Here, $\sigma_{+}=|{1}\rangle\!\langle 0|$ and $\sigma_{-}=|{0}\rangle\!\langle 1|$ are the raising and lowering operators for the qubit, respectively. The qubit ground state is $\ket{0}$ with energy $0$ and its excited state is $\ket{1}$ with energy $\Omega_{0}$ . The bath Hamiltonian, $H_{B}$ , is given by:

\displaystyle H_{B}=\sum_{k}\omega_{k}b_{k}^{\dagger}b_{k}=\sum_{k}\omega_{k}n_{k}\ .

(9)

Here, $b_{k}$ and $b_{k}^{\dagger}$ represent the annihilation and creation operators for the bosonic modes and $n_{k}$ is the number operator for mode $k$ with energy $\omega_{k}$ (we set $\hbar=1$ throughout). The qubit-cavity interaction Hamiltonian is

\displaystyle H_{SB}=\sigma_{+}\otimes B+\sigma_{-}\otimes B^{\dagger}\ ,

(10)

where $B=\sum_{k}g_{k}b_{k}$ . Here the $g_{k}$ ’s are coupling constants with dimensions of energy. Introducing complex phases in the couplings can model chirality [19].

In the interaction picture, the interaction Hamiltonian $H_{SB}$ becomes:


$\displaystyle\tilde{H}_{SB}(t)$	$\displaystyle=\sigma_{+}(t)\otimes B(t)+\sigma_{-}(t)\otimes B^{\dagger}(t)$	(11a)
$\displaystyle\sigma_{\pm}(t)$	$\displaystyle=e^{\pm i\Omega_{0}t}\sigma_{\pm}\ ,\quad B(t)=\sum_{k}e^{-i\omega_{k}t}g_{k}b_{k}\ .$	(11b)

This model is not analytically solvable in general, but it is when we make the assumption that the cavity supports at most one photon. We thus consider the initial joint system-bath state to be:

\displaystyle\ket{\phi(0)}=c_{0}(0)\ket{\psi_{0}}+c_{1}(0)\ket{\psi_{1}}+\sum_{k}\mathfrak{c}_{k}(0)\ket{\varphi_{k}}\ ,

(12)

where


$\displaystyle\ket{\psi_{0}}$	$\displaystyle=\ket{0}_{S}\otimes\ket{v}_{B}$	(13a)
$\displaystyle\ket{\psi_{1}}$	$\displaystyle=\ket{1}_{S}\otimes\ket{v}_{B}$	(13b)
$\displaystyle\ket{\varphi_{k}}$	$\displaystyle=\ket{0}_{S}\otimes\ket{k}_{B}\ .$	(13c)

Here $\ket{v}_{B}$ denotes the vacuum state of the cavity, and $\ket{k}_{B}=b_{k}^{\dagger}\ket{v}_{B}=\ket{0_{1},\cdots,0_{k-1},1_{k},0_{k+1},\cdots}$ denotes the state with one photon in mode $k$ . The subspace spanned by $\{\ket{\psi_{0}},\ket{\psi_{1}},\ket{\varphi_{k}}\}$ is referred to as the $1$ -excitation subspace, and is conserved under the Hamiltonian in Eq. 1. I.e., the joint state remains in the following form for all time $t$ :

\displaystyle\ket{\phi(t)}=c_{0}(t)\ket{\psi_{0}}+c_{1}(t)\ket{\psi_{1}}+\sum_{k}\mathfrak{c}_{k}(t)\ket{\varphi_{k}}\ ,

(14)

subject to the normalization condition:

\displaystyle|c_{0}(t)|^{2}+|c_{1}(t)|^{2}+\sum_{k}|\mathfrak{c}_{k}(t)|^{2}=1\ .

(15)

The problem remains solvable with a more general bath state assuming interaction with a continuous-mode laser field [20].

We assume that initially there are no photons in the cavity [10], hence:

\displaystyle\mathfrak{c}_{k}(0)=0\ .

(16)

We introduce a spectral density $J(\omega)$ via:

\displaystyle\sum_{k}|g_{k}|^{2}e^{-i\omega_{k}t}=\int_{0}^{\infty}d\omega J(\omega)e^{-i\omega t}\ .

(17)

The continuum of bath spectral modes is a necessary condition for irreversibility; a discrete spectrum necessarily results in recurrences.

To complete the model specification, we consider three different bath spectral densities: an impulse function centered at the cutoff frequency $\omega_{c}$ , an Ohmic function with the same cutoff frequency, and a triangular spectral density with a sharp cutoff, namely:


$\displaystyle J_{1}(\omega)$	$\displaystyle=\|g\|^{2}\delta(\omega-\omega_{c})$	(18a)
$\displaystyle J_{2}(\omega)$	$\displaystyle=\eta\omega e^{-\omega/\omega_{c}}$	(18b)
$\displaystyle J_{3}(\omega)$	$\displaystyle=\eta\omega\Theta(\omega_{c}-\omega)\ ,$	(18c)

where the Heaviside function obeys $\Theta(x)=0$ for $x<0$ and $\Theta(x)=1$ for $x\geq 0$ . In the impulse spectral density $J_{1}(\omega)$ , the bath has a singular response at the cutoff frequency, characterized by a Dirac delta function. The Ohmic spectral density $J_{2}(\omega)$ is ubiquitous in the study of the spin-boson problem [21]. The dimensionless parameter $\eta$ in the Ohmic spectral density serves as a measure of the coupling strength between the bath and the system, while the ratio of the qubit frequency to the cutoff frequency measures the ratio of the photonic energy gap, which is the energy between the ground state and the excited state, and the number of frequency modes before reaching the cutoff frequency. The triangular spectral density $J_{3}(\omega)$ is a sharp-cutoff approximation to the Ohmic spectral density, which we introduce to simplify analytical calculations. The case of a Lorentzian spectral density was studied in Ref. [2].

II.3 Exact Solution

Considering the entire qubit-cavity system as closed, its evolution is governed by the Schrödinger equation, which can be expressed as [see Appendix A]:

	$\displaystyle i\partial_{t}\ket{\phi(t)}$	$\displaystyle=\sum_{k}g_{k}\mathfrak{c}_{k}(t)e^{i(\Omega_{0}-\omega_{k})t}\ket{\psi_{1}}$
		$\displaystyle+\sum_{k}g_{k}^{*}c_{1}(t)e^{-i(\Omega_{0}-\omega_{k})t}\ket{\varphi_{k}}\ .$		(19)

Multiplying this equation by $\bra{\psi_{1}}$ or $\bra{\varphi_{k}}$ , we obtain the following set of differential equations for the amplitudes:


$\displaystyle\dot{c}_{0}(t)$	$\displaystyle=0$	(20a)
$\displaystyle\dot{c}_{1}(t)$	$\displaystyle=-i\sum_{k}g_{k}\mathfrak{c}_{k}(t)e^{i(\Omega_{0}-\omega_{k})t}$	(20b)
$\displaystyle\dot{\mathfrak{c}}_{k}(t)$	$\displaystyle=-ig_{k}^{*}c_{1}(t)e^{-i(\Omega_{0}-\omega_{k})t}$	(20c)

Integrating, we arrive at:


$\displaystyle c_{0}(t)$	$\displaystyle=c_{0}(0)$	(21a)
$\displaystyle\mathfrak{c}_{k}(t)$	$\displaystyle=-i\int_{0}^{t}dt^{\prime}g_{k}^{*}c_{1}(t^{\prime})e^{-i(\Omega_{0}-\omega_{k})t^{\prime}}$	(21b)

Let us define the memory kernel $f(t)$ as the Fourier transform of the spectral density, shifted by the qubit frequency $\Omega_{0}$ :

\displaystyle f(t)=\int_{0}^{\infty}d\omega J(\omega)e^{i(\Omega_{0}-\omega)t}\ .

(22)

Substituting Eq. 21b into Eq. 20b, we obtain:

\displaystyle\dot{c}_{1}(t)=-\int_{0}^{t}dt^{\prime}f(t-t^{\prime})c_{1}(t^{\prime})\ ,

(23)

which can be solved via a Laplace transform since the RHS is a convolution. Denoting the Laplace transform $\operatorname{Lap}$ of a general function $g(t)$ by $\hat{g}(s)$ , and recalling that $\operatorname{Lap}[\dot{g}(t)]=s\hat{g}(s)-g(0)$ , we have:

\displaystyle\hat{c}_{1}(s)=\frac{c_{1}(0)}{s+\hat{f}(s)}\ ,

(24)

and $c_{1}(t)$ is then found via the inverse Laplace transform of Eq. 24.

Next, in order to determine the interaction picture system state by tracing out the bath state, we can utilize Eq. 14 and find:

\displaystyle\tilde{\rho}(t)

\displaystyle=\mathrm{Tr}_{B}\left(|{\phi(t)}\rangle\!\langle\phi(t)|\right)=\left(\begin{array}[]{cc}1-|c_{1}|^{2}&c_{0}c_{1}^{*}(t)\\ c_{0}^{*}c_{1}(t)&|c_{1}|^{2}\end{array}\right)\ ,

(27)

so that:

\displaystyle\dot{\tilde{\rho}}(t)=\left(\begin{array}[]{cc}-\partial_{t}|c_{1}|^{2}&c_{0}\dot{c}_{1}^{*}(t)\\ c_{0}^{*}\dot{c}_{1}(t)&\partial_{t}|c_{1}|^{2}\end{array}\right)\ .

(30)

At this point, it is straightforward to verify that the dynamics are time-local,

\displaystyle\dot{\tilde{\rho}}

\displaystyle=\mathcal{K}_{S}(t)\tilde{\rho}\ ,

(31)

with a generator given by:

	$\displaystyle\mathcal{K}_{S}(t)\tilde{\rho}(t)$	$\displaystyle=-\frac{i}{2}\mathcal{S}(t)[\sigma_{+}\sigma_{-},\tilde{\rho}(t)]$
		$\displaystyle+\gamma(t)\left(\sigma_{-}\tilde{\rho}(t)\sigma_{+}-\frac{1}{2}\{\sigma_{+}\sigma_{-},\tilde{\rho}(t)\}\right)\ ,$		(32)

provided we identify:


$\displaystyle\mathcal{S}(t)$	$\displaystyle=-2\imaginary\left(\frac{\dot{c}_{1}(t)}{c_{1}(t)}\right)$	(33a)
$\displaystyle\gamma(t)$	$\displaystyle=-2\real\left(\frac{\dot{c}_{1}(t)}{c_{1}(t)}\right)$	(33b)

The rate $\gamma(t)$ can be negative, corresponding to non-Markovian dynamics according to the CP non-divisibility criterion [22].

We focus on the excited state population and the coherence, which evolve in time according to:


$\displaystyle\tilde{\rho}_{11}(t)$	$\displaystyle=\|c_{1}(t)\|^{2}=\|c_{1}(0)\|^{2}\exp{-\int_{0}^{t}\gamma(t^{\prime})dt^{\prime}}$	(34a)
$\displaystyle\tilde{\rho}_{01}(t)$	$\displaystyle=c_{0}{c}_{1}^{}(t)=c_{0}{c}_{1}^{}(0)\exp{\frac{1}{2}\int_{0}^{t}\left(i\mathcal{S}(t^{\prime})-\gamma(t^{\prime})\right)dt^{\prime}}\ .$	(34b)

Details of the derivation above can be found in Appendix A.

We next discuss the solutions for the three spectral densities of Eq. 18. In each case, we express the solution in terms of $c_{1}(t)$ or its Laplace transform.

II.4 Exact solution for three different spectral densities

II.4.1 $J_{1}=|g|^{2}\delta(\omega-\omega_{c})$

As a toy example, we consider the impulse bath spectral density, which replaces the continuum of bath modes with a single mode. Consequently, we do not expect irreversibility and indeed, the solution is oscillatory. To demonstrate this, we circumvent the Laplace transform and instead use Eq. 22 to write:

\displaystyle f(t)

\displaystyle=|g|^{2}e^{i(\Omega_{0}-\omega_{c})t}\ .

(35)

As a result, Eq. 23 becomes:

\displaystyle\dot{c}_{1}(t)

\displaystyle=-|g|^{2}\int_{0}^{t}dt^{\prime}e^{i(\Omega_{0}-\omega_{c})(t-t^{\prime})}c_{1}(t^{\prime})\ .

(36)

Differentiating both sides yields:

\displaystyle\ddot{c}_{1}(t)-i(\Omega_{0}-\omega_{c})\dot{c}_{1}(t)+|g|^{2}c_{1}(t)=0\ ,

(37)

whose solution is:

\displaystyle\frac{c_{1}(t)}{c_{1}(0)}

\displaystyle=e^{i(\Omega_{0}-\omega_{c})t/2}\left(\cos(\frac{t}{2}\delta)-i\frac{\Omega_{0}-\omega_{c}}{\delta}\sin(\frac{t}{2}\delta)\right),

(38)

where $\delta$ is a real number:

\displaystyle\delta=\sqrt{(\Omega_{0}-\omega_{c})^{2}+4g^{2}}\ .

(39)

Since the solution is perfectly periodic, any master equation approximation with a non-trivial dissipator term will deviate from this exact solution for sufficiently long times. Note that the excited state population $|c_{1}(t)|^{2}\in|c_{1}(0)|^{2}[\frac{(\Omega_{0}-\omega_{c})^{2}}{(\Omega_{0}-\omega_{c})^{2}+4g^{2}},1]$ .

The Lamb shift and decay rate are now found from Eq. 33 to be:


$\displaystyle\mathcal{S}(t)$	$\displaystyle=\left(1-\frac{(\Omega_{0}-\omega_{c})^{2}}{\delta^{2}}\right)\frac{\frac{\Omega_{0}-\omega_{c}}{\delta}}{\cot^{2}\left(\frac{t\delta}{2}\right)+\frac{(\Omega_{0}-\omega_{c})^{2}}{\delta^{2}}}$	(40a)
$\displaystyle\gamma(t)$	$\displaystyle=\left(1-\frac{(\Omega_{0}-\omega_{c})^{2}}{\delta^{2}}\right)\frac{\cot\left(\frac{t\delta}{2}\right)}{\cot^{2}\left(\frac{t\delta}{2}\right)+\frac{(\Omega_{0}-\omega_{c})^{2}}{\delta^{2}}}$	(40b)

II.4.2 $J_{2}(\omega)=\eta\omega e^{-\omega/\omega_{c}}$

Now we turn our attention to the Ohmic spectral density. The memory kernel, as given in Eq. 22, takes the form:

\displaystyle f(t)=\frac{\eta\omega_{c}^{2}e^{i\Omega_{0}t}}{(1+i\omega_{c}t)^{2}}\ .

(41)

To obtain the Laplace transform analytically, we integrate by parts and obtain:

	$\displaystyle\hat{f}(s)=$	$\displaystyle\eta\Big{[}\left(s-i\Omega_{0}\right)e^{-\frac{\Omega_{0}+is}{\omega_{c}}}\left(i\frac{\pi}{2}-\operatorname{Ci}\left(\frac{s-i\Omega_{0}}{\omega_{c}}\right)\right.$
		$\displaystyle\left.\quad-i\operatorname{Si}\left(\frac{s-i\Omega_{0}}{\omega_{c}}\right)\right)-i\omega_{c}\Big{]}\ ,$		(42)

where the sine and cosine integral functions are defined as:


$\displaystyle\operatorname{Si}(z)$	$\displaystyle\equiv\int_{0}^{z}\frac{\sin(t)}{t}dt=\pi/2-\int_{z}^{\infty}\frac{\sin(t)}{t}dt$	(43a)
$\displaystyle\operatorname{Ci}(z)$	$\displaystyle\equiv-\int_{z}^{\infty}\frac{\cos(t)}{t}dt$	(43b)
$\displaystyle\operatorname{Ei}(z)$	$\displaystyle\equiv-\int_{-z}^{\infty}\frac{e^{-t}}{t}dt\ ,$	(43c)

where for future reference we have also defined the exponential integral function.

Due to the absence of analytic expressions for the Laplace transforms of the trigonometric integral functions, we use the numerical inverse Laplace transform of Eq. 24 to obtain $c_{1}(t)$ . The Lamb shift and decay rate are then calculated numerically using Eq. 33.

II.4.3 $J_{3}(\omega)=\eta\omega\Theta(\omega_{c}-\omega)$

Finally, we compute the memory kernel for the triangular spectral density, for which we obtain:

\displaystyle f(t)=\int_{0}^{\omega_{c}}\eta\omega e^{it(\Omega_{0}-\omega)}=\eta e^{it(\Omega_{0}-\omega_{c})}\frac{1-e^{i\omega_{c}t}+i\omega_{c}t}{t^{2}}\ .

(44)

The Laplace transform is:

\displaystyle\hat{f}(s)=-i\eta\omega_{c}-\eta(s-i\Omega_{0})\ln\left(\frac{s-i\Omega_{0}}{s-i\Omega_{0}+i\omega_{c}}\right)\ ,

(45)

and once more the (numerical) inverse Laplace transform of Eq. 24 yields $c_{1}(t)$ . Finally, the Lamb shift and decay rate are again computed numerically using Eq. 33.

III Approximation methods

In this section, we compute the excited state population predictions for the damped Jaynes-Cummings model of the previous section using three different Markovian master equations and using the TCL to second and fourth orders. In each case, we first provide a brief summary of the underlying theory of the corresponding master equation, both to assist the reader who may be unfamiliar with this theory and to establish our notation.

III.1 Coarse-Grained Lindblad Equation (CG-LE)

We follow the derivation of Ref. [23]. Let us choose a dimensionless, fixed and orthogonal system operator basis for $\mathcal{B}(\mathcal{H}_{S})$ as $\{S_{\alpha}\}_{\alpha=0}^{M}$ with $S_{0}=I$ and $M=d^{2}-1$ , where $d=\dim(\mathcal{H}_{S})$ and

\displaystyle\mathrm{Tr}(S_{\alpha}^{\dagger}S_{\beta})=\frac{1}{N_{\alpha}}\delta_{\alpha\beta}\ ,

(46)

where $N_{\alpha}$ is a normalization factor. We can then always write the system-bath interaction Hamiltonian in the following form:

\displaystyle H_{SB}=\sum_{\alpha}g_{\alpha}S_{\alpha}\otimes B_{\alpha}\ ,

(47)

where $\{B_{\alpha}\}$ are dimensionless bath operators and $\{g_{\alpha}\}$ are the coupling coefficients. In the interaction picture, with $U_{S}(t)=e^{-iH_{S}t}$ and $U_{B}(t)=e^{-iH_{B}t}$ , we obtain:


$\displaystyle\tilde{H}_{SB}$	$\displaystyle=\sum_{\alpha}g_{\alpha}S_{\alpha}(t)\otimes B_{\alpha}(t)$	(48a)
$\displaystyle S_{\alpha}(t)$	$\displaystyle=U_{S}^{\dagger}(t)S_{\alpha}U_{S}(t)=\sum_{\beta}p_{\alpha\beta}(t)S_{\beta}$	(48b)
$\displaystyle B_{\alpha}(t)$	$\displaystyle=U_{B}^{\dagger}(t)B_{\alpha}U_{B}(t)=\sum_{\beta}q_{\alpha\beta}(t)B_{\beta}\ ,$	(48c)

with initial conditions $p_{\alpha\beta}(0)=q_{\alpha\beta}(0)=\delta_{\alpha\beta}$ . The dynamics of the total system+bath are described by Eqs. 4 and 5. The system state is obtained by tracing over the bath and, assuming the initial state is factorized [ $\rho_{SB}(0)=\rho(0)\otimes\rho_{B}(0)$ ], can be represented as a completely positive quantum dynamical map:

\displaystyle\tilde{\rho}(t)=\mathrm{Tr}_{B}[\tilde{\rho}_{SB}(t)]=\sum_{i=0}^{M}\tilde{K}_{i}(t)\rho(0)\tilde{K}_{i}^{\dagger}(t)

(49)

where $\{\tilde{K}_{i}\}$ are Kraus operators in the interaction picture, which are defined as:

\displaystyle\tilde{K}_{i=\{\mu\nu\}}(t)=\sqrt{\lambda_{\mu}}\bra{\nu}\tilde{U}(t)\ket{\mu}\ ,

(50)

where the initial bath state is spectrally decomposed as $\rho_{B}(0)=\sum_{\mu}\lambda_{\mu}|{\mu}\rangle\!\langle\mu|$ . Using a standard Dyson series expansion of $\tilde{U}(t)$ similar to Eq. 6, we can write:


$\displaystyle\tilde{K}_{i}(t)$	$\displaystyle=\sqrt{\lambda_{\mu}}\delta_{\mu\nu}I+\sum_{n=1}^{\infty}K_{i}^{(n)}(t)$	(51a)
	$\displaystyle=\sum_{\alpha=0}^{M}b_{i\alpha}(t)S_{\alpha}=b_{i0}I+\sum_{\alpha=1}^{M}\sum_{n=1}^{\infty}b^{(n)}_{i\alpha}(t)S_{\alpha}\ ,$	(51b)

where $K_{i}^{(n)}(t)$ is the $n$ th order term in the Dyson series, and the second line is an expansion in the system operator basis. In the weak-coupling limit ( $\max_{\alpha}g_{\alpha}t\ll 1$ ), the higher-order terms in the expansion of $K_{i}^{(n)}(t)$ become negligible, i.e., $\|K_{i}^{(n+1)}\|\sim g_{\alpha}t\|K_{i}^{(n)}\|$ . Therefore, we can approximate the exact Kraus operators by truncating the expansion to first order. This yields:


$\displaystyle\tilde{K}_{i}^{(1)}(t)$	$\displaystyle=-i\sqrt{\lambda_{\mu}}\bra{\nu}\int_{0}^{t}dt_{1}\tilde{H}_{SB}(t_{1})\ket{\mu}$	(52a)
	$\displaystyle=-i\sqrt{\lambda_{\mu}}\sum_{\alpha}g_{\alpha}\int_{0}^{t}dt_{1}S_{\alpha}(t_{1})\bra{\nu}B_{\alpha}(t_{1})\ket{\mu}$	(52b)
	$\displaystyle=-it\sqrt{\lambda_{\mu}}\sum_{\alpha\beta\gamma}g_{\alpha}S_{\beta}\bra{\nu}B_{\gamma}\ket{\mu}\Gamma_{\alpha}^{\beta\gamma}(t)\ ,$	(52c)

where:

\displaystyle\Gamma_{\alpha}^{\beta\gamma}(t)\equiv\frac{1}{t}\int_{0}^{t}dt_{1}p_{\alpha\beta}(t_{1})q_{\alpha\gamma}(t_{1})\ .

(53)

Then, to match Eqs. 52c and 51b, we have:

\displaystyle b^{(1)}_{i\alpha}(t)=-it\sqrt{\lambda_{\mu}}\sum_{\alpha^{\prime}\alpha^{\prime\prime}}g_{\alpha^{\prime}}\bra{\nu}B_{\alpha^{\prime\prime}}\ket{\mu}\Gamma^{\alpha\alpha^{\prime\prime}}_{\alpha^{\prime}}(t)\ ,

(54)

and Eq. 51 implies that $b_{i0}=\sqrt{\lambda_{\mu}}\delta_{\mu\nu}$ . Next, we can construct the process matrix $\chi(t)$ (closely related to the Choi matrix), where

\chi_{\alpha\beta}(t)=\sum_{i=\mu\nu}b_{i\alpha}(t)b^{*}_{i\beta}(t)\ ,

(55)

and truncate it to lowest order ( $n\leq 1$ ) in the Dyson expansion as follows:

$\displaystyle\chi_{00}(t)$	$\displaystyle=\sum_{\mu}\lambda_{\mu}=1$	(56)
$\displaystyle\chi^{(1)}_{\alpha 0}(t)$	$\displaystyle=-it\sum_{\alpha^{\prime}\alpha^{\prime\prime}}g_{\alpha^{\prime}}\expectationvalue{B_{\alpha^{\prime\prime}}}_{B}\Gamma_{\alpha^{\prime}}^{\alpha\alpha^{\prime\prime}}(t),~{}~{}\alpha\geq 1$	(57)
$\displaystyle\chi^{(1)}_{\alpha\beta}(t)$	$\displaystyle=t^{2}\sum_{\alpha^{\prime}\alpha^{\prime\prime}\beta^{\prime}\beta^{\prime\prime}}g_{\alpha^{\prime}}g_{\beta^{\prime}}^{*}\expectationvalue{B^{\dagger}_{\beta^{\prime\prime}}B_{\alpha^{\prime\prime}}}_{B}\times$
	$\displaystyle\Gamma^{\alpha\alpha^{\prime\prime}}_{\alpha^{\prime}}(t)\left(\Gamma^{\beta\beta^{\prime\prime}}_{\beta^{\prime}}(t)\right)^{*},~{}~{}\alpha,\beta\geq 1\ ,$	(58)

where $\expectationvalue{X}_{B}\equiv\mathrm{Tr}(\rho_{B}X)$ . Let us also define

\expectationvalue{X}_{j}\equiv\frac{1}{\tau}\int_{j\tau}^{(j+1)\tau}X(t)dt\ ,

(59)

where we call $\tau$ the coarse-graining timescale. Note that $\chi_{\alpha\beta}(0)=\delta_{\alpha 0}\delta_{\beta 0}$ . Then:

\displaystyle a_{\alpha\beta}\equiv\expectationvalue{\dot{\chi}^{(1)}_{\alpha\beta}}_{0}=\frac{1}{\tau}(\chi^{(1)}_{\alpha\beta}(\tau)-\chi_{\alpha\beta}(0))=\frac{\chi^{(1)}_{\alpha\beta}(\tau)}{\tau}

(60)

unless $\alpha=\beta=0$ , in which case we have $\expectationvalue{\dot{\chi}_{00}}_{0}=0$ .

It can be shown [23] (see also [24]) that starting from Eq. 49, substituting the various expansions above, rearranging terms, and assuming that Eq. 60 can be extended to any interval $[t,t+\tau]$ (essentially an assumption of Markovianity), that one arrives at the Lindblad equation in the interaction picture:

	$\displaystyle\dot{\tilde{\rho}}(t)$	$\displaystyle=-i[H_{LS},\tilde{\rho}(t)]$
		$\displaystyle+\sum_{\alpha,\beta=1}^{M}a_{\alpha\beta}(S_{\alpha}\tilde{\rho}(t)S_{\beta}^{\dagger}-\frac{1}{2}\{S_{\beta}^{\dagger}S_{\alpha},\tilde{\rho}(t)\})\ ,$		(61)

where the Lamb shift is given by:


$\displaystyle H_{LS}$	$\displaystyle=\frac{i}{2}\sum_{\alpha}\expectationvalue{\dot{\chi}_{\alpha 0}}S_{\alpha}-\expectationvalue{\dot{\chi}_{\alpha 0}}^{*}S_{\alpha}^{\dagger}$	(62a)
	$\displaystyle=\frac{1}{2}\sum_{\alpha}\phi_{\alpha}S_{\alpha}+\phi_{\alpha}^{*}S_{\alpha}^{\dagger}$	(62b)

with

\displaystyle\phi_{\alpha}\equiv\sum_{\alpha^{\prime}\alpha^{\prime\prime}}g_{\alpha^{\prime}}\expectationvalue{B_{\alpha^{\prime\prime}}}_{B}\Gamma^{\alpha\alpha^{\prime\prime}}_{\alpha^{\prime}}(\tau)\ ,

(63)

and the decoherence rates are:

\displaystyle a_{\alpha\beta}=\tau\sum_{\alpha^{\prime}\alpha^{\prime\prime}\beta^{\prime}\beta^{\prime\prime}}g_{\alpha^{\prime}}g_{\beta^{\prime}}^{*}\expectationvalue{B^{\dagger}_{\beta^{\prime\prime}}B_{\alpha^{\prime\prime}}}_{B}\Gamma_{\alpha^{\prime}}^{\alpha\alpha^{\prime\prime}}(\tau)\Gamma_{\beta^{\prime}}^{\beta\beta^{\prime\prime}}(\tau)^{*}\ .

(64)

The choice of the coarse-graining timescale $\tau$ is crucial. It can be understood as a free optimization parameter, constrained by the inequality

\tau_{S}\ll\tau\ll 1/\omega_{c}\ ,

(65)

where $\tau_{S}$ is the timescale over which $\tilde{\rho}(t)$ changes, which arises from the replacement of

\langle\dot{\tilde{\rho}}(t)\rangle_{0}\equiv\frac{\tilde{\rho}(\tau)-{\rho}(0)}{\tau}

(66)

by $\dot{\tilde{\rho}}(t)$ in arriving at Eq. 61.

For the spin-boson model Eq. 10, we have $S_{+}=\sigma_{+}$ , $S_{-}=\sigma_{-}$ , $B_{k}=b_{k}$ or $b_{k}^{\dagger}$ . The system-bath interaction Hamiltonian in the interaction picture is given by Eq. 11. From Eqs. 48b and 48c, we obtain:


$\displaystyle p^{\pm\pm}(t)$	$\displaystyle=e^{\pm i\Omega_{0}t}\ ,\quad p^{\pm\mp}(t)=0$	(67a)
$\displaystyle q^{\pm\pm}_{kk^{\prime\prime}}(t)$	$\displaystyle=\delta_{kk^{\prime\prime}}e^{\mp i\omega_{k}t}\ ,\quad q^{\pm\mp}_{kk^{\prime\prime}}(t)=0\ ,$	(67b)

where the $+$ or $-$ superscripts indicate that the corresponding bath operator is $b_{k}$ or $b_{k}^{\dagger}$ , respectively.

Assuming that the initial state of the bath (cavity) is the zero-temperature vacuum state $\rho_{B}(0)=|{v}\rangle\!\langle v|$ , we obtain the standard bosonic expectation values:

	$\displaystyle\langle{b_{k}^{\dagger}b_{l}}\rangle_{B}$	$\displaystyle=\langle{b_{k}^{\dagger}}\rangle_{B}=\langle{b_{k}}\rangle_{B}=\langle{b_{k}^{\dagger}b_{l}^{\dagger}}\rangle_{B}=\langle{b_{k}b_{l}}\rangle_{B}=0$
	$\displaystyle\langle b_{k}b_{l}^{\dagger}\rangle_{B}$	$\displaystyle=\delta_{kl}.$		(68)

Then, by utilizing Eq. 53, we obtain a slightly modified expression for $\Gamma$ up to first order:

\displaystyle\Gamma_{\alpha,(\alpha^{\prime}k)}^{\beta,(\beta^{\prime}k^{\prime\prime})}(t)=\frac{1}{t}\int_{0}^{t}dt_{1}p^{\alpha\beta}(t_{1})q^{\alpha^{\prime}\beta^{\prime}}_{kk^{\prime\prime}}(t_{1})\ ,

(69)

where $\alpha,\beta,\alpha^{\prime},\beta^{\prime}\in\{+,-\}$ . Using Eq. 67, we have the following non-zero $\Gamma$ ’s:


$\displaystyle\Gamma_{+,(+k)}^{+,(+k)}(t)$	$\displaystyle=\frac{e^{i(\Omega_{0}-\omega_{k})t}-1}{i(\Omega_{0}-\omega_{k})t}$	(70a)
$\displaystyle\Gamma_{-,(-k)}^{-,(-k)}(t)$	$\displaystyle=\frac{e^{-i(\Omega_{0}-\omega_{k})t}-1}{-i(\Omega_{0}-\omega_{k})t}\ .$	(70b)

We also have a slightly modified expression for $b$ by using Eq. 54:

\displaystyle b_{\mu\nu,\alpha}

\displaystyle=-it\sqrt{\lambda_{\mu}}\sum_{(\alpha^{\prime}k^{\prime})}g^{\alpha^{\prime}}_{k^{\prime}}\bra{\mu}B^{\alpha^{\prime}}_{k^{\prime}}\ket{\nu}\Gamma^{\alpha,(\alpha^{\prime}k^{\prime})}_{\alpha,(\alpha^{\prime}k^{\prime})}\ ,

(71)

where $B^{\alpha}_{k}$ represents $b_{k}$ when $\alpha=+$ and $b_{k}^{\dagger}$ when $\alpha=-$ . The Lamb shift rates are given by Eq. 63, and vanish:

\displaystyle\phi_{\alpha}=\sum_{(\alpha^{\prime}k^{\prime})}g^{\alpha^{\prime}}_{k^{\prime}}\expectationvalue{B^{\alpha^{\prime}}_{k^{\prime}}}_{B}\Gamma^{\alpha,(\alpha^{\prime}k^{\prime})}_{\alpha,(\alpha^{\prime}k^{\prime})}(\tau)=0\ ,

(72)

since the expectation values of creation and annihilation operators between vacuum states vanish, as indicated by Section III.1. This result will be seen to undermine the quality of the CG-LE and C-LE when we perform a comparison with the exact results in Section IV.

The decoherence rates are:

	$\displaystyle a_{\alpha\beta}(\tau)$	$\displaystyle=\tau\sum_{(\alpha^{\prime}k^{\prime})(\beta^{\prime}l^{\prime})}g^{\alpha^{\prime}}_{k^{\prime}}g_{l^{\prime}}^{-\beta^{\prime}}\expectationvalue{B_{l^{\prime}}^{-\beta^{\prime}}B^{\alpha^{\prime}}_{k^{\prime}}}_{B}$
		$\displaystyle\qquad\times\Gamma_{\alpha,(\alpha^{\prime}k^{\prime})}^{\alpha,(\alpha^{\prime}k^{\prime})}(\tau)\Gamma_{\beta,(\beta^{\prime}l^{\prime})}^{\beta,(\beta^{\prime}l^{\prime})}(\tau)^{*}\ .$		(73)

Using Sections III.1 and 70, we obtain:

	$\displaystyle a_{++}(\tau)$	$\displaystyle=\sum_{k}\|g_{k}\|^{2}\tau\operatorname{sinc}^{2}\left(\frac{(\Omega_{0}+\omega_{k})\tau}{2}\right)\ ,$		(74)
	$\displaystyle a_{--}(\tau)$	$\displaystyle=\sum_{k}\|g_{k}\|^{2}\tau\operatorname{sinc}^{2}\left(\frac{(\Omega_{0}-\omega_{k})\tau}{2}\right)\ ,$		(75)

\displaystyle\gamma(\tau)

\displaystyle\equiv a_{--}(\tau)=\sum_{k}|g_{k}|^{2}\tau\operatorname{sinc}^{2}\left(\frac{(\Omega_{0}-\omega_{k})\tau}{2}\right)\ ,

(76)

where $\operatorname{sinc}(x)\equiv\sin(x)/x$ . Introducing the spectral density

J(\omega)=\sum_{k}|g_{k}|^{2}\delta(\omega-\omega_{k})\ ,

(77)

we can write this as

\gamma(\tau)=\int_{0}^{\infty}d\omega J(\omega)\tau\operatorname{sinc}^{2}\left(\frac{(\Omega_{0}-\omega)\tau}{2}\right)\ .

(78)

Let us now define

\bar{\delta}(x,y)\equiv\frac{1}{2\pi}y\operatorname{sinc}^{2}\left(\frac{xy}{2}\right)\ ,\quad y\geq 0\ .

(79)

This function behaves similarly to the Dirac- $\delta$ function:


		$\displaystyle\int_{-\infty}^{\infty}\bar{\delta}(x,y)dx=1$		(80a)
		$\displaystyle\lim_{y\to\infty}\bar{\delta}(x,y)=\delta(x)\ ,$		(80b)

i.e., it is sharply peaked at $x=0$ , and the peak becomes sharper as $y$ grows. The peak width is $\sim 1/y$ . We can thus also write

\gamma(\tau)=2\pi\int_{0}^{\infty}d\omega J(\omega)\bar{\delta}\left(\Omega_{0}-\omega,\tau\right)\ .

(81)

We will show below that this representation allows us to express the RWA-LE as the $\tau\to\infty$ limit of the CG-LE and C-LE results, as expected on general grounds [18].

Finally, we obtain the interaction picture Lindblad equation as:

\displaystyle\dot{\tilde{\rho}}(t)

\displaystyle=\gamma(\tau)\left(\sigma_{-}\tilde{\rho}(t)\sigma_{+}-\frac{1}{2}\{\sigma_{+}\sigma_{-},\tilde{\rho}(t)\}\right)\ .

(82)

Taking matrix elements, we find that the populations and coherences are decoupled. Solving for the excited state population and coherence, we obtain, respectively:


$\displaystyle\tilde{\rho}_{11}(t)$	$\displaystyle=\rho_{11}(0)e^{-\gamma(\tau)t}$	(83a)
$\displaystyle\tilde{\rho}_{01}(t)$	$\displaystyle=\rho_{01}(0)e^{-\frac{1}{2}\gamma(\tau)t}\ ,$	(83b)

which is to be contrasted with the exact solution given in Eq. 34. The coarse-graining time $\tau$ can be chosen to optimize the agreement with the exact solution.

III.1.1 $J_{1}(\omega)=|g|^{2}\delta(\omega-\omega_{c})$

For the impulse spectral density we find, using Eq. 81:


$\displaystyle\gamma(\tau)$	$\displaystyle=\|g\|^{2}\tau\operatorname{sinc}^{2}\left(\frac{(\Omega_{0}-\omega_{c})\tau}{2}\right)$	(84a)
	$\displaystyle=2\pi\|g\|^{2}\bar{\delta}\left(\Omega_{0}-\omega,\tau\right)\ .$	(84b)

III.1.2 $J_{2}(\omega)=\eta\omega e^{-\omega/\omega_{c}}$

For the Ohmic spectral density we find, using Eq. 81:


$\displaystyle\gamma(\tau)$	$\displaystyle=\eta\int_{0}^{\infty}\omega e^{-\omega/\omega_{c}}\tau\operatorname{sinc}^{2}\left(\frac{(\Omega_{0}-\omega)\tau}{2}\right)d\omega$	(85a)
	$\displaystyle=2\pi\eta\int_{0}^{\infty}\omega e^{-\omega/\omega_{c}}\bar{\delta}\left(\Omega_{0}-\omega,\tau\right)d\omega$	(85b)
	$\displaystyle=\frac{\eta}{\tau}e^{-\frac{\Omega_{0}}{\omega_{c}}}\left[\left(1-\frac{\Omega_{0}}{\omega_{c}}-i\Omega_{0}\tau\right)\operatorname{Ei}\left(\frac{\Omega_{0}}{\omega_{c}}+i\Omega_{0}\tau\right)+\right.$
	$\displaystyle\left(1-\frac{\Omega_{0}}{\omega_{c}}+i\Omega_{0}\tau\right)\operatorname{Ei}\left(\frac{\Omega_{0}}{\omega_{c}}-i\Omega_{0}\tau\right)$	(85c)
	$\displaystyle\left.+2\left(\frac{\Omega_{0}}{\omega_{c}}-1\right)\operatorname{Ei}\left(\frac{\Omega_{0}}{\omega_{c}}\right)\right]+\frac{2\eta}{\tau}(1-\cos\Omega_{0}\tau)\ ,$

where the last equality is derived in Section B.1.

III.1.3 $J_{3}(\omega)=\eta\omega\Theta(\omega_{c}-\omega)$

For the triangular spectral density we find, using Eq. 81:


$\displaystyle\gamma(\tau)$	$\displaystyle=\eta\int_{0}^{\omega_{c}}\omega\tau\operatorname{sinc}^{2}\left(\frac{(\Omega_{0}-\omega)\tau}{2}\right)d\omega$	(86a)
	$\displaystyle=2\pi\eta\int_{0}^{\omega_{c}}\omega\bar{\delta}\left(\Omega_{0}-\omega,\tau\right)d\omega$	(86b)
	$\displaystyle=\frac{2\eta}{\tau}\Bigg{(}\ln\left(\frac{\omega_{c}}{\Omega_{0}}-1\right)+\operatorname{Ci}\left(\tau\Omega_{0}\right)-\operatorname{Ci}\left(\tau\left(\omega_{c}-\Omega_{0}\right)\right)$
	$\displaystyle+\frac{\omega_{c}\left(\cos\left(\tau\left(\omega_{c}-\Omega_{0}\right)\right)-1\right)}{\omega_{c}-\Omega_{0}}-\cos\left(\tau\left(\omega_{c}-\Omega_{0}\right)\right)$
	$\displaystyle+\cos\left(\tau\Omega_{0}\right)\Bigg{)}+2\eta\Omega_{0}\left(\operatorname{Si}\left(\tau\left(\omega_{c}-\Omega_{0}\right)\right)+\operatorname{Si}\left(\tau\Omega_{0}\right)\right)\ ,$	(86c)

where the last equality is derived in Section B.2.

III.2 Cumulant Lindblad Equation C-LE

This section briefly reviews an alternative derivation of the Lindblad equation, based on a cumulant expansion [18]. Similarly to the CG-LE, the C-LE approach also uses a coarse-graining time scale that can be optimized to approximate the exact result. Despite using a rather different approach to deriving the Lindblad equation, we will show that for the problem we study in this work, the C-LE ultimately results in identical expressions for the master equation and its parameters (and hence also the solution, of course) as the CG-LE.

We start by writing the system-bath interaction Hamiltonian of Eq. 1 as:

\displaystyle H_{SB}=\lambda\sum_{\alpha}S_{\alpha}\otimes B_{\alpha}

(87)

where $S_{\alpha}$ and $B_{\alpha}$ are the system and bath operators, respectively (not necessarily Hermitian), and $\lambda$ is a dimensionless parameter to be used below for a series expansion, which we eventually set equal to $1$ . Note that unlike Eq. 47, the $S_{\alpha}$ are now not a basis, and the $B_{\alpha}$ have dimensions of energy since they include the coupling constants $g_{\alpha}$ . Assuming a factorized initial condition $\rho_{SB}(0)=\rho(0)\otimes\rho_{B}(0)$ , we associate a CPTP map $\Lambda_{\lambda}$ to the reduced density matrix of Eq. 7:

\displaystyle\tilde{\rho}(t)=\Lambda_{\lambda}(t)\rho(0)\ .

(88)

This CPTP map can be related to Eq. 6 by introducing superoperators $K^{(n)}$ , which collect terms with matching power of $\lambda$ :

\displaystyle\Lambda_{\lambda}(t)=\exp\left(\sum_{n=1}^{\infty}\lambda^{n}K^{(n)}(t)\right)

(89)

This is known as the cumulant expansion. The first-order term is then:

\displaystyle K^{(1)}(t)\rho(0)=-i\int_{0}^{t}ds\mathrm{Tr}_{B}\left([\tilde{H}(s),\rho_{SB}(0)]\right)\ ,

(90)

which can be eliminated by shifting the bath operators $B_{\alpha}$ , assuming stationarity, i.e., $[H_{B},\rho_{B}(0)]=0$ (see below and Appendix C). Moving on to the second order in $\lambda$ , we have:

\displaystyle K^{(2)}(t)

\displaystyle\rho(0)=-\int_{0}^{t}ds\int_{0}^{s}ds^{\prime}\mathrm{Tr}_{B}[\tilde{H}(s),[\tilde{H}(s^{\prime}),\rho_{SB}(0)]]\ ,

(91)

where the double commutator can be rearranged using:

$\displaystyle\mathrm{Tr}_{B}$	$\displaystyle[S_{\alpha}(s)\otimes B_{\alpha}(s),[S_{\beta}(s^{\prime})\otimes B_{\beta}(s^{\prime}),\rho(0)\otimes\rho_{B}(0)]]$
$\displaystyle=$	$\displaystyle[S_{\alpha}(s),S_{\beta}(s^{\prime})\rho(0)]\mathrm{Tr}(B_{\alpha}(s)B_{\beta}(s^{\prime})\rho_{B})-$
	$\displaystyle[S_{\alpha}(s),\rho(0)S_{\beta}(s^{\prime})]\mathrm{Tr}(B_{\beta}(s^{\prime})B_{\alpha}(s)\rho_{B})$	(92)

Introducing the bath correlation function

\displaystyle\mathcal{B}_{\alpha\beta}(s,s^{\prime})\equiv\mathrm{Tr}[B_{\alpha}(s)B_{\beta}(s^{\prime})\rho_{B}]\ ,

(93)

we have:

\displaystyle\mathcal{B}_{\alpha\beta}(s,s^{\prime})=\mathcal{B}_{\alpha\beta}^{*}(s^{\prime},s)\ .

(94)

To obtain a reduced form of the second-order cumulant in Eq. 91, it is useful to define a new variable that includes the double integration of the bath correlation function as follows:

\displaystyle\mathcal{B}_{\alpha\beta\omega}(t)\equiv\int_{0}^{t}ds\int_{0}^{s}ds^{\prime}e^{i\omega(s-s^{\prime})}\mathcal{B}_{\alpha\beta}(s,s^{\prime})

(95)

We introduce two more variables that will be associated with the Lamb shift and the decoherence rate:

\displaystyle Q_{\alpha\beta\omega}(t)\equiv\frac{1}{2i}\left(\mathcal{B}_{\alpha\beta\omega}(t)-\mathcal{B}^{*}_{\alpha\beta\omega}(t)\right)=\imaginary(\mathcal{B}_{\alpha\beta\omega}(t))\ ,

(96)

and


$\displaystyle b_{\alpha\beta\omega}(t)$	$\displaystyle\equiv\int_{0}^{t}ds\int_{0}^{t}ds^{\prime}e^{i\omega(s-s^{\prime})}\mathcal{B}_{\alpha\beta}(s,s^{\prime})$	(97a)
	$\displaystyle=\mathcal{B}_{\alpha\beta\omega}(t)+\mathcal{B}^{*}_{\alpha\beta\omega}(t)$	(97b)
	$\displaystyle=2\real(\mathcal{B}_{\alpha\beta\omega}(t))\ ,$	(97c)

where the equality in the second line is shown in Appendix D.

Explicitly, in this section, $\alpha\in\{+,-\}$ , and as in Eq. 11, in the interaction picture the system operators are $S_{\pm}(t)=\sigma_{\pm}e^{\pm i\Omega_{0}t}$ and the bath operators are $B_{+}(t)=\sum_{k}g_{k}e^{-i\omega_{k}t}b_{k},~{}B_{-}(t)=\sum_{k}g_{k}^{*}e^{i\omega_{k}t}b_{k}^{\dagger}$ . The initial bath state $\rho_{B}=|{v}\rangle\!\langle v|$ and the commutation rules Section III.1 give us just one non-zero bath correlation function:


$\displaystyle\mathcal{B}_{+-}(s,s^{\prime})$	$\displaystyle=\mathrm{Tr}[\rho_{B}B_{+}(s)B_{-}(s^{\prime})]$	(98a)
	$\displaystyle=\sum_{l,l^{\prime}}g_{l}g_{l^{\prime}}^{*}e^{-i(\omega_{l}s-\omega_{l^{\prime}}s^{\prime})}\mathrm{Tr}[\rho_{B}(0)b_{l}b_{l^{\prime}}^{\dagger}]$	(98b)
$\displaystyle\mathcal{B}_{-+}(s,s^{\prime})$	$\displaystyle=\mathcal{B}_{++}(s,s^{\prime})=\mathcal{B}_{--}(s,s^{\prime})=0\ .$	(98c)

The double commutator Section III.2 can now be simplified as follows:

$\displaystyle\sum_{\alpha,\beta}\mathrm{Tr}_{B}$	$\displaystyle[S_{\alpha}(s)\otimes B_{\alpha}(s),[S_{\beta}(s^{\prime})\otimes B_{\beta}(s^{\prime}),\rho(0)\otimes\rho_{B}(0)]]$
	$\displaystyle=e^{i\Omega_{0}(s-s^{\prime})}[\sigma_{+},\sigma_{-}\rho(0)]\mathcal{B}_{+-}(s,s^{\prime})$
	$\displaystyle\quad-e^{-i\Omega_{0}(s-s^{\prime})}[\sigma_{-},\rho(0)\sigma_{+}]\mathcal{B}_{+-}(s^{\prime},s)\ .$	(99)

Using Eq. 94, we have:


	$\displaystyle K^{(2)}(t)\rho(0)$		(100a)
	$\displaystyle\quad=-\mathcal{B}_{+-,\Omega_{0}}(t)[\sigma_{+},\sigma_{-}\rho(0)]+\mathcal{B}_{+-,\Omega_{0}}^{*}(t)[\sigma_{-},\rho(0)\sigma_{+}]$
	$\displaystyle\quad=-i\imaginary(\mathcal{B}_{+-,\Omega_{0}}(t))[\sigma_{+}\sigma_{-},\rho(0)]$		(100b)
	$\displaystyle\qquad+2\real(\mathcal{B}_{+-,\Omega_{0}}(t))\left(\sigma_{-}\rho(0)\sigma_{+}-\frac{1}{2}\{\sigma_{+}\sigma_{-},\rho(0)\}\right)\ .$

Hence, the second-order cumulant takes the form:

	$\displaystyle K^{(2)}(t)\rho(0)=-iQ_{+-,\Omega_{0}}(t)[\sigma_{+}\sigma_{-},\rho(0)]$
	$\displaystyle\quad+b_{+-,\Omega_{0}}(t)\left(\sigma_{-}\rho(0)\sigma_{+}-\frac{1}{2}\{\sigma_{+}\sigma_{-},\rho(0)\}\right)\ .$		(101)

Hence, the state in the interaction picture after the CP map in Eq. 88 using a truncation up to the second order in the cumulant expansion, is:


$\displaystyle\tilde{\rho}(t)$	$\displaystyle=\Lambda_{\lambda}(t)\rho(0)\approx\exp\left(\lambda^{2}K^{(2)}(t)\right)\rho(0)$	(102a)
	$\displaystyle={\rho}(0)-i\lambda^{2}Q_{+-,\Omega_{0}}(t)[\sigma_{+}\sigma_{-},\rho(t)]$	(102b)
	$\displaystyle\quad+\lambda^{2}b_{+-,\Omega_{0}}(t)\left(\sigma_{-}\rho(t)\sigma_{+}-\frac{1}{2}\{\sigma_{+}\sigma_{-},\rho(t)\}\right)\ .$

Now we use the coarse-graining method by averaging over the coarse-graining timescale $\tau$ as in Eq. 60:


$\displaystyle\langle\dot{b}_{+-,\Omega_{0}}(t)\rangle_{0}$	$\displaystyle=\frac{b_{+-,\Omega_{0}}(\tau)}{\tau}$	(103a)
$\displaystyle\langle\dot{Q}_{+-,\Omega_{0}}(t)\rangle_{0}$	$\displaystyle=\frac{Q_{+-,\Omega_{0}}(\tau)}{\tau}\ ,$	(103b)

which, when applied to Eq. 102b, yields:

		$\displaystyle\langle\dot{\tilde{\rho}}(t)\rangle_{0}=-i\lambda^{2}\langle\dot{Q}_{+-,\Omega_{0}}(t)\rangle_{0}[\sigma_{+}\sigma_{-},\rho(0)]$		(104)
		$\displaystyle+\lambda^{2}\langle\dot{b}_{+-,\Omega_{0}}(t)\rangle_{0}\left(\sigma_{-}\rho(0)\sigma_{+}-\frac{1}{2}\{\sigma_{+}\sigma_{-},\rho(0)\}\right)\ ,$

where we also used Eq. 66. Let us now define


$\displaystyle\mathcal{S}(\tau)$	$\displaystyle\equiv 2\langle\dot{Q}_{+-,\Omega_{0}}(t)\rangle_{0}=\frac{2\imaginary(\mathcal{B}_{+-,\Omega_{0}}(\tau))}{\tau}$	(105a)
$\displaystyle\gamma(\tau)$	$\displaystyle\equiv\langle\dot{b}_{+-,\Omega_{0}}(t)\rangle_{0}=\frac{2\real(\mathcal{B}_{+-,\Omega_{0}}(\tau))}{\tau}\ ,$	(105b)

where we used Eqs. 96 and 97c.

For a general $\rho_{B}(0)$ obtained by tracing out the system from Eq. 14 we find that $[\rho_{B}(0),H_{B}]\neq 0$ (as shown in Appendix E), which means that the bath correlation function $\mathcal{B}_{\alpha\beta}(s,s^{\prime})$ is not stationary. However, it is for the vacuum bath state $\rho_{B}(0)=|{v}\rangle\!\langle v|$ which we assume to be the case throughout, so we can write

\mathcal{B}_{+-}(s,s^{\prime})=\mathcal{B}_{+-}(s-s^{\prime})=\int_{0}^{\infty}d\omega J(\omega)e^{-i\omega(s-s^{\prime})}\ .

(106)

We then obtain, using Eq. 96:


$\displaystyle S(\tau)$	$\displaystyle=\frac{1}{\tau}\imaginary\int_{0}^{\tau}ds\int_{0}^{s}ds^{\prime}e^{i\omega(s-s^{\prime})}\int_{0}^{\infty}d\omega J(\omega)e^{-i\omega(s-s^{\prime})}$	(107a)
	$\displaystyle=\frac{\tau}{2}\imaginary\int_{0}^{\infty}d\omega J(\omega)=0\ ,$	(107b)

since the spectral density is real. I.e., just like in the CG-LE case [Eq. 72], the Lamb shift vanishes.

For the decay rate we now have, using Eq. 97b:


$\displaystyle\gamma(\tau)$	$\displaystyle=\frac{1}{\tau}\int_{0}^{\infty}d\omega J(\omega)\int_{0}^{\tau}ds\int_{0}^{\tau}ds^{\prime}e^{i(\Omega_{0}-\omega)(s-s^{\prime})}$	(108a)
	$\displaystyle=\int_{0}^{\infty}d\omega J(\omega)\tau\operatorname{sinc}^{2}\left(\frac{(\Omega_{0}-\omega)\tau}{2}\right)\ ,$	(108b)

which is identical to the CG-LE result, Eq. 78.

Moreover, similar to how we arrived at Eq. 61, assuming Markovianity in the sense that Eq. 104 can be extended to any interval $[t,t+\tau]$ we again arrive at the Lindblad equation in the interaction picture, after replacing $\langle\dot{\tilde{\rho}}(t)\rangle_{0}\mapsto\dot{\tilde{\rho}}(t)$ , and setting $\lambda=1$ . The form of this Lindblad equation is identical to Eq. 82. In particular, both the excited state population and the coherence are the same as in Eq. 83. Thus, the end results of the C-LE and CG-LE are identical for the model considered in this work.

III.3 Rotating Wave Approximation Lindblad Equation (RWA-LE)

The rotating wave approximation (RWA) drops the non-secular (off-diagonal) frequency terms which appear in the C-LE [see, e.g., Eq. 98b]. This approximation is based on the idea that the terms with $\omega\neq\omega^{\prime}$ are rapidly oscillating if $t\gg|\omega-\omega^{\prime}|^{-1}$ , which thus (roughly) average to zero. Since we will assume that $t\gg\tau_{B}$ , where $\tau_{B}$ is the bath correlation time (the time over which the bath correlation function decays), the former assumption is consistent provided we also assume that the Bohr frequency differences satisfy $\min_{\omega\neq\omega^{\prime}}|\omega-\omega^{\prime}|>1/\tau_{B}$ . Combining this with the weak coupling assumption, we obtain:

\displaystyle g\ll\frac{1}{\tau_{B}}<\min_{\omega\neq\omega^{\prime}}|\omega-\omega^{\prime}|\ .

(109)

By considering the weak coupling limit, taking the bath correlation timescale as the inverse of the cutoff frequency, and considering the Bohr frequencies $\{0,\pm\Omega_{0}\}$ , Eq. 109 becomes:

\displaystyle\eta\ll 1<\Omega_{0}/\omega_{c}\ .

(110)

Furthermore, the Born approximation states that for a sufficiently large bath, the composite state factorizes:

\displaystyle\tilde{\rho}_{SB}(t)\approx\tilde{\rho}(t)\otimes\rho_{B}.

(111)

Thus, up to second order in the Dyson series, the system state evolves according to:


$\displaystyle\dot{\tilde{\rho}}$	$\displaystyle=-\mathrm{Tr}_{B}[\tilde{H}(t),\int_{0}^{t}d\tau[\tilde{H}(t-\tau),\tilde{\rho}_{SB}(t-\tau)]]$	(112a)
	$\displaystyle=-\sum_{\alpha,\beta}\mathrm{Tr}_{B}[A_{\alpha}(t)\otimes B_{\alpha}(t),\int_{0}^{t}d\tau[A_{\beta}(t-\tau)$
	$\displaystyle\qquad\otimes B_{\beta}(t-\tau),\tilde{\rho}(t-\tau)\otimes\rho_{B}],$	(112b)

where in our case $\alpha,\beta\in\{+,-\}$ . It is useful to define the stationary (single-variable) bath correlation, a special case of Eq. 93:

\displaystyle\mathcal{B}_{+-}(t,t-\tau)=\int_{0}^{\infty}d\omega J(\omega)e^{-i\omega\tau}\equiv\mathcal{B}_{+-}(\tau)\ .

(113)

Now if we assume that $t\gg\tau_{B}$ , then $\tilde{\rho}(t-\tau)\approx\tilde{\rho}(t)$ . We discuss the limitations of this approximation in Appendix F, where we show that it can lead to an unbounded error.

Expanding the double commutator in terms of the bath correlation function, we obtain:

		$\displaystyle\sum_{\alpha,\beta}\mathrm{Tr}_{B}[S_{\alpha}(t)\otimes B_{\alpha},[S_{\beta}(t-\tau)\otimes B_{\beta}(t-\tau),\tilde{\rho}(t)\otimes\rho_{B}]]$
		$\displaystyle\quad=[S_{+}(t),S_{-}(t-\tau)\tilde{\rho}(t)]\mathcal{B}_{+-}(\tau)$
		$\displaystyle\quad\quad-[S_{-}(t),\tilde{\rho}(t)S_{+}(t-\tau)]\mathcal{B}_{+-}(-\tau)\ .$		(114)

Consequently, substituting Section III.3 back into Eq. 112, we can write:

		$\displaystyle\dot{\tilde{\rho}}(t)=-\int_{0}^{t}d\tau\mathcal{B}_{+-}(\tau)e^{i\Omega_{0}\tau}[\sigma_{+},\sigma_{-}\tilde{\rho}(t)]$
		$\displaystyle\quad+\int_{0}^{t}d\tau\mathcal{B}_{+-}(-\tau)e^{-i\Omega_{0}\tau}[\sigma_{-},\tilde{\rho}(t)\sigma_{+}]\ .$		(115)

To arrive at a Lindblad form, we complete the Markovian approximation by setting the upper limit of the integral to be $\infty$ . This is justified since the bath correlation decays rapidly to zero for $t\gg 1/\omega_{c}$ . Now, let


$\displaystyle\Gamma_{\alpha\beta}(\omega)$	$\displaystyle\equiv\int_{0}^{\infty}d\tau\mathcal{B}_{\alpha\beta}(\tau)e^{i\omega\tau}$	(116a)
	$\displaystyle=\int_{0}^{\infty}d\omega^{\prime}J(\omega^{\prime})\int_{0}^{\infty}d\tau e^{i(\omega-\omega^{\prime})\tau}$	(116b)
	$\displaystyle=\pi J(\omega)+i\int_{0}^{\infty}d\omega^{\prime}\ J(\omega^{\prime})\mathcal{P}\left(\frac{1}{\omega-\omega^{\prime}}\right)\ ,$	(116c)

where we used the identity

\int_{0}^{\infty}d\tau e^{ix\tau}=\pi\delta(x)+i\mathcal{P}\left(\frac{1}{x}\right)\ ,

(117)

and where the Cauchy principal value is defined as

\mathcal{P}\left(\frac{1}{x}\right)[f]=\lim_{\epsilon\to 0}\int_{-\epsilon}^{\epsilon}\frac{f(x)}{x}dx\ ,

(118)

for smooth functions $f$ with compact support on the real line $\mathbb{R}$ .

We show in Appendix G that taking the complex conjugate of $\Gamma_{\alpha\beta}$ yields:

\displaystyle\Gamma^{*}_{{\pm\mp}}(\omega)=\int_{0}^{\infty}d\tau\mathcal{B}_{{\pm\mp}}(-\tau)e^{-i\omega\tau}\ .

(119)

This simplifies Section III.3 into Lindblad form:


$\displaystyle\dot{\tilde{\rho}}(t)$	$\displaystyle=-\Gamma_{+-}(\Omega_{0})[\sigma_{+},\sigma_{-}\tilde{\rho}]+\Gamma_{+-}^{*}(\Omega_{0})[\sigma_{-},\tilde{\rho}\sigma_{+}]$	(120a)
	$\displaystyle=-i\imaginary[\Gamma_{+-}(\Omega_{0})][\sigma_{+}\sigma_{-},\rho(t)]$	(120b)
	$\displaystyle\quad+2\real[\Gamma_{+-}(\Omega_{0})]\left(\sigma_{-}\rho(t)\sigma_{+}-\frac{1}{2}\{\sigma_{+}\sigma_{-},\rho(t)\}\right)\ .$

This result has the same form as the exact Section II.3, but with a time-independent Lamb shift and decay rate:


$\displaystyle\mathcal{S}$	$\displaystyle=2\imaginary[\Gamma_{+-}(\Omega_{0})]=2\int_{0}^{\infty}d\omega^{\prime}\ J(\omega^{\prime})\mathcal{P}\left(\frac{1}{\Omega_{0}-\omega^{\prime}}\right)$	(121a)
$\displaystyle\gamma$	$\displaystyle=2\real[\Gamma_{+-}(\Omega_{0})]=2\pi J(\Omega_{0})\ .$	(121b)

This last result is consistent with the finding that the RWA-LE is the $\tau\to\infty$ limit of the C-LE [18], since it follows from Eqs. 80b and 81 that

\lim_{\tau\to\infty}\gamma(\tau)=2\pi\int_{0}^{\infty}d\omega J(\omega)\delta(\Omega_{0}-\omega)=2\pi J(\Omega_{0})\ .

(122)

The population and coherence are given by the $\tau\to\infty$ limit of Eq. 83, i.e.,


$\displaystyle\tilde{\rho}_{11}(t)$	$\displaystyle=\rho_{11}(0)e^{-\gamma t}$	(123a)
$\displaystyle\tilde{\rho}_{01}(t)$	$\displaystyle=\rho_{01}(0)e^{-\frac{1}{2}\gamma t}\ ,$	(123b)

We are now ready to present the results for the three spectral densities.

III.3.1 $J_{1}(\omega)=|g|^{2}\delta(\omega-\omega_{c})$

For the impulse spectral density we find, using Eq. 121:


$\displaystyle\mathcal{S}$	$\displaystyle=2\|g\|^{2}\mathcal{P}\left(\frac{1}{\Omega_{0}-\omega_{c}}\right)$	(124a)
$\displaystyle\gamma$	$\displaystyle=2\pi\|g\|^{2}\delta(\Omega_{0}-\omega_{c})\ .$	(124b)

This means that the decay rate either vanishes or is singular at $\Omega_{0}=\omega_{c}$ . Hence, the RWA-LE is unsuitable for describing the model with this spectral density.

III.3.2 $J_{2}(\omega)=\eta\omega e^{-\omega/\omega_{c}}$

We can immediately write down the decay rate as $\gamma=2\pi J_{2}(\Omega_{0})$ . However, the Cauchy principal value complicates the calculation of the Lamb shift, so we use a direct method instead.

For the Ohmic spectral density, the bath correlation in Eq. 113 takes the form:

\displaystyle\mathcal{B}_{+-}(\tau)=\eta\int_{0}^{\infty}\omega e^{-\omega/\omega_{c}}e^{-i\omega\tau}d\omega=\frac{\eta\omega_{c}^{2}}{(1+i\omega_{c}\tau)^{2}}\ .

(125)

The one-sided Fourier integral in Eq. 116a becomes:


		$\displaystyle\Gamma_{+-}(\Omega_{0})=\eta\omega_{c}\int_{0}^{\infty}d(\omega_{c}\tau)\frac{e^{i\Omega_{0}\tau}}{(1+i\omega_{c}\tau)^{2}}$		(126a)
		$\displaystyle\quad=-i\eta\omega_{c}+\eta\Omega_{0}e^{-\Omega_{0}/\omega_{c}}(\pi+i\operatorname{Ei}\left(\Omega_{0}/\omega_{c}\right))\ ,$		(126b)

where we derive the second equality in Section H.1.

Thus, the Lamb shift and the decay rate are:


$\displaystyle\mathcal{S}$	$\displaystyle=2J_{2}(\Omega_{0})\operatorname{Ei}\left(\frac{\Omega_{0}}{\omega_{c}}\right)-2\eta\omega_{c}$	(127a)
$\displaystyle\gamma$	$\displaystyle=2\pi J_{2}(\Omega_{0})\ .$	(127b)

III.3.3 $J_{3}(\omega)=\eta\omega\Theta(\omega_{c}-\omega)$

We can once more immediately write down the decay rate as $\gamma=2\pi J_{3}(\Omega_{0})$ , but a direct calculation is again advantageous for arriving at the form of the Lamb shift.

For the triangular spectral density, the bath correlation in Eq. 113 takes the form:

\displaystyle\mathcal{B}_{+-}(\tau)=\eta\int_{0}^{\omega_{c}}\omega e^{-i\omega\tau}d\omega=\eta\frac{e^{-i\omega_{c}\tau}(1+i\omega_{c}\tau)-1}{\tau^{2}}

(128)

The one-sided Fourier integral in Eq. 116a becomes:

\displaystyle\Gamma_{+-}(\Omega_{0})=\eta\omega_{c}\int_{0}^{\infty}d(\omega_{c}\tau)\frac{e^{-i\omega_{c}\tau}(1+i\omega_{c}\tau)-1}{(\omega_{c}\tau)^{2}}e^{i\Omega_{0}\tau}

(129)

Thus, as we derive in Section H.2, the Lamb shift and the decay rate are:


$\displaystyle\mathcal{S}$	$\displaystyle=-2\Omega_{0}\ln\left\|\frac{\omega_{c}}{\Omega_{0}}-1\right\|-2\eta\omega_{c}$	(130a)
$\displaystyle\gamma$	$\displaystyle=2\pi J_{3}(\Omega_{0})\ .$	(130b)

Note that Eq. 110 requires $\Omega_{0}>\omega_{c}$ , but in this case it follows from Eq. 130b that $\gamma$ vanishes. This breakdown of the validity conditions of the Markov approximation, along with the issue of the potentially unbounded approximation error discussed in Appendix F, highlights that the RWA-LE has limited validity for the model we study here. Our simulation results reinforce these conclusions, as shown in Section IV.

III.4 TCL

In this section we briefly review the time-convolutionless formalism, closely following the presentation of [2] while adding a few pertinent details.

We can rewrite the Liouville equation in Eq. 3 as:

\displaystyle\frac{d\tilde{\rho}_{SB}}{dt}\equiv\lambda\mathcal{L}\tilde{\rho}_{SB}\ \ ,\qquad\mathcal{L}\equiv[\tilde{H}_{SB},\cdot]\ ,

(131)

where we introduce the Feshbach projection superoperator $\mathcal{P}$ via

\displaystyle\mathcal{P}\rho_{SB}\equiv\mathrm{Tr}_{B}(\rho)_{SB}\otimes\rho_{B}\ ,

(132)

and its orthogonal complement $\mathcal{Q}=\mathcal{I}-\mathcal{P}$ . For an arbitrary operator $A$ , we define:

\displaystyle\hat{A}\equiv\mathcal{P}A\ ,\qquad\bar{A}\equiv\mathcal{Q}A\ ,

(133)

which leads via Eq. 131 to:


$\displaystyle\partial_{t}\hat{\rho}_{SB}$	$\displaystyle=\lambda\hat{\mathcal{L}}\hat{\rho}_{SB}+\lambda\hat{\mathcal{L}}\bar{\rho}_{SB}$	(134a)
$\displaystyle\partial_{t}\bar{\rho}_{SB}$	$\displaystyle=\lambda\bar{\mathcal{L}}\hat{\rho}_{SB}+\lambda\bar{\mathcal{L}}\bar{\rho}_{SB}\ .$	(134b)

The second equation has a solution given by:

\displaystyle\bar{\rho}_{SB}(t)=\mathcal{G}(t,0)\bar{\rho}_{SB}(0)+\lambda\int_{0}^{t}\mathcal{G}(t,t^{\prime})\bar{\mathcal{L}}(t^{\prime})\bar{\rho}_{SB}(t^{\prime})dt^{\prime}\ ,

(135)

where:

\displaystyle\mathcal{G}(t,0)\equiv T_{+}e^{\lambda\int_{0}^{t}\bar{\mathcal{L}}(t^{\prime})dt^{\prime}}\ .

(136)

From Eq. 131, we obtain:

\displaystyle\tilde{\rho}_{SB}(t)=\mathcal{U}_{+}(t,t^{\prime})\tilde{\rho}_{SB}(t^{\prime})\ ,

(137)

where

\displaystyle\mathcal{U}_{+}(t,t^{\prime})=T_{+}e^{\lambda\int_{t^{\prime}}^{t}\bar{\mathcal{L}}(s)ds}\ .

(138)

We can back-propagate the system state as:

\displaystyle\tilde{\rho}_{SB}(t^{\prime})=\mathcal{U}_{-}(t^{\prime},t)\tilde{\rho}_{SB}(t^{\prime})\ ,

(139)

where $\mathcal{U}_{-}(t^{\prime},t)=[\mathcal{U}_{+}(t,t^{\prime})]^{-1}$ . When we substitute Eq. 139 back into Eq. 135, it yields:


$\displaystyle\bar{\rho}_{SB}(t)$	$\displaystyle=\mathcal{G}(t,0)\bar{\rho}_{SB}(0)+\Sigma(t)\rho_{SB}(t)$	(140a)
$\displaystyle\Sigma(t)$	$\displaystyle\equiv\lambda\int_{0}^{t}\mathcal{G}(t,t^{\prime})\bar{\mathcal{L}}(t^{\prime})\hat{\mathcal{U}}_{-}(t^{\prime},t)dt^{\prime}\ .$	(140b)

Assuming that $\mathcal{I}-\Sigma$ is invertible and using Eq. 134, we arrive at:


$\displaystyle\partial_{t}\hat{\rho}_{SB}(t)$	$\displaystyle=\mathcal{J}(t)\bar{\rho}_{SB}(0)+\mathcal{K}(t)\hat{\rho}_{SB}(t)$	(141a)
$\displaystyle\mathcal{J}(t)$	$\displaystyle\equiv\lambda\hat{\mathcal{L}}(t)[\mathcal{I}-\Sigma]^{-1}\mathcal{G}(t,0)\mathcal{Q}$	(141b)
$\displaystyle\mathcal{K}(t)$	$\displaystyle\equiv\lambda\hat{\mathcal{L}}[\mathcal{I}-\Sigma]^{-1}\mathcal{P}\ .$	(141c)

If the inhomogeneity $\mathcal{J}(t)$ vanishes (as is the case for a factorized initial system-bath state), then the resulting master equation is time-local and known as the time-convolutionless (TCL) master equation. By expressing $\mathcal{I}-\Sigma$ as a geometric series, we obtain:

\displaystyle[\mathcal{I}-\Sigma]^{-1}=\sum_{n=0}^{\infty}\Sigma^{n}(t)\ .

(142)

This allows us to write the TCL generator as an expansion in powers of $\lambda$ :

\displaystyle\mathcal{K}(t)=\lambda\hat{\mathcal{L}}(t)\left(\sum_{n=0}^{\infty}\Sigma^{n}(t)\right)\mathcal{P}=\sum_{n=1}^{\infty}\lambda^{n}\mathcal{K}_{n}(t)\ .

(143)

For Gaussian baths, like the one considered here, the odd-order terms vanish, i.e., $\mathcal{K}_{2n+1}=0$ , which follows from the vanishing multi-time bath correlation function utilizing Wick’s theorem [25, 26]. In general, the $n$ ’th order term is given by [27, 28]:

\displaystyle\mathcal{K}_{n}(t)=\int_{0}^{t}\!\!\!dt_{1}\int_{0}^{t}\!\!\!dt_{2}\cdots\int_{0}^{t_{n-2}}dt_{n-1}\langle\mathcal{L}(t)\mathcal{L}(t_{1})...\mathcal{L}(t_{n-1})\rangle_{\text{oc}}\ ,

(144)

where the ordered cumulants are defined as:

	$\displaystyle\langle\mathcal{L}(t)$	$\displaystyle\mathcal{L}(t_{1})...\mathcal{L}(t_{n-1})\rangle_{\text{oc}}=$
		$\displaystyle\sum(-1)^{q}\mathcal{P}\mathcal{L}(t)\cdots\mathcal{L}(t_{i})\mathcal{P}\mathcal{L}(t_{j})\cdots\mathcal{L}(t_{m})\mathcal{P}\ ,$		(145)

where the sum is over all possible arrangements of $q$ $\mathcal{P}$ ’s and $n$ $\mathcal{L}$ ’s such that there is at least one $\mathcal{L}$ between $\mathcal{P}$ ’s and there is a time ordering $t\geq...\geq t_{i}\geq t_{j}...\geq t_{m}$ .

Since Eq. 31 is already time-local, we can relate it to the TCL generator $\mathcal{K}(t)$ via [16]:

\displaystyle\mathcal{K}_{S}(t)\rho(t)=\mathrm{Tr}_{B}[\mathcal{K}(t)(\rho(t)\otimes\rho_{B})]\ .

(146)

This leads exactly to the ansatz given in Section II.3, i.e., the TCL master equation is given by

	$\displaystyle\dot{\tilde{\rho}}$	$\displaystyle=\mathcal{K}_{S}(t)\tilde{\rho}=-\frac{i}{2}\mathcal{S}(t)[\sigma_{+}\sigma_{-},\tilde{\rho}(t)]$
		$\displaystyle\qquad+\gamma(t)\left(\sigma_{-}\tilde{\rho}(t)\sigma_{+}-\frac{1}{2}\{\sigma_{+}\sigma_{-},\tilde{\rho}(t)\}\right)\ .$		(147)

The solution is given by Eq. 34, to the desired order in perturbation theory, i.e.:


$\displaystyle\tilde{\rho}_{11}(t)$	$\displaystyle=\rho_{11}(0)\exp{-\int_{0}^{t}\gamma^{(n)}(t^{\prime})dt^{\prime}}$	(148a)
$\displaystyle\tilde{\rho}_{01}(t)$	$\displaystyle=\rho_{01}(0)\exp{\frac{1}{2}\int_{0}^{t}\left(i\mathcal{S}^{(n)}(t^{\prime})-\gamma^{(n)}(t^{\prime})\right)dt^{\prime}}\ ,$	(148b)

where $n=2$ for TCL2, $n=4$ for TCL4, etc.

To find the perturbative expansion of the Lamb shift and the decay rate, we observe that $\sigma_{+}$ is an eigenoperator of the generator:

\displaystyle\mathcal{K}_{S}(t)\sigma_{+}=-\frac{1}{2}(\gamma(t)+iS(t))\sigma_{+}\ .

(149)

Using the superoperator $\mathcal{L}(t)=-i[H_{SB}(t),\cdot]$ , we can verify [29] that $\sigma_{+}\otimes\rho_{B}$ is an eigenoperator of $\mathcal{L}(t)\mathcal{L}(t_{1})$ with eigenvalue $-f(t-t_{1})$ , where $f(t)$ is defined in Eq. 22. Substituting Eq. 143 into Eq. 149, we obtain:

	$\displaystyle\mathcal{K}_{S}(t)\sigma_{+}=\sum_{n=1}^{\infty}(-\lambda^{2})^{n}\int_{t_{0}}^{t}dt_{1}\int_{t_{0}}^{t_{1}}dt_{2}\cdots\int_{t_{0}}^{t_{2n-2}}dt_{2n-1}$
	$\displaystyle\quad\langle f(t-t_{1})f(t_{2}-t_{3})\cdots f(t_{2n-2}-t_{2n-1})\rangle_{\text{oc}}\sigma_{+}\ .$		(150)

Here, the terms in the summation follow the rule of considering all possible arrangements of the memory kernel $f(t_{i}-t_{j})$ . It is important to ensure that the times are properly ordered as $t_{m}\leq...\leq t_{i}\leq t_{j}\leq...\leq t$ . Consequently, we can expand $\gamma(t)$ and $S(t)$ as:

\displaystyle\gamma(t)=\sum_{n=1}^{\infty}\lambda^{2n}\gamma_{2n}(t)\ ,\quad\mathcal{S}(t)=\sum_{n=1}^{\infty}\lambda^{2n}\mathcal{S}_{2n}(t)\ .

(151)

In particular, if we define the function

\displaystyle Z(t,t^{\prime})\equiv\int_{0}^{t}dt_{1}f(t^{\prime}-t_{1})

we obtain, as shown in Ref. [16], the following expressions


$\displaystyle\gamma_{2}(t)+iS_{2}(t)$	$\displaystyle=2\int_{0}^{t}dt_{1}f(t-t_{1})=2\lim_{t^{\prime}\rightarrow t}Z(t,t^{\prime})$	(152a)
$\displaystyle\gamma_{4}(t)+iS_{4}(t)$	$\displaystyle=2\int_{0}^{t}dt_{1}\int_{0}^{t_{1}}dt_{2}\int_{0}^{t_{2}}dt_{3}$	(152b)
	$\displaystyle\!\!\!\!\!\left[f(t-t_{2})f(t_{1}-t_{3})+f(t-t_{3})f(t_{1}-t_{2})\right]\ .$

Thus,


		$\displaystyle\mathcal{S}_{2}(t)[\gamma_{2}(t)]=2\imaginary[\real]\lim_{t^{\prime}\rightarrow t}Z(t,t^{\prime})$		(153a)
		$\displaystyle\mathcal{S}_{4}(t)[\gamma_{4}(t)]=2\imaginary[\real]\int_{0}^{t}dt_{1}\int_{0}^{t_{1}}dt_{2}\left(f(t-t_{2})Z(t_{1},t_{2})\right.$
		$\displaystyle\qquad\left.+f(t_{1}-t_{2})Z(t,t_{2})\right)\ ,$		(153b)

where the imaginary and real parts are taken for $\mathcal{S}$ and $\gamma$ , respectively. The TCL2 and TCL4 cases correspond, respectively, to the following substitutions into Eq. 148:


		$\displaystyle\text{TCL2:}\begin{cases}\mathcal{S}^{(2)}(t)=\mathcal{S}_{2}(t)\\ \gamma^{(2)}(t)=\gamma_{2}(t)\end{cases}$		(154a)
		$\displaystyle\text{TCL4:}\begin{cases}\mathcal{S}^{(4)}(t)=\mathcal{S}_{2}(t)+\mathcal{S}_{4}(t)\\ g^{(4)}(t)=\gamma_{2}(t)+\gamma_{4}(t)\end{cases}$		(154b)

This follows from Eq. 151 with $\lambda=1$ (recall that $\lambda$ is a formal expansion parameter).

Next, we consider our three spectral densities.

III.4.1 $J_{1}=|g|^{2}\delta(\omega-\omega_{c})$

We use the integral of the shifted memory kernel from Eq. 35:

\displaystyle Z(t,t_{1})=\frac{2|g|^{2}\sin\frac{(\Omega_{0}-\omega_{c})t_{1}}{2}}{\Omega_{0}-\omega_{c}}e^{i\frac{(\Omega_{0}-\omega_{c})(2t-t_{1})}{2}}\ .

(155)

The second-order Lamb shift and decay rate are, using Eq. 153a:


$\displaystyle\mathcal{S}_{2}(t)$	$\displaystyle=2\|g\|^{2}\frac{1-\cos((\Omega_{0}-\omega_{c})t)}{\Omega_{0}-\omega_{c}}$	(156a)
$\displaystyle\gamma_{2}(t)$	$\displaystyle=2\|g\|^{2}t\operatorname{sinc}[(\Omega_{0}-\omega_{c})t)]\ .$	(156b)

The fourth-order Lamb shift and decay rate are, using Eq. 153:


$\displaystyle\mathcal{S}_{4}(t)$	$\displaystyle=4\|g\|^{4}\frac{-\sin^{2}(\Omega_{0}-\omega_{c})t+(\Omega_{0}-\omega_{c})t\sin(\Omega_{0}-\omega_{c})t}{(\Omega_{0}-\omega_{c})^{3}}$	(157a)
$\displaystyle\gamma_{4}(t)$	$\displaystyle=2\|g\|^{4}\frac{2(\Omega_{0}-\omega_{c})t\cos(\Omega_{0}-\omega_{c})t-\sin 2(\Omega_{0}-\omega_{c})t}{(\Omega_{0}-\omega_{c})^{3}}$	(157b)

III.4.2 $J_{2}(\omega)=\eta\omega e^{-\omega/\omega_{c}}$

Again we use the integral of the shifted memory kernel from Eq. 41:

$\displaystyle Z(t,t_{1})$	$\displaystyle=\frac{i\eta\omega_{c}e^{i\Omega_{0}t}}{1+i\omega_{c}t}-\frac{i\eta\omega_{c}e^{i\Omega_{0}(t-t_{1})}}{1+i\omega_{c}(t-t_{1})}$
	$\displaystyle\quad+iJ(\Omega_{0})\Big{(}\operatorname{Ei}\left(\frac{\Omega_{0}}{\omega_{c}}(1+i\omega_{c}(t-t_{1}))\right)$
	$\displaystyle\qquad-\operatorname{Ei}\left(\frac{\Omega_{0}}{\omega_{c}}(1+i\omega_{c}t)\right)\Big{)}\ .$	(158)

The second-order Lamb shift and decay rate are, using Eq. 153a:


$\displaystyle\mathcal{S}_{2}(t)$	$\displaystyle=-2\eta\omega_{c}+2\eta\omega_{c}\frac{\cos(\Omega_{0}t)+\omega_{c}t\sin(\Omega_{0}t)}{1+(\omega_{c}t)^{2}}$
	$\displaystyle-2J_{2}(\Omega_{0})\left(\real\left[\operatorname{Ei}\left(\frac{\Omega_{0}}{\omega_{c}}(1+i\omega_{c}t)\right)\right]-\operatorname{Ei}\left(\frac{\Omega_{0}}{\omega_{c}}\right)\right)$	(159a)
$\displaystyle\gamma_{2}(t)$	$\displaystyle=2\eta\omega_{c}\frac{\omega_{c}t\cos(\Omega_{0}t)-\sin(\Omega_{0}t)}{1+(\omega_{c}t)^{2}}$
	$\displaystyle\quad+2J_{2}(\Omega_{0})\imaginary\left[\operatorname{Ei}\left(\frac{\Omega_{0}}{\omega_{c}}(1+i\omega_{c}t)\right)\right]\ .$	(159b)

The Markovian decay rate is obtained in the long time limit. Using $\lim_{x\rightarrow\infty}\operatorname{Ei}(ix+1)=i\pi$ , we obtain:

\displaystyle\gamma_{M}=\lim_{t\rightarrow\infty}\gamma_{2}(t)=2\pi\eta\Omega_{0}e^{-\frac{\Omega_{0}}{\omega_{c}}}=2\pi J(\Omega_{0})\ ,

(160)

in agreement with Eq. 127b.

For the fourth-order decay rate, Eq. 153 does not admit a closed-form solution and needs to be evaluated numerically.

III.4.3 $J_{3}(\omega)=\eta\omega\Theta(\omega_{c}-\omega)$

We use the shifted integral of the memory kernel in Eq. 44:

$\displaystyle Z(t,t_{1})$	$\displaystyle=\frac{\eta}{t(t-t_{1})}\left(e^{i(\Omega_{0}-\omega_{c})(t-t_{1})}t-e^{i\Omega_{0}(t-t_{1})}t\right.$
	$\displaystyle\left.+(e^{i\Omega_{0}t}-e^{i(\Omega_{0}-\omega_{c})t})(t-t_{1})\right)$	(161)
	$\displaystyle+i\eta\Omega_{0}\left[\operatorname{Ei}(i(\Omega_{0}-\omega_{c})t)-\operatorname{Ei}(i\Omega_{0}t)\right.$
	$\displaystyle\left.-\operatorname{Ei}(i(\Omega_{0}-\omega_{c})(t-t_{1}))+\operatorname{Ei}(i\Omega_{0}(t-t_{1}))\right]\ .$

Thus:


$\displaystyle\mathcal{S}_{2}(t)$	$\displaystyle=-2i\eta\Big{(}\frac{\sin(\Omega_{0}t)-\sin((\Omega_{0}-\omega_{c})t)}{t}$
	$\displaystyle\quad-\Omega_{0}\ln\|1-\frac{\omega_{c}}{\Omega_{0}}\|-\omega_{c}$
	$\displaystyle\quad+\real(\Omega_{0}\left(\operatorname{Ei}(i(\Omega_{0}-\omega_{c})t)-\operatorname{Ei}(i\Omega_{0}t))\right)\Big{)}$	(162a)
$\displaystyle\gamma_{2}(t)$	$\displaystyle=2\eta\Big{(}\frac{\cos(\Omega_{0}t)-\cos((\Omega_{0}-\omega_{c})t)}{t}$
	$\displaystyle\quad-\pi\Omega_{0}\Theta(\omega_{c}-\Omega_{0})$
	$\displaystyle\quad+\Omega_{0}\imaginary\left[\operatorname{Ei}(i\Omega_{0}t)-\operatorname{Ei}(i(\Omega_{0}-\omega_{c})t)\right]\Big{)}\ .$	(162b)

We recover the Markov approximation in the large $t$ limit:

\displaystyle\lim_{t\rightarrow\infty}\gamma_{2}(t)=2\pi\eta J_{3}(\Omega_{0})\ .

(163)

While the TCL result Eq. 162 does not entirely match the exact result, it does describe an oscillatory behavior similar to the exact solution, which is entirely absent in the Markovian limit. A similar phenomenon has been described in the context of non-equilibrium dynamics in Ref. [30].

Moreover, recall that the Markovian rate vanishes when $\Omega_{0}>\omega_{c}$ [Eq. 130b], resulting in the absence of decay. In contrast, when $\Omega_{0}>\omega_{c}$ we find:

\displaystyle\int_{0}^{\infty}\gamma_{2}(t^{\prime})dt^{\prime}=2\eta\left(\frac{\omega_{c}}{\Omega_{0}-\omega_{c}}+\ln\left(\frac{\Omega_{0}-\omega_{c}}{\Omega_{0}}\right)\right)\ ,

(164)

so that the asymptotic limit for the population in the TCL2 approximation (which coincides with the Redfield equation), $e^{-\int_{0}^{\infty}\gamma_{2}(t^{\prime})dt^{\prime}}$ , is non-zero. This qualitatively recovers the non-zero decay property of the exact solution.

For the fourth-order decay rate, Eq. 153 again does not admit a closed-form solution and needs to be evaluated numerically.

		exact	C-LE/CG-LE	RWA-LE	TCL2	TCL4
population $\rho_{11}$		Eq. 34a	Eq. 83a	Eq. 123a	Eqs. 148a and 154a	Eqs. 148a and 154b
coherence $\rho_{01}$		Eq. 34b	Eq. 83b	Eq. 123b	Eqs. 148b and 154a	Eqs. 148b and 154b
Lamb shift $\mathcal{S}$	$J_{1}$	Eq. 40a	$0$	Eq. 124a	Eq. 156a	Eq. 157a
	$J_{2}$		$0$	Eq. 127a	Eq. 159
	$J_{3}$		$0$	Eq. 130a	Eq. 162
decay rate $\gamma$	$J_{1}$	Eq. 40b	Eq. 84	Eq. 124b	Eq. 156b	Eq. 157b
	$J_{2}$		Eq. 85	Eq. 127b	Eq. 159
	$J_{3}$		Eq. 86	Eq. 130b	Eq. 162

Table 1: Analytical results. The equations listed in the table are explicit analytical results. The results for empty cells are obtained numerically.

IV Analysis

In this section, we compare the exact solution for the excited state population $\rho_{11}(t)$ and the coherence to the results of the various approximation schemes we described above. We also compare the exact Lamb shift and decay rate $\gamma(t)$ to the corresponding quantities predicted by these approximation schemes. We do this for all three spectral densities where possible, for the CG-LE/C-LE, RWA-LE, TCL2 (Redfield), and TCL4. These results are summarized with corresponding equation references in Table 1.

IV.1 Computation of the optimal coarse graining time in CG-LE

Before presenting the comparison to the exact results, we first explain our methodology for optimizing the coarse-graining time $\tau$ within the framework of the CG-LE, since in the ensuing comparison we use the optimal $\tau$ values. Recall that the coarse-graining time needs to satisfy the condition $\omega_{c}\tau\ll 1$ [Eq. 65].

To determine the appropriate coarse-graining timescale $\tau$ , we minimize a metric that quantifies the deviation from the exact solution. We employ the integrated trace-norm distance norm [18] as our chosen metric, defined as:

\displaystyle\mathcal{D}_{[0,T]}\equiv\frac{1}{2T}\int_{0}^{T}dt\|\rho_{\text{exact}}(t)-\rho_{\text{approx}}(t)\|_{1}\ ,

(165)

where $\|M\|_{1}\equiv\mathrm{Tr}\sqrt{M^{\dagger}M}$ . Since $M=\rho_{\text{exact}}-\rho_{\text{approx}}$ is Hermitian and traceless, it can be written as $\left(\begin{array}[]{cc}a&b\\ b^{*}&-a\\ \end{array}\right)$ in the qubit case, so that $\mathrm{Tr}\sqrt{M^{2}}=2\sqrt{a^{2}+|b|^{2}}$ , and we can simplify the integrand as follows:

	$\displaystyle\\|\rho_{\text{exact}}-\rho_{\text{approx}}\\|_{1}$		(166)
	$\displaystyle=2\sqrt{(\rho_{11,\text{exact}}-\rho_{11,\text{approx}})^{2}+\|\rho_{01,\text{exact}}-\rho_{01,\text{approx}}\|^{2}}\ .$

In Fig. 1(a), we illustrate the distance $\mathcal{D}$ between the exact solution and both the CG-LE and the RWA-LE as a function of the dimensionless quantity $\omega_{c}\tau$ , where we vary the coarse-graining time $\tau$ . In this example, we use the parameters $\eta=1$ and $\Omega_{0}/\omega_{c}=1$ , and the numerical integration is performed up to a total time of $\omega_{c}T=100$ . This upper limit is justified by the observation that the distance is minimized at relatively short times. For example, in Fig. 1(a) the minimum occurs at $\omega_{c}\tau=0.50139$ , thus satisfying Eq. 65. Notably, the coarse-graining solution outperforms the RWA-LE solution, consistent with the findings reported in [18].

In Fig. 1(b), we present the minimum coarse-graining time as a function of the qubit frequency $\Omega_{0}/\omega_{c}$ for $\eta=1$ across the three different spectral models. The coarse-graining times used in the later figures are indicated by the red dots in the plot. We selected $\Omega_{0}/\omega_{c}\in[0.15,2]$ to account for a range of small to large qubit frequencies.

The most notable conclusion from Fig. 1(b) is that condition $\omega_{c}\tau\ll 1$ cannot be satisfied for the impulse ( $J_{1}$ ) and triangular ( $J_{3}$ ) spectral densities within the range of $\Omega_{0}/\omega_{c}$ values shown. However, it is satisfied for the Ohmic spectral density ( $J_{2}$ ), the most physically relevant of the three. We thus expect the CG-LE to perform poorly for $J_{1}$ and $J_{3}$ , but to perform relatively well for $J_{2}$ . These expectations are borne out in our results below.

Note further that for the Ohmic spectral density, the coarse-graining timescale increases as the qubit frequency decreases, tending to the RWA-LE solution in line with Eq. 121. Conversely, for the spectral density $J_{1}=|g|^{2}\delta(\omega-\omega_{c})$ , characterized by strong non-Markovian behavior and oscillations that prevent the fitting of a Markovian exponential decay, the coarse-graining time diverges at $\omega-\omega_{c}$ to attempt to fit the exact solution.

We remark that to ensure the validity of RWA-LE, the weak coupling limit Eq. 110 needs to be satisfied, a point we comment on in more detail below.

IV.2 Exact solution vs TCL and Markov approximations for the Ohmic spectral density

We now present the results of a comparison between the exact solution and the different master equations for the Ohmic spectral density, the most physically interesting of the three densities. Our general expectations are that the TCL approximations will capture some of the aspects of the non-Markovian dynamics and will thus outperform the two Markovian master equations, and that the CG-LE will outperform the RWA-LE due to the ability to optimize the coarse-graining timescale $\tau$ in the former. This expectation depends on the validity condition $\omega_{c}\tau\ll 1$ not being violated, as explained in the previous section. Indeed, we find that this holds for all the examples we consider in this section, at least in the sense that we find $\omega_{c}\tau<1$ in all cases.

For the initial condition, we consider a simple scenario where $c_{1}(0)=c_{0}(0)=\frac{1}{2}$ and $c_{k}(0)=0$ . In this setup, the system is initially in the state $\ket{+}=\frac{1}{\sqrt{2}}(\ket{0}+\ket{1})$ , while the bath is initially in the vacuum state $\ket{v}$ of the cavity. To determine the amplitude for the exact solution, from which we obtain the population $\rho_{11}(t)=|c_{11}(t)|^{2}$ of the excited state $\ket{1}$ , we need to apply a numerical inverse Laplace transform to Eq. 24, where the Laplace transform of the memory kernel $\hat{f}(s)$ is given by Section II.4.2. Irrespective of the specific positive values of $\eta$ , $\omega_{c}$ , and $\Omega_{0}$ , the roots of Eq. 24 are found to be complex. This means that the amplitude $c_{1}(t)$ is oscillatory. This trend persists even in cases of weak coupling, resulting in a damped oscillation of the amplitude.

IV.2.1 General observations

Figs. 2, 3, 4 and 5 offer a comparison of the temporal dynamics against the exact solution across four distinct approximations: TCL2, TCL4, CG-LE ( $=$ C-LE), and RWA-LE. The comparison metrics comprise the trace-norm distance $\|\rho_{\text{exact}}-\rho_{\text{approx}}\|_{1}$ , decay rate $\gamma$ , Lamb shift $S$ , population $\rho_{11}$ , and coherence $\rho_{01}$ , all with the Ohmic spectral density $J_{2}$ , but with different values of the coupling $\eta$ and qubit frequency ratio $\Omega_{0}/\omega_{c}$ . The equations plotted for each curve are listed in Table 1.

During the initial time intervals ( $\omega_{c}t\approx 1$ ), the TCL2 and TCL4 solutions generally exhibit much closer agreement with the exact solution than the Markovian CG-LE and RWA-LE, as evidenced especially by the trace-norm distance curves. The plots representing population and the real part of the coherence reveal that the TCL approximations aptly capture the Zeno effect observed in the exact solution – characterized by a gradual concave decay at short times. This stands in contrast to the monotonic exponential decay exhibited by the Markovian approaches.

Note that, as mentioned above in the discussion of Eq. 121, the TCL2 decay rate $\gamma$ gradually approaches the asymptotic behavior of the RWA-LE decay rate. Concerning the CG-LE, we note that it exhibits a smaller trace-norm distance from the exact solution compared to the RWA-LE, as anticipated in Ref. [18]. This is due to the aforementioned ability to optimize the coarse-graining time. Both Markovian approximations exhibit a constant decay rate and Lamb shift, but in the case of the CG-LE, the Lamb shift $S$ is zero, resulting in the coherence being purely real.

IV.2.2 Weak and strong coupling

Recall that the RWA-LE needs to satisfy Eq. 110, in particular $\eta\ll 1$ . In practice, we will use both weak ( $\eta=0.5$ ) and strong ( $\eta=1$ ) coupling. We thus expect the RWA-LE to be a relatively poor approximation in the later case, at least for relatively short evolutions.

Figs. 2 and 3 compare the results for weak and strong coupling. Using a qubit frequency equal to the cutoff frequency $\Omega_{0}/\omega_{c}=1$ , Fig. 2 shows the strong coupling results, while Fig. 3 shows the weak coupling results.

After rising first, the exact decay rate [Eq. 33b] for strong coupling ( $\eta=1$ ) in Fig. 2(b) exhibits a negative trend and continues to exhibit non-monotonic behavior for even longer evolution times. Intervals during which $\gamma$ turns negative correspond to non-Markovian dynamics, leading to an increase in population, as evidenced in Fig. 2(d). Consequently, during intermediate times, the approximation methods struggle to accurately match the exact solution, which displays a remarkably slow population decay.

In contrast, for weak coupling ( $\eta=1/2$ ), the decay rate in Fig. 3(b) remains positive, causing the population and coherence to approach zero more rapidly. This aligns closely with the behavior exhibited by the approximation methods, as demonstrated in Figs. 3(d), 3(e) and 3(f). This trend is also affirmed by the trace-norm distance plots in Figs. 2(a) and 3(a), where weaker coupling results in a smaller trace-norm distance difference between the approximation methods and the exact solution.

IV.2.3 Weak and strong qubit frequency

Fixing the coupling at $\eta=1$ [technically at the upper limit of Eq. 110], we compare the cases where the qubit frequency $\Omega_{0}$ is either less than or greater than the cutoff frequency $\omega_{c}$ , as shown in Figs. 4 and 5, respectively. Recall that Eq. 110 also imposes the validity condition $\Omega_{0}/\omega_{c}>1$ on the RWA-LE, so we expect better agreement in Fig. 5, as in indeed the case.

In more detail, from the trace-norm distance plots in Figs. 4(a) and 5(a), we observe that all the approximations exhibit closer agreement with the exact result when the qubit frequency is relatively high, especially at longer evolution times. This trend is also noticeable in the population (Figs. 4(d) and 5(d)) and coherence decay (Figs. 4(e), 5(e), 4(f) and 5(f)). For low qubit frequencies, as in Fig. 4(b), the exact decay rate $\gamma$ exhibits non-Markovian negative phases. Correspondingly, the population in Fig. 4(d) increases initially and then stabilizes at a non-zero value. The approximations fail to match the exact solution’s behavior for extended periods, with the TCL approximation being effective primarily during short times ( $\omega_{c}t<1.5$ ). This behavior is characteristic of strong non-Markovian behavior, which is expected in models with non-flat spectral densities such as the ones considered here.

Conversely, for higher qubit frequencies in Fig. 5(b), the behavior of the TCL decay rate resembles that of the exact solution. Even the exact Lamb shift in Fig. 5(c) closely aligns with the approximation methods, with the notable exception of CG-LE (where the Lamb shift is zero). Consequently, in this scenario, the approximations reasonably reproduce the exact behavior, as demonstrated in Figs. 5(d), 5(e) and 5(f). The TCL4 approximation is particularly good according to all six of our metrics.

IV.3 Exact solution vs TCL and Markov approximations for the impulse spectral density

For the impulse spectral density $J_{1}(\omega)=|g|^{2}\delta(\omega-\omega_{c})$ , the exact solution for the population $|c_{1}(t)|^{2}$ of the excited population [see Eq. 38] is purely periodic. Consequently, no Markovian approximation with exponential decay can accurately capture the exact solution. This renders the RWA-LE inadequate for the description, as remarked in the discussion of Eq. 124b. However, since the CG-LE retains the coarse-graining time-scale as an optimization parameter, we use it to understand how close a Markovian approximation can still be to the exact solution in this rather extreme non-Markovian case, even though we have already concluded that the validity condition of the CG-LE cannot be satisfied [recall the discussion in Section IV.1].

Figs. 6 and 7 depict the time evolution of the same metrics as in the previous section but for the impulse spectral density $J_{1}$ with two different qubit frequencies: $\Omega_{0}/\omega_{c}=1/2$ and $\Omega_{0}/\omega_{c}=2$ . For this spectral density, we have closed analytical expressions for the TCL Lamb shift [see Eqs. 156a and 157a] and decay rate [see Eqs. 156b and 157b]. In general, we observe that as expected, the CG-LE is a poor approximation to both the short-time and longer-time oscillatory behavior of the decay rate, population, and coherence.

In contrast, the TCL2 and TCL4 approximations in Eqs. 156b and 157b are fairly accurate for short-time dynamics. However, the populations and coherences under TCL4 deviate rapidly from the exact solution for longer times, as demonstrated in Figs. 6(d), 7(d), 6(e), 6(f), 7(e) and 7(f). The deviation is accentuated by the recurrent intervals of negative decay rates in the exact solution and TCL approximations. This is especially noticeable in TCL4, where the oscillations of the decay rates increase with time, taking on higher negative values, as seen in Figs. 6(b) and 7(b). Consequently, TCL4 develops an instability reflected in the diverging oscillation frequency of the coherence seen in Figs. 6(e) and 6(f). An explanation of the breakdown of the TCL approximation is given in Appendix I.

When comparing the low qubit frequency case ( $\Omega_{0}/\omega_{c}=1/2$ ) in Figs. 6(d), 6(e) and 6(f) with the high qubit frequency case ( $\Omega_{0}/\omega_{c}=2$ ) in Figs. 7(d), 7(e) and 7(f), a distinction emerges. In the former, TCL2 exhibits a decay to zero similar to CG-LE, while in the latter, TCL2 mirrors the oscillatory behavior of the exact solution. This observation underscores the higher accuracy of the TCL approximation as the qubit frequency increases. However, while TCL4 is a better approximation at short times $\omega_{c}t<1.5$ , it is worse than TCL2 at long times for the impulse spectral density.

IV.4 Exact solution vs TCL and Markov approximations for the triangular spectral density

The triangular spectral density $J_{3}(\omega)=\eta\omega\Theta(\omega_{c}-\omega)$ is intermediate between the impulse density $J_{1}(\omega)$ and the Ohmic density $J_{2}(\omega)$ because it matches $J_{2}(\omega)$ for low frequencies and drops to zero for frequencies exceeding the cutoff. In Figs. 8 and 9, we perform a comparison for $\Omega_{0}/\omega_{c}=1/2$ and $\Omega_{0}/\omega_{c}=1.8$ . Recall that here too, the validity condition of the CG-LE cannot be satisfied, as discussed in Section IV.1.

Similar to the impulse case $J_{1}$ , for $J_{3}$ , we observe a strong non-Markovian oscillatory behavior in the exact solution’s population and coherence, in Figs. 8(d), 8(e), 8(f), 9(d), 9(e) and 9(f), respectively. The decay rate (as seen in Figs. 8(b) and 9(b)) oscillates between positive and negative values.

As previously discussed, the TCL approximation is capable of capturing the quantum Zeno effect, which leads to a better agreement with the exact solution than the Markovian approximation does for short evolution times. However, a significant difference emerges in the exact decay rate between low and high qubit frequencies, as shown in Figs. 8(b) and 9(b). In the latter case, the decay rate oscillates and converges to zero (equivalent to the rate of the Markov approximations). Conversely, in the former case, the decay rate exhibits discontinuous behavior, rendering the approximation methods unsuitable for fitting the exact solution. In general, when $\Omega_{0}<\omega_{c}$ , the discontinuous behavior of the decay rate leads to oscillatory population dynamics, where the population first decays to zero, then revives, and subsequently decays again. However, for $\Omega_{0}>\omega_{c}$ , the decay rate converges to zero, resulting in non-decaying population dynamics, similar to what was observed for $J_{1}$ in Fig. 7(d).

Furthermore, as depicted in Fig. 8(c), the exact Lamb shift $S$ diverges when $\Omega_{0}<\omega_{c}$ in contrast to the zero Lamb shift in the CG-LE, the constant non-zero $S$ in the RWA-LE, and TCL2 (which converges to RWA). However, TCL4 diverges in an attempt to match the exact $S$ . In contrast, for the Lamb shift $S$ when $\Omega_{0}>\omega_{c}$ (as shown in Fig. 9(c)), the TCL approximations align with the oscillatory behavior of the exact solution. Even the constant RWA-LE approximation closely resembles the asymptotic behavior. However, the CG solution maintains a zero Lamb shift.

Moving on to the population and coherence in $\Omega_{0}<\omega_{c}$ , as illustrated in Figs. 8(d), 8(e) and 8(f), we observe that the approximations rapidly decay to zero, matching the exact solution at short times. However, the oscillatory behavior at later times remains beyond the reach of the approximations, similar to what was observed for $J_{1}$ in Fig. 6(d).

In contrast, for the population and coherence at $\Omega_{0}>\omega_{c}$ , as depicted in Figs. 9(d), 9(e) and 9(f), the TCL approximation exhibits an oscillatory and non-decaying behavior, akin to the exact solution, matching the Zeno effect at short times. TCL4, in particular, matches the exact solution well up to $\omega_{c}t\approx 2$ . The CG-LE approximation performs much better than RWA-LE for short times in describing the excited state population, which remains constant in the latter due to its zero decay rate [as can be seen from Eq. 130b]. The two are comparable in describing the coherence. For long times, the Markov approximations, which display monotonic decay, fail to capture the exact behavior. CG-LE, which aims to minimize the trace-norm distance with the exact solution, decays slowly to achieve its objective.

In summary, for the triangular spectral density, when $\Omega_{0}>\omega_{c}$ , the non-Markovian oscillatory behavior is captured qualitatively by the TCL approximation, as opposed to CG or RWA. In contrast, when $\Omega_{0}<\omega_{c}$ , all approximations exhibit rapid decay and fail to accurately capture the exact solution.

V Summary and conclusions

In this work, we studied the dynamics of a qubit in a cavity interacting with bosonic baths described by three different spectral densities: impulse, Ohmic, and triangular. The model is exactly solvable within the single-excitation subspace, and this allowed us to perform a comprehensive comparison to a number of different master equations, both Markovian (CG-LE, C-LE, and RWA-LE) and non-Markovian (TCL2 and TCL4). Of the three spectral densities, the Ohmic model is the most physically relevant, and the other two were introduced primarily to allow us to reach closed-form analytical results, as well as to study extreme non-Markovian dynamics.

For weak coupling and a large qubit frequency, the Ohmic case leads to a quasi-exponential decay, characteristic of the Markovian limit. In this regime, therefore, we found that the purely Markovian master equations are able to approximate the exact solution rather well. Outside of this regime, in particular also for the impulse and triangular spectral densities, the Markovian master equations perform poorly, but we found that TCL is still relatively accurate, in particular for short evolution times. The TCL is also able to capture distinctly non-Markovian features such as bath-induced population and coherence oscillations.

Within the regime of validity of the Markovian master equations, we found the CG-LE and C-LE to be better approximations than the standard RWA-LE. This was achieved by optimizing the coarse-graining time to minimize the difference between the CG-LE and the exact solution.

Overall, this work shows that low-order quantum master equations can be accurate when operated in their guaranteed regime of validity (short evolution times, in particular), but significant caution must be exercised in trusting their predictions outside of these regimes, as they can dramatically deviate from the exact dynamics. The TCL approach stands out as significantly more accurate than the Markovian master equations, even when the latter are fine-tuned via the optimization of a free parameter such as the coarse-graining time. The standard Lindblad equation based on the rotating wave approximation is particularly suspect.

Future research may wish to address—within the context of the same analytically solvable model as studied here, or similar analytically solvable models—the accuracy of TCL at higher orders, as well as a variety of other master equations, such as the phenomenological post-Markovian master equation [31, 32], Floquet-based master equations for periodic driving [33, 34], time-dependent Markovian quantum master equations [35, 36, 37, 6, 38], the universal Lindblad equation [39], the geometric-arithmetic master equation [40, 41], regularized cumulant-based master equations [42], as well as various adiabatic master equations [43, 44, 45]. Another interesting direction worth considering is to incorporate more general initial states, such as a hermal Gibbs state, a coherent state [20] or a squeezed initial bath state.

Acknowledgements.

This research was supported by the ARO MURI grant W911NF-22-S-0007, by the National Science Foundation the Quantum Leap Big Idea under Grant No. OMA-1936388, and by the Defense Advanced Research Projects Agency (DARPA) under Contract No. HR001122C0063. We thank Prof. Todd Brun for useful discussions and Dawei Zhong for providing the analysis given in Appendix E.

Appendix A Detailed exact solution

Considering the Hamiltonian in Eq. 1, the general solution for the $1$ -excitation subspace is given by Eq. 14. This solution can be expressed as a linear combination of the basis eigenvectors $\{\ket{\psi}_{0},\ket{\psi}_{1},\ket{\varphi_{k}}\}$ , which are tensor products of the qubit system basis $\{\ket{0},\ket{1}\}$ and the $1$ -excitation bath basis $\{\ket{v},\ket{k}\}_{k\in 1,2,..}$ , as shown in Eq. 13. Here, $\ket{v}$ denotes the vacuum state with no photons, and $\ket{k}=b_{k}^{\dagger}\ket{v}$ represents the state with one photon in mode $k$ . The coefficients of the linear combination satisfy the normalization condition given in Eq. 15. For the evolution of the closed system within the $1$ -excitation subspace, we utilize the Schrödinger equation with the interaction Hamiltonian in Eq. 11. Using


$\displaystyle\sigma_{-}(t)\ket{0}$	$\displaystyle=\sigma_{+}(t)\ket{1}=B(t)\ket{v}=0$	(167a)
$\displaystyle\sigma_{+}(t)\ket{0}$	$\displaystyle=e^{i\Omega_{0}t}\ket{1}\ ,\quad\sigma_{-}(t)\ket{1}=e^{-i\Omega_{0}t}\ket{0}$	(167b)
$\displaystyle B(t)\ket{k}$	$\displaystyle=g_{k}e^{-i\omega_{k}t}\ket{v}\ ,\quad B^{\dagger}(t)\ket{v}=\sum_{k}g_{k}^{*}e^{i\omega_{k}t}\ket{k}\ ,$	(167c)

this leads to Section II.3:


$\displaystyle i\dot{\ket{\phi}}$	$\displaystyle=\dot{c}_{0}(t)\ket{\psi_{0}}+\dot{c}_{1}(t)\ket{\psi_{1}}+\sum_{k}\dot{\mathfrak{c}}_{k}(t)\ket{\varphi_{k}}$	(168a)
	$\displaystyle=\lambda\tilde{H}_{SB}\ket{\phi(t)}$	(168b)
	$\displaystyle=\lambda[\sigma_{+}(t)\otimes B(t)+\sigma_{-}(t)\otimes B^{\dagger}(t)]\times$
	$\displaystyle\Big{(}c_{0}(t)\ket{\psi_{0}}+c_{1}(t)\ket{\psi_{1}}+\sum_{k}\dot{\mathfrak{c}}_{k}(t)\ket{\varphi_{k}}\Big{)}$	(168c)
	$\displaystyle=\lambda\Big{(}\sigma_{+}(t)\ket{0}\otimes B(t)\sum_{k}\mathfrak{c}_{k}(t)\ket{k}+$
	$\displaystyle\qquad\qquad c_{1}(t)\sigma_{-}(t)\ket{1}\otimes B^{\dagger}(t)\ket{v}\Big{)}$	(168d)
	$\displaystyle=\lambda\Big{(}\sum_{k}g_{k}\mathfrak{c}_{k}(t)e^{i(\Omega_{0}-\omega_{k})t}\ket{\psi_{1}}+$
	$\displaystyle\qquad\qquad\sum_{k}g_{k}^{*}c_{1}(t)e^{-i(\Omega_{0}-\omega_{k})t}\ket{\varphi_{k}}\Big{)}$	(168e)

We can obtain the joint system-bath density matrix using Eq. 14:

		$\displaystyle\rho_{SB}(t)=\left(\|c_{0}\|^{2}\|{0}\rangle\!\langle 0\|+c_{0}c_{1}^{}\|{0}\rangle\!\langle 1\|+c_{1}c_{0}^{}\|{1}\rangle\!\langle 0\|\right.$
		$\displaystyle\left.+\|c_{1}\|^{2}\|{1}\rangle\!\langle 1\|\right)\otimes\|{v}\rangle\!\langle v\|+\|{0}\rangle\!\langle 0\|\otimes\left(\sum_{j,k}\mathfrak{c}_{j}\mathfrak{c}_{k}^{*}\|{j}\rangle\!\langle k\|\right.$
		$\displaystyle\left.+\sum_{k}(c_{0}\mathfrak{c}_{k}^{}\|{v}\rangle\!\langle k\|+c_{0}^{}\mathfrak{c}_{k}\|{k}\rangle\!\langle v\|)\right)+\sum_{k}c_{1}\mathfrak{c}_{k}^{*}\|{1}\rangle\!\langle 0\|$
		$\displaystyle\otimes\|{v}\rangle\!\langle k\|+\sum_{k}\mathfrak{c}_{k}c_{1}^{*}\|{0}\rangle\!\langle 1\|\otimes\|{k}\rangle\!\langle v\|\ .$		(169)

Taking the partial trace over the bath, we obtain the system state in the matrix form shown in Eq. 27. Its time derivative is given in Eq. 30. If we write the system density matrix in terms of the populations $\rho_{00},\rho_{11}$ and its coherences $\rho_{01},\rho_{10}$ , we have:


$\displaystyle\rho_{11}(t)$	$\displaystyle=\|c_{1}(t)\|^{2}=1-\rho_{00}(t)$	(170a)
$\displaystyle\rho_{01}(t)$	$\displaystyle=c_{0}\dot{c}^{}_{1}(t)=\rho_{10}^{}(t)\ .$	(170b)

Explicitly substituting the system density matrix into the ansatz Section II.3 along with Eq. 31, we obtain:

\displaystyle\mathcal{K}_{S}(t)\rho=\left(\begin{array}[]{cc}\gamma(t)\rho_{11}&\frac{1}{2}(i\mathcal{S}(t)-\gamma(t))\rho_{01}\\ \frac{1}{2}(-i\mathcal{S}(t)-\gamma(t))\rho_{10}&-\gamma(t)\rho_{11}\end{array}\right)

(173)

Comparing this with Eq. 30 gives us the following differential equations:


$\displaystyle\partial_{t}\|c_{1}(t)\|^{2}$	$\displaystyle=-\gamma(t)\|c_{1}(t)\|^{2}$	(174a)
$\displaystyle c_{0}^{*}\dot{c}_{1}(t)$	$\displaystyle=-\frac{1}{2}\left(i\mathcal{S}(t)+\gamma(t)\right)c_{0}^{*}c_{1}(t)\ .$	(174b)

From the first equation, we obtain the population evolution Eq. 34a, and from the second one, we can take its real and imaginary parts to obtain Eq. 33.

Appendix B Simplification of the decay rate expressions for the CG-LE

B.1 Ohmic spectral density

Starting from Eq. 85 we have:


$\displaystyle\gamma(\tau)$	$\displaystyle=\int_{0}^{\infty}\eta\omega e^{-\omega/\omega_{c}}\tau\operatorname{sinc}^{2}\left(\frac{(\Omega_{0}-\omega)\tau}{2}\right)d\omega$	(175a)
	$\displaystyle=\frac{2}{\tau}\eta\int_{-\Omega_{0}}^{\infty}d\nu\frac{(\nu+\Omega_{0})e^{-(\nu+\Omega_{0})/\omega_{c}}}{\nu^{2}}(1-\cos\nu\tau)\ ,$	(175b)

where we used the change of variables $\nu=\omega-\Omega_{0}$ . This integral can be split into two parts. By using the exponential integral function Eq. 43c we have:

	$\displaystyle\int_{-\Omega_{0}}^{\infty}d\nu\frac{e^{-\nu/\omega_{c}}}{\nu}(1-\cos\nu\tau)=-\operatorname{Ei}\left(\frac{\Omega_{0}}{\omega_{c}}\right)$
	$\displaystyle\quad+\frac{1}{2}\left[\operatorname{Ei}\left(\frac{\Omega_{0}}{\omega_{c}}+i\Omega_{0}\tau\right)+\operatorname{Ei}\left(\frac{\Omega_{0}}{\omega_{c}}-i\Omega_{0}\tau\right)\right]\ .$		(176)

For the second integral, we can utilize integration by parts. Setting $u=e^{-\nu/\omega_{c}}(1-\cos\nu\tau)$ and $dv=\nu^{-2}d\nu$ , we have:

	$\displaystyle\int_{-\Omega_{0}}^{\infty}udv=-\frac{e^{\Omega_{0}/\omega_{c}}}{\Omega_{0}}(1-\cos\Omega_{0}\tau)$
	$\displaystyle+\int_{-\Omega_{0}}^{\infty}d\nu\frac{e^{-\nu/\omega_{c}}}{\nu}\left(\tau\sin\nu\tau-\frac{1}{\omega_{c}}(1-\cos\nu\tau)\right)\ .$		(177)

Since

	$\displaystyle\int_{-\Omega_{0}}^{\infty}d\nu\frac{e^{-\nu/\omega_{c}}}{\nu}\tau\sin\nu\tau=$		(178)
	$\displaystyle-\frac{i}{2}\tau\left[\operatorname{Ei}\left(\frac{\Omega_{0}}{\omega_{c}}+i\Omega_{0}\tau\right)-\operatorname{Ei}\left(\frac{\Omega_{0}}{\omega_{c}}-i\Omega_{0}\tau\right)\right]\ ,$

we can combine all the terms, and arrive at the final expression for $\gamma(\tau)$ given in Eq. 85.

B.2 Triangular spectral density

Starting from Eq. 86 we have:


$\displaystyle\gamma(\tau)$	$\displaystyle=\int_{0}^{\omega_{c}}\eta\omega\tau\operatorname{sinc}^{2}\left(\frac{(\Omega_{0}-\omega)\tau}{2}\right)d\omega$	(179a)
	$\displaystyle=\frac{2}{\tau}\int_{-\Omega_{0}}^{\omega_{c}-\Omega_{0}}d\nu\frac{\eta(\nu+\Omega_{0})}{\nu^{2}}(1-\cos\nu\tau)\ .$	(179b)

where we again used the change of variables $\nu=\omega-\Omega_{0}$ . Now for $\Omega_{0}<\omega_{c}$ , we have:

	$\displaystyle\int_{-\Omega_{0}}^{\omega_{c}-\Omega_{0}}d\nu\frac{1-\cos\nu\tau}{\nu}=$
	$\displaystyle\ \ \ \ln\left(\frac{\omega_{c}-\Omega_{0}}{\Omega_{0}}\right)-\operatorname{Ci}((\omega_{c}-\Omega_{0})\tau)+\operatorname{Ci}(\Omega_{0}\tau)\ ,$		(180)

where the sine and cosine integral functions are given in Eqs. 43a and 43b, respectively. The second term can be obtained by using integration by parts:

	$\displaystyle\int_{-\Omega_{0}}^{\omega_{c}-\Omega_{0}}d\nu\frac{1-\cos\nu\tau}{\nu^{2}}$	$\displaystyle=-\frac{1-\cos\nu\tau}{\nu}\bigg{\|}_{-\Omega_{0}}^{\omega_{c}-\Omega_{0}}$
		$\displaystyle+\tau\int_{-\Omega_{0}}^{\omega_{c}-\Omega_{0}}\frac{\sin(\nu\tau)}{\nu}\ .$		(181)

The last term is the sine integral function. Combining, we arrive at the total rate as given by Eq. 86.

Appendix C Why the first order cumulant $K^{(1)}$ in the C-LE can be made to vanish

Let $\langle{B}\rangle=\mathrm{Tr}(\rho_{B}B)$ and define a new, shifted bath operator:

B^{\prime}\equiv B-\langle{B}\rangle I_{B}\ .

(182)

Its expectation value vanishes:

\langle{B^{\prime}}\rangle=\langle{B}\rangle-\langle{B}\rangle\langle{I_{B}}\rangle=0\ .

(183)

The corresponding bath interaction picture operator is:

	$\displaystyle B^{\prime}(t)$	$\displaystyle=U_{B}^{\dagger}(t)B^{\prime}U_{B}(t)=U_{B}^{\dagger}(t)(B-\langle{B}\rangle I_{B})U_{B}(t)$
		$\displaystyle=B(t)-\langle{B}\rangle I_{B}\ ,$		(184)

where, as before, $U_{B}(t)=e^{-iH_{B}t}$ . Assuming stationarity $[H_{B},\rho_{B}(0)]=0$ immediately implies $[U_{B}(t),\rho_{B}]=0$ . In this case also $\langle{B^{\prime}(t)}\rangle=0$ , since then:


$\displaystyle\langle{B^{\prime}(t)}\rangle$	$\displaystyle=\mathrm{Tr}[U_{B}^{\dagger}(t)B^{\prime}U_{B}(t)\rho_{B}]$	(185a)
	$\displaystyle=\mathrm{Tr}[B^{\prime}U_{B}(t)\rho_{B}U_{B}^{\dagger}(t)]=\langle{B^{\prime}}\rangle=0\ .$	(185b)

Let

H^{\prime}_{SB}\equiv A\otimes B^{\prime}\ ,\quad\Delta H_{S}^{\prime}\equiv\langle{B}\rangle A\ ,\quad H^{\prime}_{S}\equiv H_{S}+\Delta H^{\prime}_{S}\ .

(186)

Correspondingly, the system-bath interaction can be written as

H_{SB}=H_{SB}^{\prime}+A\otimes(B-B^{\prime})=H_{SB}^{\prime}+\Delta H_{S}^{\prime}\otimes I_{B}\ .

(187)

Thus, we can write the full Hamiltonian as:


$\displaystyle H$	$\displaystyle=H_{S}\otimes I_{B}+H_{SB}+I_{S}\otimes H_{B}$	(188a)
	$\displaystyle=H^{\prime}_{0}+H^{\prime}_{SB}\ ,\quad H^{\prime}_{0}\equiv H^{\prime}_{S}\otimes I_{B}+I_{S}\otimes H_{B}\ ,$	(188b)

where $H^{\prime}_{0}$ defines a new, shifted interaction picture Hamiltonian.

Correspondingly, in this new interaction picture $\tilde{H}^{\prime}(t)=U_{0}^{\prime{\dagger}}(t)H^{\prime}_{SB}U^{\prime}_{0}(t)=A(t)\otimes B^{\prime}(t)$ , and we find, using Eq. 185:

	$\displaystyle\mathrm{Tr}_{B}[\tilde{H}^{\prime}(t),\rho_{SB}(0)]$	$\displaystyle=\mathrm{Tr}_{B}[A(t)\otimes B^{\prime}(t),\rho_{S}(0)\otimes\rho_{B}(0)]$
		$\displaystyle=\langle{B^{\prime}(t)}\rangle[A(t),\rho_{S}(0)]=0\ .$		(189)

Therefore, $K^{\prime(1)}(t)\rho(0)=0$ , with $K^{\prime(1)}$ defined within the shifted interaction picture and with the modified system-bath interaction $H^{\prime}_{SB}$ .

The extension to the case when $H_{SB}$ has the general form $H_{SB}=\sum_{\alpha}A_{\alpha}\otimes B_{\alpha}$ is immediate; in this case $B^{\prime}_{\alpha}=B_{\alpha}-\langle{B_{\alpha}}\rangle I_{B}$ and:

H_{S}^{\prime}=\sum_{\alpha}\langle{B_{\alpha}}\rangle A_{\alpha}\ ,\quad H^{\prime}_{SB}=\sum_{\alpha}A_{\alpha}\otimes B^{\prime}_{\alpha}\ .

(190)

Now, for our bath operators $B_{+}(t)=\sum_{k}g_{k}e^{i\omega_{k}t}b_{k}$ and $B_{-}=B_{+}^{*}$ and the bath state $\rho_{B}=|{v}\rangle\!\langle v|$ we have that

\langle{B_{+}(t)}\rangle=\mathrm{Tr}[\rho_{B}B_{+}(t)]=\sum_{k}g_{k}e^{i\omega_{k}t}\langle{b_{k}}\rangle=0\ ,

(191)

and analogously $\langle{B_{-}(t)}\rangle=0$ , both arising from the fact that the annihilation and creation operators average to zero [Section III.1]. As a result, in our case, in fact $B^{\prime}(t)=B(t)$ .

Appendix D Proof of Eq. 97b

It is useful to relate $\mathcal{B}_{\alpha\beta\omega}(t)$ with $b_{\alpha\beta\omega}(t)$ :

	$\displaystyle\mathcal{B}_{\alpha\beta\omega}(t)=\int_{0}^{t}ds\int_{0}^{s}ds^{\prime}e^{i(\omega(s-s^{\prime})}\mathcal{B}_{\alpha\beta}(s,s^{\prime})$		(192)
	$\displaystyle=\left[\int_{0}^{t}ds\int_{0}^{t}ds^{\prime}-\int_{0}^{t}ds\int_{s}^{t}ds^{\prime}\right]e^{i(\omega^{\prime}s-\omega s^{\prime})}\mathcal{B}_{\alpha\beta}(s,s^{\prime})$
	$\displaystyle=\left[\int_{0}^{t}ds\int_{0}^{t}ds^{\prime}-\int_{0}^{t}ds^{\prime}\int_{0}^{s^{\prime}}ds\right]e^{i\omega(s-s^{\prime})}\mathcal{B}_{\alpha\beta}(s,s^{\prime})$
	$\displaystyle=b_{\alpha\beta\omega}(t)-\int_{0}^{t}ds\int_{0}^{s}ds^{\prime}e^{-i\omega(s^{\prime}-s)}\mathcal{B}_{\alpha\beta}(s^{\prime},s)$
	$\displaystyle=b_{\alpha\beta\omega}(t)-\mathcal{B}^{*}_{\alpha\beta\omega}(t)\ ,$		(193)

where in the last line we used Eq. 94.

Appendix E $[H_{B},\rho_{B}]\neq 0$

The bath state can be obtained via a partial trace of the system from the state Eq. 14, whose explicit pure density matrix is given in Appendix A:

	$\displaystyle\rho_{B}(t)=\mathrm{Tr}_{S}\|{\phi(t)}\rangle\!\langle\phi(t)\|=(\|c_{0}\|^{2}+\|c_{1}\|^{2})\|{v}\rangle\!\langle v\|$		(194)
	$\displaystyle+\sum_{k}\mathfrak{c}_{k}(t)c_{0}^{}\|{k}\rangle\!\langle v\|+\sum_{k}\mathfrak{c}_{k}^{}(t)c_{0}\|{v}\rangle\!\langle k\|+\sum_{k,k^{\prime}}\mathfrak{c}_{k}\mathfrak{c}^{*}_{k^{\prime}}\|{k}\rangle\!\langle k^{\prime}\|\ .$

Now using the bath Hamiltonian in Eq. 9 and using the identities $n_{j}\ket{v}=0$ and $n_{j}\ket{k}=\delta_{jk}\ket{j}$ , we have:


$\displaystyle H_{B}\rho_{B}$	$\displaystyle=\sum_{j,k}\omega_{j}\mathfrak{c}_{j}(t)\mathfrak{c}_{k}^{}(t)\|{j}\rangle\!\langle k\|+\sum_{j}\omega_{j}c_{0}^{}\mathfrak{c}_{j}(t)\|{j}\rangle\!\langle v\|$	(195a)
	$\displaystyle=(\rho_{B}H_{B})^{\dagger}\neq\rho_{B}H_{B}$	(195b)
	$\displaystyle=\sum_{j,k}\omega_{j}\mathfrak{c}_{k}(t)\mathfrak{c}_{j}^{}(t)\|{k}\rangle\!\langle j\|+\sum_{j}\omega_{j}c_{0}\mathfrak{c}_{j}^{}(t)\|{v}\rangle\!\langle j\|.$	(195c)

In the particular case where $\rho_{B}(0)=|{v}\rangle\!\langle v|$ , or $c_{0}=0$ and $\mathfrak{c}_{k}(0)=0$ , it trivially follows that $[H_{B},\rho_{B}(0)]=0$ .

Appendix F Inadequacy of the Markov approximation

In the main text, we showed that we can reduce Eq. 112 to Eq. 120. However, this may lead to an unbounded approximation error, as we now show in detail.

Before the Markov approximation, Eq. 112 contains terms of the following form:

\displaystyle\int_{0}^{t}d\tau\mathcal{B}_{\alpha\beta}(\pm\tau)e^{\pm i\Omega_{0}\tau}\sigma_{\alpha}\sigma_{\beta}\tilde{\rho}(t-\tau)\ ,

(196)

where $\alpha,\beta\in\{+,-\}$ . The Markov approximation replaces the latter with

\displaystyle\int_{0}^{\infty}d\tau\mathcal{B}_{\alpha\beta}(\pm\tau)e^{\pm i\Omega_{0}\tau}\sigma_{\alpha}\sigma_{\beta}\tilde{\rho}(t)\ .

(197)

Therefore the approximation error is the difference between these two quantities, which we write as

\delta=\|\Delta_{1}+\Delta_{2}\|\leq\|\Delta_{1}\|+\|\Delta_{2}\|\ ,

(198)

where $\|\cdot\|$ represents the operator norm and


$\displaystyle\Delta_{1}$	$\displaystyle\equiv\int_{0}^{\infty}d\tau\mathcal{B}_{\alpha\beta}(\pm\tau)e^{\pm i\Omega_{0}\tau}\sigma_{\alpha}\sigma_{\beta}(\tilde{\rho}(t)-\tilde{\rho}(t-\tau))$	(199a)
$\displaystyle\Delta_{2}$	$\displaystyle\equiv\int_{t}^{\infty}d\tau\mathcal{B}_{\alpha\beta}(\pm\tau)e^{\pm i\Omega_{0}\tau}\sigma_{\alpha}\sigma_{\beta}\tilde{\rho}(t-\tau)\ .$	(199b)

Let $\|\cdot\|_{1}$ denote the trace norm and observe that $\|AB\|\leq\|A\|\|B\|_{1}$ for any pair of operators $A$ and $B$ .

For the Ohmic spectral density $J_{2}(\omega)$ we have, using Eq. 125:


$\displaystyle\\|\Delta_{2}\\|$	$\displaystyle\leq\int_{t}^{\infty}d\tau\|\mathcal{B}_{\alpha\beta}(\pm\tau)\|\,\\|\sigma_{\alpha}\\|\,\\|\sigma_{\beta}\\|\,\\|\tilde{\rho}(t-\tau)\\|_{1}$	(200a)
	$\displaystyle\leq\int_{t}^{\infty}d\tau\|\mathcal{B}_{+-}(\pm\tau)\|=\eta\omega_{c}\int_{t}^{\infty}\frac{d(\omega_{c}\tau)}{1+(\omega_{c}\tau)^{2}}$	(200b)
	$\displaystyle=\left(\frac{\pi}{2}-\arctan(\omega_{c}\tau)\right)\eta\omega_{c}\ .$	(200c)

This quantity goes to zero for $\omega_{c}\tau\gg 1$ as required. On the other hand, the error term $\Delta_{1}$ is unbounded. First, by the mean value theorem, there is a point $t^{\prime}\in[t-\tau,t]$ such that:

\displaystyle\|\tilde{\rho}(t)-\tilde{\rho}(t-\tau)\|\leq\tau\sup_{t^{\prime}\in[t-\tau,t]}\|\dot{\tilde{\rho}}(t^{\prime})\|\ ,

(201)

so that:

\displaystyle\|\Delta_{1}\|\leq\int_{0}^{\infty}d\tau\,\tau|B_{\alpha\beta}|\sup_{t^{\prime}\in[t-\tau,t]}\|\dot{\tilde{\rho}}(t^{\prime})\|\ .

(202)

We can bound $\|\dot{\tilde{\rho}}(t^{\prime})\|$ using the state evolution in Section III.3, where we undo the Markovian approximation by replacing $\tilde{\rho}(t)$ with $\tilde{\rho}(t-\tau)$ [i.e., returning to Eq. 112]:


$\displaystyle\\|\dot{\tilde{\rho}}(t)$	$\displaystyle\\|\leq\int_{0}^{t}d\tau\|B_{+-}(\tau)\|\,\\|[\sigma_{+},\sigma_{-}\tilde{\rho}(t-\tau)]\\|$
	$\displaystyle\quad+\int_{0}^{t}d\tau\|B_{-+}(-\tau)\|\,\\|[\sigma_{+},\tilde{\rho}(t-\tau)\sigma_{-}]\\|$	(203a)
	$\displaystyle\leq 4\int_{0}^{t}d\tau\|B_{+-}(\tau)\|\,\\|\sigma_{+}\\|\,\\|\sigma_{-}\\|\\|\tilde{\rho}(t-\tau)\\|_{1}$	(203b)
	$\displaystyle\leq 4\int_{0}^{t}d\tau\|B_{+-}(\tau)\|\leq 4\int_{0}^{\infty}d\tau\|B_{+-}(\tau)\|$	(203c)
	$\displaystyle=2\pi\eta\omega_{c}\ ,$	(203d)

where in the last line we used Eq. 200c evaluated at $t=0$ .

While the integral

\displaystyle\int_{0}^{\infty}\tau^{n}|\mathcal{B}_{+-}(\tau)|d\tau=\frac{\pi}{2}\eta\omega_{c}^{1-n}\sec(n\pi/2)

(204)

converges for $|n|<1$ , it does not for $|n|\geq 1$ . Indeed, for the Ohmic spectral density we have:

\displaystyle\int_{0}^{\infty}\tau|\mathcal{B}_{+-}(\tau)|d\tau=\lim_{\tau\rightarrow\infty}\frac{\eta}{2}\ln[1+(\omega_{c}\tau)^{2}]\ ,

(205)

which diverges. While this diverging upper bound does not prove that the error $\delta$ itself diverges, it does suggest that this is indeed the case and that hence the approximation of replacing $\rho(t-\tau)$ by $\rho(t)$ is inaccurate. Indeed, the exact solution exhibits excited state population oscillations instead of purely Markovian exponential decay for all values of the parameters of the Ohmic density.

A similar situation arises for the triangular spectral density $J_{3}(\omega)$ . We find, numerically, that the integral $\int_{0}^{\infty}\tau^{n}|\mathcal{B}_{+-}(\tau)|d\tau$ [recall Eq. 128] diverges for $n>0$ , and therefore the Markov approximation is inadequate.

For a rigorous error bound, see Ref. [6].

Appendix G Proof of Eq. 119

Recall that $\Gamma_{\alpha\beta}(\omega)\equiv\int_{0}^{\infty}d\tau\mathcal{B}_{\alpha\beta}(\tau)e^{i\omega\tau}$ and $H_{SB}=\sigma_{+}\otimes B_{+}+\sigma_{-}\otimes B_{-}$ and $B_{+}^{\dagger}=B_{-}$ . Thus, if $\alpha\neq\beta$ , $B_{\alpha}=B^{\dagger}_{\beta}$ . It follows that


$\displaystyle\mathcal{B}_{\beta\alpha}^{*}(\tau)$	$\displaystyle=\mathrm{Tr}\left[\left(\rho_{B}B_{\beta}(\tau)B_{\alpha}\right)^{\dagger}\right]$	(206a)
	$\displaystyle=\mathrm{Tr}\left[\rho_{B}U_{B}(\tau)B^{\dagger}_{\alpha}U_{B}^{\dagger}(\tau)B^{\dagger}_{\beta}\right]$	(206b)
	$\displaystyle=\mathrm{Tr}\left[\rho_{B}U_{B}(\tau)B_{\beta}U_{B}^{\dagger}(\tau)B_{\alpha}\right]=\mathcal{B}_{\beta\alpha}(-\tau)\ .$	(206c)

Hence, $\mathcal{B}_{\pm\mp}^{*}(\tau)=\mathcal{B}_{\pm\mp}(-\tau)$ [as well as $\mathcal{B}^{*}_{\pm\pm}(\tau)=\mathcal{B}_{\mp\mp}(-\tau)$ , though we don’t use this result], which yields $\Gamma^{*}_{{\pm\mp}}(\omega)=\int_{0}^{\infty}d\tau\mathcal{B}_{{\pm\mp}}(-\tau)e^{-i\omega\tau}$ .

Appendix H Simplification of the Lamb shift expressions for the RWA-LE

H.1 $J_{2}(\omega)=\eta\omega e^{-\omega/\omega_{c}}$

Here we derive Eq. 126b. Our starting point is Eq. 126a, which we write as


		$\displaystyle\Gamma_{+-}(\Omega_{0})/(\eta\omega_{c})=\int_{0}^{\infty}d(\omega_{c}\tau)\frac{e^{i\Omega_{0}\tau}}{(1+i\omega_{c}\tau)^{2}}$		(207a)
		$\displaystyle\quad=\int_{0}^{\infty}\frac{e^{i\alpha x}dx}{(1+ix)^{2}}$		(207b)
		$\displaystyle\quad=\int_{-\infty}^{\infty}\frac{e^{i\alpha x}dx}{(1+ix)^{2}}-\int_{-\infty}^{0}\frac{e^{i\alpha x}dx}{(1+ix)^{2}}\ ,$		(207c)

where $\alpha=\Omega_{0}/\omega_{c}$ and $x=\omega_{c}\tau$ .

The complex exponential integral is defined as follows [46]:

\displaystyle\text{E}_{1}(z)\equiv\int_{z}^{\infty}\frac{e^{-u}}{u}du=-\int_{-\infty}^{-z}\frac{e^{u}}{u}du\ ,

(208)

where $z=x+iy$ and $|\arg(z)|\leq\pi/2$ . The following property holds for $y>0$ :

\displaystyle-\text{E}_{1}(-y)=\operatorname{Ei}(y)+i\pi\ ,

(209)

where the real exponential integral $\operatorname{Ei}$ was defined in Eq. 43c.

The first integral on the right hand side of Eq. 207c can be computed via the residue theorem:

\displaystyle\int_{-\infty}^{\infty}\frac{e^{i\alpha x}dx}{(1+ix)^{2}}=-\int_{-\infty}^{\infty}\frac{e^{i\alpha x}dx}{(x-i)^{2}}=-2\pi i(i\alpha e^{-\alpha})\ ,

(210)

since the second order pole of the analytic function is at $x=i$ , so that the residue is:

\displaystyle\text{Res}=\frac{d}{dx}e^{i\alpha x}\bigg{|}_{x=i}=i\alpha e^{-\alpha}\ .

(211)

For the second integral on the right side of Eq. 207c we use a change of variable $z=x-i$ and integrate by parts:


	$\displaystyle\int_{-\infty}^{0}\frac{e^{i\alpha x}dx}{(1+ix)^{2}}=-\int_{-\infty}^{0}\frac{e^{i\alpha x}dx}{(x-i)^{2}}=-e^{-\alpha}\int_{-\infty}^{-i}\frac{e^{i\alpha z}dz}{z^{2}}$		(212a)
	$\displaystyle=-e^{-\alpha}\left[-\frac{e^{i\alpha z}}{z}\bigg{\|}_{-\infty}^{-i}+i\alpha\int_{-\infty}^{-i}\frac{e^{i\alpha z}}{z}dz\right]$		(212b)
	$\displaystyle=i-i\alpha e^{-\alpha}\int_{-\infty}^{-i}\frac{e^{i\alpha z}}{z}dz=i-i\alpha e^{-\alpha}\int_{-\infty}^{\alpha}\frac{e^{u}}{u}du\ .$		(212c)

Using the complex exponential integral Eq. 208 alongside Eq. 209, we have:


$\displaystyle\int_{-\infty}^{0}\frac{e^{i\alpha x}dx}{(1+ix)^{2}}$	$\displaystyle=i-i\alpha e^{-\alpha}(-\text{E}_{1}(-\alpha))$	(213a)
	$\displaystyle=\pi\alpha e^{-\alpha}+i-i\alpha e^{-\alpha}\operatorname{Ei}(\alpha)\ .$	(213b)

Hence, the integral of interest, Eq. 207b, becomes

	$\displaystyle\int_{0}^{\infty}\frac{e^{i\alpha x}dx}{(1+ix)^{2}}$	$\displaystyle=\pi\alpha e^{-\alpha}-i+i\alpha e^{-\alpha}\operatorname{Ei}(\alpha)$
		$\displaystyle=-i+\alpha e^{-\alpha}(\pi+i\operatorname{Ei}(\alpha))\ .$		(214)

Collecting these results we obtain Eq. 126b.

H.2 $J_{3}(\omega)=\eta\omega\Theta(\omega_{c}-\omega)$

Here we derive Eq. 130. Our starting point is Eq. 129, which we write as

\displaystyle\Gamma_{+-}(\Omega_{0})/(\eta\omega_{c})=\int_{0}^{\infty}dx\frac{e^{-ix}(1+ix)-1}{x^{2}}e^{i\alpha x}\ ,

(215)

where, again, $\alpha=\Omega_{0}/\omega_{c}$ and $x=\omega_{c}\tau$ .

We start from the following indefinite integral, solved via integration by parts:

\displaystyle\int\frac{e^{iux}}{x^{2}}dx=-\frac{e^{iux}}{x}+iu\int\frac{e^{iux}}{x}dx\ .

(216)

Hence the integral in Eq. 215 is:


	$\displaystyle=\int_{0}^{\infty}\frac{e^{i(\alpha-1)x}}{x^{2}}dx-\int_{0}^{\infty}\frac{e^{i\alpha x}}{x^{2}}dx+i\int_{0}^{\infty}\frac{e^{i(\alpha-1)x}}{x}dx$		(217a)
	$\displaystyle=\left.\frac{e^{i\alpha x}(1-e^{-ix})}{x}\right\|_{0}^{\infty}-i\alpha\int_{0}^{\infty}\frac{e^{i\alpha x}(1-e^{-ix})}{x}dx$		(217b)
	$\displaystyle=-i+2\alpha\int_{0}^{\infty}e^{i(2\alpha-1)u}\operatorname{sinc}(u)du$		(217c)

To compute the last integral, we consider the real and imaginary parts separately. The real part is:


	$\displaystyle\int_{0}^{\infty}\cos(2\alpha-1)u\frac{\sin u}{u}du$		(218a)
	$\displaystyle=\begin{cases}\int_{0}^{\infty}\frac{\sin 2\alpha u+\sin(2-2\alpha)u}{2u}du=\frac{\pi}{2}&0<\alpha<1\\ \int_{0}^{\infty}\frac{\sin 2\alpha u-\sin(2\alpha-2)u}{2u}du=0&\alpha>1\end{cases}$		(218b)

The imaginary part is:


	$\displaystyle\int_{0}^{\infty}\sin(2\alpha-1)u\frac{\sin u}{u}du$		(219a)
	$\displaystyle\quad=\int_{0}^{\infty}\frac{\cos(2\alpha-2)u-\cos 2\alpha u}{2u}du$		(219b)
	$\displaystyle\quad=\begin{cases}-\frac{1}{2}\ln\left(\frac{1}{\alpha}-1\right)&0<\alpha<1\\ -\frac{1}{2}\ln\left(1-\frac{1}{\alpha}\right)&\alpha>1\end{cases}\ .$		(219c)

Therefore the integral in Eq. 215 is

\displaystyle\frac{\Gamma_{+-}(\Omega_{0})}{\eta\omega_{c}}=\begin{cases}\pi\alpha-i\left(1+\alpha\ln\left(\frac{1}{\alpha}-1\right)\right)&0<\alpha<1\\ -i\left(1+\alpha\ln\left(1-\frac{1}{\alpha}\right)\right)&\alpha>1\end{cases}\ .

(220)

Collecting these results and using Eq. 121 we obtain Eq. 130.

Appendix I Breakdown of the TCL approximation

For strong coupling or a small gap $\Omega_{0}<\omega_{c}$ , as in Figs. 6 and 8, we observe that the TCL approximation is accurate for short times but fails to capture the subsequent oscillatory behavior. For these cases the operator $I-\Sigma$ in Eq. 141 is not invertible. For example, in the case of the impulse spectral density $J_{1}$ , there is a common time $t_{0}$ where the excited state population $\rho_{11}(t_{0})=|c_{1}(t_{0})|^{2}=0$ independently from the initial condition. This can be seen directly from Eq. 38 or from Eq. 40b when $\gamma$ diverges. The time $t_{0}$ is the minimum time $t_{n}$ such that

t_{n}=\frac{2}{|\delta|}\left(\arctan\frac{|\delta|}{|\Omega_{0}-\omega_{c}|}+n\pi\right),

(221)

where $n$ is and integer and $\delta$ is given in Eq. 39. Since TCL is a time-local master equation, it is impossible to invert the evolution for $t\geq t_{0}$ back to its initial condition. This implies that the TCL gives inconsistent result when the exact solution for the population vanishes.

A similar argument can be given for the triangular spectral density $J_{3}$ as shown in Fig. 8. In contrast, this argument is not valid for the infinite range Ohmic spectral density $J_{2}$ , where the population decreases but only vanishes as $t\to\infty$ .

References

Alicki and Lendi [2007] R. Alicki and K. Lendi, Quantum Dynamical Semigroups and Applications (Springer Science & Business Media, 2007).
Breuer and Petruccione [2002a] H.-P. Breuer and F. Petruccione, The Theory of Open Quantum Systems (Oxford University Press, Oxford, 2002).
Rivas and Huelga [2012] A. Rivas and S. F. Huelga, Open Quantum Systems: An Introduction, SpringerBriefs in Physics (Springer-Verlag, Berlin Heidelberg, 2012).
Lindblad [1976] G. Lindblad, On the generators of quantum dynamical semigroups, Comm. Math. Phys. 48, 119 (1976).
Gorini et al. [1976] V. Gorini, A. Kossakowski, and E. C. G. Sudarshan, Completely positive dynamical semigroups of ${N}$ -level systems, J. Math. Phys. 17, 821 (1976).
Mozgunov and Lidar [2020] E. Mozgunov and D. Lidar, Completely positive master equation for arbitrary driving and small level spacing, Quantum 4, 227 (2020).
Krovi et al. [2007] H. Krovi, O. Oreshkov, M. Ryazanov, and L. D., Phys. Rev. A. 76, 052117 (2007).
Jing and Wu [2018] J. Jing and L.-A. Wu, Decoherence and control of a qubit in spin baths: an exact master equation study, Scientific Reports 8, 1471 (2018).
Jaynes and Cummings [1963] E. T. Jaynes and F. W. Cummings, Comparison of quantum and semiclassical radiation theories with application to the beam maser, Proceedings of the IEEE 51, 89 (1963).
Garraway [1997] B. M. Garraway, Nonperturbative decay of an atomic system in a cavity, Physical Review A 55, 2290 (1997).
Carmichael [1993] H. Carmichael, An open systems approach to quantum optics, Lecture notes in physics No. m18 (Springer-Verlag, Berlin, 1993).
Koch et al. [2007a] J. Koch, T. M. Yu, J. Gambetta, A. A. Houck, D. I. Schuster, J. Majer, A. Blais, M. H. Devoret, S. M. Girvin, and R. J. Schoelkopf, Charge-insensitive qubit design derived from the cooper pair box, Phys. Rev. A 76, 042319 (2007a).
Leppäkangas et al. [2018] J. Leppäkangas, J. Braumüller, M. Hauck, J.-M. Reiner, I. Schwenk, S. Zanker, L. Fritz, A. Ustinov, M. Weides, and M. Marthaler, Phys. Rev. A 97, 052321 (2018).
Stehli et al. [2023] A. Stehli, J. D. Brehm, T. Wolz, A. Schneider, H. Rotzinger, M. Weides, and A. V. Ustinov, Quantum emulation of the transient dynamics in the multistate landau-zener model, npj Quantum Information 9, 61 (2023).
Breuer et al. [1999] H.-P. Breuer, B. Kappler, and F. Petruccione, Stochastic wave-function method for non-markovian quantum master equations, Physical Review A 59, 1633 (1999).
Vacchini and Breuer [2010] B. Vacchini and H.-P. Breuer, Exact master equations for the non-markovian decay of a qubit, Physical Review A 81, 042103 (2010).
Koch et al. [2007b] J. Koch, T. M. Yu, J. Gambetta, A. A. Houck, D. I. Schuster, J. Majer, A. Blais, M. H. Devoret, S. M. Girvin, and R. J. Schoelkopf, Charge-insensitive qubit design derived from the Cooper pair box, Physical Review A 76, 042319 (2007b).
Majenz et al. [2013] C. Majenz, T. Albash, H.-P. Breuer, and D. A. Lidar, Coarse graining can beat the rotating-wave approximation in quantum markovian master equations, Phys. Rev. A 88, 012103 (2013).
Roccati and Cilluffo [2024] F. Roccati and D. Cilluffo, Controlling markovianity with chiral giant atoms (2024), arXiv:2402.15556 [quant-ph] .
Ahmadi et al. [2024] B. Ahmadi, R. R. Rodríguez, R. Alicki, and M. Horodecki, Approximation scheme and non-hermitian renormalization for the description of atom-field-system evolution, Phys. Rev. A 109, 012408 (2024).
Leggett et al. [1987] A. J. Leggett, S. Chakravarty, A. T. Dorsey, M. P. A. Fisher, A. Garg, and W. Zwerger, Dynamics of the dissipative two-state system, Rev. Mod. Phys. 59, 1 (1987).
Rivas et al. [2010] Á. Rivas, S. F. Huelga, and M. B. Plenio, Entanglement and non-markovianity of quantum evolutions, Physical review letters 105, 050403 (2010).
Lidar et al. [2001] D. A. Lidar, Z. Bihary, and K. Whaley, From completely positive maps to the quantum Markovian semigroup master equation, Chem. Phys. 268, 35 (2001).
Bacon et al. [1999] D. Bacon, D. A. Lidar, and K. B. Whaley, Robustness of decoherence-free subspaces for quantum computation, Phys. Rev. A 60, 1944 (1999).
Gaudin [1960] M. Gaudin, Une démonstration simplifiée du théorème de wick en mécanique statistique, Nuclear Physics 15, 89 (1960).
Fogedby [2022] H. C. Fogedby, Field-theoretical approach to open quantum systems and the Lindblad equation, Physical Review A 106, 022205 (2022).
Van Kampen [1974a] N. G. Van Kampen, A cumulant expansion for stochastic linear differential equations. i, Physica 74, 215 (1974a).
Van Kampen [1974b] N. G. Van Kampen, A cumulant expansion for stochastic linear differential equations. ii, Physica 74, 239 (1974b).
Breuer and Petruccione [2002b] H. Breuer and F. Petruccione, The Theory of Open Quantum Systems (Oxford University Press, 2002).
Ban [2010] M. Ban, Nonequilibrium dynamics of the dispersive Jaynes–Cummings model by a non-Markovian quantum master equation, Journal of Physics A: Mathematical and Theoretical 43, 335305 (2010).
Shabani and Lidar [2005] A. Shabani and D. A. Lidar, Completely positive post-markovian master equation via a measurement approach, Physical Review A 71, 020101 (2005).
Zhang et al. [2022] H. Zhang, B. Pokharel, E. M. Levenson-Falk, and D. Lidar, Predicting non-markovian superconducting-qubit dynamics from tomographic reconstruction, Physical Review Applied 17, 054018 (2022).
Alicki et al. [2006] R. Alicki, D. A. Lidar, and P. Zanardi, Internal consistency of fault-tolerant quantum error correction in light of rigorous derivations of the quantum Markovian limit, Phys. Rev. A 73, 052311 (2006).
Hartmann et al. [2017] M. Hartmann, D. Poletti, M. Ivanchenko, S. Denisov, and P. Hänggi, Asymptotic floquet states of open quantum systems: the role of interaction, New Journal of Physics 19, 083011 (2017).
Haddadfarshi et al. [2015] F. Haddadfarshi, J. Cui, and F. Mintert, Completely positive approximate solutions of driven open quantum systems, Physical Review Letters 114, 130402 (2015).
Yamaguchi et al. [2017] M. Yamaguchi, T. Yuge, and T. Ogawa, Markovian quantum master equation beyond adiabatic regime, Physical Review E 95, 012136 (2017).
Dann et al. [2018] R. Dann, A. Levy, and R. Kosloff, Time-dependent markovian quantum master equation, Physical Review A 98, 052129 (2018).
Meglio et al. [2023] G. D. Meglio, M. B. Plenio, and S. F. Huelga, Time dependent markovian master equation beyond the adiabatic limit (2023), arXiv:2304.06166 [quant-ph] .
Nathan and Rudner [2020] F. Nathan and M. S. Rudner, Universal lindblad equation for open quantum systems, Physical Review B 102, 115109 (2020).
Davidović [2020] D. Davidović, Completely Positive, Simple, and Possibly Highly Accurate Approximation of the Redfield Equation, Quantum 4, 326 (2020).
Davidović [2022] D. Davidović, Geometric-arithmetic master equation in large and fast open quantum systems, Journal of Physics A: Mathematical and Theoretical 55, 455301 (2022).
Winczewski et al. [2021] M. Winczewski, A. Mandarino, M. Horodecki, and R. Alicki, Bypassing the intermediate times dilemma for open quantum system (2021), arXiv:2106.05776 [quant-ph] .
Albash et al. [2012] T. Albash, S. Boixo, D. A. Lidar, and P. Zanardi, Quantum adiabatic markovian master equations, New J. of Phys. 14, 123016 (2012).
Yip et al. [2018] K. W. Yip, T. Albash, and D. A. Lidar, Quantum trajectories for time-dependent adiabatic master equations, Physical Review A 97, 022116 (2018).
Campos Venuti and Lidar [2018] L. Campos Venuti and D. A. Lidar, Error reduction in quantum annealing using boundary cancellation: Only the end matters, Phys. Rev. A 98, 022315 (2018).
Corrington [1961] M. S. Corrington, Applications of the Complex Exponential Integral, Mathematics of Computation 15, 1 (1961).