Pseudomode approach and vibronic non-Markovian phenomena in light harvesting complexes

TL;DR

This paper explores the pseudomode approach's theoretical foundations and its application to non-Markovian vibronic phenomena in light-harvesting complexes, with extensions to finite-temperature scenarios.

Contribution

It connects the pseudomode approach with well-known quantum models and applies it to describe non-Markovian effects in biological complexes.

Findings

Connection established between pseudomode approach and Friedrichs/Jaynes-Cummings models

Application to Fenna-Matthews-Olson complexes with experimental estimations

Discussion on extending the approach to finite temperatures

Abstract

The pseudomode approach is discussed with an emphasis to Gorini-Kossakowski-Sudarshan-Lindblad form of this approach. The connection of the pseudomode approach with solutions of both the Friedrichs model and Jaynes-Cummings model with dissipation at zero temperature is shown. The obtained results are applied to non-Markovian phenomena description in the Fenna-Matthews-Olson complexes, estimations based on experimental data are done. A deformation approach to generalization of the pseudomode approach to the finite-temperature case is discussed.

Equations

$$\rho =R\oplus 0+(1-\mathrm{Tr}R)|0⟩⟨0|,$$

$$\frac{d}{dt}R(t)=-i{H}_{\mathrm{eff}}R(t)+iR(t)H_{\mathrm{eff}}^{†},$$

$$R(0)=R_0,$$

$$H_{\mathrm{eff}}=H-\frac{i}{2}\sum _{l=1}^n{\gamma }_l|l⟩⟨l|,$$

$$\frac{d}{dt}\mathrm{Tr}R(t)=\mathrm{Tr}\frac{H_{\mathrm{eff}}-H_{\mathrm{eff}}^{†}}{i}R(t),$$

$$\frac{d}{dt}\mathrm{Tr}R(t)=-\sum _l{\gamma }_l⟨l|R(t)|l⟩,$$

$$\frac{d}{dt}\rho (t)=-i[H\oplus 0,\rho (t)]+\sum _l\left(L_l\rho (t)L_l^{†}-\frac{1}{2}{L}_l^{†}{L}_l\rho (t)-\frac{1}{2}\rho (t)L_l^{†}{L}_l\right),$$

$$\rho (0)={\rho }_0,$$

$$\frac{d}{dt}\rho (t)=\frac{d}{dt}R(t)\oplus 0-\frac{d}{dt}\mathrm{Tr}R(t)|0⟩⟨0|.$$

$$\frac{d}{dt}R(t)\oplus 0$$

$$=-i[H,R(t)]\oplus 0-\frac{1}{2}\sum _{l=1}^n{\gamma }_l|l⟩⟨l|R(t)\oplus 0-\frac{1}{2}\sum _{l=1}^n{\gamma }_{l}R(t)|l⟩⟨l|\oplus 0=$$

$$=-i[H\oplus 0,R(t)\oplus 0]-\frac{1}{2}\sum _{l=1}^n{\gamma }_l|l⟩⟨0||0⟩⟨l|R(t)\oplus 0-\frac{1}{2}\sum _{l=1}^n{\gamma }_{l}R(t)\oplus 0|l⟩⟨0||0⟩⟨l|=$$

$$=-i[H\oplus 0,\rho (t)]-\frac{1}{2}\sum _l\left(L_l^{†}{L}_l\rho (t)+\rho (t)L_l^{†}{L}_l\right),$$

$$\frac{d}{dt}\mathrm{Tr}R(t)|0⟩⟨0|=-\sum _l{\gamma }_l⟨l|R(t)|l⟩|0⟩⟨0|=-\sum _l{\gamma }_l|0⟩⟨l|R(t)\oplus 0|l⟩⟨0|=-\sum _l{L}_l\rho (t)L_l^{†}.$$

$$\frac{d}{dt}|\tilde{\psi }(t)⟩=-i{H}_{\mathrm{eff}}|\tilde{\psi }(t)⟩,$$

$$|l⟩$$

$$|0⟩$$

$$\hat{\rho }=\sum _{l,k=0}^n{\rho }_{lk}|\hat{l}⟩⟨\hat{k}|.$$

$$\frac{d}{dt}\hat{\rho }(t)=-i[\hat{H},\hat{\rho }(t)]+\sum _{l=1}^n{\gamma }_l\left(a_l\hat{\rho }(t)a_l^{†}-\frac{1}{2}{a}_l^{†}{a}_l\hat{\rho }(t)-\frac{1}{2}\hat{\rho }(t)a_l^{†}{a}_l\right),\hat{H}=\sum _{l,k=1}^n{H}_{lk}{a}_l^{†}{a}_k.$$

$$\frac{d}{dt}\hat{\rho }(t)=-i[{\omega }_1{\sigma }^{†}\sigma +{\omega }_2{b}^{†}b+g({\sigma }^{†}b+\sigma b^{†}),\hat{\rho }(t)]+\gamma \left(b\hat{\rho }(t)b^{†}-\frac{1}{2}{b}^{†}b\hat{\rho }(t)-\frac{1}{2}\hat{\rho }(t)b^{†}b\right).$$

