Calculation of coherences in Foerster and modified Redfield theories of excitation energy transfer

TL;DR

This paper extends Foerster and modified Redfield theories to include the dynamics of quantum coherences, aligning well with exact numerical methods and sometimes outperforming standard Redfield results.

Contribution

It derives formulas for coherences within these theories, overcoming a key limitation and improving accuracy in modeling excitation energy transfer.

Findings

Formulas for coherences match hierarchical equations of motion results.

Modified Redfield coherence formulas can outperform standard Redfield.

Enhanced theoretical framework improves understanding of energy transfer dynamics.

Abstract

Foerster and modified Redfield theories play one of the central roles in the description of excitation energy transfer in molecular systems. However, in the present state, these theories describe only the dynamics of populations of local electronic excitations or delocalized exciton eigenstates, respectively, i.e., the diagonal elements of the density matrix in the corresponding representation. They do not give prescription for propagating the off-diagonal elements of the density matrix (coherences). This is commonly accepted as a limitation of these theories. Here we derive formulas for the dynamics of the coherences in the framework of Foerster and modified Redfield theories and, thus, remove this limitation. These formulas provide excellent correspondence with numerically exact calculations according to the hierarchical equations of motion. Also we show that, even within the range of applicability of the standard Redfield theory, the formulas for coherences derived in the framework of the modified Redfield theory provide, in some cases, more precise results.

Figures

Equations

$$H=H^{\mathrm{el}}+H^{\mathrm{Coul}}+H^{\mathrm{ph}}+H^{\text{el-ph}}+H^{\mathrm{reorg}},$$

$$H^{\mathrm{el}}=\sum _{n=1}^N|n⟩{\epsilon }_n^0⟨n|,$$

$$H^{\mathrm{Coul}}=\sum _{n=1}^N\sum _{m>n}^N(J_{nm}|n⟩⟨m|+\text{h.c.}),$$

$$H^{\mathrm{ph}}=\sum _{n=1}^N{H}_n^{\mathrm{ph}},$$

$$H_n^{\mathrm{ph}}=\sum _i\left(\frac{p_{ni}^2}{2M_{ni}}+\frac{1}{2}{M}_{ni}{\omega }_{ni}^2{q}_{ni}^2\right),$$

$$H^{\text{el-ph}}=\sum _{n=1}^N|n⟩u_n⟨n|,u_n=\sum _i{M}_{ni}{\omega }_{ni}^2{d}_{ni}{q}_{ni},$$

$$H^{\mathrm{reorg}}=\sum _{n=1}^N|n⟩{\lambda }_n⟨n|,{\lambda }_n=\frac{1}{2}\sum _i{M}_{ni}{\omega }_{ni}^2{d}_{ni}^2,$$

$$\dot{\rho }(t)=-i[H,\rho (t)].$$

$$H=H_0+H^{\prime },$$

$$H_0=\sum _{\alpha }|\alpha ⟩H_{\alpha }⟨\alpha |,H_{\alpha }=⟨\alpha ⟩H\alpha ,$$

$$H^{\prime }=\sum _{\alpha \ne \beta }|\alpha ⟩H_{\alpha \beta }^{\prime }⟨\beta |,H_{\alpha \beta }^{\prime }=⟨\alpha ⟩H\beta .$$

$${\dot{\rho }}^{(I)}(t)=-i[V(t),{\rho }^{(I)}(t)]\equiv -i\mathcal{L}(t){\rho }^{(I)}(t),$$

$$\mathcal{P}A=\sum _{\alpha }{\mathrm{Tr}}_B(A_{\alpha \alpha })|\alpha ⟩⟨\alpha |{\rho }_{\alpha },A_{\alpha \alpha }=⟨\alpha ⟩A\alpha .$$

$$[\mathcal{P},{\mathcal{L}}_0]=0,\mathcal{P}\mathcal{L}(t)\mathcal{P}=0,$$

$$\rho (t)=\mathcal{P}\rho (t)+\mathcal{Q}\rho (t).$$

$$\mathcal{P}\dot{\rho }(t)$$

$$\mathcal{Q}\dot{\rho }(t)$$

$$\begin{gathered} \mathcal{Q}\rho (t)={\mathrm{T}}_{+}\exp \left\{-i{\int }_0^t\mathcal{Q}\mathcal{L}(\tau )d\tau \right\}\mathcal{Q}\rho (0) \\ -i{\int }_0^{t}d\tau {\mathrm{T}}_{+}\exp \left\{-i{\int }_{\tau }^t\mathcal{Q}\mathcal{L}({\tau }^{\prime })d{\tau }^{\prime }\right\}\mathcal{L}(\tau )\mathcal{P}\rho (\tau ), \end{gathered}$$

$$\begin{gathered} {\mathrm{T}}_{+}\exp \left\{-i{\int }_0^{t}f(\tau )d\tau \right\} \\ =1+\sum _{n=1}^{\infty }(-i{)}^n{\int }_0^{t}d{\tau }_1{\int }_0^{{\tau }_1}d{\tau }_2\dots {\int }_0^{{\tau }_{n-1}}d{\tau }_n \\ f({\tau }_1)f({\tau }_2)\cdots f({\tau }_n) \end{gathered}$$

$$\mathcal{Q}\rho (t)=-i{\int }_0^t\mathcal{L}(\tau )\mathcal{P}\rho (\tau )d\tau .$$

$$\mathcal{P}\dot{\rho }(t)=-{\int }_0^{t}d\tau \mathcal{P}\mathcal{L}(\tau )\mathcal{L}(0)\mathcal{P}\rho (t-\tau ).$$

$$\mathcal{P}\dot{\rho }(t)=\left\{-{\int }_0^{\infty }d\tau \mathcal{P}\mathcal{L}(\tau )\mathcal{L}(0)\right\}\mathcal{P}\rho (t).$$

$$\mathcal{P}\rho (t)=\sum _{\alpha }{p}_{\alpha }(t)|\alpha ⟩⟨\alpha |{\rho }_{\alpha },$$

$$p_{\alpha }(t)={\mathrm{Tr}}_B⟨\alpha ⟩\rho (t)\alpha ,$$

$${\dot{p}}_{\alpha }(t)=\sum _{\beta \ne \alpha }[K_{\alpha \beta }{p}_{\beta }(t)-K_{\beta \alpha }{p}_{\alpha }(t)],$$

$$K_{\alpha \beta }=2\Re {\int }_0^{\infty }d\tau \mathrm{Tr}[e^{i{H}_{\beta }\tau }{H}_{\beta \alpha }^{\prime }{e}^{-i{H}_{\alpha }\tau }{H}_{\alpha \beta }^{\prime }{\rho }_{\beta }].$$

$$\rho (0)=\sum _{\alpha }{p}_{\alpha }(0)|\alpha ⟩⟨\alpha |{\rho }_{\alpha },$$

$$\rho (0)=\sum _{\alpha }{p}_{\alpha }(0)|\alpha ⟩⟨\alpha |{\rho }_g,$$

$$\mathcal{P}A=({\mathrm{Tr}}_{B}A){\rho }_g,H^{\prime }=H^{\text{el-ph}}$$

$$K_{nm}=2|J_{nm}{|}^2\Re {\int }_0^{\infty }d\tau e^{i{H}_m\tau }{e}^{-i{H}_n\tau }{\rho }_m,$$

$$K_{nm}=2|J_{nm}{|}^2\Re {\int }_0^{\infty }d\tau F_m^{\ast }(\tau )A_n(\tau ),$$

$$F_m(\tau )$$

$$A_n(\tau )$$

$${\rho }_n^{\mathrm{ph},\mathrm{g}}$$

Full text

I Introduction

