Non-Markovian evolution of multi-level system interacting with several reservoirs. Exact and approximate

TL;DR

This paper presents an exactly solvable model for multi-level quantum systems interacting with multiple reservoirs at zero temperature, comparing exact and approximate decay and decoherence rates.

Contribution

It introduces a solvable model and systematically compares exact solutions with common approximate master equations for multi-reservoir interactions.

Findings

Exact and approximate rates are quantitatively compared.

Parameter space is classified based on discrepancies between exact and approximate rates.

Insights into the validity of approximate methods in different regimes.

Abstract

An exactly solvable model for the multi-level system interacting with several reservoirs at zero temperatures is presented. Population decay rates and decoherence rates predicted by exact solution and several approximate master equations, which are widespread in physical literature, are compared. The space of parameters is classified with respect to different inequities between the exact and approximate rates.

Figures

Equations

$$\mathcal{H}\equiv (\mathbb{C}\oplus {\mathbb{C}}^N)\otimes \underset{i=1}{\overset{N}{⨂}}{\mathfrak{F}}_b({\mathcal{L}}^2(\mathbb{R})).$$

$${\hat{H}}_S=\sum _i{\epsilon }_i|i⟩⟨i|+\sum _{i\ne j}{J}_{ij}|i⟩⟨j|=0\oplus H_S,i,j=1,\dots ,N,$$

$${\hat{H}}_B=\sum _{i=1}^N\int {\omega }_k{b}_{k,i}^{†}{b}_{k,i}dk.$$

$${\hat{H}}_I=\sum _i\int \left(g_k^{\ast }|0⟩⟨i|\otimes b_{k,i}^{†}+g_k|i⟩⟨0|\otimes b_{k,i}\right)dk,$$

$$\frac{d}{dt}|\Psi (t)⟩=-i\hat{H}|\Psi (t)⟩,$$

$$|\Psi (0)⟩=(|\psi (0)⟩+{\psi }_0(0)|0⟩)\otimes |\Omega ⟩,⟨0|\psi (0)⟩=0,$$

$$|{\psi }_F⟩=|\psi ⟩\oplus \sum _{i=1}^N\int dk{\psi }_{k,i}|k,i⟩,$$

$$|{\hat{\psi }}_F⟩=|\psi ⟩\otimes |\Omega ⟩+|0⟩\otimes \sum _{i=1}^N\int dk{\psi }_{k,i}{b}_{k,i}^{†}|\Omega ⟩.$$

$$|\Psi (t)⟩={\psi }_0(0)|0⟩\otimes |\Omega ⟩+|{\hat{\psi }}_F(t)⟩,$$

$$\frac{d}{dt}|{\psi }_F(t)⟩=-i{H}_F|{\psi }_F(t)⟩,|{\psi }_F(0)⟩=|\psi (0)⟩\oplus 0,$$

$$H_B=\sum _{i=1}^N\int {\omega }_k|k,i⟩⟨k,i|dk,H_I=\sum _{i=1}^N\int \left(g_k^{\ast }|k,i⟩⟨i|+g_k|i⟩⟨k,i|\right)dk.$$

$$\hat{N}=\sum _{i=1}^N|i⟩⟨i|\otimes I+I\otimes \sum _{i=1}^N\int b_{k,i}^{†}{b}_{k,i}dk,$$

$$G(t)=\int |g_k{|}^2{e}^{-i{\omega }_{k}t}dk$$

$$\frac{d}{dt}|\psi (t)⟩=-i{H}_S|\psi (t)⟩-{\int }_0^{t}dsG(t-s)|\psi (s)⟩$$