$$\frac{d}{dt}\rho (t)=-i[{\omega }_1|1⟩⟨1|+{\omega }_2|2⟩⟨2|+g(|1⟩⟨2|+|2⟩⟨1|),\rho (t)]+$$

$$+\gamma \left(|0⟩⟨2|\rho (t)|2⟩⟨0|-\frac{1}{2}|2⟩⟨2|\rho (t)-\frac{1}{2}\rho (t)|2⟩⟨2|\right)$$

$${\mathrm{Tr}}_{{\mathcal{X}}_i,i\in I}\hat{\rho }={\hat{\rho }}_{\overline{I}}+\left({\rho }_{00}+\sum _{l\in I}{\rho }_{ll}\right)|\hat{0}⟩⟨\hat{0}|,{\hat{\rho }}_{\overline{I}}=\sum _{l,k\in \overline{I}}{\rho }_{lk}|\hat{l}⟩⟨\hat{k}|,$$

$${\mathrm{Tr}}_I\rho ={\rho }_{\overline{I}}+\left({\rho }_{00}+\sum _{l\in I}{\rho }_{ll}\right)|0⟩⟨0|,{\rho }_{\overline{I}}=\sum _{l,k\in \overline{I}}{\rho }_{lk}|l⟩⟨k|.$$

$$H_F=\int {\omega }_k|k⟩⟨k|dk+{\omega }_1|1⟩⟨1|+\int \left(g_k^{\ast }|k⟩⟨1|+g_k|1⟩⟨k|\right)dk.$$

$${\hat{H}}_{SB}=\int {\omega }_{k}I\otimes b_k^{†}{b}_{k}dk+{\omega }_1|1⟩⟨1|\otimes I+\int \left(g_k^{\ast }|0⟩⟨1|\otimes b_k^{†}+g_k|1⟩⟨0|\otimes b_k\right)dk.$$

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

$$|\psi (0)⟩={\psi }_1(0)|1⟩,$$

$$|\psi (t)⟩={\psi }_1(t)|1⟩+\int dk{\psi }_k(t)|k⟩,$$

$$\frac{d{\psi }_1(t)}{dt}=-i{\omega }_1{\psi }_1(t)-{\int }_0^{t}d{t}^{\prime }G(t-t^{\prime }){\psi }_1(t^{\prime }),$$

$${\psi }_k(t)=-i{g}_k^{\ast }{\int }_0^{t}d\tau e^{i{\omega }_k(\tau -t)}{\psi }_1(\tau ).$$

$$\frac{d}{dt}{\psi }_1(t)$$

$$\frac{d}{dt}{\psi }_k(t)$$

Full text

Pseudomode approach and vibronic non-Markovian phenomena in light harvesting complexes111The research was supported by RSF (project No. 17-71-20154 ).

**A.E. Teretenkov222Steklov Mathematical Institute of Russian Academy of Sciences, Moscow, Russia; Lomonosov Moscow State University, Moscow, Russia

E-mail:taemsu@mail.ru**

The pseudomode approach is discussed with an emphasis to Gorini-Kossakowski-Sudarshan- Lindblad form of this approach. The connection of the pseudomode approach with solutions of both the Friedrichs model and Jaynes- Cummings model with dissipation at zero temperature is shown. The obtained results are applied to non-Markovian phenomena description in the Fenna-Matthews-Olson complexes, estimations based on experimental data are done. A deformation approach to generalization of the pseudomode approach to the finite-temperature case is discussed.

1 Introduction

We use the pseudomode approach developed in [1, 2, 3, 4, 5] to describe non-Markovian evolution caused by strong coupling of a system and a reservoir. We apply the approach to the theory of non-Markovian exciton transport phenomena in light harvesting complexes. Vibrational degrees of freedom play the role of the reservoir in such systems. And states arising as a result of the strong coupling between excitons and these vibrational degrees of freedom are called vibrons. Thus, we refer to non-Markovian phenomena due to strong coupling as vibronic ones.

The pseudomode approach allows one to describe the evolution of a system which is strongly coupled with a reservoir if the following conditions are met (more specific and strict statements could be found below in the main text of the article):

The system has a finite number of levels. 2. 2.

The reservoir is at zero temperature. 3. 3.

Interaction with the reservoir is described by a spectral density with a finite number of poles in the lower half-plane.

Then this approach allows one to describe evolution of the system in terms of solutions for a finite number of linear ordinary differential equations. Similar to [2, 3] we are especially interested if it is possible to represent these equations as a master equation with a Gorini-Kossakowski- Sudarshan- Lindblad (GKSL) generator [6, 7].

In contrast with [1, 2, 3, 4, 5] we separate three subjects, which, in our opinion, are logically independent, in the pseudomode approach. These three subjects correspond to three sections 2-4 of our article. Such a presentation allows us to clarify the relation between the pseudomode approach and the solution of the Jaynes- Cummings model with dissipation at zero temperature [8] as well as the Friedrichs model [9].

