Higher-order corrections to the Redfield equation with respect to the system-bath coupling based on the hierarchical equations of motion

TL;DR

This paper develops a general method to derive higher-order corrections to the Redfield quantum master equation using hierarchical equations of motion, improving its accuracy for stronger system-bath couplings.

Contribution

The authors introduce a novel approach based on HEOM to systematically obtain higher-order corrections to the Redfield equation, extending its applicability.

Findings

Derived conditions for the validity of the Redfield equation.

Established criteria for the secular approximation in open quantum systems.

Provided a framework for more accurate modeling of system-bath interactions.

Abstract

The Redfield equation describes the dynamics of a quantum system weakly coupled to one or more reservoirs and is widely used in theory of open quantum system. However, the assumption of weak system-reservoir coupling is often not fully adequate and higher-order corrections to the Redfield equation with respect to the system-bath coupling is required. Here we propose a general method of derivation of higher-order corrections to the Redfield quantum master equation based on the hierarchical equations of motion (HEOM). Also we derive conditions of validity of the Redfield equation as well as the additional secular approximation for it.

Equations

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

$$H^{\mathrm{el}}=\sum _{j=1}^N|j⟩{\epsilon }_j⟨j|+\sum _{j=1}^N\sum _{k>j}^N(J_{jk}|j⟩⟨k|+\text{h.c.}),$$

$$H^{\mathrm{ph}}=\sum _{j=1}^N{H}_j^{\mathrm{ph}},H_j^{\mathrm{ph}}=\sum _i\left(\frac{p_{ji}^2}{2M_{ji}}+\frac{1}{2}{M}_{ji}{\omega }_{ji}^2{q}_{ji}^2\right),$$

$$H^{\text{el-ph}}=\sum _{j=1}^N{V}_j\otimes u_j,V_j=|j⟩⟨j|,u_j=\sum _i{M}_{ji}{\omega }_{ji}^2{d}_{ji}{q}_{ji},$$

$${\dot{\rho }}^{\mathrm{tot}}(t)=-i[H^{\text{el-ph}}(t),{\rho }^{\mathrm{tot}}(t)],$$

$${\mathcal{J}}_j(\omega )=\frac{\pi }{2}\sum _i{M}_{ji}{\omega }_{ji}{d}_{ji}^2[\delta (\omega -{\omega }_{ji})-\delta (\omega +{\omega }_{ji})].$$

$$C_j(t)=\mathrm{Tr}\{u_j(t)u_j{\rho }_j^{\mathrm{ph}}\}$$

$$\begin{array}{rl}{C}_j(t)& ={\int }_{-\infty }^{+\infty }{\omega }^2{\mathcal{J}}_j(\omega )[n_{\mathrm{BE}}(\omega )+1]e^{-i\omega t}d\omega \\ & ={\int }_0^{\infty }{\omega }^2{\mathcal{J}}_j(\omega )\coth \frac{\beta \omega }{2}\cos \omega td\omega -i{\int }_0^{\infty }{\omega }^2{\mathcal{J}}_j(\omega )\sin \omega td\omega \end{array}$$

$${\omega }^2{\mathcal{J}}_j(\omega )=2{\lambda }_j\frac{\omega {\gamma }_j}{{\omega }^2+{\gamma }_j^2},$$

$$C_j(t)={\lambda }_j{\gamma }_j\left(\frac{2}{\beta {\gamma }_j}-i\right)e^{-{\gamma }_{j}t}$$

$$\dot{\rho }(t)=\sum _{j=1}^N{\Phi }_j(t){\sigma }_{{\mathbf{e}}_j}(t),$$

$${\dot{\sigma }}_{\mathbf{n}}(t)=-\sum _{j=1}^N{n}_j{\gamma }_j{\sigma }_{\mathbf{n}}(t)+\sum _{j=1}^N\left[{\Phi }_j(t){\sigma }_{\mathbf{n}+{\mathbf{e}}_j}(t)+{\lambda }_j{n}_j{\Theta }_j(t){\sigma }_{\mathbf{n}-{\mathbf{e}}_j}(t)\right].$$

$${\Theta }_j(t)=i{\gamma }_j\left(\frac{2}{\beta {\gamma }_j}[V_j(t),\cdot ]-i\{V_j(t),\cdot \}\right).$$

$$c_j={\lambda }_j{\gamma }_j\sqrt{{\left(\frac{2}{\beta {\gamma }_j}\right)}^2+1}.$$

$${\tilde{\sigma }}_{\mathbf{n}}(t)={\left(\prod _{j=1}^N{\lambda }_j^{n_j}{n}_j!\right)}^{-1}{\sigma }_{\mathbf{n}}(t).$$