Mathematical description of excitation energy transfer (EET) in molecular systemsMay and Kühn (2011); Valkunas, Abramavicius, and Mančal (2013); Renger, May, and Kühn (2001) is an important branch of theory of open quantum systems and chemical physics. This field of research finds applications in quantum technologies and biological systems. An important example of the latter is light-harvesting pigment-protein complexes in photosynthetic systems of plants and some bacteria.Mohseni et al. (2014) These complexes are responsible for the absorption of light and the transfer of the excitation energy into photosynthetic reaction centers. Recently, experimental evidences of long-lived quantum coherences in EET in light-harvesting complexes were obtained.Engel et al. (2007); Collini et al. (2010). This rises questions about the origin of these long-lived coherences and their possible role in EET. Note that EET in light-harvesting complexes is highly efficient: usually, more than 95% of the absorbed quanta of energy are transferred to reaction centers.Mohseni et al. (2014)

From the viewpoint of theory of open quantum systems,May and Kühn (2011); Valkunas, Abramavicius, and Mančal (2013); Breuer and Petruccione (2002); Accardi, Lu, and Volovich (2002); Accardi and Kozyrev (2002); Rivas and Huelga (2012) the electronic degrees of freedom of molecules constitute “a system”, which is coupled to “a bath” consisting of the vibrational degrees of freedom and the environment of the molecules. Quantum master equations for either the whole reduced density matrix of the system or only its diagonal part in some basis, is a widely used tool of theory of open quantum system.

There are three main approaches (or theories) leading to Markovian quantum master equations. If the dipole couplings between different molecules can be treated as small parameters, then the corresponding perturbation theory leads to Förster resonance energy transfer theory.Förster (1946, 1948) On the other hand, treating the system-bath coupling as a small perturbation parameter leads to theory of weak-coupling limit, or Redfield theory (which is often referred to as standard Redfield theory).Redfield (1965); Davies (1974); Accardi, Lu, and Volovich (2002); Accardi and Kozyrev (2002) Alternatively, treating only the off-diagonal part of the system-bath interaction Hamiltonian (in the exciton representation) leads to the modified Redfield theory.Zhang et al. (1998); Yang and Fleming (2002).

All three approaches have shown to be valid for the description of EET in molecular systems in certain regimes and provide useful insights into our understanding of these phenomena. In Ref. Seibt and Mančal, 2017 it is stressed that, despite of theoretical development of more advanced models (for example, with distinguished vibrational modesKolli et al. (2012); Plenio, Almeida, and Huelga (2013); Novoderezhkin et al. (2017)) and non-Markovian master equations (for example, the hierarchical equations of motionsIshizaki and Fleming (2009a) or the polaron-representation master equationJang et al. (2008); Kolli, Nazir, and Olaya-Castro (2011)), these three theories will continue to play an important role for the description of EET for years to come.

Unfortunately, unlike the standard Redfield theory, Förster and the modified Redfield theories describe the dynamics of only the diagonal part of the reduced density matrix of the system (populations) in the local electronic basis and exciton eigenstate basis, respectively. They do not describe the dynamics for the off-diagonal elements of the reduced density matrix (coherences). This is commonly accepted as a limitation of these theories.Yang and Fleming (2002); Ishizaki and Fleming (2009b); Novoderezhkin and van Grondelle (2010); Valkunas, Abramavicius, and Mančal (2013); Seibt and Mančal (2017) In Ref. Yang and Fleming, 2002 it is stressed that the population transfer rates are not enough for the description of optical experiments and the off-diagonal elements are needed.

Here we derive formulas for the dynamics of coherences in the framework of Förster and modified Redfield theories and, thus, remove this limitation. Our formulas provide good correspondence with the numerically exact method of hierarchical equations of motion (HEOM).Ishizaki and Fleming (2009a) Also we show that modified Redfield theory is at least as good as standard Redfield theory in predicting the exciton coherences and local site populations and, in some cases (even within the range of validity of the standard Redfield equation), provides more accurate results.

Note that, in Ref. Hwang-Fu, Chen, and Cheng, 2015, a phenomenological equation for the dynamics of coherences within modified Redfield theory is proposed. According to this equation, the coherences exponentially decay to zero. In contrast, we derive formulas for coherences using a rigorous method. Namely, we use the Zwanzig projector (super)operator method.Zwanzig (1960) Usually it is used for derivation of quantum master equations. But, in fact, it can be used for more general purposes as well. We show that the coherences do not decay to zero, but tend to some stationary non-zero values. This is confirmed by the HEOM method as well.

The paper is organized as follows. The theoretical background is reviewed in Sec. II. It includes the Hamiltonian (Sec. II.1), a general scheme of Förster and modified Redfield theories in terms of the Zwanzig projection operator formalism (Sec. II.2) and the application of this general scheme to the particular cases of Förster (Sec. II.3) and modified Redfield (Sec. II.4) theories.

The new formulas for coherences are derived in Sec. III. Namely, subsection III.1 is devoted to the case of the equilibrium initial state of the bath and no initial coherences. In subsection III.2, the case of a non-equilibrium initial state of the bath is considered. In subsection III.3, these results are generalized to the case of an arbitrary initial density matrix, i.e., with initial coherences. Every subsection in Sec. III includes formulas in the framework of the general scheme and their application to Förster and modified Redfield theories. In Sec. IV, we provide examples of calculations of coherences according to the derived formulas and compare them with the numerically exact HEOM method.

After that, in Sec. V, we discuss several issues. The developed approach allows to calculate the evolution of the whole electronic density matrix (not just its diagonal elements) for an arbitrary initial electronic density matrix (i.e., also not necessarily diagonal). In subsection V.1 we discuss the properties of the corresponding dynamical map $\Lambda_{t}$, which maps the initial electronic density matrix $\sigma(0)$ to the evolved one $\sigma(t)$. In particular, we argue that the dynamics of the whole density matrix can be treated as non-Markovian, despite the fact that the dynamics of the populations alone is Markovian. In subsection V.2 we discuss in which sense the Förster mechanism of EET can be said to be “incoherent hopping”. Finally, in subsection V.3, we give analytical rough estimates of the magnitude of electronic coherences in the Förster approach and discuss the range of validity of this approach. In particular, we show that the Förster approach may be adequate even if the intersite Coulombic coupling is much larger than the reorganization energy.

II Theoretical background

II.1 Hamiltonian

The Hamiltonian describing EET processes in molecular aggregates is as follows (the so called Frenkel exciton Hamiltonian):Yang and Fleming (2002); Mohseni et al. (2014); Seibt and Mančal (2017)

[TABLE]

where “h.c.” stands for Hermitian conjugate. Here the terms $H^{\rm el}$ and $H^{\rm Coul}$ represent the electronic (system) part: $N$ is the number of monomers (e.g., individual molecules) in the aggregate, $\ket{n}$ represents the excited electronic state of the $n$th monomer (with all other monomers being in the ground state), $\varepsilon^{0}_{n}$ is the electronic excitation energy of the $n$th monomer, $J_{nm}$ is the dipole Coulombic coupling constant between the $n$th and $m$th monomers. These coupling constants are responsible for EET between the monomers. Each monomer $n$ is coupled to its own phononic bath consisting of harmonic oscillators, with $q_{ni}$ and $p_{ni}$ being the position and momentum operators of the $i$th phonon mode of the corresponding bath. The parameters $M_{ni}$ and $\omega_{ni}$ are the mass and frequency of the corresponding mode, and $d_{ni}$ is the displacement of the equilibrium configuration of the mode between the ground and excited electronic states of the monomer. These displacements $d_{ni}$ play the role of coupling constants between the system (electronic degrees of freedom) and the bath (phononic degrees of freedom). The system-bath interaction manifests itself in the interaction Hamiltonian $H^{\text{el-ph}}$ and also in renormalization of the electronic excitation energies described by the term $H^{\rm reorg}$, with $\lambda_{n}$ being the reorganization energy of the $n$th monomer.