In section 2 we define a map which allows one to construct an equation with a GKSL generator by evolution with non-Hermitian Hamiltonian and non-increasing normalization. Let us note that such a map was discussed in [2]. But GKSL equations obtained there were equations in the infinite-dimensional Hilbert space. Hence their solution is much more complex then the solution of initial equations with a non-Hermitian Hamiltonian. Moreover, we weaken conditions from [2] on the non-Hermitian Hamiltonian.

In section 3 we start with the GKSL equation from the previous section and construct a GKSL equation in a larger (possibly infinite-dimensional) Hilbert space. We present the former GKSL equation which is a restriction of the latter one on 0-particle and 1-particle subspaces. Actually this is how the infinite-dimensional GKSL equations mentioned above arise in [2]. They have a form of equations which describes the Jaynes- Cummings model with dissipation at zero temperature. However, the solution of these equations from [8, formula (5.2)] is just based on the restriction of the initial model on 0-particle and 1-particle subspaces. The results of section 3 allow us to generalize these solution on the case of multiple two-level atoms being strongly coupled with multiple bosonic modes.

In section 4 we present the pseudomode approach just as the method to solve the Friedrichs model which reduces this model to evolution with non-Hermitian Hamiltonian. And only after that we use the results of sections 2 and 3 to show how this approach allows one to solve a non-Markovian problem for the spin-boson model in the rotating wave approximation in terms of finite-dimensional GKSL equations, which was the original purpose to introduce this approach in [10].

In section 5 we apply the pseudomode approach to describe non-Markovian decay of the two-level system in resonance with the frequency of the spectral density peak similar to paragraph 10.1.2 from [11]. However, we describe the transition into Markovian limit in terms of Van Hove-Bogolyubov scaling [12] for the reduced density matrix rather than time-convolutionless projection operator method from [11].

In section 6 we discuss the application of the obtained results to describe the non-Markovian phenomena in the Fenna-Matthews-Olson (FMO) light harvesting complexes at cryogenic temperatures. We show that the presence of a strong coupling with a reservoir leads to long-living coherences and populations oscillations in non-degenerate systems. Such a behavior was discovered experimentally in [13].

In section 7 we develop a deformation approach to description of non-Markovian phenomena at a finite temperature. Namely, we heuristically write a generator in the non-Markovian case. But then we show that the reduced density matrix in the Van Hove-Bogolyubov limit satisfies a Markovian master equation and thus one could regard our generator as a deformation of the Markovian one. Moreover, in the zero temperature limit our generator coincides with the exact generator from the previous section.

In the conclusion section we summarize the main results of our work and discuss several directions for further studies.

2 Normalization reconstruction

Definition 1**.**

Let $R$ be such a nuclear self-adjoint non-negatively definite operator in the Hilbert space $\mathcal{H}$ that $\operatorname{Tr}R\leqslant 1$ (hereinafter we call this operator non-normalized density matrix). Let us map this operator into the operator $\rho$ in $\mathcal{H}\oplus\mathbb{C}$ by formula

[TABLE]

where $|0\rangle$ is the basis vector in $\mathcal{H}\oplus\mathbb{C}$ corresponding to the subspace $\mathbb{C}$. We call such a map normalization reconstruction.

The operator $\rho$ is obvious to be a density matrix (a self-adjoint non-negatively definite operator which trace equals one) in $\mathcal{H}\oplus\mathbb{C}$. Let us note that if the matrix $R$ is already ”well-normalized”, i.e. $\operatorname{Tr}R=1$, then, nevertheless, we obtain the new matrix $\rho=R\oplus 0$ in the ”larger” Hilbert space $\mathcal{H}\oplus\mathbb{C}$ according to our definition. It could seem to be unnatural but we will see below that our definition is useful. If $R$ is projector on the vector $|\psi\rangle$ we define that the pure state $|\psi\rangle\oplus 0$ is a normalization reconstruction of the pure state $|\psi\rangle$.

Let us consider a finite-dimensional Hilbert space $\mathbb{C}^{n}$ and a matrix $H_{\rm eff}\in\mathbb{C}^{n\times n}$. Let us pose a Cauchy problem for the matrix-valued function $R(t):\mathbb{R}_{+}\rightarrow\mathbb{C}^{n\times n}$

[TABLE]

[TABLE]

where $R_{0}$ is a non-normalized density matrix and † is a Hermitian conjugation.

Definition 2**.**

We call equation (2) a von Neumann equation with the non-Hermitian Hamiltonian $H_{\rm eff}$.

Definition 3**.**

Let the solution $R(t)$ of Cauchy problem (2)-(3) satisfy $\operatorname{Tr}R(t)\leqslant 0$ for arbitrary initial non-normalized density matrix $R_{0}$0, then we call equation (2) von Neumann equation with the non-Hermitian Hamiltonian $H_{\rm eff}$ and non-increasing normalization.

Lemma 1**.**

Equation (2) is a von Neumann equation with a non-Hermitian Hamiltonian and non-increasing normalization if and only if $H_{\rm eff}$ could be presented in the following form

[TABLE]