$$\dot{\rho }(t)={\nu }^2\sum _{j=1}^N{\lambda }_j{\Phi }_j(t){\sigma }_{{\mathbf{e}}_j}(t),$$

$${\dot{\sigma }}_{\mathbf{n}}(t)=-\sum _{j=1}^N{n}_j{\gamma }_j{\sigma }_{\mathbf{n}}(t)+\sum _{j=1}^N\left[{\nu }^2(n_j+1){\lambda }_j{\Phi }_j(t){\sigma }_{\mathbf{n}+{\mathbf{e}}_j}(t)+{\Theta }_j(t){\sigma }_{\mathbf{n}-{\mathbf{e}}_j}(t)\right],$$

$$\mathcal{N}\gg \frac{{\omega }_{\max }}{\min ({\gamma }_1,\dots ,{\gamma }_N)},$$

$${\dot{\sigma }}_{\mathbf{n}}(t)=-\sum _{j=1}^N{n}_j{\gamma }_j{\sigma }_{\mathbf{n}}(t)+\sum _{j=1}^N{\Theta }_j(t){\sigma }_{\mathbf{n}-{\mathbf{e}}_j}(t),$$

$${\sigma }_{\mathbf{n}}(t)=\sum _{j=1}^N{\int }_0^{t}ds{e}^{-\left({\sum }_{k=1}^N{\gamma }_k{n}_k\right)s}{\Theta }_j(t-s){\sigma }_{\mathbf{n}-{\mathbf{e}}_j}(t-s).$$

$${\sigma }_{\mathbf{n}}(t)\cong \sum _{j=1}^N\left\{{\int }_0^{t}ds{e}^{-\left({\sum }_{k=1}^N{\gamma }_k{n}_k\right)s}{\Theta }_j(t-s)\right\}{\sigma }_{\mathbf{n}-{\mathbf{e}}_j}(t),$$

$$\begin{array}{rl}{\dot{\sigma }}_{\mathbf{n}}(t)& \cong -\sum _{j=1}^N{n}_j{\gamma }_j{\sigma }_{\mathbf{n}}(t)+\sum _{j=1}^N{\Theta }_j(t){\sigma }_{\mathbf{n}-{\mathbf{e}}_j}(t)\\ & +\sum _{j,k=1}^N{\nu }^2(n_j+1){\lambda }_j{\Phi }_j(t)\left\{{\int }_0^{t}ds{e}^{-\left({\sum }_{l=1}^N{\gamma }_l{n}_l+{\gamma }_j\right)s}{\Theta }_k(t-s)\right\}{\sigma }_{\mathbf{n}+{\mathbf{e}}_j-{\mathbf{e}}_k}(t).\end{array}$$

$$\dot{\rho }(t)\cong \sum _{j=1}^N{\nu }^2{\lambda }_j{\int }_0^{t}ds{\Phi }_j(t){\Theta }_j(t-s)e^{-{\gamma }_{j}s}\rho (t).$$

$${\sigma }_{\mathbf{n}}(t)=\sum _{m=0}^{\infty }{\nu }^{2m}{\sigma }_{\mathbf{n}}^{(m)}(t).$$

$$\dot{\rho }(t)=\sum _{m=1}^{\infty }{\nu }^{2m}\sum _{j=1}^N{\lambda }_j{\Phi }_j(t){\sigma }_{{\mathbf{e}}_j}^{(m-1)}(t),$$

$${\dot{\sigma }}_{\mathbf{n}}^{(m)}(t)=-\sum _{j=1}^N{n}_j{\gamma }_j{\sigma }_{\mathbf{n}}^{(m)}(t)+\sum _{j=1}^N\left[(n_j+1){\lambda }_j{\Phi }_j(t){\sigma }_{\mathbf{n}+{\mathbf{e}}_j}^{(m-1)}(t)+{\Theta }_j(t){\sigma }_{\mathbf{n}-{\mathbf{e}}_j}^{(m)}(t)\right],$$

$${\sigma }_{\mathbf{n}}^{(m)}(t)=\sum _{j=1}^N{\int }_0^{t}ds\left[(n_j+1){\lambda }_j{\Phi }_j(s){\sigma }_{\mathbf{n}+{\mathbf{e}}_j}^{(m-1)}(s)+{\Theta }_j(s){\sigma }_{\mathbf{n}-{\mathbf{e}}_j}^{(m)}(s)\right]e^{-\left({\sum }_{k=1}^N{\gamma }_k{n}_k\right)(t-s)}.$$