The dynamics of the density matrix $\rho(t)$ of both electronic and phononic degrees of freedom is given by the von Neumann equation

[TABLE]

II.2 General scheme of Förster and modified Redfield theories

Since the bath part includes an infinite number of degrees of freedom, even numerical solution of equation (2) is challenging. Perturbation theory and derivation of quantum master equations for a finite number of degrees of freedom is one of the main tools of theory of open quantum systems. The Zwanzig projection superoperator method is widely used for derivation of quantum master equations. It is based on distinguishing between “slow” and “fast” degrees of freedom. Here we describe a particular application of this method. We follow Refs. Yang and Fleming, 2002 and Seibt and Mančal, 2017. Let $\{\ket{\alpha}\}$ be an orthonormal basis of the system such that the off-diagonal part of the Hamiltonian $\braket{\alpha}{H}{\beta}$, $\alpha\neq\beta$, can be treated as a small perturbation. Then we can divide the Hamiltonian into two parts (reference and perturbation Hamiltonians):

[TABLE]

where

[TABLE]

We assume the the evolution generated by the unperturbed part $H_{0}$ can be treated exactly. In the interaction representation with respect to $H_{0}$, we have

[TABLE]

where $\rho^{(I)}(t)=e^{iH_{0}t}\rho(t)e^{-iH_{0}t}$, $\mathcal{L}(t)=[H^{\prime}(t),\cdot]$, $H^{\prime}(t)=e^{iH_{0}t}H^{\prime}e^{-iH_{0}t}.$ In the following, we will always work in the interaction representation and will omit the superscript $(I)$.

Let us introduce a projection (super)operator $\mathcal{P}$ acting on an arbitrary operator $A$ on the whole space of electronic and phononic degrees of freedom as follows:

[TABLE]

Here $\operatorname{Tr}_{B}$ is the partial trace over the bath (phonons), and $\rho_{\alpha}=e^{-\beta H_{\alpha}}/\operatorname{Tr}e^{-\beta H_{\alpha}}$ is the equilibrium state of the bath corresponding to the system state $\ket{\alpha}$ (here $\beta$ is the inverse temperature, not to be confused with the index $\beta$ in Eq. (4)). The two additional properties for such projection operator are satisfied:

[TABLE]

where $\mathcal{L}_{0}=[H_{0},\cdot]$.

Let us also denote $\mathcal{Q}=1-\mathcal{P}$. The density operator $\rho(t)$, thus, can be divided into two parts:

[TABLE]

Inserting the resolution of identity $\mathcal{P}+\mathcal{Q}$ in front of $\rho(t)$ in the right-hand side of Eq. (5) yields the system of equations

[TABLE]

where properties (7) were used.

If we treat $\mathcal{P}\rho(t)$ as a known function, then a formal solution of equation (9b) for $\mathcal{Q}\rho(t)$ is:

[TABLE]

where

[TABLE]

is the chronological exponential. If we now substitute the formal solution (10) into Eq. (9a), we will obtain the so called Nakajima–Zwanzig equation for $\mathcal{P}\rho(t)$. It is formally exact and equivalent to the von Neumann equation (5). But we are going to derive an approximate master equation for $\mathcal{P}\rho$ using the perturbation theory with respect to $H^{\prime}$.

The operator $\mathcal{P}$ projects onto the subspace of slow degrees of freedom: as we see from Eq. (9a), $\mathcal{P}\rho$ is changed slowly if $H^{\prime}$ is small. It is assumed that the free (unperturbed) dynamics generated by $H_{0}$ moves $\mathcal{Q}\rho$ to zero: $e^{-iH_{0}t}\mathcal{Q}\rho e^{iH_{0}t}\to 0$ as $t\to\infty$ in some sense. This condition (or its analogue) is required for the perturbation theory to be valid on arbitrary large times and, explicitly or implicitly, is used in derivations of master equations. But establishing a precise mathematical formulation of this condition and a rigorous proof of it is still desired. In our case, this condition means that the bath relaxes to the state $\rho_{\alpha}$ provided that the system is in the state $\ket{\alpha}$, and the off-diagonal elements $\braket{\alpha}{\rho}{\beta}$, $\alpha\neq\beta$, disappear under the free dynamics. This is a fast relaxation process.

As we see from Eq. (9b), the perturbation Hamiltonian $H^{\prime}$ slowly moves $\mathcal{Q}\rho(t)$ away from zero. The smallness of the perturbation with respect to the rate of the bath relaxation assures the smallness of $\mathcal{Q}\rho(t)$ for all times. This allows to use the perturbation theory with respect to $H^{\prime}$ and leave only the lowest-order terms in Eq. (10) for arbitrarily large times.

Often it is said that $\mathcal{P}$ and $\mathcal{Q}$ project onto the subspaces of “relevant” and “irrelevant” degrees of freedom, but this is not always true: for example, here $\mathcal{Q}\rho$ includes the coherences $\braket{\alpha}{\rho}{\beta}$, $\alpha\neq\beta$, which are often of interest as well. From the mathematical point of view it is better to say about “slow” and “fast” degrees of freedom. The separation between the time scales of the slow dynamics induced by $H^{\prime}$ and the fast relaxation dynamics induced by $H_{0}$ assures the validity of perturbation theory. We discuss this in more detail in Sec. V.3.

Let, for simplicity, $\mathcal{Q}\rho(0)=0$. From Eq. (10), in the first order with respect to $\mathcal{L}(t)$, we have

[TABLE]

The substitution of this formula into Eq. (9a) yields

[TABLE]

Equation (13) is a non-Markovian quantum master equation for $\mathcal{P}\rho$. Here the non-Markovianity (see, e.g., Refs. Breuer and Petruccione, 2002 and Rivas and Huelga, 2012 for more details) manifests itself as a dependence of the right-hand side not only on the current state $\mathcal{P}\rho(t)$ but also on the past states: $\mathcal{P}\rho(t-\tau)$.

Under some conditions, we can pass to a Markovian quantum master equation. The operator $\mathcal{P}$ includes the partial trace over the bath, which contains an infinite number of degrees of freedom. So, the integrand in Eq. (13) reflects the fast bath relaxation processes and is expected to decay to zero much faster than $\mathcal{P}\rho$ evolves. In this case we can replace $\mathcal{P}\rho(t-\tau)$ by $\mathcal{P}\rho(t)$ in the right-hand side of Eq. (13) and extend the upper limit of integration to infinity there. This is the so called Born–Markov approximation,Breuer and Petruccione (2002); Valkunas, Abramavicius, and Mančal (2013) which leads to the Markovian quantum master equation

[TABLE]

The substitution of Eq. (4b) and the expression

[TABLE]

where

[TABLE]

into Eq. (14) yields the following master equations for $p_{\alpha}(t)$:

[TABLE]

where the transfer rates $K_{\alpha\beta}$ are given by

[TABLE]

Note that the condition $\mathcal{Q}\rho(0)=0$ means that the initial state has the form

[TABLE]

i.e., there are no initial coherences (off-diagonal parts of the density matrix) and the bath is in equilibrium (adjusted to the corresponding states of the system). A generalization of this scheme to a non-equilibrium case, where the initial state has the form

[TABLE]