where $H=H^{\dagger}$, $|l\rangle,l=1,\ldots,n$ is a certain orthonormal basis in $\mathbb{C}^{n}$: $\langle l|l^{\prime}\rangle=\delta_{ll^{\prime}}$ and $\gamma_{l}$ are non-negative real numbers: $\gamma_{l}\geqslant 0$.

Proof**.**

1. Let (2) be a solution of the von Neumann equation with a non-Hermitian Hamiltonian and non-increasing normalization. Let us apply the trace to both sides of equation (2)

[TABLE]

then $\frac{d}{dt}\operatorname{Tr}R(t)|_{t=0}\leqslant 0$ for any initial condition is equivalent to $\frac{H_{\rm eff}-H_{\rm eff}^{\dagger}}{i}\leqslant 0$ as one could assume $R(0)=|v\rangle\langle v|$ . Then non-negatively definite Hermitian matrix could be represented in the form $\frac{H_{\rm eff}-H_{\rm eff}^{\dagger}}{i}=-\sum_{l}\gamma_{l}|l\rangle\langle l|$, where $\gamma_{l}\geqslant 0$ are opposite to eigenvalues of this matrix and $|l\rangle$ are eigenvectors of the matrix forming an orthonormal basis. Besides, let us define $H=\frac{1}{2}(H_{\rm eff}+H_{\rm eff}^{\dagger})$. Thus, we have (4).

2. Inversely, if we assume (4), then it follows from (2) that

[TABLE]

which ensures non-increasing normalization for an arbitrary matrix $R(t)\geqslant 0$. The fact that the condition $R(t)\geqslant 0$ is conserved by equation (2) follows from the explicit solution $R(t)=e^{-iH_{\rm eff}t}R(0)e^{iH_{\rm eff}^{\dagger}t}$. Indeed, if $R(0)\geqslant 0$, then $\langle v|R(t)|v\rangle=\langle v^{\prime}|R(0)|v^{\prime}\rangle\geqslant 0$, where $|v^{\prime}\rangle=e^{iH_{\rm eff}^{\dagger}t}|v\rangle$. ∎

Propostion 1**.**

Let $R_{t}$ be a solution of Cauchy problem (2)-(3) for the von Neumann equation with the non-Hermitian Hamiltonian $H_{\rm eff}$ and non-increasing normalization. Then the density matrix $\rho_{t}$ obtained by normalization reconstruction (1) is a solution of the Cauchy problem for the GKSL equation in the Lindblad form

[TABLE]

[TABLE]

where $L_{l}=\sqrt{\gamma_{l}}|0\rangle\langle l|$, $\rho_{0}$ is obtained by normalization reconstruction of $R_{0}$ and $H$, $\gamma_{l}$, $|l\rangle$ are defined in lemma 1.

Proof**.**

The matrix $\rho(t)$ has form (1) with respect to definition 1. Hence, we have

[TABLE]

Taking into account (2) and lemma 1 we obtain

[TABLE]

where the orthogonality condition $\langle 0|l\rangle=0,l\neq 0$ and the definition $L_{l}=\sqrt{\gamma_{l}}|0\rangle\langle l|$ are used. Taking into account equation (5) we have

[TABLE]

By substituting these expression into (7) we obtain (6). ∎

Just as a usual von Neumann equation in case of a pure initial state could be reduced to a Schroedinger equation, the solution of the Cauchy problem (2)-(3) with the initial condition $R_{0}=|\tilde{\psi}_{0}\rangle\langle\tilde{\psi}_{0}|$ could be represented in the form $R(t)=|\tilde{\psi}(t)\rangle\langle\tilde{\psi}(t)|$, where $|\tilde{\psi}(t)\rangle$ is a solution of the equation

[TABLE]

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

Definition 4**.**

Let us call equation (8), where $H_{\rm eff}$ is defined by formula (4), a Schroedinger equation with the non-Hermitian Hamiltonian and non-increasing normalization.

It is equation (8) that was considered in [2]. In this section we have introduced its natural generalization (2).

3 One-particle second quantization

Let us consider a Hilbert space of the form $\mathbb{C}^{n}\oplus\mathbb{C}$. Let $\mathcal{X}_{i}$ be auxiliary Hilbert spaces with $\dim\mathcal{X}_{i}\geqslant 2$. Let us introduce the linear injection $\hat{}:\mathbb{C}^{n}\oplus\mathbb{C}\rightarrow\otimes_{i=1}^{n}\mathcal{X}_{i}$ which is defined on the basis by the following rule

[TABLE]

We call such a map a one-particle second quantization. The term ”one-particle” refers to the fact that we consider only the states in $\otimes_{i=1}^{n}\mathcal{X}_{i}$ which correspond to the image of $\mathbb{C}^{n}\oplus\mathbb{C}$ for such an injective map. In a similar way let us map the density matrices $\rho=\sum_{l,k=0}^{n}\rho_{lk}|l\rangle\langle k|$ into the matrices

[TABLE]

Propostion 2**.**