$$\begin{gathered} {\sigma }_{\mathbf{n}}^{(m)}(t)=\sum _{\begin{array}{c}\mathrm{paths}\\ (\mathbf{n};m)\to (0;0)\end{array}}{\int }_0^{t}d{s}_1{\int }_0^{s_1}d{s}_2\dots {\int }_0^{s_{2m+|\mathbf{n}|-1}}d{s}_{2m+|\mathbf{n}|} \\ \left[\prod _{l=1}^{2m+|\mathbf{n}|}{\mathrm{Op}}_l(s_l)\right]e^{-\left({\sum }_{k=1}^N{\gamma }_k{n}_k\right)t}\rho (s_{2m+|\mathbf{n}|}), \end{gathered}$$

$${\mathbb{Z}}_{+}^{N+1}=\{({\tilde{n}}_1,\dots ,{\tilde{n}}_N;\tilde{m})|{\tilde{n}}_i\ge 0\text{ for all }i,\tilde{m}\ge 0\}$$

$${\mathrm{Op}}_l(s)=n_j^{(l)}{\lambda }_j{\Phi }_j(s)e^{-{\gamma }_{j}s}\equiv n_j^{(l)}{\lambda }_j{\tilde{\Phi }}_j(s);$$

$${\int }_0^{t}d{s}_{2m+|\mathbf{n}|}{\int }_{s_{2m+|\mathbf{n}|}}^{t}d{s}_1{\int }_{s_{2m+|\mathbf{n}|}}^{s_1}d{s}_2\dots {\int }_{s_{2m+|\mathbf{n}|}}^{s_{2m+|\mathbf{n}|-2}}d{s}_{2m+|\mathbf{n}|-1}$$

$$\dot{\rho }(t)=\sum _{m=1}^{\infty }{\nu }^{2m}{\int }_0^{t}ds{\mathcal{K}}^{(2m)}(t,s)\rho (s),$$

$${\mathcal{K}}^{(2)}(t,s)=\sum _{j=1}^N{\lambda }_j{\tilde{\Phi }}_j(t){\tilde{\Theta }}_j(s),$$

Full text

1 Introduction

Since the experimental observations of quantum coherent effects in excitation energy transfer (EET) in photosynthetic light-harvesting complexes [1, 2], theoretical description of quantum EET in biological systems attracts much attention (see, e.g., [3]). This field of research offers also certain challenges for mathematical research.

One of the mathematical challenges is derivation of quantum master equation governing the dynamics of EET. An EET process in molecular systems can be described from the viewpoint of theory of open quantum systems [4, 5, 6]. In this case, the electronic degrees of freedom of molecules constitute “a system”, which is coupled to “a bath” (or “a reservoir”) consisting of the vibrational degrees of freedom and the environment of the molecules. Quantum master equations for the reduced density matrix of the system (i.e., of the electronic degrees of freedom) is a widely used tool of theory of open quantum system. There are well-known rigorous derivations of Markovian quantum master equations for the cases of weak system-bath coupling [7, 8, 9] and low density of the particles in the bath [10, 11, 12] as well as a general form of a generator of a quantum dynamical semigroup – the so called Gorini–Kossakowski–Sudarshan–Lindblad (GKSL) form [13, 14]. Another well-known limiting case is the weak dipole coupling between electronic excitations of molecules. This leads to Förster theory of EET [15, 16]. The corresponding EET dynamics is also Markovian. Note that, in the physical literature, weak system-bath coupling theory is also referred to as Redfield theory [17].

Unfortunately, the mentioned limiting regimes leading to Markovian dynamics are not satisfied in light-harvesting complexes. Hence, non-Markovian quantum master equations are required. A possible approach is the derivation of higher-order (non-Markovian) corrections to the known Markovian master equations. For example, the Redfield quantum master equation corresponds to the second-order perturbation theory with respect to the system-bath coupling. General methods using the projection operator technique or cumulant expansion method allow to derive the corrections of an arbitrary order within a given perturbation theory [4]. In [18], explicit expressions of the fourth-order corrections to the Redfield equation were obtained.

In this message, we propose another way of derivation of explicit expressions for higher-order corrections to the Redfield equation. This approach is based on the hierarchical equations of motion (HEOM) [19]. the HEOM is a widely used approach to the non-Markovian dynamics of open quantum systems. The application of this method to the description of EET in molecular systems was presented in [20, 21]. This method is numerically exact, but computationally expensive, hence, approximate equations based on the HEOM are desired. Note also that the HEOM is suitable only for special types of the spectral density (a function which specifies the system-bath interaction), hence, the same is true for approximate schemes derived from the HEOM.