where $\rho_{g}=e^{-\beta H^{\rm ph}}/\operatorname{Tr}e^{-\beta H^{\rm ph}}$ is given in Ref. Seibt and Mančal, 2017. However, if the bath relaxation is much faster than the population transfer, then the non-equilibrium corrections are negligible and we can still use Eq. (18). For simplicity, we will assume this case and use Eq. (18) for population transfer rates.

As we will see in the following subsections, particular choices of the basis $\{\ket{\alpha}\}$ (and the corresponding perturbation Hamiltonians $H^{\prime}$) lead to Förster and modified Redfield theories. Already here, within the general scheme, we see that master equation (17) describes only the dynamics of populations $p_{\alpha}(t)$, or, in other words, of the diagonal part of the reduced density operator of the system $\sigma(t)=\operatorname{Tr}_{B}\rho(t)$ in the basis $\{\ket{\alpha}\}$. The off-diagonal part of it (coherences between different states) are dropped. The calculation of coherences in the framework of this scheme is the subject of the present paper.

Finally, let us note that the choice

[TABLE]

followed by an analogous derivationMay and Kühn (2011); Mohseni et al. (2014); Valkunas, Abramavicius, and Mančal (2013); Rivas and Huelga (2012); Davies (1974) leads to standard Redfield (or weak-coupling limit) theory and a Markovian quantum master equation for the whole reduced density matrix of the system $\sigma(t)$.

II.3 Förster theory

The Förster theory corresponds to the case where the intersite coupling constants $J_{nm}$ are small. We choose $\{\ket{\alpha}\}=\{\ket{n}\}$, i.e., the local electronic basis. Correspondingly, $H^{\prime}=H^{\rm Coul}$, $H_{n}=\braket{n}{H}{n}=\varepsilon_{n}^{0}+\lambda_{n}+H^{\rm ph}+u_{n}$. The application of the general formula (18) to this particular case leads to

[TABLE]

where $\rho_{m}=e^{-\beta H_{m}}/\operatorname{Tr}e^{-\beta H_{m}}$. Using the cumulant expansion method (which is approximate in general, but is exact if the reservoir consists of harmonic oscillators,Mukamel (1995) which is our case), we arrive at:Yang and Fleming (2002)

[TABLE]

[TABLE]

[TABLE]

where

[TABLE]

and

[TABLE]

is the lineshape function,Mukamel (1995)

[TABLE]

is the correlation function of the $n$th bath, $u_{n}(t)=e^{iH^{\rm ph}_{n}t}u_{n}e^{-iH^{\rm ph}_{n}t}$.

II.4 Modified Redfield theory

Denote $E_{k}$ and $\ket{k}=\sum_{n=1}^{N}\varphi_{kn}\ket{n}$ the eigenvalues and the corresponding eigenstates of the system Hamiltonian $H^{\rm el}+H^{\rm reorg}+H^{\rm Coul}$ (the spectrum is assumed to be non-degenerate). Taking $\{\ket{\alpha}\}=\{\ket{k}\}$, i.e., the exciton basis, and

[TABLE]

where

[TABLE]

$a_{kk^{\prime}}(n)=\overline{\varphi_{kn}}\varphi_{k^{\prime}n}$, corresponds to modified Redfield theory. In other words, this corresponds to the case when the coupling between different delocalized excitations can be treated as small. The application of the general formula (18) and the cumulant expansion method to this particular case leads toYang and Fleming (2002)

[TABLE]

where $\rho_{k}=e^{-\beta H_{k}}/\operatorname{Tr}e^{-\beta H_{k}}$, $H_{k}=\braket{k}{H}{k}$,

[TABLE]

[TABLE]

[TABLE]

and $E_{k}^{0}=E_{k}-\lambda_{k}$.

III Dynamics of coherences

III.1 The case of equilibrium initial state of the bath

Quantum master equation allows to describe the dynamics of the “slow” degrees of freedom $\mathcal{P}\rho(t)$. In the case of the Förster and modified Redfield approaches, this corresponds to the populations $p_{n}(t)$ and $p_{k}(t)$, respectively. However, the Zwanzig projection operator formalism, namely, formula (10) allows also to describe the dynamics of the “fast” degrees of freedom as well. A problem is that $\mathcal{Q}\rho(t)$ is infinite-dimensional, so, expression (10) cannot be calculated directly. However, we may be interesting in a particular observable $A$. In this paper we are interested in coherences, i.e., the off-diagonal elements of the density matrix, and we take

[TABLE]

Denote $A(t)=e^{iH_{0}t}Ae^{-iH_{0}t}$, then

[TABLE]

where we have denoted $\sigma(t)=\operatorname{Tr}_{B}\rho(t)$ the reduced density operator of the system. In this subsection we consider the case $\mathcal{Q}\rho(0)=0$, i.e., the bath is initially in equilibrium (see Eq. (19)). We have

[TABLE]

where we have used Eq. (12) and $\operatorname{Tr}A\mathcal{P}\rho=\braket{\beta}{\mathcal{P}\rho}{\alpha}=0$. For other choices of $A$ (different from Eq. (38)), the term $\operatorname{Tr}A\mathcal{P}\rho$ may be non-zero, but it can be treated as known since $\mathcal{P}\rho(t)$ is given by the solution of quantum master equation (13) or (14).

Here we assume that the Born–Markov approximation, which allows to turn Eq. (13) into Eq. (14), is valid. Then, again, the time dependence of the term $A(t)\mathcal{L}(t-\tau)$ in Eq. (41) represents fast bath relaxation processes, while $\mathcal{P}\rho(t)$ is evolved more slowly. So, we can replace $\mathcal{P}\rho(t-\tau)$ by $\mathcal{P}\rho(t)$:

[TABLE]

Since $p_{\alpha}(t)$ and $p_{\beta}(t)$ are known from the solution of the master equation, formula (42) allows to calculate the dynamics of the coherences.

Remark 1*.*

Note that, in contrast to Eq. (14), here we do not extend the upper limit of integration to infinity. This approximation works well for large $t$, but introduce a significant error for small $t$ and, in some cases, moderate $t$. Indeed, from Eq. (42) we see that $\sigma_{\beta\alpha}(0)=0$, which is in agreement with the initial condition $\mathcal{Q}\rho(0)=0$. However, the extension of the upper limit of integration to infinity would give $\sigma_{\beta\alpha}(0)\neq 0$. In contrast, the extension of the upper limit of integration to infinity in Eq. (13) produce the error in the time derivative of $\mathcal{P}\rho(t)$, but not in $\mathcal{P}\rho(t)$ directly. Since $\mathcal{P}\rho(t)$ is changed slowly, the error in $\mathcal{P}\dot{\rho}(t)$ for small $t$ does not produce a significant error in $\mathcal{P}\rho(t)$. Thus, on small times, the coherences evolve faster than the populations. Also one can notice that the discussed error has the second order of smallness (with respect to $H^{\prime}$) for $\mathcal{P}\dot{\rho}(t)$, but only the first order of smallness for $\sigma_{\beta\alpha}(t)$.

For large times $t$, the integrals in Eq. (42) saturate to some constant values, so, the time integration can be extended to infinity, and the dynamics of the coherences are driven by the populations with the constant coefficients:

[TABLE]

where

[TABLE]

However, as we will observe in Sec. V.1, the coherences may oscillate around the mean values (43) for a long time comparable with the EET time scale. To describe these oscillations, we shall use more general formula (42).

From Eq. (43) we see that the coherences in Förster and modified Redfield theories do not decay to zero (like in the aforementioned phenomenological approach of Ref. Hwang-Fu, Chen, and Cheng, 2015), but tend to some constant values dependent on the stationary populations. The decay of the coherences under the free (unperturbed) dynamics generated by the Hamiltonian $H_{0}$ is compensated by their “pumping” from the populations induced by $H^{\prime}$.