$$\{\begin{array}{l}\frac{d}{dt}P|{\psi }_F(t)⟩=-iP{H}_{F}P|{\psi }_F(t)⟩-iP{H}_{F}Q|{\psi }_F(t)⟩,\\ \frac{d}{dt}Q|{\psi }_F(t)⟩=-iQ{H}_{F}P|{\psi }_F(t)⟩-iQ{H}_{F}Q|{\psi }_F(t)⟩.\end{array}$$

$$Q|{\psi }_F(t)⟩=e^{-iQ{H}_{F}Qt}Q|{\psi }_F(0)⟩-i{\int }_0^{t}ds{e}^{-iQ{H}_{F}Q(t-s)}Q{H}_{F}P|{\psi }_F(s)⟩.$$

$$\frac{d}{dt}P|{\psi }_F(t)⟩=-iP{H}_{F}P|{\psi }_F(t)⟩$$

$$-P{H}_{F}Q{\int }_0^{t}ds{e}^{-iQ{H}_{F}Q(t-s)}Q{H}_{F}P|{\psi }_F(s)⟩.$$

$$P|{\psi }_F(t)⟩=|\psi (t)⟩,P{H}_{F}P=H_S,Q|{\psi }_F(0)⟩=0,QHQ=H_B,$$

$$Q{H}_{F}P=\sum _{i=1}^N\int g_k^{\ast }|k,i⟩⟨i|dk,P{H}_{F}Q=\sum _{i=1}^N{g}_k|i⟩⟨k,i|dk,$$

$$P{H}_{F}Q{e}^{-iQ{H}_{F}Qt}Q{H}_{F}P=\sum _{i,j=1}^N\int dkd{k}^{\prime }{g}_k{g}_{k^{\prime }}^{\ast }|i⟩⟨k,i|e^{-i{H}_{B}t}|k^{\prime },j⟩⟨j|=$$

$$=\sum _{i,j=1}^N\int dkd{k}^{\prime }{g}_k{g}_{k^{\prime }}^{\ast }|i⟩e^{-i{\omega }_{k}t}{\delta }_{ij}\delta (k-k^{\prime })⟨j|=I_N\int dk{e}^{-i{\omega }_{k}t}|g_k{|}^2=G(t)I_N.$$

$$\frac{d}{dt}|{\psi }_I(t)⟩=-{\int }_0^{t}dsG(t-s)e^{i{H}_S(t-s)}|{\psi }_I(s)⟩$$

$${\rho }_S(t)\equiv {\mathrm{tr}}_R|\Psi (t)⟩⟨\Psi (t)|,$$

$${\rho }_S(t)=0\oplus |\psi (t)⟩⟨\psi (t)|+{\psi }_0(0{)}^{\ast }|\psi (t)⟩⟨0|+{\psi }_0(0)|0⟩⟨\psi (t)|+(1-||\psi (t)|{|}^2)|0⟩⟨0|$$

$${\rho }_{SI}(t)=0\oplus |{\psi }_I(t)⟩⟨{\psi }_I(t)|+{\psi }_0^{\ast }(0)|{\psi }_I(t)⟩⟨0|+{\psi }_0(0)|0⟩⟨{\psi }_I(t)|+(1-||{\psi }_I(t)|{|}^2)|0⟩⟨0|.$$

$$G(t)=g^2{e}^{-\frac{\gamma }{2}|t|-i\epsilon t},\gamma >0,g>0,\epsilon \in \mathbb{R},$$

$$G(t)={\int }_{-\infty }^{+\infty }\frac{d\omega }{2\pi }{e}^{-i\omega t}\mathcal{J}(\omega ),\mathcal{J}(\omega )={\int }_{-\infty }^{+\infty }G(t)e^{i\omega t}dt$$

$$\mathcal{J}(\omega )=\frac{\gamma g^2}{{\left(\frac{\gamma }{2}\right)}^2+(\omega -\epsilon {)}^2},$$

$$|\phi (t)⟩\equiv -ig{\int }_0^{t}ds{e}^{-\left(\frac{\gamma }{2}+i\epsilon \right)(t-s)}|\psi (s)⟩,$$

$$\frac{d}{dt}|\tilde{\psi }(t)⟩=-i{H}_{\mathrm{eff}}|\tilde{\psi }(t)⟩,H_{\mathrm{eff}}=\left(\begin{array}{cc}{H}_S& g{I}_N\\ g{I}_N& \left(\epsilon -i\frac{\gamma }{2}\right)I_N\end{array}\right)$$

$$\frac{d}{dt}|\psi (t)⟩=-i{H}_S|\psi (t)⟩-ig|\phi (t)⟩$$

Full text

Non-Markovian evolution of multi-level system interacting with several reservoirs. Exact and approximate111The research was supported by RSF (project No. 17-71-20154).

**A.E. Teretenkov222Steklov Mathematical Institute of Russian Academy of Sciences, ul. Gubkina 8, Moscow 119991, Russia; Lomonosov Moscow State University, Leninskie Gory 1, Moscow 119991, Russia

E-mail:taemsu@mail.ru**

An exactly solvable model for the multi-level system interacting with several reservoirs at zero temperatures is presented. Population decay rates and decoherence rates predicted by exact solution and several approximate master equations, which are widespread in physical literature, are compared. The space of parameters is classified with respect to different inequities between the exact and approximate rates.

1 Introduction

In this work we consider the evolution of the multi-level system each level of which interacts with its own bosonic reservoir at zero temperature. For the simplicity we assume that the reservoirs are similar and the coupling of each level to its reservoir is also similar (the explicit mathematical model is described in Section 2). The main aim of the work is the comparison of the exact evolution of the reduced density matrix of the system (obtained by the pseudomode method, see Section 3) with the approximate evolution defined by master equations which are widespread in the physical literature. Namely, we consider the Nakajima-Zwanzig equation in the Born approximation, the non-Markovian Redfield equation and the Markovian Redfield equation. Usually in physical literature equations for the reduced density matrix are derived in the following way. The initial Liouville-von Neumann (linear differential) equation is reduced to the Nakajima-Zwanzig (linear integro-differential) equation [1, 2] (or some equivalent equation obtained by exclusion of reservoirs degrees of freedom). The Nakajima–Zwanzig is an exact equation for the reduced density matrix. Then the following four assumptions are subsequently done:

Born approximation [3, p. 131], [4, p. 7]. This equation is also integro-differential, which allows one to consider it as a non-Markovian one [5, 6]. Sometimes this approximation is also called the Redfield approximation or the second-order approximation [7, p. 249], [8, Subsec. 11.2], [9]. 2. 2.

Assumption that the reduced density matrix inside the integral could be taken at the same time as outside, which leads to the non-Markovian Redfield equation [10], [3, p. 132], [11]. Actually, this approximation is very close to Markovianity but only on the long times. 3. 3.

The full Markovian approximation, which leads to the Markovian Redfield equation [12, p. 141], [3, p. 132] , [4, p. 7], . 4. 4.

Secular approximation [3, p. 132], [12, p. 145], which leads to Markovian equations with the Gorini–Kossakowski–Sudarshan–Lindblad (GKSL) generator [14, 13].

At the same time only the equation which obtained after all four assumptions has mathematically strict justification [15, 16]. And its derivation goes back to [17, 18] and is based on the van Hove–Bogolyubov scaling [16, Sec. 1.8]. Moreover, only the last equation guarantees positivity and even complete positivity [14, 13]. At the same time the Redfield equation can violate the positivity, which indeed could be fixed by applying the slippage-operators to the initial conditions [19, 20, 21] or by some other methods [22]. On the other hand, in physical calculations the Redfield equation without secular approximation is frequently used [24, 23]. Some advantage of the Redfield equation in comparison with GKSL equations of the secular approximation are discussed in [25, 26]. For the damped oscillator the Redfield equation can both be translation invariant and have canonical equilibrium state [27] rather than GKSL equation [28]. If the aim is to take into account the non-Markovian effects then only Born approximation is usually done [9]. Thus, it is suggested in physical literature that the equations of the first one or two approximations could be more accurate than the GKSL equation which needs all four approximations. So it is natural to study the exactly solvable model to compare its prediction with the approximate ones.

In Section 2 we present the initial problem for the system and the reservoirs. It is important to mention that we introduce the Hamiltonian which is already in the rotating wave approximation form. Hence we leave the discussion of this approximation by itself out of the range of our study. At the same time, if one derives GKSL equations by the stochastic limit approach, then it is not necessary to assume this approximation and its a corollary of the van Hove–Bogolyubov scaling [16].

In Section 3 we present the pseudomode approach and obtain the exact evolution of the reduced density matrix. This approach was developed in [29, 30, 31, 32, 33, 34]. In [35] we have shown that the Friedrichs model [36] naturally arises as an intermediate step in this approach. In a quite general form the description of the non-Markovian evolution in the Friedrichs approximation was discussed in [37].

In Sections 4 and 5 we present the Nakajima-Zwanzig equation in Born approximation and the Redfield equation (both non-Markovian and Markovian) for our special case, accordingly, and solve them. In our case the Markovian Redfield equation appears to be identical to the one with the secular approximation.

Finally, in Section 6 we compare the solutions for exact and approximate equations. In Conclusions we summarize our results and suggest the directions for the further studies.

2 Schroedinger equation for system and reservoirs

We consider the evolution in the Hilbert space

[TABLE]

Here $\mathbb{C}\oplus\mathbb{C}^{N}$ is a $(N+1)$-dimensional Hilbert space with a pointed one-dimensional subspace which corresponds to the degrees of freedom of the $(N+1)$-level system. Let $|i\rangle,i=0,1,\ldots,\ N$ be an orthonormal basis in such a space and $|0\rangle$ correspond (be collinear) to the pointed subspace. $\mathfrak{F}_{b}(\mathcal{L}^{2}(\mathbb{R}))$ are bosonic Fock spaces which describe the reservoirs. There is one reservoir for each excited level of the system. Let $|\Omega\rangle$ be a vacuum vector for the reservoirs. Let us also introduce the creation and annihilation operators which satisfy the canonical commutation relations: $[b_{k,i},b_{k^{\prime},j}^{\dagger}]=\delta_{ij}\delta(k-k^{\prime})$, $[b_{k,i},b_{k^{\prime},j}]=0$, $b_{k,i}|\Omega\rangle=0$.

We consider the system Hamiltonian in the general form

[TABLE]

without assumption that $\hat{H}_{S}$ is diagonalized in the basis $|i\rangle$. From the physical point of view $|i\rangle$ plays the role of local basis [38, 39]. The only restriction is that it is non-zero only in the excited subspace i.e. in the subspace which is an orthogonal complement to the pointed one-dimensional subspace. The restriction of $\hat{H}_{S}$ to the excited subspace is denoted by $H_{S}$ in formula (1).

The reservoir Hamiltonian is a sum of identical Hamiltonians of free bosonic fields (with the identical dispersion relation $\omega_{k}$)

[TABLE]

The interaction is described by the following Hamiltonian

[TABLE]

i.e. each level interacts only with its own reservoir and the functions $g_{k}$ (called the form factors [40]) are the same for all reservoirs. Let us note that from the physical point of view such a Hamiltonian assumes that the dipole approximation (the only terms which are linear in creation and annihilation operators involved in the Hamiltonian) and the rotating wave approximation (the terms of the form $|i\rangle\langle 0|\otimes b_{k,i}^{\dagger}$ are absent) are justified. The validity of the latter one in a general case is controversial [41, 42], but we do not discuss this question in our study.

We consider the Schroedinger equation

[TABLE]

with the Hamiltonian $\hat{H}=\hat{H}_{S}\otimes I+I\otimes\hat{H}_{B}+\hat{H}_{I}$ and the initial condition

[TABLE]

i.e. we assume the initial condition to be completely factorized as in [35]. Thus, non-Markovian effects related to non-factorized initial state are also out of the range of our study. A quite general approach for such effects was suggested in [43]. We also assume that the initial state of the system is pure and the reservoirs are in the vacuum states, i.e. in their Gibbs states at zero temperature.

3 Pseudomode approach and exact evolution

First of all let us show that 1-particle restriction of the Hamiltonian $\hat{H}$ is related to a generalized Friedrichs model. Let us introduce an injective map $\hat{}:h\rightarrow\mathcal{H},$ where $h=\mathbb{C}^{N}\oplus\bigoplus_{i=1}^{N}\mathcal{L}^{2}(\mathbb{R})$ which for any $|\psi_{F}\rangle\in h$ of the form

[TABLE]

where $\psi_{k,i}$ are from the Schwartz space $\mathcal{S}(\mathbb{R})$ and $\int dk\cdot|k,i\rangle$ is a Fourier transform for each $i=1,\ldots,m$, defines $|\hat{\psi}_{F}\rangle\in\mathcal{H}$ by the formula

[TABLE]

Such a map we called one-particle second quantization in [35], it is a special realization of the idea to consider non-composite systems as composite ones suggested in [44, 45, 46, 47].

Theorem 1**.**

The solution of the Cauchy problem (4)-(5) has the form

[TABLE]

where $|\psi_{F}(t)\rangle$ is a solution of the Cauchy problem for the Friedrichs model:

[TABLE]

where $H_{F}=H_{S}\oplus H_{B}+H_{I}$, $H_{S}$ is defined by formula (1) and

[TABLE]

The proof of this theorem is based on the direct substitution. But the deep reason for the preservation of 0-particle and 1-particle subspaces consists in the presence of the integral of motion

[TABLE]

which is nothing else but the total number of particles in reservoir and excitations in the system.

This N-level generalized Friedrichs model is close to the one considered in [48], but there was only one reservoir coupled to all the excited states of the system.

Now we are going to obtain the reduced evolution but for a state vector rather than a density matrix as it is usually done for the master equation derivation. Let us define the projection $P$ on the linear subspace $\mathbb{C}^{N}$ in the Hilbert space $h$ and the projection $Q=I_{N}-P$ on the orthogonal complement to this subspace.

Theorem 2**.**

Let the integral

[TABLE]

converge for all $t\in\mathbb{R}_{+}$ and define the continuous function $G(t)$, then $|\psi(t)\rangle=P|\psi_{F}(t)\rangle$ satisfies the integro-differential equation

[TABLE]

with the initial condition $|\psi(t)\rangle|_{t=0}=|\psi(0)\rangle$, where $H_{S}$ is defined by formula (1).

Proof.

Let us prove this theorem by the projection approach to emphasize the closeness of this approach to the Nakajima-Zwanzig projection one. Let us represent equation (7) as the system

[TABLE]

Let us solve the second equation as a linear differential with respect to $Q|\psi_{F}(t)\rangle$ considering $-iQH_{F}P|\psi_{F}(t)\rangle$ as an inhomogeneity

[TABLE]

Substituting into the first one we obtain

[TABLE]

Taking into account

[TABLE]

we have

[TABLE]

Substituting into (3) we obtain (10).

∎

Let us note that the solution of equation (10) with the initial condition $|\psi(t)\rangle|_{t=0}=|\psi(0)\rangle$ and a continuous function $G(t)$ exists and is unique [58, Sec. 2.1]

Thus, equation (10) is an analog of the Nakajima-Zwanzig equation but for the state vector. Such approach is called Feshbach projection approach [49, 50]. Its application to open systems and its generalization describing completely positive evolution of a reduced density matrix could be found in [51].

In fact, we have already used the pseudomode approach in the interaction picture in [35], but we have not dwellt on it explicitly. So let us introduce here the equation in the interaction picture in the explicit form.

Corollary 1**.**

The vector $|\psi_{I}(t)\rangle\equiv e^{iH_{S}t}|\psi(t)\rangle$ satisfies the equation

[TABLE]

with the initial condition$|\psi_{I}(t)\rangle|_{t=0}=|\psi(0)\rangle$.

It is important that $|\psi(t)\rangle$ and $|\psi_{I}(t)\rangle$ unambiguously define the evolution of a reduced density matrix or of a reduced density matrix in the interaction picture, accordingly.

Lemma 1**.**

The reduced density matrix

[TABLE]

where $\mathrm{tr}_{R}\;$ is a partial trace with respect to the space $\bigotimes\limits_{i=1}^{N}\mathfrak{F}_{b}(\mathcal{L}^{2}(\mathbb{R}))$ and $|\Psi(t)\rangle$ is defined by formula (6), has the form

[TABLE]

as well as the reduced density matrix in the interaction picture $\rho_{SI}(t)\equiv e^{i\hat{H}_{S}t}\rho_{S}(t)e^{-i\hat{H}_{S}t}$ has the form

[TABLE]

This lemma could be proved by direct calculation (see also discussion in [35, Sec. 4]).

In this paper we focus on the case, when $G(t)$ has the special form

[TABLE]

although the approach presented below could be immediately generalized to the case of the combination of exponentials $G(t)=\sum_{j}g_{j}^{2}e^{-\frac{\gamma_{j}}{2}|t|-i\varepsilon_{j}t}$ (see [35]). In physics usually the spectral density $\mathcal{J}(\omega)$ [31] defined by

[TABLE]

rather than the function $G(t)$ is given from the experimental data. In our case, when $G(t)$ has form (14), we have

[TABLE]

i.e. nothing else but Lorentzian spectral density.

Theorem 3**.**

Let $|\psi(t)\rangle$ be a solution of (10) with the initial condition $|\psi(t)\rangle|_{t=0}=|\psi(0)\rangle$ in the case, when $G(t)$ has form (14), then $|\tilde{\psi}(t)\rangle\equiv|\psi(t)\rangle\oplus|\varphi(t)\rangle\in\mathbb{C}^{N}\oplus\mathbb{C}^{N}$, where

[TABLE]

satisfies the Schroedinger equation with a non-Hermitian (dissipative) Hamiltonian

[TABLE]

with the initial condition $|\tilde{\psi}(t)\rangle=|\psi(0)\rangle\oplus 0$.

Proof.

Substituting (14) in (10) and taking into account (16) we obtain

[TABLE]

and differentiating (16) with respect to $t$ we have

[TABLE]

By combining these two ordinary differential equations into one for the vector $|\tilde{\psi}(t)\rangle=|\psi(t)\rangle\oplus|\varphi(t)\rangle$ we obtain (17). ∎

Actually, the non-Hermitian part of $H_{\rm eff}$ is very close to the optical potential concept from the nuclear physics [52, 53].

Corollary 2**.**

The vector $|\tilde{\psi}_{I}(t)\rangle\equiv|\psi_{I}(t)\rangle\oplus e^{iH_{S}t}|\varphi(t)\rangle$ satisfies the equation

[TABLE]

with the initial condition $|\tilde{\psi}_{I}(t)\rangle=|\psi(0)\rangle\oplus 0$.

Equation (18) could be solved in the global basis, i.e. the eigenbasis of $H_{S}$ [38, 39]. We will numerate the global basis by the Greek letters:

[TABLE]

We also define the energy detunings

[TABLE]

Corollary 3**.**

Let us decompose $|\psi_{I}(t)\rangle=\sum_{\alpha=1}^{N}\psi_{\alpha}(t)|\alpha\rangle,\oplus e^{iH_{S}t}|\varphi(t)\rangle=\sum_{\alpha=1}^{N}\varphi_{\alpha}(t)|\alpha\rangle$, then the coefficients of the decomposition satisfy the systems of equations

[TABLE]

The eigenvalues of the matrix $-iH_{I,\rm eff}$ equal

[TABLE]

We are interested in the population decay rates and the decoherence rates. They are defined by the real parts of the eigenvalues

[TABLE]

The population decay rates in the global basis are characterized by the minimal absolute value of these eigenvalues

[TABLE]

as well as the decoherence rates are characterized by

[TABLE]

At the end of this section let us concentrate on the fact that equation (17) is equivalent to a GKSL equation for a density matrix in the $(2N+1)$-dimensional space.

Theorem 4**.**

The $(2N+1)$ by $(2N+1)$ matrix

[TABLE]

satisfies the GKSL equation

[TABLE]

where $|0\rangle,|1\rangle,\ldots,|N\rangle,|\tilde{1}\rangle,\ldots,|\tilde{N}\rangle$ is the basis of $\mathbb{C}^{2N+1}$.

The proof of this theorem could be found in [35, Proposition 1]. Thus, we have built up the GKSL description of the non-Markovian evolution of the reduced density matrix at the cost of introduction of additional degrees of freedom (see recent discussion in [35], [54]). This GKSL equation (similar to the Friedrichs model) could be considered as a one-particle restriction of models in larger Hilbert spaces. For example, multi-level generalizations of the Jaynes-Cummings model with dissipation [55] could be considered in such a way.

4 Nakajima-Zwanzig equation in Born approximation

The Nakajima-Zwanzig equation in the Born approximation [3, p. 131], [56, p.9] at zero temperature (the state of the reservoir is $|\Omega\rangle\langle\Omega|$)

[TABLE]

where $\hat{H}_{I}(t)$ is the iteration Hamiltonian in the interaction picture, i.e.

[TABLE]

In our case $\hat{H}_{S},\hat{H}_{B},\hat{H}_{I}$ are defined by formulae (1), (2) and (3). This equation could be obtained in the second order of the perturbation theory with respect to coupling constant between a system and a reservoir [3, Subsec. 9.1.1], [57].

Lemma 2**.**

The interaction Hamiltonian (23) has the form

[TABLE]

where

[TABLE]

Moreover, $b_{\alpha,i}(t)|\Omega\rangle=0$ and the following commutation relations are held

[TABLE]

where $G(t)$ is defined by (9).

Proof.

Let us represent the interaction Hamiltonian in the global basis

[TABLE]

Taking into account $e^{i\hat{H}_{S}t}|\alpha\rangle\langle 0|e^{-i\hat{H}_{S}t}=e^{iE_{\alpha}t}|\alpha\rangle\langle 0|$ and $e^{i\hat{H}_{B}t}b_{k,i}e^{-i\hat{H}_{B}t}=e^{-i\omega_{k}t}b_{k,i}$ we obtain

[TABLE]

Thus, we obtain (24). $b_{\alpha,i}(t)|\Omega\rangle=0$ and the second of commutation relations (26) follows from definition (25) of $b_{\alpha,i}(t)$ which is a linear combination of annihilation operators. Let us calculate

[TABLE]

Thus, we obtain the first of commutation relations (26). ∎

Lemma 3**.**

In our case, i.e., when $\hat{H}_{I}(t)$ is defined by formula (24), equation (22) takes the form

[TABLE]

where the projection

[TABLE]

is introduced.

Proof.

Let us expand the following expression

[TABLE]

Thus, only two terms could be calculated, the other two could be obtained by interchange of $t$ and $s$. By lemma 2 we obtain

[TABLE]

and

[TABLE]

where $[\Pi,\hat{H}_{S}]=[0\oplus I_{N},0\oplus H_{S}]=0$ is used. By substituting the obtained expressions into (4) and then into (22) we obtain (3). ∎

Equation (3) has a convolution form and, hence, could be solved by means of Laplace transform [58, p. 30]. These convolution form is preserved in the Schroedinger picture:

[TABLE]

Lemma 4**.**

The solution of equation (3) could be represented in the form

[TABLE]

where the vector-valued function $|\psi_{I}(t)\rangle$ satisfies equitation (12) with the initial condition $|\psi_{I}(t)\rangle|_{t=0}=|\psi(0)\rangle$, the matrix-valued function $\sigma(t)$ satisfies the equitation

[TABLE]

and the scalar function $\rho_{00}(t)$ could be defined from the normalization condition $\rho_{00}(t)=1-\mathrm{Tr}\;\sigma(t)$.

Thus, without additional assumptions about function $G(t)$ we obtain that the coherences between excited and ground state of the system (defined by the vector $|\psi_{I}(t)\rangle$) coincide for exact (13) and approximate (30) solutions.

Now let us consider the special case, when the function $G(t)$ is defined by formula (14). Similar to the solution of equation (12) for $|\psi_{I}(t)\rangle$ (see theorem 3) one could use a similar technique to solve equation (31). This technique is also very close to auxiliary density matrices method developed in [59, 60].

Lemma 5**.**

Let $\sigma(t)$ satisfy equation (31), $G(t)$ have form (14) and

[TABLE]

then $\sigma(t)$ and $X(t)$ satisfy the system of linear matrix equations with constant coefficients

[TABLE]

As well as for the case of equations (18) it is natural to solve these system in the global basis.

Corollary 4**.**

If one decomposes $\sigma(t)$ and $X(t)$ in the global basis $\sigma(t)=\sum_{\alpha\beta}\sigma_{\alpha\beta}(t)|\alpha\rangle\langle\beta|$, $X(t)=\sum_{\alpha\beta}X_{\alpha\beta}(t)|\alpha\rangle\langle\beta|$, then the coefficients of the decomposition satisfy the systems

[TABLE]

So the solution of system (32) is reduced to the calculation of $3\times 3$-matrix exponential. Hence, the evolution is defined by the eigenvalues of such a matrix which could be expressed as zeros of the characteristic polynomial

[TABLE]

Namely,

[TABLE]

In particular the population decay rate in the global basis is defined by zeros of the characteristic polynomial

[TABLE]

with real coefficients. Namely,

[TABLE]

5 Redfield equation

In physical literature [3, p. 132], [12, p. 141] the following equation is frequently used:

[TABLE]

which is called the (non-Markovian) Redfield master equation [10]. In (36) we have taken into account zero temperature of the reservoirs. This equation differs from equation (22) by the fact that $\rho_{SI}$ inside the integral is a function of $t$ rather than $s$. Thus, equation (36) has a time-local generator but this generator is time-dependent. Moreover, at $t=0$ the integral in the left hand side of (22) vanishes, which leads to zero population and coherence decay rates at $t=0$.

After the full Markovian approximation [3, p. 132], [12, p. 145] Redfield equation (36) takes the form

[TABLE]

which we call the Markovian Redfield master equation.

Lemma 6**.**

In our case, i.e., when $\hat{H}_{I}(t)$ is defined by (24), equation (36) takes the form

[TABLE]

and equation (37) takes the form

[TABLE]

where the projection $\Pi$ is defined by formula (28).

The proof of this lemma is analogous to the proof of lemma 3.

Corollary 5**.**

Let $G(t)$ be defined by (14), $\gamma>0$ and

[TABLE]

then equation (6) takes the form

[TABLE]

and equation (6) takes the form

[TABLE]

where

[TABLE]

Proof.

Actually, in general, (6) could be represented in form (41) with

[TABLE]

In our case, when $G(t)$ is defined by (14), integration with respect to $s$ leads to (40). In the case of (42) one would integrate to infinity and obtain (42). ∎

In particular, equation (6) has a time-independent generator, i.e. non-secular terms do not appear and it is already in the secular approximation. Actually, this is a result of the rotating wave approximation in initial interaction Hamiltonian (3). One could show [3, p. 136] that in general the Redfield equation in the secular approximation has the GKSL form, but let us present the GKSL form for special case (42) explicitly:

[TABLE]

where the braces denote the anticommutator $\{A,B\}\equiv AB+BA$.

Let us note that similarly to $\hat{H}_{S}$ the matrix $\hat{Y}(t)$ defined by (40) is supported on the subspace $\mathbb{C}^{N}$ of $\mathbb{C}\oplus\mathbb{C}^{N}$, i.e. $\hat{Y}(t)=0\oplus Y(t)$.

Corollary 6**.**

Equation (41) has a solution of form (30), where the vector-valued function $|\psi_{I}(t)\rangle$ is the solution of the Cauchy problem

[TABLE]

and the matrix-valued function $\sigma(t)$ is a solution of the Cauchy problem

[TABLE]

Analogously, equation (42) has a solution of form (30), where the vector-valued function $|\psi_{I}(t)\rangle$ is a solution of the Cauchy problem

[TABLE]

and the matrix-valued function $\sigma(t)$ is a solution of the Cauchy problem

[TABLE]

Corollary 7**.**

If one decomposes $|\psi(t)\rangle$ and $\sigma(t)$ in the global basis $|\psi_{I}(t)\rangle=\sum_{\alpha}\psi_{\alpha}(t)|\alpha\rangle,\sigma(t)=\sum_{\alpha\beta}\sigma_{\alpha\beta}(t)|\alpha\rangle\langle\beta|$, then Cauchy problems (43), (44) take the form

[TABLE]

[TABLE]

accordingly, and Cauchy problems (45), (46) take the form

[TABLE]

[TABLE]

accordingly, where

[TABLE]

and the function $\mathcal{J}(\omega)$ is defined by formula (15).

Similar to (20) it is natural to define the decoherence rates between excited and ground states of the system by

[TABLE]

and they are time-dependent. The population decay rates and the decoherence rates should be defined as follows

[TABLE]

At the large time $t\rightarrow+\infty$ the rates predicted by the Redfield equation approach the rates predicted by the Redfield equation in the secular approximation

[TABLE]

[TABLE]

Although, for the purposes of our article just the population decay rates and the decoherence rates are needed, one could obtain the explicit expression for the populations and coherences. Therefore let us state the following lemma.

Lemma 7**.**

The solutions of Cauchy problems (47) and (48) have the form

[TABLE]

[TABLE]

accordingly.

To prove that one should directly integrate (47) and (48). Let us note that this integration with respect to time could be regarded as an analog of time-deformation [61], but it depends on the global basis.

6 Comparison

First of all let us compare decoherence rates between excited and ground states for different equations mentioned in the previous sections.

Theorem 5**.**

Let $\eta_{\alpha 0}^{\rm exact}(\gamma,\Delta E_{\alpha},g)$, $\eta_{\alpha 0}^{\rm B}(\gamma,\Delta E_{\alpha},g)$ and $\eta_{\alpha 0}^{\rm GKSL}(\gamma,\Delta E_{\alpha},g)$ be the real-valued functions defined for $\gamma>0,g>0,\Delta E_{\alpha}\in\mathbb{R}$ by (21),(33) and (50), accordingly, then

[TABLE]

[TABLE]

where

[TABLE]

The proof of (52) follows immediately from lemma 4 and the discussion after that. The proof of (53) is based on the direct comparison of (21) and (50) for the ground-excited case.

Corollary 8**.**

Let $\eta_{\alpha 0}^{\rm R}(\gamma,\Delta E_{\alpha},g,t^{*})$ be defined by (49) for $\gamma>0,g>0,\Delta E_{\alpha}\in\mathbb{R},t\geqslant 0$ and $\eta_{\alpha 0}^{\rm exact}(\gamma,\Delta E_{\alpha},g)$ be defined by (21) for $\gamma>0,g>0,\Delta E_{\alpha}\in\mathbb{R}$. Let $F_{1}(u,v)$ be defined by (54). If $F_{1}\left(\frac{\Delta E_{\alpha}^{2}}{g^{2}},\frac{\gamma^{2}}{g^{2}}\right)>0$, then there exists $t^{*}\in\mathbb{R}$ (may be not unique) such that $\eta_{\alpha 0}^{\rm R}(\gamma,\Delta E_{\alpha},g,t^{*})=\eta_{\alpha 0}^{\rm exact}(\gamma,\Delta E_{\alpha},g)$.

The proof follows from the fact that $\eta_{\alpha 0}^{\rm R}(t)$ run over all the points from $[0,\eta_{\alpha 0}^{\rm GKSL}]$ (but may be do not only them) and theorem 5.

Now we are going to compare the population decay rates of excited states for different equations mentioned in previous sections. For that we need a technical lemma based on Routh-Hurwitz criterion.

Lemma 8**.**

Let $x\in\mathbb{R}$, then the roots of the cubic polynomial

[TABLE]

lie in the half-plane $\mathrm{Re}\;\lambda<x$ if and only if

[TABLE]

Proof.

Let us substitute $\lambda=\tilde{\lambda}+x$, where $x\in\mathbb{R}$, into (55)

[TABLE]

Applying the Routh-Hurwitz criterion [62, p. 194] (in the particular case of a cubic equation it is also called the Vyshnegradsky criterion [63, p. 37]) to this polynomial in $\tilde{\lambda}$ we obtain: $\mathrm{Re}\;\tilde{\lambda}<0$, i.e. $\mathrm{Re}\;\lambda<x$, if and only if (56) are satisfied. ∎

Theorem 6**.**

Let $\eta_{\alpha}^{\rm exact}(\gamma,\Delta E_{\alpha},g)$, $\eta_{\alpha}^{\rm B}(\gamma,\Delta E_{\alpha},g)$ and $\eta_{\alpha}^{\rm GKSL}(\gamma,\Delta E_{\alpha},g)$ be the real-valued functions defined for $\gamma>0,g>0,\Delta E_{\alpha}\in\mathbb{R}$ by (20),(35) and (51), accordingly, then

[TABLE]

where

[TABLE]

[TABLE]

[TABLE]

and $F_{1}(u,v)$ is defined by formula (54).

Proof.

1. To compare $\eta_{\alpha}^{\rm exact}$ and $\eta_{\alpha}^{\rm B}$ let us apply lemma 8 to the equation $f_{\alpha}(\lambda)=0$, where $f_{\alpha}(\lambda)$ is defined by formula (34), i.e.

[TABLE]

and let $x=\mathrm{Re}\;\lambda_{\alpha,\pm}^{\rm exact}$ which is defined by (19). The third of inequalities (56) appears to be satisfied without further assumptions and the first one follows from the second one. Thus, only the second inequality has to be satisfied. It takes the form

[TABLE]

where the function $F_{2}(u,v)$ is defined by formula (6).

2. The comparison of $\eta_{\alpha}^{\rm exact}=2\eta_{\alpha 0}^{\rm exact}$ and $\eta_{\alpha}^{\rm GKSL}=2\eta_{\alpha 0}^{\rm GKSL}$ is equivalent to the comparison of $\eta_{\alpha 0}^{\rm exact}$ and $\eta_{\alpha 0}^{\rm GKSL}$, which was done in theorem 5.

3. The comparison of $\eta_{\alpha}^{\rm GKSL}$ and $\eta_{\alpha}^{\rm B}$ is similar to the first one. The only distinction is that one should assume $x=-\eta_{\alpha 0}^{\rm GKSL}$, then functions (6) and (59) occur from lemma 8. ∎

The results of theorems 5 and 6 are presented in figure 1. The following feature should be noted.

Corollary 9**.**

There are only two points of the half-plane $\gamma/g>0,\Delta E_{\alpha}/g\in\mathbb{R}$ such that

[TABLE]

These points are approximately $\Delta E_{\alpha}/g\approx\pm 0.55$, $\gamma/g\approx 3.55$.

If one wants to have simpler conditions for the comparison of $\eta_{\alpha}^{\rm exact}$ and $\eta_{\alpha}^{\rm B}$ than defined by (6), then one could use the following proposition.

Corollary 10**.**

If $\gamma>\sqrt{64/3}g\approx 4.62g$ or $|\Delta E_{\alpha}|>\sqrt{\frac{1}{3}\left(-4+\sqrt[3]{44-3\sqrt{177}}+\sqrt[3]{44+3\sqrt{177}}\right)}g\approx 0.81g$, then $\eta_{\alpha}^{\rm exact}(\gamma,\Delta E_{\alpha},g)<\eta_{\alpha}^{\rm B}(\gamma,\Delta E_{\alpha},g)$. If $\gamma\leqslant 3\sqrt{2g^{2}-4\Delta E_{\alpha}^{2}}$, then $\eta_{\alpha}^{\rm exact}(\gamma,\Delta E_{\alpha},g)>\eta_{\alpha}^{\rm B}(\gamma,\Delta E_{\alpha},g)$.

Figure 1 shows, where one population decay rates are greater than other ones, but it does not show how close they are. So we have done it numerically in figure 2. One could see that $\eta_{\alpha}^{\rm B}$ is close to $\eta_{\alpha}^{\rm exact}$ not only near the curve $\eta_{\alpha}^{\rm B}=\eta_{\alpha}^{\rm exact}$ but also for scientifically small $g$. The region, where $\eta_{\alpha}^{\rm B}$ is close to $\eta_{\alpha}^{\rm GKSL}$, is also wide and not concentrated only near the curve $\eta_{\alpha}^{\rm B}=\eta_{\alpha}^{\rm GKSL}$. Interestingly the second region is not even close to be a subset of the first one, so there are areas of parameters when GKSL equations could provide the better fit than the non-Markovian Nakajima-Zwanzig equation in the Born approximation for the populations decay rate which is experimentally observable. But for a sufficiently narrow peak in the spectral density (small $\gamma$) the non-Markovian Nakajima-Zwanzig equation in the Born approximation reproduces the population decay better than the GKSL one.

7 Conclusions

We have considered the model of the multi-level system interacting with several local reservoirs at zero temperature. We have compared the population decay rates and the decoherence rates for exact solution and several approximate master equations: the Nakajima-Zwanzig equation in the Born approximation, the Redfield equation. It was shown:

Both the initial model and the approximate master equations are exactly solvable in the global basis. 2. 2.

The Nakajima-Zwanzig equation in the Born approximation gives an exact result for the coherences between the excited states and the ground states (without additional assumptions for the spectral density of the reservoir). 3. 3.

The conditions for all possible inequalities between excited-ground decoherence rates and population rates in the global basis for the exact, Born and Markovian Redfield cases are fully characterized by theorems 5 and 6. 4. 4.

Both numerically and analytically we have shown that there exist the cases when the Markovian GKSL equation reproduces the population decay better than the non-Markovian Nakajima-Zwanzig equation in the Born approximation, but this is not the case for the sufficiently narrow spectral density.

In our opinion the following directions for the further studies could be fruitful.

Application to the real systems. As in [35] in this study we were inspired by vibronic non-Markovian phenomena in light harvesting complexes. The approach described here could be applied to the one-exciton models [65, 64] of the Fenna-Matthews-Olson complexes at cryogenic temperatures. For them the non-Markovian phenomena were experimentally observed [66], which leads to the sufficient interest in the quantum phenomena in photosynthetic systems [67, 68, 69, 70, 71]. 2. 2.

Finite-temperature analysis. The fact that (8) is an integral of motion for our system allows one to separate the equations with fixed number of particles. So may be the exact finite-temperature solutions could be obtained on this way. 3. 3.

Multiple Lorentzian and non-Lorentzian generalization of the results described. Multiple Lorentzian peaks case for the spectral density could be considered in a straight forward way by the methods from [29, 30, 31, 32, 33, 34, 35]. Non-Lorentzian case could be dealt with by general Laplace transform methods, but we think that the approach from [72, 73] could provide more physical insight.

8 Acknowledgments

The author thanks A. S. Trushechkin for sufficient help in setting the main goals of this study and fruitful discussion at all the steps of the study. The author thanks I. V. Volovich, S. V. Kozyrev, B. O. Volkov and A. I. Mikhailov for the useful discussion of the problems considered in the work.