Another important issue we address in this message is the range of validity of the Redfield equation and the additional secular approximation for it in terms of physical parameters. The mathematical derivations of master equations use formal limits, which are not suitable for physics since they do not give an answer whether the corresponding limiting case can be applied to a particular physical system or not. The conditions involving the physical parameters are required. In the physical literature, some heuristic conditions of applicability of the Redfield equation in terms of physical parameters are commonly used. An example is the condition that the dipole couplings between electronic degrees of freedom must be much larger than the bath relaxation rates [6]. However, as can be seen from numerical experiments. the Redfield equation works fine even when the dipole couplings are small provided that the system-bath couplings are also small with respect to the bath relaxation rates. Moreover, in this case, the Redfield equation can be simplified using the so called local approach [22, 23, 24].

There are debates concerning the secular approximation for the Redfield equation [5, 25, 26]. The secular approximation involves discarding of the highly-oscillating terms from the Redfield equation. This allows to express the equation in the GKSL form, which guarantees the preservation of positivity. But, in many cases, this approximation is inadequate and looses important effects. Note that, recently, other approximations to the Redfield equation leading to equations of the GKSL form were proposed [27, 28]. In some cases, the aforementioned local approach also leading to a GKSL equation can be used.

In this message we propose a general method of derivation of higher-order corrections to the Redfield quantum master equation based on the HEOM. In Sec. 2, we give the Frenkel exciton Hamiltonian and the HEOM. In Sec. 3, we derive the general formula of higher-order corrections to the Redfield equation. Finally, in Sec. 4, we obtain conditions of validity of the Redfield equation and the secular approximation for it using the derived corrections.

2 Theoretical background

The Hamiltonian describing EET processes in molecular aggregates is as follows (the so called Frenkel exciton Hamiltonian) [3, 20, 21]:

[TABLE]

where “h.c.” stands for Hermitian conjugate. Also we assume $\hbar=1$. Here $H^{\rm el}$ is the electronic (system) Hamiltonian: $N$ is the number of monomers (e.g., individual molecules) in the aggregate, $\ket{j}$ represents the excited electronic state of the $j$th site (molecule) with all other sites being in the ground state, $\varepsilon_{j}$ is the electronic excitation energy of the $j$th site, $J_{jk}$ is the dipole Coulombic coupling constant between the $j$th and $k$th sites. These coupling constants are responsible for EET between the sites. Each site $j$ is coupled to its own phononic bath consisting of harmonic oscillators, with $q_{ji}$ and $p_{ji}$ being the position and momentum operators of the $i$th phonon mode of the corresponding bath. The parameters $M_{ji}$ and $\omega_{ji}$ are the mass and frequency of the corresponding mode, and $d_{ji}$ is the displacement of the equilibrium configuration of the mode between the ground and excited electronic states of the site. These displacements $d_{ji}$ 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}}$.

We consider the initial state (density operator) of both electronic and phononic degrees of freedom of the form $\rho^{\rm tot}(0)=\rho(0)\otimes\rho^{\rm ph}$, where $\rho(0)$ is the initial electronic density matrix and $\rho^{\rm ph}=e^{-\beta H^{\rm ph}}/\operatorname{Tr}e^{-\beta H^{\rm ph}}$ is the equilibrium phononic state. The density operator in the interaction representation $\rho^{\rm tot}(t)=U_{I}(t)\rho^{\rm tot}(0)U_{I}(t)^{\dagger}$, $U_{I}(t)=e^{iH_{0}t}e^{iHt}$, $H_{0}=H^{\rm el}+H^{\rm ph}$, satisfies the von Neumann equation

[TABLE]

where $H^{\text{el-ph}}(t)=e^{iH_{0}t}H^{\text{el-ph}}e^{-iH_{0}t}$, $[A,B]=AB-BA$, and (in the following) $\{A,B\}=AB+BA$.

The system-bath interaction is often specified in terms of the spectral density. The spectral density function for the bath coupled to the $j$th site defined as

[TABLE]

Then, the correlation function of the $j$th bath

[TABLE]

where $u_{j}(t)=e^{iH^{\rm ph}_{j}t}u_{j}e^{-iH^{\rm ph}_{j}t}$, can be expressed as

[TABLE]

Here $n_{\rm BE}(\omega)=(e^{\beta\omega}-1)^{-1}$ is the Bose–Einstein distribution, $\beta$ is the inverse temperature of the baths.

Consider the Drude–Lorentz spectral density:

[TABLE]

where $\lambda_{j}$ and $\gamma_{j}$ are the reorganization energy and Debye (cutoff) frequency of the corresponding bath. 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.

Following [20, 21], we impose a high-temperature condition $\beta\gamma_{j}\ll 1$ for all $j$. Then we can approximate $\coth(\beta\omega/2)$ by $2/(\beta\omega)$ in Eq. (5) and obtain the expression for the correlation function:

[TABLE]

for $t\geq 0$ and $C_{j}(-t)=C^{*}_{j}(t)$.