Now we apply the general formula (42) to Förster and modified Redfield theories using the cumulant expansion technique. In the Förster case $A=\ket{n}\bra{m}$, we obtain

[TABLE]

In the modified Redfield case $A=\ket{k}\bra{k^{\prime}}$, we obtain

[TABLE]

III.2 The case of non-equilibrium initial state of the bath

Now we extend the results of the previous subsection to the case of a non-equilibrium initial state of the bath. Let the initial state of the molecular system be given by Eq. (20), i.e., there is still no initial coherences, but the bath is in a non-equilibrium state. The bath state $\rho_{g}$ describes the bath equilibrium state when no electronic excitation is present. Right after the excitation, the bath is still in the state $\rho_{g}$ (the Condon approximation), which is non-equilibrium for an excited electronic state. Then

[TABLE]

and the first-order approximation of Eq. (10) (with respect to $\mathcal{L}(t)$) is

[TABLE]

The substitution of Eq. (48) into Eq. (III.1) gives

[TABLE]

The first term has been already calculated in Eq. (42). Consider the second term:

[TABLE]

where $\rho_{g}^{\alpha}(t)=e^{-iH_{\alpha}t}\rho_{g}e^{iH_{\alpha}t}$. So, $\sigma_{\beta\alpha}(t)$ has the form

[TABLE]

where $\sigma_{\beta\alpha}^{\rm eq}(t)$ is the expression for the equilibrium case given by Eq. (42) and $\sigma_{\beta\alpha}^{\rm noneq}(t)$ is a non-equilibrium correction:

[TABLE]

The non-equilibrium terms are significantly non-zero only on the time scale of the bath relaxation. Since the populations evolve more slowly, we can replace $p_{\alpha}(0)$ and $p_{\beta}(0)$ by $p_{\alpha}(t)$ and $p_{\beta}(t)$ in Eq. (52) and obtain

[TABLE]

However, in the calculations below, this additional approximation introduces small but noticeable error and, so, we will use formulas (51)–(52).

Remark 2*.*

The substitution of Eq. (48) into Eq. (9a) yields Eq. (14) with a non-equilibrium correction, which is calculated in Ref. Seibt and Mančal, 2017. Here, for simplicity, we assume that the bath relaxes much faster than the population transfer occurs and this correction is negligible. Then $p_{\alpha}(t)$ and $p_{\beta}(t)$ are well described by master equation (17) with the equilibrium transfer rates (18).

But even in this case, the non-equilibrium corrections for coherences (52) may be essential because, on small times, coherences evolve faster than the populations, see Remark 1. But in general, solutions of the master equation with the non-equilibrium (time-dependent) transfer rates from Ref. Seibt and Mančal, 2017 also can be used for $p_{\alpha}(t)$ and $p_{\beta}(t)$ in Eqs. (42), (51), and (53). Also note that the non-equilibrium generalization of the Förster transfer rates was firstly derived in Ref. Jang, Jung, and Silbey, 2002, but a compact expression involving the lineshape functions $g_{n}(t)$ was obtained in Ref. Seibt and Mančal, 2017.

The traces involving $\rho_{g}$ in Eqs. (52) and (53) for Förster and modified Redfield theories appear also during the derivation of the non-equilibrium transfer rates for the populations and were also calculated in Ref. Seibt and Mančal, 2017. For Förster theory, Eq. (52) takes the form

[TABLE]

where

[TABLE]

Note that $F_{m}(\tau,t)$ can be replaced by $F_{m}(\tau)$ in Eq. (54) for large $t$ because $\lim_{\tau\to\infty}\dot{g}_{m}(\tau)=-\lambda_{m}$Zhang et al. (1998) and the integral expressions are significantly non-zero only for small $\tau$. So, the non-equilibrium correction vanishes for large times.

For modified Redfield theory, Eq. (52) takes the form

[TABLE]

where

[TABLE]

[TABLE]

[TABLE]

Note that $(\Im\dot{g})_{kk^{\prime}k^{\prime}k^{\prime}}(t)=\Im\dot{g}_{kk^{\prime}k^{\prime}k^{\prime}}(t)$ whenever $a_{kk^{\prime}}(n)$ is real for all $n$.

III.3 The case of initial coherences in the system

Finally, we consider the case when the initial state has non-zero off-diagonal elements. We will consider the initial state of the form

[TABLE]

where $\sigma(0)$ is an arbitrary density operator of the electronic degrees of freedom. We divide $\mathcal{Q}\rho(0)$ into the diagonal and off-diagonal parts:

[TABLE]

In both Eqs. (9a) and (III.1) only the off-diagonal part of $\mathcal{Q}\rho(t)$ matters. We have

[TABLE]

Note that the influence of $\mathcal{Q}\rho(0)^{\rm diag}$ (of the last term in Eq. (62)) has been already taken into account in the previous subsection. Here we will study the influence of only the off-diagonal part $\mathcal{Q}\rho(0)^{\text{off-diag}}$.

Firstly, the initial coherence infer the population transfer on small times. As we discussed above, the substitution of the last term of Eq. (62) into Eq. (9a) produces the non-equilibrium correction to the population transfer rates which is assumed to be negligible. The substitution of the first two terms of the right-hand side of Eq. (62) into Eq. (9a) yields Eq. (14) with the following additional term in the right-hand side:

[TABLE]

where

[TABLE]

Consequently, the master equation (17) is modified to

[TABLE]

Now let us consider the influence of initial coherences on their further dynamics. The substitution of Eq. (62) into Eq. (III.1) gives $\sigma_{\beta\alpha}(t)$ as a sum of three terms:

[TABLE]

where the basic term $\sigma_{\beta\alpha}^{\rm eq}(t)$ is given by Eq. (42), the non-equilibrium correction $\sigma_{\beta\alpha}^{\rm noneq}(t)$ is given by Eq. (52), and the new term (a correction caused by the initial coherence) $\sigma_{\beta\alpha}^{\rm init}(t)$ is given by

[TABLE]

Now let us apply Eqs. (65) and (67) to the particular theories. In Förster theory, the traces in Eqs. (65) and (67) can be easily calculated, which results in

[TABLE]

and

[TABLE]

In the modified Redfield theory, the trace in Eq. (65) is a particular case of the traces calculated in Ref. Seibt and Mančal, 2017, while the trace in Eq. (67) can be calculated by a slight modification of the derivations there. It turns out that

[TABLE]

where

[TABLE]

where

[TABLE]

Remark 3*.*

The following question about self-consistency of our calculations may be asked. Let $\mathcal{Q}\rho(0)=0$. At some instant $t_{0}>0$, coherences are non-zero and given by Eq. (42). If we apply the calculations of this subsection for the case of initial coherences, do we reproduce formulas (17), (18), and (42) for an arbitrary time $t>t_{0}$? An important point we should notice is that $\rho(t_{0})$ differs from the state given by Eq. (60). Namely,

[TABLE]

So, the calculations of this subsection should be modified for $\mathcal{Q}\rho(t_{0})$ given by Eqs. (73). If we do this, then, since $\mathcal{Q}\rho(t)$ can be expressed in two ways as

[TABLE]

we will reproduce formulas (17), (18), and (42). Informally speaking, the coherences and the state of the bath in Eqs. (73) are not arbitrary, but adjusted with each other in a special way.

IV Calculations

In this section we will compare the calculations of the coherences according to the derived formulas with the calculations according to the numerically exact HEOM (hierarchical equations of motion) method.Ishizaki and Fleming (2009a) The lineshape functions $g_{n}(t)$ can be expressed asMukamel (1995)

[TABLE]

where $\mathcal{J}_{n}(\omega)$ is the spectral density function for the bath coupled to the $n$th site defined as