Let $a_{l}$ be arbitrary operators such that $a_{l}^{\dagger}|\hat{0}\rangle=|\hat{l}\rangle$ and $a_{l}|\hat{l}\rangle=|\hat{0}\rangle$. Let $\rho_{t}$ satisfy equation (6) with coefficients defined in proposition 1. Then $\hat{\rho}_{t}$ defined by formula (10) satisfies the equation

[TABLE]

The proof of this proposition is based on direct substitution.

There are three usual choices for operations $a_{l}$ in the physical applications: bosonic annihilation operators $b_{l}$ satisfying canonical commutation relations $[b_{l},b_{k}^{\dagger}]=\delta_{lk},[b_{l},b_{k}]=0$, fermionic annihilation operators $f_{l}$ satisfying canonical anticommutation relations $\{f_{l},f_{k}^{\dagger}\}=\delta_{lk},\{f_{l},f_{k}\}=0$ or two-level annihilation operators $\sigma_{l}$ satisfying $\{\sigma_{l},\sigma_{l}^{\dagger}\}=I,\{\sigma_{l},\sigma_{l}\}=0,[\sigma_{l},\sigma_{k}]=0,[\sigma_{l},\sigma_{k}^{\dagger}]=0,l\neq k$. It is important to note that it is not necessary to choose the operators of the same type for all the indices. For example, one could choose some of $a_{l}$ to be bosonic and some of them to be two-level. This is the case of interaction of two-level atoms and modes of a bosonic (electromagnetic, phonon etc.) field. The operators of different kind are assumed to commute with each other.

In particular, both the Jaynes- Cummings model (introduced in [14], bibliographical review could be found in [15]) and its dissipative zero-temperature generalization [8, 10, 16, 17] have form (11). Namely, the Jaynes- Cummings model with dissipation at zero temperature is defined by the equation

[TABLE]

Proposition 2 allows one to reduce the problem of finding one-particle solutions of this equation to the solution of the following equation:

[TABLE]

which has form (11).

Let us also note that if we assume all the operators $a_{l}$ in equation (11) to be either only bosonic ones or only fermionic ones then we obtain a particular case of the GKSL equation with a quadratic generator [18, 19, 20, 21, 22]. Gaussian solutions could be obtained for such equations. However, in this case one could also obtain one-particle solutions as an example of the simplest non-Gaussian solutions for GKSL editions with quadratic generators.

The most important feature of the map (9) is that it allows one to consider the initial system as a composite one. In particular, one could calculate partial traces with respect to subspace $\mathcal{X}_{i}$, where $i\in I$ and $I$ is a subset of the set of numbers from $1$ to $n$. Moreover, if the density matrix is obtained by formula (10), then the particular trace has very simple expression:

Propostion 3**.**

Let $\hat{\rho}=\sum_{l,k=0}^{n}\rho_{lk}|\hat{l}\rangle\langle\hat{k}|$, then the partial trace with respect to spaces indexed by $I$ could be calculated by the formula:

[TABLE]

where $\overline{I}=\{0,\ldots,n\}\setminus I$.

One could prove this proposition by direct calculation. Let us emphasize that according to our definitions [math] is always contained in $\overline{I}$ rather than $I$.

Thus, if one traces out the subspaces indexed by $I$, then the result is just a matrix without entries corresponding to the indices from $I$ and with vacuum state population changed to its initial value plus the sum of populations corresponding to the indices from $I$. The latter one guarantees normalization conservation.

Moreover, this proposition motivates the concept of the trace with respect to a set of indices for initial $(n+1)\times(n+1)$ matrices $\rho$ by formula

[TABLE]

In general, the embedding of the density matrices of so-called non-composite systems (see [26]) into spaces which have a tensor product form appeared to be quite fruitful for the derivation of an entropy inequalities. The non-composite system is a system described by a Hilbert space which is not of the tensor product form or at least this form is not given explicitly. General discussion (far beyond our particular one-particle case) and more complete bibliography on this topic could by found in [23, 24, 25]. In regard to light harvesting complexes the trace with respect to an index set was discussed in [27] for mutual information calculation. The results of this work were discussed in [28] from the holography point of view.

4 Friedrichs model and pseudomode approach

The Friedrichs model was introduced in [9] and has been widely discussed in scientific papers in connection with the perturbation theory for continuous spectrum [29], non-integrate large Poincare’s systems and spontaneous time-symmetry breaking [31, 30, 32, 33], Gamow vectors in the rigged Hilbert spaces [34, 36, 35]. Moreover, this model has generalizations which allow one to obtain non-Markovian equations which are second order differential equations with respect to time [37]. We rather focus on the features of the Friedrichs model which are related to pseudomode theory developed in [1, 2, 3, 4, 5].

The Friedrichs model is defined by the Hamiltonian in $\mathbb{C}\oplus L^{2}(\mathbb{R})$:

[TABLE]

The spin-boson model in the rotating wave approximation is one of possible second quantizations according to proposition 2. It is defined by the Hamiltonian in $\mathbb{C}^{2}\otimes\mathfrak{F}_{b}(L^{2}(\mathbb{R}))$

[TABLE]