Then, in [21], it is shown that the dynamics of the reduced density operator of the system $\rho(t)=\operatorname{Tr}_{R}\rho^{\rm tot}(t)$ satisfies the following hierarchical equations of motion (HEOM):

[TABLE]

Here $\mathbf{n}=(n_{1},\ldots,n_{N})\in\mathbb{Z}_{+}^{N}$ (i.e., all $n_{j}$ are non-negative integers), $\mathbf{e}_{j}=(0,\ldots,0,1,0,\ldots,0)$, where 1 on the $j$th position, $\Phi_{j}(t)$ and $\Theta_{j}(t)$ are the superoperators defined as $\Phi_{j}(t)=i[V_{j}(t),\,\cdot\,]$ and

[TABLE]

The operators $\sigma_{\mathbf{n}}$ act in the same Hilbert space as $\rho$ (namely, in $\mathbb{C}^{N}$) and are referred to as auxiliary density operators (ADOs). Note that Eq. (8a) can be expressed as a particular case of Eq. (8b) if we put $\rho(t)=\sigma_{\mathbf{0}}(t)$, where $\mathbf{0}=(0,\ldots,0)$. Also we put by definition $\sigma_{\mathbf{n}}(t)\equiv 0$ if $\mathbf{n}$ has negative elements. The initial conditions for the ADOs are $\sigma_{\mathbf{n}}(0)=0$ for all $\mathbf{n}$ with positive elements.

Suppose that $\rho(t)$ is of order one, then, from Eq. (8b), the magnitude of $\sigma_{\mathbf{n}}$ is proportional to $\prod_{j=1}^{N}c_{j}^{n_{j}}$, where

[TABLE]

Hence, the magnitude of $\sigma_{\mathbf{n}}$ indefinitely increases for large as $|\mathbf{n}|\equiv n_{1}+\ldots n_{N}$ if $c_{j}>1$ for all $j$, and the last condition cannot be excluded. In [29], a rescaling of ADOs was proposed. For our purposes, it will be convenient to introduce a slightly different rescaling, namely:

[TABLE]

Then the magnitudes of $\tilde{\sigma}_{\mathbf{n}}(t)$ are bounded and tend to zero as $|\mathbf{n}|\to\infty$. Also, since we will study the limiting case of small system-bath couplings, we substitute all $d_{ji}$ in Eq. (1d) by $\nu d_{ji}$, where $\nu$ is a formal small dimensionless parameter. This is equivalent to the replacement of all $\lambda_{j}$ by $\nu^{2}\lambda_{j}$, see Eqs. (3) and (6). Then the hierarchy (8) is rewritten as

[TABLE]

where we have removed the tildes from the ADOs since, in the following, we will consider only rescaled ADOs.

There are several ways of truncation of the infinite hierarchy of equations (10). One commonly accepted way is setting $\sigma_{\mathbf{n}}(t)\equiv 0$ for $|\mathbf{n}|$ larger some threshold $\mathcal{N}$. According to [20, 21], this threshold value should satisfy

[TABLE]

where $\omega_{\max}$ is the largest difference between the eigenvalues of $H^{\rm el}$. Then, for $|\mathbf{n}|=\mathcal{N}$, we have

[TABLE]

or,

[TABLE]

If we assume that $\sigma_{\mathbf{n}-\mathbf{e}_{j}}(t-s)$ evolves much slower than $e^{-\left(\sum_{k=1}^{N}\gamma_{k}n_{k}\right)s}\,\Theta_{j}(t-s)$ decays (with the increase of $s$), then we can perform the following approximation:

[TABLE]

and, hence, for $|\mathbf{n}|=\mathcal{N}-1$,

[TABLE]

So, $\sigma_{\mathbf{n}}(t)$ for $|\mathbf{n}|=\mathcal{N}$ do not have to be stored in the computer memory and $|\mathbf{n}|=\mathcal{N}-1$ is actually the last level of the hierarchy. Equations (11) and (12) for the last level of the hierarchy represent two truncation schemes. The first scheme was considered, for example, in [20, 21], the second one was considered in [31, 30, 29, 32].

In particular, Eq. (12) for $|\mathbf{n}|=\mathbf{0}$, i.e., for $\sigma_{\mathbf{0}}(t)\equiv\rho(t)$, is the Redfield equation (this is a well-known fact and the details will be given in Sec. 4):

[TABLE]

This corresponds to the second-order perturbation with respect to $\nu$. Extending the upper limit of integration in Eq. (13) to infinity corresponds to the Markovian approximation. Correspondingly, the hierarchies (10) with the truncation (12) for $|\mathbf{n}|\leq 1$ provides higher-order corrections to the Redfield equation with respect to the small parameter $\nu$. In the next section, we propose higher-order corrections to the Redfield equation not involving ADOs, i.e. in the form of closed single equations for $\rho(t)$ (like the Redfield equation itself and the fourth-order correction in [18]).