[TABLE]

It can be expressed in terms of the correlation function $C_{n}(t)$ (see Eq. (28)) as

[TABLE]

We will consider the dimer case $N=2$ and assume that the spectral density function for both molecules has the Drude–Lorentz form:

[TABLE]

where $\lambda$ and $\gamma$ are the reorganization energy and Debye frequency, which are identical for both molecules. The reorganization energy characterizes the strength of the coupling of electronic and phononic degrees of freedom, and the Debye frequency characterizes a time scale of fluctuation of the electronic energy and dissipation of the phonon reorganization energy.

Also we impose a high-temperature condition $\beta\gamma<1$, where $\beta=1/k_{\rm B}T$ is the inverse temperature of the baths, $T$ is the temperature and $k_{\rm B}$ is the Boltzmann’s constant. Then we can approximate $\coth(\beta\omega/2)$ by $2/(\beta\omega)$ in Eq. (75) and obtain the expression for the lineshape function (the same for both monomers):

[TABLE]

In all examples we put $T=300~{}\rm{K}$ and $\gamma^{-1}=100~{}\rm{fs}$ ($\gamma=53.08~{}\rm{cm}^{-1}$). Then $\beta\gamma\approx 0.24$. Also in all examples, $\varepsilon^{0}_{1}-\varepsilon^{0}_{2}=100~{}\rm{cm}^{-1}$. Since $J_{12}$ is the only intersite dipole coupling, we will denote it as $J$. We will vary the parameters $J$ and $\lambda$.

To distinguish between the matrix elements of the electronic density operator in the local basis $\sigma_{mn}$ and in the exciton basis $\sigma_{k^{\prime}k}$ when the indexes are concrete numbers, we will use the following denotations: the matrix elements in the local electronic basis will be denoted as $\sigma_{mn}$, while those in the exciton basis will be denoted as $\sigma^{\rm ex}_{k^{\prime}k}$. The same rule will be applied to the Dirac notations and the bath equilibrium states: $\ket{k}$ and $\rho_{k}$ will be redenoted as $\ket{k}_{\rm ex}$ and $\rho_{k}^{\rm ex}$.

IV.1 Förster theory

Let us consider the parameters $J=10~{}\rm{cm}^{-1}$ and $\lambda=100~{}\rm{cm}^{-1}$, and the initial state

[TABLE]

(a particular case of Eq. (20)). The results of the calculations of the electronic coherence $\sigma_{21}(t)$ within Förster theory according to the derived formula (45) with and without the non-equilibrium correction (54) in comparison with HEOM is presented on Fig. 1. We see an excellent agreement of the calculations according to Eqs. (45) and (54) with the numerically exact result. As for Eq. (45) without the non-equilibrium correction, it gives the same results for large times, but fails to describe the coherence (the real part of it) on small times.

The calculations according to Eq. (45) without the non-equilibrium correction can be compared with the calculations according to HEOM with the equilibrium initial state

[TABLE]

(a particular case of the state (19)) is given on Fig. 2. We see an excellent agreement.

The present choice of parameters corresponds to the case of slow nuclear motion.May and Kühn (2011); Mukamel (1995) Let us choose the parameters corresponding to fast nuclear motion: $J=10~{}\rm{cm}^{-1}$, $\lambda=1~{}\rm{cm}^{-1}$. Then the non-equilibrium corrections are negligible as shown on Fig. 3. Also we see that Förster theory provides good results, despite the regime $\lambda\ll J$, while Förster theory is believed to be valid only in the inverse case $J\ll\lambda$. Förster theory gives here good results also for the populations $\sigma_{11}(t)$ and $\sigma_{22}(t)$: the error does not exceed 3.5% for all times. We discuss the range of validity of Förster theory in Sec. V.3.

Let us return to the case $J=10~{}\rm{cm}^{-1}$ and $\lambda=100~{}\rm{cm}^{-1}$, when the non-equilibrium correction is essential, and consider the case of initial electronic coherences. Namely, let us consider the case when the initial electronic state is not a local (electronic) excitation, but an exciton eigenstate:

[TABLE]

Using the derived formulas, we can calculate the dynamics of the whole density matrix (not only its diagonal part in the local basis) in the framework of Förster theory. Let us calculate, in the framework of Förster theory, the exciton population $\sigma_{11}^{\rm ex}(t)$ as well as the exciton coherence $\sigma_{21}^{\rm ex}(t)$ and compare the results with the calculations within the modified Redfield theory, which is traditionally used for the description of exciton population transfer. An exciton population is a linear combination of the local site populations and intersite (electronic) coherences. Since we use the non-equilibrium corrections for the electronic coherences, we use also non-equilibrium, time-dependent generalizations for population transfer rates $K_{\alpha\beta}$ in Eq. (17) for both modified Redfield and Förster theories,Seibt and Mančal (2017) for the honest comparison between the two theories. This is the only place in our work where we have used the non-equilibrium, time-dependent population transfer rates.

The results of such comparison are presented on Fig. 4 (for the exciton populations) and Fig. 5 (for the exciton coherence). We see that both Förster and modified Redfield theories are in good agreement with the numerically exact HEOM calculations, but Förster theory is, nevertheless, in better agreement.

IV.2 Modified Redfield theory

In this subsection we test the derived formulas for exciton coherences within modified Redfield theory. Examples for the cases $\lambda\ll J\ll\gamma$ and $J\ll\lambda\ll\gamma$ was already considered on Fig. 3 and Figs. 4 and 5. In both cases, Förster theory provides good results and the modified Redfield theory can be compared with it. Consider now the case $\lambda\ll\gamma\ll J$, which is certainly beyond the range of validity of the Förster theory, but is in the range of validity of standard and modified Redfield theories (see Sec. V.3): $J=100~{}\rm{cm}^{-1}$ and $\lambda=2~{}\rm{cm}^{-1}$. Consider the initial state as an exciton eigenstate (82). The results of the application of Eq. (46) (without the non-equilibrium correction) for the exciton coherence $\sigma_{21}^{\rm ex}(t)$ in comparison with the calculations according to the standard non-secular Redfield master equation and HEOM are presented on Fig. 6. We see that Eq. (46) in the framework of the modified Redfield theory is in significantly better agreement with the numerically exact calculations than the standard Redfield master equation. Note that the secular Redfield master equation predicts zero exciton coherence whenever the initial exciton coherence is zero.

Using the derived formulas, we can calculate the whole density matrix (not only its diagonal part in the exciton basis) in the framework of the modified Redfield theory. Let us calculate, in the framework of the modified Redfield theory, not the exciton populations, but the local site population $\sigma_{11}(t)$, for the same initial state and parameters. The results are presented on Fig. 7. Again, we see that the calculations according to the modified Redfield theory gives significantly better results than both secular and non-secular standard Redfield master equations.

Let us now consider the initial state as a local electronic excitation (80), which has non-zero exciton coherences. The calculation of the coherence $\sigma_{21}^{\rm ex}$ is presented in Fig. 8. Now all three approaches: HEOM, the modified Redfield approach with formula (46), and the secular and non-secular standard Redfield master equations give approximately the same results.

V Discussion

V.1 Properties of the dynamical map

The developed approach allows to calculate the evolution of the whole electronic density matrix (not just its diagonal elements) for an arbitrary initial electronic density matrix (i.e., also not necessarily diagonal). If we consider the initial system-bath state of the form (60), then the evolution of the electronic density matrix according to the master equation (65) and formula (66) can be described as an action of a dynamical mapBreuer and Petruccione (2002); Rivas and Huelga (2012) $\Lambda_{t}$:

[TABLE]