This is the model being solved in [1, 2, 3, 4, 5]. However, as the temperature of the reservoir was assumed equal to zero, then the one-particle solutions of the Schroedinger equation with Hamiltonian (15) were considered. Hence, the problem was actually reduced to the solution of the Schroedinger equation with the Hamiltonian $H_{F}$, even though it was not mentioned in these works.

Let us consider the Schroedinger equation with such a Hamiltonian

[TABLE]

with the initial condition

[TABLE]

then the following proposition is held.

Propostion 4**.**

The Cauchy problem (14), (16), (17) has the solution

[TABLE]

where $\psi_{1}(t)$ satisfies the equation

[TABLE]

where $G(t)=\int dke^{-i\omega_{k}t}|g_{k}|^{2}$ and $\psi_{k}(t)$ could be expressed it terms of $\psi_{1}(t)$ by the formula

[TABLE]

Actually the proof of this proposition could be found in [2]. The main idea consists in the substitution of (18) into (16). As a result, we obtain the following equations

[TABLE]

We integrate the latter equation with initial condition (17), i.e. $\psi_{k}(0)=0$, and obtain (20). Then we substitute (20) into the former one and obtain (19).

The spectral density $J(\omega)$ is usually defined in physical applications rather than the function $G(t)$. These functions are related by the Fourier transform:

[TABLE]

The psedomode approach is essentially a method for solving integro-differential equation (14) in the case when $G(t)$ is a finite linear span of exponents:

[TABLE]

where $\gamma_{l}>0$, $\omega_{l}\in\mathbb{R}$ and $g_{l}$ are complex in the general case. The spectral densities with a finite number of poles in the lower half-plane lead to such a form of the function $G(t)$ (see [2]). In particular, the spectral density

[TABLE]

leads to (21), i.e. the spectral density has the form of a sum of a finite number of peaks.

Propostion 5**.**

Let $\psi_{1}(t)$ be a solution of equation (19) with $G(t)$ of form (21). Let us consider the Hilbert space $\mathbb{C}^{n+1}$ with basis $|1\rangle$, $|\tilde{l}\rangle$. Let us define the non-Hermitian Hamiltonian $H_{\rm eff}\in\mathbb{C}^{(n+1)\times(n+1)}$ by the formula

[TABLE]

then the vector

[TABLE]

where

[TABLE]

satisfies the equation

[TABLE]

Proof**.**

Taking into account definition (24) we represent equation (19) in the following form

[TABLE]

Differentiating (24) we have

[TABLE]

The obtained equations coincide with equation (25) in the component-wise representation. ∎

In work [2] it was discussed if the equation with non-Hermitian Hamiltonian of form (22) can be represented in a GKSL form. However, the general method for such a representation was obtained only for the real $g_{l}$ case. At the same time, proposition 1 allows one to do it in the more general case when the matrix $-i(H_{\rm eff}-H_{\rm eff}^{*})$ is non-positive definite. Nevertheless, GKSL form (6) could be defined simpler if $g_{l}$ are real, as $-i(H_{\rm eff}-H_{\rm eff}^{*})$ has already been diagonalized and, hence, $L_{l}=\sqrt{\gamma_{l}}|0\rangle\langle l|$ and $H=\omega_{1}|1\rangle\langle 1|+\sum_{l=1}^{n}\left(\omega_{l}|\tilde{l}\rangle\langle\tilde{l}|+g_{l}(|\tilde{l}\rangle\langle 1|+|1\rangle\langle\tilde{l}|)\right)$.

Similar to proposition 2 the solution of the Schroedinger equation with Hamiltonian (15) with the initial condition $|\hat{\psi}(0)\rangle=\psi_{0}(0)|1\rangle+\psi_{1}(0)|1\rangle$ has the form $|\hat{\psi}(0)\rangle=\psi_{0}(0)|0\rangle\otimes|0\rangle+\psi_{1}(t)|1\rangle\otimes|0\rangle+\int dk\psi_{k}(t)|0\rangle\otimes b_{k}^{\dagger}|0\rangle$, where $\psi_{1}(t)$ and $\psi_{k}(t)$ are defined by (19) and (20). According to proposition 3 we calculate the reduced density matrix $\hat{\rho}_{S}(t)\equiv\operatorname{Tr}_{\mathfrak{F}_{b}(L^{2}(\mathbb{R}))}|\hat{\psi}(t)\rangle\langle\hat{\psi}(t)|$ by formula (12):

[TABLE]

If the density matrix $\rho$ (which is defined by (1) and satisfies GKSL equation (6)) is a result of the normalization recovery map for the non-normalized density matrix $|\tilde{\psi}(t)\rangle\langle\tilde{\psi}(t)|$, where $|\tilde{\psi}(t)\rangle$ is defined formula (23), then the partial trace of $\rho$ with respect to indices $\tilde{l}$ takes the form

[TABLE]

by formula (13).