3 Higher-order quantum master equations

Expand the ADOs $\sigma_{\mathbf{n}}(t)$, $\mathbf{n}\neq\mathbf{0}$, in the formal series with respect to the small parameter $\nu$:

[TABLE]

The substitution of this expansion into the hierarchy (10) and equating the expressions with equal orders of $\nu$ in both sides of Eq. (10b) to each other, gives

[TABLE]

$m=0,1,2,\ldots,$ with the agreement $\sigma^{(m)}_{\mathbf{n}}(t)\equiv 0$ for $m=-1$. The $2M$th-order quantum master equation corresponds to the neglection of the terms for $m>M$ in Eq. (14a) and substitution of the functions $\sigma^{(m)}_{\mathbf{e}_{j}}(t)$ for $0\leq l\leq M-1$ and $j=1,\ldots,N$ by their explicit expressions via $\rho(t)$.

Rewrite Eq. (14b) in the form

[TABLE]

Iterating this expression, we obtain

[TABLE]

where summation is over admissible paths on the lattice

[TABLE]

starting from $(\mathbf{n};m)$ and ending with $(\mathbf{0};0)$. A path is a sequence of points $(\mathbf{n}^{(l)};m^{(l)})$, $l=0,\ldots,L$, in $\mathbb{Z}_{+}^{N+1}$. Only the following steps from $(\mathbf{n}^{(l-1)};m^{(l-1)})$ to $(\mathbf{n}^{(l)};m^{(l)})$ are admissible:

1. (i)

$\tilde{\mathbf{n}}^{(l)}=\tilde{\mathbf{n}}^{(l-1)}+\mathbf{e}_{j}$ for some $j$ and $\tilde{m}^{(l)}=\tilde{m}^{(l-1)}-1$. Then, in (15),

[TABLE] 2. (ii)

$\tilde{\mathbf{n}}^{(l)}=\tilde{\mathbf{n}}^{(l-1)}-\mathbf{e}_{j}$ for some $j$ and $\tilde{m}^{(l)}=\tilde{m}^{(l-1)}$. Then ${\rm Op}_{l}(s)=\Theta_{j}(s)e^{\gamma_{j}s}\equiv\tilde{\Theta}_{j}(s)$.

Also transitions to the points $(\mathbf{0};\tilde{m})$ for $\tilde{m}>0$ are forbidden.

The end point of the path is $(\mathbf{0};0)$. To complete the path starting from the point $(\mathbf{n};m)$, we have to do overall $m$ steps of type (i), which reduces the last coordinate by one. But the “cost” of each such step is the increase of $|\tilde{\mathbf{n}}|$ by one. Hence, we also have to do overall $m+|\mathbf{n}|$ steps of type (ii) in order to reduce $|\tilde{\mathbf{n}}|$ to zero. Thus, the length of the path from $(\mathbf{n};m)$ to $(\mathbf{0};0)$ is $L=2m+|\mathbf{n}|$.

The order of integration in Eq. (15) can be changed to

[TABLE]

The substitution of Eq. (15) with order of integration (16) to (14a) gives the formal series

[TABLE]

where

[TABLE]

$m\geq 2$, with the agreement $s_{0}=t$ and $s_{2m-1}=s$. The paths in (18b) start from the point $(\mathbf{0};m)$. The admissible steps are given above. Transitions to the points $(\mathbf{0};\tilde{m})$ for $0<\tilde{m}<m$ are still forbidden.

As we said before, the master equation of order $2M$ corresponds to keeping the first $M$ terms in (17) and neglecting all the remaining terms.

Equation (17) contains convolution with respect to time and, thus, is time-nonlocal. Let us derive a time-local (convolutionless) equation, see [4, 5, 6] for general theory. For this aim, we can express

[TABLE]

Iteration of this equation and the substitution of the result to Eq. (17) gives

[TABLE]

where

[TABLE]

Again, the time-local master equation of order $2M$ corresponds to keeping the first $M$ terms in (17) and neglecting all the remaining terms.

Formulas (17)–(18) for the time-nonlocal corrections to the Redfield equation and (19)–(20) for the time-local (convolutionless) ones are the main result of this section.

Consider the examples. For $m=1$, there are exactly $N$ admissible paths from $(\mathbf{0};1)$ to $(\mathbf{0};0)$:

[TABLE]

where the operators ${\rm Op}$ corresponding to the steps (see Eq. (15)) are indicated above the arrows. Hence, for $M=2$, the Redfield equation (13) is restored.