This map is trace-preserving because the master equation for populations (65) is obviously trace-preserving. It does not preserve positivity when $H^{\prime}_{\beta\alpha}$ are large, i.e., when we are not in the range of validity of the corresponding perturbation theory. It is an open question whether $\Lambda_{t}$ preserves positivity for small enough $H^{\prime}_{\beta\alpha}$.

A further question is whether the dynamical map $\Lambda_{t}$ is Markovian or not. If we use the semigroup property $\Lambda_{t+s}(\sigma)=\Lambda_{s}(\Lambda_{t}(\sigma))$ as the definition of Markovianity,Breuer and Petruccione (2002); Rivas and Huelga (2012) then the dynamics is non-Markovian. From Eqs. (42), (52), (67), and (65), we see that the dynamics is non-Markovian, bath-dependent.

From the other side, consider large times, when the integrals in Eq. (42) saturate and corrections (52) and (65) caused by initial conditions become negligible. Then the population transfer is described by the Markovian master equation (17) and the coherences are completely defined by the populations according to Eq. (43). This dynamics can be treated as Markovian, in a sense that the knowledge of the present state of the system (without the knowledge of the time instant $t$) is sufficient to predict further dynamics.

For example, as we see on Fig. 2, the “Markovian” dynamics begins after very short initial time period, which is much shorter than the population transfer time scale. Since

[TABLE]

the population transfer time scale is equal to $(K_{12}+K_{21})^{-1}\approx 3800~{}\rm{fs}$ for the parameters on Fig. 2.

However, in general, the time when the integrals in Eq. (42) saturate and, thus, Eq. (42) can be replaced by Eq. (43) may be comparable with the population relaxation time. On Fig. 9, we compare the calculation of the exciton coherence within modified Redfield theory according to Eq. (46) and according to Eq. (43) for the parameters and initial state as on Fig. 6. We see that the decay of oscillations takes place on a time scale which is comparable with the population transfer time scale (the latter is approximately 3500 fs). Thus, one can speak about non-Markovian dynamics of the whole density matrix, despite the fact that the master equation for populations (17) alone is Markovian.

Note also that, for the considered parameters, these oscillations have little influence on the population transfer because, as we noticed above, $\mathcal{Q}\rho(t)$ influence not directly $\mathcal{P}\rho(t)$, but the time derivative $\mathcal{P}\dot{\rho}(t)$, see Eq. (9a). If the frequency of oscillations of the right-hand side of Eq. (9a) is much larger than the population transfer rate, then these oscillations can be neglected (a kind of secular approximation for populations). This assures the Markovian nature of the population transfer.

V.2 On relation between populations and coherences in Förster and modified Redfield EET mechanisms

The Förster mechanism of EET is often said to be incoherent hopping.Jang et al. (2008); May and Kühn (2011); Ishizaki and Fleming (2009a) The presented results allow to discuss what does this statement exactly mean. As we have shown, this means neither that there are no electronic coherences in the system nor that the influence of coherences on the populations transfer is negligible. Eq. (9a) means precisely that the population transfer is directly driven by the coherences: by the definitions of $H^{\prime}$ and $\mathcal{P}$, only the off-diagonal part of $\mathcal{Q}\rho(t)$ matters in the right-hand side of Eq. (9a).

As we have shown, “incoherent hopping” means that either, after a short initial time (relative to the EET time scale), coherences can be described by Eq. (43) and, thus, are completely determined by the populations, or coherences rapidly oscillate around the values given by Eq. (43). The influence of these fast oscillations on the population dynamics can be neglected and, so, in both cases we obtain a closed Markovian master equation for populations.

In other words, “incoherent hopping” actually means that the feedback loop from the populations to coherences and back to the time derivatives of the populations has no time delay. The same is true for the modified Redfield theory concerning the exciton populations and coherences.

V.3 Analytical estimates of the magnitude of coherences and the range of validity of Förster theory

It is believed that the range of validity of Förster theory is $J\ll\lambda$: the Coulomb intersite couplings are less than the system-bath couplings (expressed in the reorganization energy $\lambda$).May and Kühn (2011); Kolli, Nazir, and Olaya-Castro (2011); Seibt and Mančal (2017) But we see on Fig. 3 that Förster theory may work well even in the inverse case $\lambda\ll J$. To understand this, we discuss here the range of validity of Förster theory. A rigorous analysis of the approximations made in Sec. II.2 (for both theories) should be a subject for a separate research. Here we restrict ourselves to some mathematical intuitions. Also, during the analysis, we will obtain rough analytical estimates of electronic coherences in the Förster approach: formulas (86) below for the regime of slow nuclear motion and (87) for the regime of slow nuclear motion.

To arrive at master equation (17) with rates (18) from the exact equations (9), we have done the following approximations: (i) second-order perturbation theory in (10), (ii) replacing $\mathcal{P}\rho(t-\tau)$ by $\mathcal{P}\rho(t)$ in Eq. (13), and (iii) extending the upper limit of integration to infinity in Eq. (13). To derive the formulas for the coherences, we have used approximations (i) and (ii). Let us discuss each of these approximations.

As we already discussed in Sec. II.2, we can leave only the lowest-order terms in Eq. (10) whenever $\mathcal{Q}\rho(t)\ll 1$ for all times. In particular, this is satisfied if the perturbation Hamiltonian $H^{\prime}$ (in Förster theory it is proportional to $J$) is much smaller than the the relaxation rate of $\mathcal{Q}\rho$ induced by $H_{0}$.

So, to establish the range of validity of Förster theory, we should discuss the relaxation of $\mathcal{Q}\rho$. We can observe that the relaxation of $\mathcal{Q}\rho^{\rm diag}$ and $\mathcal{Q}\rho^{\text{off-diag}}$ (see Eqs. (61)) takes place with different rates. For the Drude–Lorentz spectral density (78), the relaxation rate of $\mathcal{Q}\rho^{\rm diag}$ can be associated with the Debye frequency $\gamma$. This can be seen, for example, from the comparison of the non-equilibrium emission function (55) with the equilibrium one (24) and the expression for the lineshape function (79): $\Im[g_{m}(t)-g_{m}(t-\tau)]-(-\lambda_{m})\sim e^{-\gamma t}$. So, the condition $J\ll\gamma$ is necessary for the validity of Förster theory.

From Eqs. (43), (45), (24), and (25), we can roughly estimate the magnitude of $\mathcal{Q}\rho^{\text{off-diag}}$ as

[TABLE]

where $\Delta\varepsilon^{0}$ is the characteristic difference between the electronic excitation energies of two monomers. Since we are interested here only in rough estimates, we do not consider the dependence of the quantities on the individual monomers. From Eq. (85), we see that $\mathcal{Q}\rho^{\text{off-diag}}$ is small whenever $\lambda$ is large (since $g(\tau)$ is proportional to $\lambda$) or $\Delta\varepsilon^{0}-2\lambda$ is large (since fast oscillations of the integrand decrease the value of the integral).

Consider this in more detail: namely, consider two limiting cases of fast and slow nuclear motionMay and Kühn (2011); Mukamel (1995) depending on the dimensionless parameter $\kappa=\sqrt{\beta\gamma^{2}/2\lambda}$.

The case $\kappa\ll 1$ corresponds to slow nuclear motion. In this limit we can approximate $e^{-\gamma t}\approx 1-\gamma t+(\gamma t)^{2}/2$ in the expression (79) for $g(t)$, neglect the imaginary part of $g(\tau)$ and obtain $g(\tau)\approx\lambda\tau^{2}/\beta$. Then

[TABLE]