Thus, we consider the spin-boson model in the rotating wave approximation with the reservoir which is initially in the vacuum state (at zero temperature) and $G(t)$ has form (21). Then the reduced density matrix for this model is equal to the reduced density matrix $\hat{\rho}_{S}(t)$, i.e. the trace of the matrix $\rho(t)$ with respect to the indices, where $\rho(t)$ is the solution of GKSL equation (6). Hence, the states $|\tilde{l}\rangle$ were called pseudomodes in [1], because the unitary evolution of the system and the reservoir with infinite-number degrees of freedom is equivalent to GKSL evolution in the finite-dimensional Hilbert space from the point of view of reduced density matrix.

5 Non-Markovian resonance decay

Let us consider the case when there is only one term in formula (21), i.e. $G(t)=g^{2}e^{-\left(\frac{\gamma}{2}+i\omega_{1}\right)t}$. Moreover, we assume $g$ to be real and the resonance condition $\omega_{0}=\omega_{1}$ to be satisfied. Then formula (22) takes the form

[TABLE]

As the Markovian limit is considered below and it is usually expressed in the interaction representation, so we also transform our Hamiltonian into the interaction representation

[TABLE]

and we write equation (25) in the interaction representation:

[TABLE]

The solution in the case of initial condition $\psi_{1}(0)=1,\psi_{\tilde{1}}(0)=0$ has the form

[TABLE]

where $\Delta=\frac{1}{4}\sqrt{\gamma^{2}-16g^{2}}$. The solution after taking the trace:

[TABLE]

In our last paragraph we have mainly reproduced in our notation the results of paragraph 10.1.2 from [11]. However, the transition to the Markovian limit was expressed there in terms of time-dependent population transfer rate which occurs in the time-convolutionless projection operator method [38, 39], [11, Section 9.2]. But we use Van Hove-Bogolyubov scaling as it was done in the stochastic limit (see [12], particularly the discussion in sections 1.8 and 1.27). Let us consider this solution as a function of time and the coupling constant $\rho_{S}(t,g)$. Then the transition to the Markovian limit is described by the formula

[TABLE]

Thus, the decay rate in the Markovian limit $\gamma_{0}=\frac{4g^{2}}{\gamma}$.

In the strong coupling limit the oscillating dynamics occurs

[TABLE]

Thus, the transition from the Markovian regime to the strongly non-Markovian one corresponds to the transition from the relaxation regime to the oscillation one.

Proposition 1 allows us to write GKSL equation (6) for the density matrix $\rho_{t}$ with the coefficients

[TABLE]

Namely, we have the equation

[TABLE]

Let us define $D_{lk}(\rho)=|k\rangle\langle l|\rho|l\rangle\langle k|-\frac{1}{2}|l\rangle\langle l|\rho-\frac{1}{2}\rho|l\rangle\langle l|$, $h_{kl}(\rho)=-i[|k\rangle\langle l|+|l\rangle\langle k|,\rho]$. This equation in terms of Markovian decay rate $\gamma_{0}$ takes the form:

[TABLE]

Propostion 6**.**

Let $\rho(t,g,\gamma_{0})$ be a solution of the Cauchy problem

[TABLE]

with the initial condition

[TABLE]

Then $\rho_{M}(t,\gamma_{0})\equiv\lim\limits_{\lambda\rightarrow 0}\rho_{S}\left(\frac{t}{\lambda^{2}},\lambda g,\lambda^{2}\gamma_{0}\right)$ satisfies the equation

[TABLE]

The proof of this proposition is based on direct calculation of the reduced density matrix which takes form

[TABLE]

Then, as $\lambda\rightarrow 0$, the limit takes the form

[TABLE]

Let us note that besides the time-convolutionless projection operator method mentioned above the non-Markovian resonance decay was studied by means of master equations with a time-non-local generator in [40].

6 Non-Markovian phenomena in FMO complexes

The following features were discovered for FMO complexes at cryogenic temperature $\beta^{-1}=77\,\text{K}$:

Long-living coherences (about 1 ps), while the Markovian models predicted coherence lifetime not exceeding 250 fs [13]. 2. 2.

Oscillation of the populations in the case of absence of degeneracy, which also does not coincide with the Markovian dynamics (see the discussion in [57, p. 154]).

These discoveries gave rise to the sufficient interest in the quantum phenomena in photosynthetic systems [41, 42, 43, 44, 45].

In spite of the fact that realistic models consider the FMO complex in the one-exciton approximation (see [47, 27] for the details) as a 7-level system [27], we follow [48, 49] and consider an essentially simplified dimer model of the FMO complex. We also consider the one-exciton approximation and assume the states $|1\rangle$ and $|0\rangle$ from the previous section to correspond to the so-called global basis (see the discussion in [46]).