Consider the fourth-order corrections: $m=2$. The admissible paths from $(\mathbf{0};2)$ to $(\mathbf{0};0)$ are:

[TABLE]

Hence,

[TABLE]

For $N=1$ (single bath), formula (21) is reduced to a particular case (for the Drude–Lorentz spectral density (6) and high-temperature approximation) of the fourth-order kernel derived in [18]. Formula (21) for $N>1$ provides a generalization of such kernel for multibath case. Note that, in contrast to $\mathcal{K}^{(2)}$, the fourth-order kernel $\mathcal{K}^{(4)}$ is not reduced to a sum of independent contributions from each bath: the terms with $j\neq k$ describe the interaction between the baths via the system.

Also

[TABLE]

Analogously, the sixth-order kernel can be derived:

[TABLE]

[TABLE]

4 On the range of validity of the Redfield equation

In this section we discuss the range of validity of the Redfield equation (13) as well as the secular approximation for it. At first, let us express it in a more convenient way. Denote $\operatorname{spec}H^{\rm el}$ the spectrum of $H^{\rm el}$ and

[TABLE]

the set of Bohr frequencies, i.e., differences between the eigenvalues of $H^{\rm el}$. Note that if $\omega\in\Omega$, then also $-\omega\in\Omega$. For $\varepsilon\in\operatorname{spec}H^{\rm el}$, denote $P_{\varepsilon}$ the projector onto the eigensubspace corresponding to $\varepsilon$. Also, for $\omega\in\Omega$, denote

[TABLE]

Then

[TABLE]

the same equalities are satisfied for $\Theta_{j}$ as well.

If we express the terms in (13) as

[TABLE]

then, using $V_{j}^{(-\omega)}=V_{j}^{(\omega){\dagger}}$, where ${\dagger}$ denotes the Hermitian conjugation, the Redfield equation can be rewritten as

[TABLE]

where h.c. stands for Hermitian conjugate terms, and

[TABLE]

For $t\gg\gamma_{j}^{-1}$, we can take the limit $t\to\infty$ in Eq. (27) and equation (26) is reduced to an autonomous matrix differential equation:

[TABLE]

where

[TABLE]

If the time scale of evolution of $\rho(t)$ is much larger than $\max_{j}\gamma_{j}^{-1}$, then we can replace $\Gamma_{j}(\omega,t)$ by $\Gamma(\omega)$ in (26) for all times: on small times, where the finiteness of the limit of integration in Eq. (27) is essential, $\rho(t)$ does not succeed to evolve significantly and the introduced error is small. In other words, equation (28) can be used for the description of evolution of $\rho(t)$ for all $t\geq 0$.

The solution of Eq. (28) can be regarded as an action of a semigroup $\{\Lambda_{t}\}_{t\geq 0}$: $\rho(t)=\Lambda_{t}\rho(0)$. The semigroup property $\Lambda_{s}\Lambda_{t}=\Lambda_{s+t}$ represents the Markovian property of the dynamics. Equation (28) is referred to as the Markovian Redfield equation. Accordingly, equation (26) is often referred to as the non-Markovian Redfield equation, but it differs from Eq. (28) only on small times (of order $\max_{j}\gamma_{j}^{-1}$).

Unfortunately, the semigroup generated by Eq. (28), generally speaking, does not preserve the positivity: the condition of positive semidefiniteness of $\rho(t)$ may be not satisfied even if this condition is satisfied for $\rho(0)$. To fix this drawback, a further approximation is commonly made. Equation (28) contains oscillating terms proportional to $e^{i(\omega^{\prime}-\omega)t}$ for $\omega^{\prime}\neq\omega$. If the time scale of evolution of $\rho(t)$ is much larger than $\max|\omega^{\prime}-\omega|^{-1}$, where the maximum is taken over $\omega,\omega^{\prime}\in\Omega$, $\omega\neq\omega^{\prime}$, then the terms with $\omega\neq\omega^{\prime}$ can be neglected as rapidly oscillating. This approximation is called the secular approximation. Then Eq. (28) is reduced to

[TABLE]

where

[TABLE]

is the Lamb shift Hamiltonian. Equation (30) has the Gorini–Kossakowski–Sudarshan–Lindblad (GKSL) form and, hence, preserves the positivity ([13, 14], see also [4, 5]). Equation (30) is referred to as the secular Redfield equation.

Now we discuss the range of validity of the Redfield equation (13) (or (26)) as well as its Markovian and secular versions (28) and (30) in terms of physical parameters rather than formal dimensionless parameter $\nu$. So, we put $\nu=1$ since this is a formal parameter and all information about the system-bath interaction is given in the physical parameters like $\lambda_{j}$, $\gamma_{j}$, etc.