where $F(x)=e^{-x^{2}}\int_{0}^{x}e^{\tau^{2}}d\tau$ is the Dawson’s integral.Temme (2010) For large $x$, $F(x)$ decays to zero as $(2x)^{-1}$. Since we already have the condition $J\ll\gamma$ and are considering the case $\kappa\ll 1$, we have $J\sqrt{\beta/2\lambda}\ll 1$, which means that $\mathcal{Q}\rho^{\text{off-diag}}$ is small.

Consider the inverse limiting case $\kappa\gg 1$ (fast nuclear motion). In this case we neglect the terms $e^{-\gamma t}$ and $-1$ in $g(t)$ and obtain

[TABLE]

Hence $\mathcal{Q}\rho^{\text{off-diag}}$ is small whenever $J$ is much smaller than the maximum of $4\lambda/\beta\gamma$ and $\Delta\varepsilon^{0}$. For example, Fig. 3 corresponds to the case $J\gg\lambda$, but $J\ll\gamma$ and $J\ll\Delta\varepsilon^{0}$, which makes Förster theory adequate.

Consider now approximation (ii): the replacement of $\mathcal{P}\rho(t-\tau)$ by $\mathcal{P}\rho(t)$ in Eq. (13). $\mathcal{P}\rho(t)$ evolves with the rate $K_{nm}+K_{mn}$ (see Eq. (84)). Consider again the case of slow nuclear motion first. As we see from Eq. (86), the characteristic time of decay of the integrand in Eq. (13) is $\sqrt{\beta/4\lambda}$. As we concluded above, $J\ll\sqrt{2\lambda/\beta}$, so,

[TABLE]

Thus, the time scale of population transfer $(K_{nm}+K_{mn})^{-1}$ is much larger than the time scale of $\mathcal{Q}\rho^{\text{\rm off-diag}}$ decay, which justifies approximation (ii) for the case of slow nuclear motion.

Consider the case of fast nuclear motion. In this case,

[TABLE]

As we concluded above, $J^{2}$ must be much smaller than the denominator, hence, the population transfer rate $K_{nm}+K_{mn}$ is much smaller that the decay rate $4\lambda/\beta\gamma$ of the integrand in Eq. (13) (Eq. (87)), which justifies approximation (ii) for the case of fast nuclear motion.

Approximation (iii) is based on the same assumption as approximation (ii). But in some cases the extension of the upper limit of integration in Eq. (13) to infinity produces a noticeable error on small times. In this case we can leave the upper limit of integration equal to $t$.

Note that in our analysis we used simplified formula (43) for coherences because, as we discussed in the end of subsection V.1, the oscillations around the mean values given by the simplified formula does not affect much the population transfer.

Summarizing, the range of validity of the Förster approach is: $J\ll\gamma$ and, in the case of fast nuclear motion, $J\ll\max\{4\lambda/\beta\gamma,\Delta\varepsilon^{0}\}$.

A similar analysis for modified Redfield theory is slightly more complicated due to pre-exponential factors in Eqs. (31) and (46) containing time derivatives of the lineshape functions, but also can be performed. In particular, $H^{\prime}\ll\gamma$ is a necessary condition for the validity of this approach. The main (well-known) characteristic feature of modified Redfield theory is that $H^{\prime}\sim\lambda$ if the excitons are highly delocalized, but $H^{\prime}\sim\lambda J/\Delta\varepsilon$ for localized excitons (for $J\ll\Delta\varepsilon$). Hence, the range of validity of modified Redfield theory includes the range of validity of standard Redfield theory (small $\lambda$) and is also adequate when large $\lambda$ is compensated by large exciton delocaization due to static disorder $\Delta\varepsilon^{0}$. A serious limitation of the modified Redfield approach is the case of not static, but dynamic localization of excitons due to strong interaction with phonons (polaron effect).Novoderezhkin and van Grondelle (2010, 2017)

Also the range of validity of modified Redfield theory intersects with the range of validity of Förster theory. As we can see from Fig. 3 and Figs. 4 and 5, the modified Redfield approach provides slightly more accurate results than Förster theory for the regime $\lambda\ll J\ll\gamma$ (even if Förster theory is also adequate), but Förster theory provides slightly better results for the regime $J\ll\lambda$. Also, from Figs. 6–8 we can conclude that modified Redfield theory is at least as good as standard Redfield theory in the description of exciton coherences and, in some cases (even within the range of validity of the standard Redfield equation), provides more accurate results.

VI Conclusions

The result of the present paper is that we are not restricted by calculations of only the diagonal elements (in either local or exciton basis) of the density operator when we use Förster or modified Redfield theory. The developed approach allows to calculate the evolution of the whole electronic density matrix for an arbitrary initial electronic density matrix (i.e., also not necessarily diagonal).

The general formula for coherence is Eq. (66) with the particular terms given in Eqs. (42), (52), and (67). The applications of these formulas to Förster theory are Eqs. (45), (54), and (69). The applications of these formulas to modified Redfield theory are Eqs. (46), (56), and (71). The modifications of the master equations for populations in the case of initial coherences are given in Eqs. (65) (the general formula), (68) (for Förster theory), and (70) (for modified Redfield theory). These formulas are rigorously derived using the Zwanzig projection operator method and provide good agreement with the numerically exact calculations according to HEOM.

Moreover, we have shown that the derived formulas for exciton coherences and local site populations within the modified Redfield theory provide at least as accurate results as the standard Redfield equation and, in some cases (even within the range of validity of the standard Redfield equation), provides more accurate results (see Figs. 6 and 7). Hence, the modified Redfield approach can be used for the investigation of the role of coherences in EET.

Extending the presented formalism to more complicated systems, for example, to EET between weakly connected clusters of molecules (which is often described by the multichromophoric generalization of Förster theoryNovoderezhkin and Razjivin (1996); Sumi (1999); Mukai, Abe, and Sumi (1999); Scholes and Fleming (2000); Jang, Newton, and Silbey (2004) and combined modified-Redfield–Förster approachNovoderezhkin and van Grondelle (2010); Jang et al. (2003); Novoderezhkin and van Grondelle (2017)), systems with degenerate spectrum,Aref’eva, Volovich, and Kozyrev (2015); Volovich and Kozyrev (2016); Kozyrev (2018) systems with common vibrational modes of different monomers,May and Kühn (2011); Jang, Jung, and Silbey (2002); Kolli et al. (2012); Novoderezhkin et al. (2017) etc. remains an open problem.

It would be interesting to apply the developed formalism to biological systems and, in particular, to coherent light harvesting in photosynthetic complexes. In Ref. Kolli et al., 2012, the role of quantized vibrations in EET processes in light-harvesting complexes by cryptophyte algae was studied using the polaron transformation with the subsequent weak system-bath coupling approach. Note that this approach also faces the the same problem with the description of electronic coherences: we need $\mathcal{Q}\rho(t)$ degrees of freedom for them, and essentially the same method (based on formulas (III.1) and (10)) was used. The authors of Ref. Kolli et al., 2012 point out a characteristic property of the considered light-harvesting complexes: highly localized exciton eigenstates (due to static disorder). As is well known and discussed in Sec. V.3, this case fall into the range of validity for the modified Redfield approach,Yang and Fleming (2002); Seibt and Mančal (2017); Novoderezhkin and van Grondelle (2010) so, probably, this approach also can be applied for studying these systems. As we discussed in Sec. V.1, the description of dynamics of the whole density matrix within the modified Redfield approach also (like the polaron-representation Markovian master equation) captures some non-Markovian effects which are believed to be noticeable in light-harvesting complexes.

Acknowledgements.

I am grateful to D. Abramavičius, T. Mančal, S. Rajabi, T. Renger, J. Seibt, and A. E. Teretenkov for fruitful discussions and comments. This work was supported by the Russian Science Foundation (project 17-71-20154).