Data from [49] suggest $\omega_{0}=202\,\text{cm}^{-1},\beta^{-1}\approx 53\,\text{cm}^{-1}$, then $n=\frac{1}{e^{\beta\omega_{0}}-1}\approx 0.02$, hence one could use zero-temperature approximation $n\approx 0$. Moreover, if we assume the Huang-Rhys factor $S=0.02$ as in [49], then the coupling constant $g=\sqrt{S}\omega_{0}\approx 29\,\text{cm}^{-1}$. The typical transition time [13] predicted by the Markovian models equals 250 fs, which corresponds the Markovian coherence decay rate $\frac{\gamma_{0}}{2}=133\,\text{cm}^{-1}$ . Let us evaluate $\frac{\gamma}{4}=g^{2}/\gamma_{0}\approx 3.2\,\text{cm}^{-1}$, i.e. $\frac{4}{\gamma}$ corresponds to 10.4 ps. Thus, our model predicts the sufficient increase of coherence lifetime, which was observed experimentally. Although our prediction could overestimate the increase as we assume an ideal resonance between reservoir modes and the transition and as we consider an oversimplified model for the FMO complex. Let us note that the condition for the transition to the oscillatory regime is held $\gamma^{2}-16g^{2}<0$. Moreover, our model predicts the frequency of population oscillations $|\Delta|\approx 29\,\text{cm}^{-1}$. The presence of such oscillations does not coincide with the Markovian models in the case of absence of degeneracy but was observed experimentally [13]. Let us emphasize that the presence of degeneracy also allows one to explain such oscillations by the presence of the dark states [50, 51].

Let us also note that the data of two-dimensional echo spectroscopy were actually observed in the experiment [52, 53], namely the dependence of diagonal peaks and cross-peaks on so-called population time was observed. But if we consider four diagrams [53, p. 178] making contribution to the cross-peaks, then two of them turn to be zero in the one-exciton approximation, one of them makes population-time-independent contribution at zero temperature. And the last one is just defined by coherences dynamics. Two diagrams contribute to diagonal peaks, one of them is also population-time-independent at zero temperature and the other is defined by the population dynamics. Thus, the relation between the spectroscopy data and reduced density matrix dynamics is direct in our case (for example it is not the case for [49]).

The pseudomode approach uses the low-temperature approximation. As $n=\frac{1}{e^{\beta\omega_{0}}-1}$ contains only the product of the reverse temperature and the frequency, then it could be used for the description of high-frequency peaks considered in [54]. However, in the next section we offer a possible way to remove this restriction.

7 Finite temperature case and Markov chain deformation

Proposition 6 from section 5 allows one two consider non-Markovianity as a ”deformation” of the initial Markov chain changing the decay to the interaction with a mediating level and a subsequent decay. An important feature of such a deformation is that the initial Markovian dynamics for the reduced density matrix is reproduced as $\lambda\rightarrow 0$. So let us extend such an idea to the finite-temperature case. Namely, let us consider the evolution with the generator

[TABLE]

In many aspects this generator is similar to the generator from [55] for the three-level system interacting with several reservoirs and a coherent field.

The next proposition generalizes proposition 6 guaranteeing transition to the Markovian equation without Lamb and Stark shifts for the reduced density matrix (see [11, p. 154], [56]) as $\lambda\rightarrow 0$.

Propostion 7**.**

Let $\rho(t,g,\gamma_{0},n)$ be a solution of the Cauchy problem

[TABLE]

Then $\rho_{M}(t,\gamma_{0},n)\equiv\lim\limits_{\lambda\rightarrow 0}\rho_{S}\left(\frac{t}{\lambda^{2}},\lambda g,\lambda^{2}\gamma_{0},n\right)$ satisfies the equation

[TABLE]

The proof the this proposition is similar to the proof of proposition 6. But the expression for $\rho_{S}(t,g,\gamma_{0},n)$ (a fortiori for $\rho(t,g,\gamma_{0},n)$) is much more complex, so we show here only the expression for

[TABLE]

which corresponds to the solution of the Markovian equation.

This proposition allows one to use generator (26) for description of the non-Markovian phenomena in the light harvesting complexes at finite (in particular, room) temperature.

8 Conclusion

The non-Markovian resonance decay at zero temperature was described by means of the pseudomode approach. The obtained results explain experimentally observed phenomena in the FMO complexes. Unlike [49] we use much simpler equations which at the same time allow one to examine the role of non-Markovian phenomena in arising long-living coherences and population oscillations in a non-degenerate system. The deformation approach was suggested. It allows one to obtain equations at finite temperature.

In our opinion, the following directions for the further studies should be mentioned:

To study the non-Markovian phenomena in more detailed models of the FMO complex as well as other light harvesting complexes at cryogenic temperatures and at room temperatures but for high frequency peaks in spectral density. 2. 2.

To apply the heuristic equations from section 7 for description of the non-Markovian phenomena at room temperature. 3. 3.

To compare the results of calculations by the heuristic equations and by other methods (time-convolutionless projection operator method, hierarchical equations of motion etc.) for the simplest systems. 4. 4.

To expand the theory of Markov chain deformation, in particular, to study how the deformation of one transition in a quantum Markov chain influences the deformation of other transitions. 5. 5.

To attempt to justify the deformation approach strictly beyond zero temperature approximation.

Thus, the results of this work could serve as a basis for the further studies.

9 Acknowledgments

The author thanks I. V. Volovich, S. V. Kozyrev, A. I. Mikhailov and A. S. Trushechkin for fruitful discussion of the problems considered in the work.