The time scale of evolution of $\rho(t)$ in Eqs. (26)–(30) is, roughly, $\min_{j,\omega}|\Gamma_{j}(\omega)|^{-1}$, where

[TABLE]

Note that we work in the high-temperature approximation $\beta\gamma_{j}<1$ and, thus, expression (31) can be simplified to $|\Gamma_{j}(\omega)|=2\lambda_{j}\gamma_{j}/\beta\sqrt{\gamma_{j}^{2}+\omega^{2}}$. However, we keep the more general expression (31) trying to extrapolate the results to the low-temperature case as well. According to the previous discussion, the range of validity of the Markovian approximation is $\max_{j,\omega}|\Gamma_{j}(\omega)|\ll\min_{j}\gamma_{j}$, or

[TABLE]

Also the range of validity of the secular approximation is

[TABLE]

It remains to establish the range of validity of the second-order approximation

[TABLE]

which was the starting point for all equations of this section. To estimate its range of validity, let us analyze the fourth-order corrections:

[TABLE]

If the fourth-order corrections are small, we can conclude that the second-order approximation is adequate. Of course, this is not a rigorous proof, but rough estimates. The second term in the square brackets (the double integral with $\mathcal{K}^{(2)}$) represents the evolution of $\rho$ from $s$ to $t$ in the second-order approximation. The influence of the second term in the square brackets is negligible if the evolution of $\rho$ in the second-order approximation is much slower than the decay of $\mathcal{K}^{(2)}(t,s)$. Since the $\mathcal{K}^{(2)}(t,s)$ decays as $\exp[-\min_{j}\gamma_{j}(t-s)]$, and the rate of evolution of $\rho$ is given by $\max_{j,\omega}|\Gamma(\omega)|$, we again arrive at the condition (32).

Consider not the last term of (34):

[TABLE]

For rough estimation of this term, consider the case of large $t$, namely, $t\gg\max_{j}\gamma_{j}^{-1}$. Then, using the representation (25), we obtain

[TABLE]

The magnitude of each term is

[TABLE]

We can estimate the last factor in (35) from above by putting $\omega_{4}+\omega_{3}=0$ and $\omega_{4}+\omega_{3}-\omega_{2}=0$. Then, expression (35) is much smaller than the magnitude of the second-order contribution (31) (for $\omega=\omega_{4}$) whenever

[TABLE]

We can see that this is a more stringent version of condition (32).

Thus, condition (36) is a (rough) sufficient condition for the validity of the Redfield equation (both Markovian and non-Markovian versions), and (33) is a sufficient condition of the secular approximation. This two inequalities are main results of this section.

5 Conclusions

We have derived explicit general formulas for the corrections to the Redfield equation of arbitrary order with respect to the system-bath coupling based on the HEOM: formulas for the time-nonlocal equation (17)–(18) and time-local (convolutionless) one (19)–(20). For the considered particular case (see below), the fourth-order corrections coincide with the corrections previously derived in the framework of the projection operator method [18]. For high orders, the explicit expressions of the equations become cumbersome and, hence, not practically useful. However, the existence of a general scheme which allows to obtain closed equations for the density matrix (without auxiliary density operators) of arbitrary order seems to be of some interest.

In this work, we adopted the HEOM for the special case of the Drude–Lorentz spectral density (6) and the high-temperature approximation. Also we considered a particular system-bath Hamiltonian. In this case, the bath correlation function is an exponent (7). In general, the HEOM is applicable to the case when the correlation function is a sum of exponents. The proposed method of derivation of the corrections to the Redfield equation works for this case as well. Theoretically, the exponential decomposition for the bath correlation functions can be applied for the general case, but it may be computationally inefficient.

Recently, a HEOM based on an alternative a decomposition scheme for the bath correlation functions was proposed [33]. Namely, the decomposition scheme in [33] is based on the Chebyshev polynomials and Bessel functions. The corresponding HEOM is referred to as the C-HEOM (Chebyshev-HEOM). It is suitable for the cases where the exponential decomposition is inefficient, but has its own shortcomings. It would be interesting to derive corrections to the Redfield equations based on the C-HEOM.

Also we have derived conditions of the validity of the Redfield equation (Ineq. (36)) as well as the secular approximation (Ineq. (33)) for it in terms of physical parameters. The derivation of Ineq. (36) is still heuristic, but has more solid basis that the heuristic conditions commonly used in the physical literature. Rigorous derivation of such conditions is still required. Note that, in [34], the validity of the Redfield equation and the secular approximation is studied for the model which is exactly solvable by the pseudomode method [35]. A heuristic derivation of the range of validity of the Förster and modified Redfield theories is proposed in [36].

Acknowledgments. This work was supported by the Russian Science Foundation (project 17-71-20154).