Quantum corrections of the truncated Wigner approximation applied to an exciton transport model

TL;DR

This paper develops a modified path integral approach for exciton transport in open quantum systems, enabling quantum corrections beyond the truncated Wigner approximation and demonstrating its application to a donor-acceptor model.

Contribution

It introduces a new approximation method that extends the truncated Wigner approximation for exciton transport, based on a modified path integral formalism.

Findings

The new method captures quantum fluctuations more accurately.

Application to a donor-acceptor model shows improved results.

The approach provides a clearer understanding of quantum effects in exciton transport.

Abstract

We modify the path integral representation of exciton transport in open quantum systems such that an exact description of the quantum fluctuations around the classical evolution of the system is possible. As a consequence, the time evolution of the system observables is obtained by calculating the average of a stochastic difference equation which is weighted with a product of pseudo-probability density functions. From the exact equation of motion one can clearly identify the terms that are also present if we apply the truncated Wigner approximation. This description of the problem is used as a basis for the derivation of a new approximation, whose validity goes beyond the truncated Wigner approximation. To demonstrate this we apply the formalism to a donor-acceptor transport model.

Equations

$$\text{allowdisplaybreaks}\hat{H}$$

$${\hat{H}}_S$$

$${\hat{H}}_{B,SB}$$

$${\hat{\rho }}_{tot}$$

$${\hat{\rho }}_B^{eq}$$

$$Z$$

$$i{\mathcal{S}}_S=$$

$${\mathbf{a}}_n(t)=$$

$$i{\mathcal{S}}_{B,SB}^n=$$

$$+\frac{i}{\pi }{\int }_0^{t}d\tau {\int }_0^{\tau }d{\tau }^{\prime }{n}_n^{\times }(\tau )n_n^o({\tau }^{\prime }){\mathcal{D}}_n(\tau -{\tau }^{\prime }),$$

$$n_n^{\times }=$$

$$n_n^o=$$

$${\mathcal{F}}_n(t)$$

$${\mathcal{D}}_n(t)$$

$$J_n(\omega )$$

$$\displaystyle{\rm tr}_{S}\big{[}\hat{O}\hat{\rho}(t)\big{]}$$

$$D[\mathbf{\psi },\mathbf{\eta }]$$

$$a_n^{f/b}$$

$$i{\mathcal{S}}_S$$

$$+{\mathbf{\eta }}^{\ast T}(\tau )(i{\partial }_{\tau }-h)\mathbf{\psi }(\tau )],$$

$$\mathbf{f}(t)$$

$${\int }_0^{t}d\tau {\mathbf{\psi }}^{\ast T}(\tau ){\partial }_{\tau }\mathbf{\eta }(\tau )=$$

$$\displaystyle\hskip 36.98857pt-\int^{t}_{0}d\tau\bm{\eta}^{T}(\tau)\partial_{\tau}\bm{\psi}^{*}(\tau)+\bm{\psi}^{*T}(\tau)\bm{\eta}(\tau)\Big{|}^{t}_{0},$$

$${\rho }_{\mathcal{W}}({\mathbf{\psi }}^{\ast },\mathbf{\psi })=\prod _{n=1}^{\mathcal{N}}\int \frac{d\Re {\eta }_{n}d\Im {\eta }_n}{2{\pi }^2}$$

$$\times e^{-\frac{1}{4}|{\eta }_n{|}^2-|{\psi }_n{|}^2+\frac{1}{2}({\psi }_n^{\ast }{\eta }_n-{\eta }_n^{\ast }{\psi }_n)}⟨\mathbf{\psi }+\frac{1}{2}\mathbf{\eta }|{\hat{\rho }}_S|\mathbf{\psi }-\frac{1}{2}\mathbf{\eta }⟩,$$

$${\mathcal{O}}_{\mathcal{W}}(\mathbf{\psi },{\mathbf{\psi }}^{\ast })=\prod _{n=1}^{\mathcal{N}}\int \frac{d\Re {\eta }_{n}d\Im {\eta }_n}{2\pi }$$

$$\times e^{-\frac{1}{2}|{\eta }_n{|}^2}⟨\mathbf{\psi }-\frac{1}{2}\mathbf{\eta }|\hat{O}|\mathbf{\psi }+\frac{1}{2}\mathbf{\eta }⟩,$$

$${\rho }_{\mathcal{W}}({\mathbf{\psi }}^{\ast },\mathbf{\psi })$$

$$\times {\rho }_{\mathcal{W}}^{abs}({\mathbf{\psi }}^{\ast },\mathbf{\psi }),$$

$${\rho }_{\mathcal{W}}^{abs}({\mathbf{\psi }}^{\ast },\mathbf{\psi })$$

$$N({\rho }_{\mathcal{W}})$$

$$\mathrm{Sign}({\rho }_{\mathcal{W}},{\mathbf{\psi }}^{\ast },\mathbf{\psi })$$

$${\int }_{-\infty }^{\infty }\frac{d\phi }{\pi }{e}^{i2\phi f}$$

$$d\mathbf{\psi }(t)$$

$${\mathcal{O}}_{\mathcal{W}}({\mathbf{\psi }}^{\ast }(t),\mathbf{\psi }(t))N({\rho }_{\mathcal{W}})\mathrm{Sign}({\rho }_{\mathcal{W}},{\mathbf{\psi }}_0^{\ast },{\mathbf{\psi }}_0).$$

$$\frac{1}{\pi }\mathcal{D}(t)$$

$$\frac{1}{\pi }\mathcal{F}(t)$$

$$\displaystyle\exp\big{[}iS^{n}_{B,SB}\big{]}$$

$$\displaystyle\times\exp\big{[}iS^{n,TL}_{B,SB}\big{]}\rho^{n}_{\Phi}(\bm{\phi}^{n}_{0}),$$

$$i{S}_{B,SB}^{n,TL}$$

$$+\Delta t{\mathbf{v}}^n{n}_{n,\tau }^o+$$

$$\int D[{\mathbf{\phi }}^n]$$

$$\int D[{\mathbf{\varphi }}^n]$$

$$\int D[{\mathbf{x}}^n]$$

Full text

exciton transport, truncated Wigner approximation

I Introduction

The exciton transport in large networks that are coupled to some environment is still a topic that attracts attention because of the experimental results from the last decade. Ultrafast non-linear spectroscopy used to probe the energy transfer in FMO complexes has revealed strong vibrational coherences of multiple pigments, that persist at much larger times scales than initially expected Brixner et al. (2005); Engel et al. (2007); Fujihashi et al. (2015). This has opened a discussion on the question whether or not these coherences have a significant impact on the efficiency of transport processes.

To study theoretically possible quantum effects that are responsible for efficient transport properties of a given network, one has to pay special attention to the way how the effect of the environment on the system is modeled. The chosen description has to be simple enough to lead to an efficient and fast scheme for the simulation of different realizations of the network and it has to be able to take into account the most important effects of the environment on the system. For example, one can construct a time local equation for the system density matrix where the effects of the environment, like dephasing, are described by Lindblad operators. To take into account the non-Markovian effects arising from the interaction between the system and the environment we have to model the latter as an infinitely large set of harmonic oscillators, that are linearly coupled to the system. These new degrees of freedom can be integrated out exactly from the problem, which makes the effective action of the system non-quadratic in the system fields and also time non-local. This description of the problem requires the use of advanced numerical methods like the Hierarchical equations of motion (HEOM) technique Tanimura and Kubo (1989); Tanimura (1990); Ishizaki and Tanimura (2005); Tanimura (2006); Zhu et al. (2013); Moix and Cao (2013), the quasi-adiabatic propagator path integral method Makri and Makarov (1995a, b); Thorwart et al. (2000), the multilayer multiconfiguration time-dependent Hartree approach Wang et al. (2006), the density matrix renormalization group method Chin et al. (2010) and other methods Mühlbacher and Kleinekathöfer (2012); Segal et al. (2010); Weiss et al. (2008), which in most cases are numerically expensive.

We are interested in approximate methods, where one can find a good compromise between the accuracy of the results and the applied numerical effort. One of the methods, whose complexity scales linearly with the size of the network and thus can be preferred for the case of working with large networks, is the truncated Wigner approximation (TWA). In this approximation one neglects some of the non-quadratic terms in the effective action. This approximation is often applied to describe the dynamics of ultracold atomic gases Polkovnikov (2003a); Isella and Ruostekoski (2005, 2006); Scott et al. (2006); Martin and Ruostekoski (2010); Norrie et al. (2006); Polkovnikov and Wang (2004); Barnett et al. (2011) and to explain some experimental results Gross et al. (2011). Its validity for closed systems at short time scales was studied in Polkovnikov (2003b). Its ability to reproduce the non-Gaussian statistics of the single mode anharmonic oscillator was shown in Corney and Olsen (2015), and its applicability to the calculation of multitime correlation functions was studied in Berg et al. (2009). Relying on these results the application of the approach to open quantum systems was also made in Ivanov et al. (2013), assuming that the effect of the environment, expressed mainly in the form of additional noise, will wash out all quantum scattering processes which occur at long times.

In the following, a different two-step approximation that extends beyond the TWA, will be derived for the special case of working with a Frenkel exciton Hamiltonian. In the first step, similarly to the TWA, some of the variables in the corresponding path integral representation of the problem will be analytically integrated out, the integration of the remaining set of variables in the path integral will be equivalent to solving an equation of motion of these variables. We will do the analytic integration without neglecting any terms from the action, but the price that we have to pay is the introduction of a set of pseudo-probability density functions in the path integral representation of the problem. From the obtained equations of motion we can clearly identify the contributions, that were also present in the TWA as well as the new terms. Since this representation of the problem is numerically expensive, in the second step we adopt a numerically applicable approximation that extends beyond the TWA for the case of preparing the system initially in the single exciton manifold. We will refer to it as the corrected truncated Wigner approximation (CTWA). The ability of both approximations to reproduce the exact system dynamics will be studied and we will see that both approximations reproduce exactly the short time behavior of the system. In the limit of weak couplings the CTWA can also extend this range to longer time scales.

The paper is structured as follows. In Sec. II.1 we introduce the system of interest and its path integral representation, in Sec. II.2 the Wigner formalism is briefly described, in Sec. II.3 a mapping procedure is presented, which makes the action local in time at the cost of introducing integrals over new variables. This step is necessary if we want to integrate out exactly the fields of the path integral representation of the system. In Sec. II.4 the TWA is briefly explained and the corresponding stochastic equation for the semi-classical description of our problem is derived. It will be used later for comparison with the equation obtained from CTWA. An exact stochastic equation, which takes all terms into account, that have been neglected in the TWA is derived in Sec. III and a numerically applicable approximation of this equation is derived in Sec. IV. The accuracy of the new approximation is tested in Sec. V for a donor-acceptor transport model. The paper concludes with a summary in Sec. VI.

II Theory

II.1 Model system

We consider a system (S) composed of $\mathcal{N}$ sites, each of them being linearly coupled to a bath (B) of harmonic oscillators. The corresponding Hamiltonian is given by:

[TABLE]

where $\hat{a}^{\dagger}_{n},\hat{b}^{\dagger}_{nk}$ create an excitation at the $n$-th site of energy $h_{nn}$ or in the $k$-th bath mode of energy $\omega_{nk}$ that is coupled to the $n$-th site. The initial state of the total system is of the form:

[TABLE]

where $\hat{\rho}_{S}$ and $\hat{\rho}^{eq}_{B}$ are the density matrix of the system and the equilibrium density matrix of the bath, $\beta=1/(k_{B}T)$, $T$ is the temperature and $N$ is a normalization constant. We assume that the initial state of the system $\hat{\rho}_{S}$ lies in the single exciton subspace. Since the Hamiltonian (1) conserves the total number of excitations it reduces to the Frenkel exciton Hamiltonian Yang and Fleming (2002). Since we work in the single exciton subspace we have the freedom to choose the operators $\hat{a}^{\dagger}_{j},\hat{a}_{j}$ to be bosonic or fermionic. In the following we assume that they are bosonic.

We are interested in the path integral representation of the Keldysh partition function of the problem Kamenev (2011). After analytically integrating out the bath degrees of freedom the expression has the following form:

[TABLE]

where $i\mathcal{S}_{S}$ origins from the system Hamiltonian $\hat{H}_{S}$ and $i\mathcal{S}^{n}_{B,SB}$ describes the effect of the $n$-th bath on the system. They are given by Ishizaki and Fleming (2009); Tanimura (2012); Tanimura and Kubo (1989); Tanimura (1990); Ishizaki and Tanimura (2005); Tanimura (2006):

[TABLE]

where $a^{f}_{n}/a^{b}_{n}$ are the fields lying on the forward/backward part of the Keldysh contour and $\sigma_{z}$ is the Pauli $z$ matrix. The noise ($\mathcal{F}$) and dissipation ($\mathcal{D}$) kernels are given by:

[TABLE]

where $J_{n}(\omega)$ is the spectral density.

II.2 Wigner formalism

We can express the expectation value of some system operator $\hat{O}$ in the following path integral form:

[TABLE]

where $\psi_{n}/\eta_{n}$ are the quantum/classical fields which are obtained from the variable transformation

[TABLE]

and $\mathcal{S}=\mathcal{S}_{S}+\sum_{n}\mathcal{S}^{n}_{B,SB}$. The $i\mathcal{S}^{n}_{B,SB}$ term is defined in (II.1) where $n^{\times}_{n}=\psi^{*}_{n}\eta_{n}+\eta^{*}_{n}\psi_{n}$, $n^{o}_{n}=2|\psi_{n}|^{2}+\frac{1}{2}|\eta_{n}|^{2}$ and $i\mathcal{S}_{S}$ is given by:

[TABLE]

where the elements $h_{nn^{\prime}}$ of the matrix $h$ are given in (2). Eq. (18) can be obtained from $i\mathcal{S}_{S}$, defined in (7), by applying the variable transformation given in (17) and then using the identity

[TABLE]

and neglecting the boundary terms. The Wigner distribution function $\rho_{\mathcal{W}}$ and Weyl symbol of the operator $\mathcal{O}_{\mathcal{W}}$ are defined as follows:

[TABLE]

where the operator $\hat{O}$ is in its normal ordered form and $|\bm{\psi}\pm\frac{1}{2}\bm{\eta}\rangle$ is a coherent state with the property $\hat{a}_{n}|\bm{\psi}\pm\frac{1}{2}\bm{\eta}\rangle=(\psi_{n}\pm\frac{1}{2}\eta_{n})|\bm{\psi}\pm\frac{1}{2}\bm{\eta}\rangle$ and $\langle\bm{\psi}|\bm{\psi}^{\prime}\rangle=\exp[\sum^{\mathcal{N}}_{n=1}\psi^{*}_{n}\psi^{\prime}_{n}]$. The same expression can be obtained if we represent the operator $\hat{O}$ in terms of symmetrized polynomials of $\hat{a}^{\dagger}_{n}$ and of $\hat{a}_{n}$ and then carry out the replacement $(\hat{a}^{\dagger}_{n},\hat{a}_{n})\rightarrow(\psi^{*}_{n},\psi_{n})$. The Wigner transform $\rho_{\mathcal{W}}$ of the density matrix is, in general, a pseudo-probability density function, i.e. it is normalized to one but it can take negative values. For the practical calculation of the path integral one can replace it with the following expression:

[TABLE]

In this case $\rho^{abs}_{\mathcal{W}}$ is a real probability density function. The derivation of Eqs. (15), (16), (18), (21) and (22) is shown in Polkovnikov (2003b).

If we neglect the contribution $i\mathcal{S}_{B,SB}$ from the environment we can integrate out the $\eta$-variables by the use of the equation

[TABLE]

for $\varphi\in\{\Re\eta,\Im\eta\}$. If we keep the $\psi_{n,0}$ variables fixed, the set of Dirac-delta function will define a unique path for the time evolution of the $\psi_{n,\tau}$ ($\tau>0$) variables. The path is described by the following equation:

[TABLE]

A single contribution to the observable ${\rm tr}[\hat{O}\hat{\rho}(t)]$ is given by:

[TABLE]

To calculate the expectation value of ${\rm tr}[\hat{O}\hat{\rho}(t)]$ we have to sample $\psi_{n,0}$ from $\rho^{abs}_{\mathcal{W}}$ and calculate the mean value of (31). To take into account the effect of the environment, i.e. of $i\mathcal{S}_{B,SB}$, we have to apply to it additional transformations which will be explained in the next subsection.

II.3 New mapping of the time non-local parts of the action to time local expressions

We will map the time non-local part of the action to a time local expression at the cost of introducing integrals over an additional set of variables which is a special realization of the idea proposed in Stockburger and Grabert (2002). To apply this mapping we assume that the noise ($\mathcal{F}$) and dissipation kernel ($\mathcal{D}$) of the action can always be represented as a sum of exponentially decaying functions:

[TABLE]

where all $\alpha^{D}_{l},\alpha^{F}_{l}$ constants are non-zero. It is important that the set of exponentially decaying functions in $\mathcal{D}$ is entirely included in the corresponding set of exponentially decaying functions in $\mathcal{F}$. This condition is automatically fulfilled if we assume that the spectral density is of the form $J(\omega)=\Theta(\omega)J^{\prime}(\omega)$ with $J^{\prime}(\omega)$ an odd function and with simple poles which do not lie on the real axis. Then by use of the residual theorem one can show that the set of poles in the upper half plane of $J^{\prime}(\omega)$ is equal to $\{-i\lambda_{l}\}_{l\in L}$ and the set of poles of $\coth(\omega/(2T))$ is equal to $\{-i\lambda_{l}\}_{l\in\tilde{L}}$.

The mapping of $i\mathcal{S}^{n}_{B,SB}$ given in (II.1) to a time-local action $i\mathcal{S}^{n,TL}_{B,SB}$ is defined as follows (we will denote the time arguments of the variables as subscripts):

[TABLE]

The integer $\mathcal{M}_{n}$ is equal or larger than the number of elements in $L\cup\tilde{L}$ and the $A^{n}$-matrix is chosen such that the set of its eigenvalues coincides with $\{\lambda_{l}\}_{l\in L\cup\tilde{L}}$ (some eigenvalues may appear more than once). The initial distribution function of the $\phi$-fields is given by:

[TABLE]

The covariance matrix $\Sigma^{n}$ is well defined since all eigenvalues of $A$ have a negative real part. The exact values of $\bm{v}^{n}$, $\bm{\varepsilon}^{n}$, $B^{n}$, $A^{n}$ are determined after integrating out the $\bm{x}^{n}$, $\bm{\phi}^{n}$, $\bm{\varphi}^{n}$ variables in (34) and comparing the result with (II.1). Since the number of free variables is larger than the number of equations that have to be fulfilled, we have some freedom in the choice of $\bm{v}^{n}$, $\bm{\varepsilon}^{n}$, $B^{n}$, $A^{n}$. All details about the proof of (34), (35) are given in the Appendix A, where we have assumed that the algebraric multiplicity of every eigenvalue of the $A^{n}$ matrix is equal to one. In this case the values of $\bm{v}^{n}$, $\bm{\varepsilon}^{n}$, $B^{n}$, $A^{n}$ are chosen such that the equation

[TABLE]

is fulfilled. In addition, the constants have to be chosen such that the part of the expression

[TABLE]

that depends on the difference $(t-t^{\prime})$ is equal to $\frac{1}{\pi}\mathcal{F}_{n}(t-t^{\prime})$. In the following we will also add the constraint, that $v^{n}_{m}=0$ if $b^{n}_{m}=0$, which is not a necessary condition for (34) to hold, but it will be needed to show (60) and (61).

The ansatz that we have used in (34), (35) origins from the idea that different Hamiltonian operators of the environment can produce the same effective action $i\mathcal{S}_{B,SB}$ after integrating out the environmental degrees of freedom. One of the best known examples is the replacement of the Hamiltonian of each of the $\mathcal{N}$ environments in (3) with the Hamiltonian of a single harmonic oscillator (nuclear mode) that is coupled to one of the sites of the system and also to a Markovian bath, which is also composed of non-interacting harmonic oscillators Garg et al. (1985). If we integrate out the bath degrees of freedom from the path integral expression of the problem but keep the nuclear modes, then the contribution of the environment plus the nuclear mode will have the same form as the right-hand side of (34). The $\phi$ ($\varphi$) variables will correspond to the real or imaginary part of the classical (quantum) variables of the nuclear mode. Similarly to the $\eta$ variables in (18), the $\varphi$ variables can be integrated out analytically by the use of (29) and the set of Dirac-delta functions will define a path for the time evolution of $\bm{\phi}^{n}_{t}$. The $x$-variables will be integrated by the use of Monte-Carlo methods, which effectively will make the equations for the $\phi$-variables stochastic. The noise in these equations origins from the Markovianity of the bath to which the nuclear mode was coupled. Since the dimension of the new $\bm{\varphi}^{n}_{t}$, $\bm{\phi}^{n}_{t}$, $\bm{x}^{n}_{t}$ in (35) can be larger than two we can just assume that each site of the system is coupled to more than one nuclear modes, where the nuclear modes can be coupled to the same or to different Markovian environments.

Finally, we have to mention that $\rho^{n}_{\Phi}(\bm{\phi})$ can be interpreted as the probability density function, where all nuclear modes are in equilibrium with the rest of the environment for the case that they are not coupled to the system. To demonstrate this we can decouple the nuclear modes from the system in the interval $[-t_{R},0]$ ($-t_{R}<0$) by setting $\bm{\varepsilon}^{n}=\bm{v}^{n}=0$ in (35) and let them evolve until they reach a steady state at $t=0$. The equation of motion for the $\phi^{n}_{m,\tau}$ variables which is derived after integrating out the $\varphi^{n}_{m,\tau}$ variables by using (29), is given by:

[TABLE]

The solution of this equation at $t=0$ is given by:

[TABLE]

where in the second line we have taken the limit $-t_{R}\rightarrow-\infty$. In the first (second) line the sums over $\tau$ is from $-t_{R}$ ($-\infty$) to [math] by taking steps of length ${\scriptstyle\Delta}t$. By interpreting every $x^{n}_{m,t}$ as a realization of a normally distributed random variable $X^{n}_{m,t}\sim\mathcal{N}(0,1)$ it follows that the steady state solution is also a random variable with zero mean and variance equal to

[TABLE]

where we have used that $\langle\bm{X}^{n}_{t}\bm{X}^{nT}_{t^{\prime}}\rangle=\delta_{t,t^{\prime}}\mathbbm{1}$. The last expression is equal to the definition of the covariance matrix of the multivariate probability density function $\rho^{n}_{\Phi}$ in (40), which proves our statement.

II.4 Truncated Wigner approximation

We will conclude this section by giving the equations of motion for the $\psi_{n,\tau}$, $\phi^{n}_{m,\tau}$-variables, which are obtained by applying the TWA. For the current problem this approximation is equivalent to neglecting the $\frac{1}{2}|\eta_{n}|^{2}$-terms in $n^{o}_{n}$, defined in (11), which appear in every $\mathcal{S}^{n,TL}_{B,SB}$ contribution to the action. We can integrate out the $\eta_{n}$, $\varphi^{n}_{m,\tau}$ terms by the use of (29). If we keep the variables $x^{n}_{m,\tau}$, $\phi^{n}_{m,0}$, $\psi_{n,0}$ fixed, the set of Dirac-delta functions defines a unique path for the time evolution of the $\psi_{n,\tau}$ and $\phi^{n}_{m,\tau}$ variables, which is given by:

[TABLE]

This is the Euler-Maruyama discretization of the following set of stochastic equations:

[TABLE]

where the set $\{W^{n}_{m}(t)\}$ ($n\in\{1\ldots\mathcal{N}\}$, $m\in\{1\ldots\mathcal{M}_{n}\}$) is a set of independent Wiener processes, which have the property $dW^{n}_{m}(t)dW^{n^{\prime}}_{m^{\prime}}(t)=dt\delta_{nn^{\prime}}\delta_{mm^{\prime}}$. To calculate the expectation value of ${\rm tr}[\hat{O}\hat{\rho}(t)]$ we have to sample the variables $x^{n}_{m,\tau}$ from a normal distribution $\mathcal{N}(0,1)$ and the variables $\phi^{n}_{m,0}$, $\psi_{n,0}$ from $\rho^{n}_{\Phi}$, $\rho^{abs}_{\mathcal{W}}$, defined in (40), (24), and calculate the mean value of (31).

III Stochastic equation with all quantum corrections

In this section we will analytically integrate all $\varphi$- and $\eta$-variables from (15) without neglecting any terms in the action. The idea is to make all terms in $i\mathcal{S}^{n}_{B,SB}$ ($n\in\{1\ldots\mathcal{N}\}$) to be linear in $(\varphi,\Re\eta,\Im\eta)$ and also completely imaginary, which will allow us to integrate them out via (29). To do this we will introduce first a set of probability density functions $f_{R|X}(r|x),f_{\Theta}(\theta)$:

[TABLE]

where $J_{0}(x)$ is the Bessel function of first kind. They are defined such that the following equation is true for every $n\in\{1,\ldots\mathcal{N}\}$:

[TABLE]

The sum in the second line of (61) is taken only over those $m\in\{1,\ldots\mathcal{M}_{n}\}$, where $b^{n}_{m}>0$. The new differentials and variables are defined as follows:

[TABLE]

The proof of (60) is given in the Appendix B. Here we will only point out that the right-hand side of (60) differs from the TWA, applied to the left-hand side of the same equation, only by the ${\scriptstyle\Delta}\chi$-term in $i\tilde{S}^{n,TL}_{B,SB}$. If we neglect the ${\scriptstyle\Delta}\chi$ term in the second line of (61), then we can integrate out the $\theta^{n}_{m,t}$, $r^{n}_{m,t}$ variables in the path integral on the right-hand side of (60) since nothing will depend on them. In addition, via the transformation $x^{n}_{m,\tau}\rightarrow-{\rm Sign}(v^{n}_{m})x^{n}_{m,\tau}$ at the right-hand side of (60), the $D[\bm{x}^{n}]$ term does not change, but the $\sqrt{{\scriptstyle\Delta}t}\tilde{B}^{n}\bm{x}^{n}_{\tau}$ term in (61) transforms to $\sqrt{{\scriptstyle\Delta}t}B^{n}\bm{x}^{n}_{\tau}$ which proves our statement.

We have to remind, that $f_{R|X}$ is a pseudo-probability density function, as shown in Fig. 1, and it can be replaced by the following expression:

[TABLE]

With the formula (29) we can integrate out all $\varphi^{n}_{m,\tau}$, $\Re\eta_{n,\tau}$, $\Im\eta_{n,\tau}$ variables from the path integral representation of the expectation value of the operator $\hat{O}$. If we keep the variables $x^{n}_{m,\tau}$, $r^{n}_{m,\tau}$, $\theta^{n}_{m,\tau}$, $\phi^{n}_{m,0}$, $\psi_{n,0}$ fixed, the set of Dirac-delta functions defines the following equation for the time evolution of $\bm{\psi}_{\tau}$, $\bm{\phi}^{n}_{\tau}$:

[TABLE]

where only the $n$-th element of the $\bm{\kappa}^{nm}$ vector is non-zero. The last sum in (72) is taken only over those $(n,m)$, where $b^{n}_{m}>0$. The equations (47), (48), obtained from the TWA, differ from (72), (73) only by the $B^{n}$ matrix, which is replaced by $\tilde{B}^{n}$ and by the absence of the term proportional to ${\scriptstyle\Delta}\chi$. To calculate the expectation value of the operator $\hat{O}$ we sample the variables $x^{n}_{m,\tau}$, $\theta^{n}_{m,\tau}$, $r^{n}_{m,\tau}$, $\phi^{n}_{m,0}$, $\psi_{n,0}$ from the probability density functions $f_{X}$, $f_{\Theta}$, $f^{abs}_{R|X}$, $\rho^{n}_{\Phi}$, $\rho^{abs}_{\mathcal{W}}$ given in (58), (55), (67), (40), (24) and calculate the mean value of:

[TABLE]

Since all normalisation constants $N(f_{R|X},x^{n}_{m,\tau})$ are always larger than one, it follows that their product will grow exponentially in time. This growth has to be compensated by the exponential decay of $\mathcal{O}_{\mathcal{W}}$ and/or by the alternating sign of the product of all sign-functions of every sampling of the random variables. In both cases this will require an exponential growth of the number of trajectories to obtain a good approximation of the expectation value. An additional problem comes from the ${\scriptstyle\Delta}\chi$-term in (72), because it is proportional to $({\scriptstyle\Delta}t)^{1/4}$ and this requires the use of very small time steps to obtain an accurate trajectory. Although being practically inapplicable, this new representation of the problem is a good starting point for the derivation of corrections that go beyond the truncated Wigner approximation.

IV Corrections of the TWA

To obtain an approximation that can be applied to the calculation of practical problems, we have to find a way to eliminate the pseudo-probability density function $f_{R|X}$ from the path integral representation of ${\rm tr}[\hat{O}\hat{\rho}(t)]$. To do this we take into account the fact, that we are only interested in the time evolution of the different elements of the density matrix. Since the system lies in the single exciton subspace, the density matrix and its elements are given by:

[TABLE]

The Weyl symbol of the $\hat{a}^{\dagger}_{n^{\prime}}\hat{a}_{n}$ operator is equal to:

[TABLE]

An equation of motion for $\psi_{n}(t)\psi^{*}_{n^{\prime}}(t)$ can be obtained if we take the equation for $\bm{\psi}$ from (72) and derive from it the equation for $\bm{\psi}^{*T}$. Then we take the product of both equations such that on the left-hand side of the equation we obtain the matrix $\bm{\psi}\bm{\psi}^{*T}$. In the new equation we replace all polynomials in ${\scriptstyle\Delta}\chi$ with their average over $\theta$ and $r$. Up to the fourth moment the only non-vanishing expectation values are

[TABLE]

The last result can be explained by the properties of the $f_{\Theta}$-function and the phase $\exp[i\theta^{n}_{m,\tau}]$, which leaves only those expectation values to be non-zero, where every ${\scriptstyle\Delta}\chi$ is multiplied with its complex conjugate. In this case the expression does not depend on the $\theta$-variable. After taking the average over all $r^{n}_{m,\tau}$, $\theta^{n}_{m,\tau}$ variables, the equation of motion for the matrix $\bm{\psi}\bm{\psi}^{*T}$ has the following form:

[TABLE]

If we take the limit ${\scriptstyle\Delta}t\rightarrow 0$ and neglect the terms proportional to ${\scriptstyle\Delta}t^{s}$ ($s>1$), the equation is equal to the Euler-Maruyama discretization of the following stochastic equation:

[TABLE]

where $\bm{W}^{n}(t)$ is defined in (52). A single realization of the observable is obtained by sampling $\phi^{n}_{m,0}$, $\psi_{n,0}$ from $\rho^{n}_{\Phi}$, $\rho^{abs}_{\mathcal{W}}$ and calculating the mean value of (31). An alternative way to derive (82) and the possibility to calculate multitime correlation functions within this approximation is given in the Appendix D.

If we neglect the second line of (82), we obtain again the TWA, since $(\bm{\psi}\bm{\psi}^{*T})(t)$ has a solution of the form

[TABLE]

which factorizes into a product of the solutions for $\bm{\psi}$ and $\bm{\psi}^{*}$ from (50).

If we use the solution of (83), which is equal to (120) after replacing $n^{o}(\tau)$ with $2(\bm{\psi}\bm{\psi}^{*T})_{nn}(\tau)$, in the definition of $\tilde{h}(t)=h-{\rm diag}[\bm{\varepsilon}^{1T}\bm{\phi}^{1}(t),\ldots\bm{\varepsilon}^{\mathcal{N}T}\bm{\phi}^{\mathcal{N}}(t)]$, then (82) transforms into a differential equation with a memory kernel. This memory kernel plays a crucial role for the effect of the new ${\scriptstyle\Delta}\chi$-dependent term on the evolution of the system. If the memory term is absent, then the noise generated from the ${\scriptstyle\Delta}\chi$-terms will not have any effect on the $\bm{\psi}\bm{\psi}^{*T}$ matrix on average.

In the end we mention that (82) preserves the trace of the $\bm{\psi}\bm{\psi}^{*T}$-matrix on average (and also the trace of $\rho=\langle\bm{\psi}\bm{\psi}^{\dagger}\rangle-\mathbbm{1}/2$). Since $\bm{\kappa}^{nm}\bm{\kappa}^{nm\dagger}$ is a diagonal matrix one can show by induction that

[TABLE]

which is zero on average.

V Examples

We consider the example of having a system of $\mathcal{N}=2$ sites, each of them being linearly coupled to an independent bath of harmonic oscillators. The corresponding spectral densities will be the same and are given by:

[TABLE]

where $\Theta$ is the Heaviside step function. Since $J^{\prime}(\omega)$ is odd, we can use the residual theorem to calculate the noise and dissipation kernels of the action $i\mathcal{S}^{n}_{B,SB}$ $(n=1,2)$:

[TABLE]

The $(C(\Omega t),S(\Omega t))$ functions are replaced with their upper/lower definition, when we use the upper/lower sign of $\pm,\mp$ in the definition of $J^{\prime}(\omega)$. If we use the lower definition of $(C(t),S(t))$, we have the additional constraint, that $\gamma>\Omega$. Depending on the situation we will use the $a^{\prime}_{1}$, $a^{\prime}_{2}$ instead of $a_{1}$, $a_{2}$. Both pairs of variables are given by:

[TABLE]

and are related to the reorganization energy as follows:

[TABLE]

For simplicity we approximate ${\coth}(\frac{\omega}{2T})\approx\frac{2T}{\omega}$ and neglect all sums over the $l$-index in (89), (93), (94), which is justified at high temperatures. In this case both functions $\mathcal{F}$ and $\mathcal{D}$ can be decomposed as a sum of two exponentially decaying functions with the exponents

[TABLE]

The time non-local part of the action $i\mathcal{S}^{n}_{B,SB}$, generated from every one of the reservoirs, can be mapped to a time local action in the same way as described in II.3. Additional information about the choice of $A^{n},\bm{\varepsilon}^{n},\bm{v}^{n},B^{n}$ used in (34),(35) can be found in the Appendix C.

We will assume that both spectral densities are equal and we will consider the two cases, where $a^{\prime}_{1}=0$ or $\Omega=0$. Both spectral densities now have the following form:

[TABLE]

All parameters of the system will be given in units of the difference between the energy levels of the donor and the acceptor $\Delta=h_{11}-h_{22}$. The coupling between both sites is set to $h_{12}/\Delta=0.4$ and the temperature of the reservoirs is set to $T/\Delta=2$, unless it is not mentioned explicitly. We assume, that the excitation is initially localized at the first site, which produces the following Wigner distribution function:

[TABLE]

The pseudo probability density function $\rho^{1}_{\mathcal{W}}$ can be represented as a product of $\rho^{1,abs}_{\mathcal{W}}$, ${\rm Sign}(\rho^{1}_{\mathcal{W}},\psi^{*}_{1},\psi_{1})$, $N(\rho^{1}_{\mathcal{W}})$ in the same way as explained in Sec. II.2.

Information about the quality of the proposed approximation can be obtained from Fig. 2,3 where we have compared our results with those from the TWA and with the results obtained by the use of the HEOM method which will be referred to as the exact results. In all figures we plot the population $P_{1}$ of the first site. To understand why our approximation works better in some regimes, we will look at the elements of ${\rm diag}[\bm{\varepsilon}^{1T}\bm{\phi}^{1}(t),\bm{\varepsilon}^{2T}\bm{\phi}^{2}(t)]$, that are contained in $\tilde{h}(t)$, defined in (49). We can replace $\bm{\phi}^{n}(t)$ with the solution of (83), which is obtained from (120) after replacing $n^{o}(\tau)$ with $2|\psi(\tau)|^{2}$. It follows that the diagonal matrix plays a role of a memory term, whose elements are given by ($n$-index is omitted):

[TABLE]

where $a^{\prime}$, $\mu$ and $g(t)$ are given by

[TABLE]

The second term of (108) is obtained by making use of the fact that the vectors $\bm{v},\bm{\varepsilon}$ are defined such that (42) is fulfilled. The equivalence $\mathcal{D}(t)/\pi=a^{\prime}g(t)$ can be seen by a direct substitution of (95), (88) into (107). The proof of $\bm{\varepsilon}^{T}e^{At}Bd\bm{W}(t)=(a^{\prime})^{1/2}g(t)\mu dW_{1}(t)$ is given in the Appendix C.

We are mainly interested in the contributions to (108) containing $g(t)$ since they are responsible for the impact of the $\bm{\kappa}\bm{\kappa}^{*T}$-dependent terms on the time evolution of the system. From the property $\int^{\infty}_{0}g(\tau)d\tau=1$ it follows, that for the two different spectral densities given in (102), the strength of the memory kernel $g(t)$ is the same but the way how the previous values of $\psi(t)$, $W(t)$ are taken into account is different.

The common feature of all cases, where the approximation is accurate at large time scales, is that the weight of the $g(t)$ function is uniformly distributed over large time scales and it does not change sign. The measure, that defines if a time interval is sufficiently large, is determined by the time scale $\tau_{S}=1/h_{max}$, where $h_{max}$ is the maximum absolute value of the matrix representation of the commutator $H^{\times}_{S}:X\mapsto[h,X]$. For our system Hamiltonian it follows that $\tau_{S}=1/\Delta$.

If we consider the case, where $\Omega=0$ for the spectral densities of both environments, then the condition, that $g(t)$ is uniformly distributed over a large time interval is equivalent to $\tau_{S}\ll 1/\gamma\Leftrightarrow\Delta\gg\gamma$. This can be seen in Fig. 2, where an increase of $\gamma/\Delta$ from $0.1$ to $1$ decreases the quality of the approximation at large time scales ($t\Delta=100$).

If we consider the case $a^{\prime}_{1}=0$ and $(C,S)=(\cos,\sin)$ for the spectral functions of both environments, then the condition for $g(t)$ is equivalent to $h_{max}>\gamma>\Omega$, where the second inequality comes from the restriction that $g(t)$ does not have to change sign over the $\tau_{S}$-timescale. The change of the quality of the approximation by changing $\gamma,\Omega$ of both spectral densities can be seen in Fig.3. In the first (second) column of the figure we see how an increase of $\gamma$ ($\Omega$) leads to a decrease of the quality of the approximation because the condition $h_{max}>\gamma$ ($\gamma>\Omega$) is violated.

The common feature of all figures is that our approximation reproduces more accurate results at large time scales than the TWA does, as long as the condition for the form of the memory kernel $g(t)$ is fulfilled and the system bath coupling $a^{\prime}$ is sufficiently small. At larger couplings as well as for parameters of the environment that violate the condition for the memory kernel the quality of the approximation decreases as it is shown in the second column of Fig. 2. Similarly to the TWA, the time range of validity of our approximation increases by an increase of the temperature.

To test the new approximation at low temperatures we have considered again the case $a^{\prime}_{1}=0$, $(C,S)=(\cos,\sin)$. The time evolution of the population of the donor site is plotted in Fig. 4 for the same six cases as those considered in Fig. 3 but with temperature decreased from $T/\Delta=2$ to $T/\Delta=0.2$. We see that the CTWA still behaves better than the TWA but, as expected, both approximations fail to describe the short time behavior of the system accurately. The main reason is that the relative weight between the temperature dependent noise kernel $\mathcal{F}(t)$ of the action that can always be taken exactly into account by the use of stochastic methods, and the temperature independent dissipation kernel $\mathcal{D}(t)$ whose effects are always approximated in the TWA and CTWA, decrease. In this case the noise produced by $\mathcal{F}(t)$ is not strong enough to overcompensate the effects from the dissipation kernel after applying the CTWA.

VI Summary

In this paper we have derived an alternative method to describe the exciton transport in open quantum systems. Instead of trying to derive an exact equation of motion for the reduced density matrix of the system, $\hat{\rho}_{S}$, we have used the fact that $\hat{\rho}_{S}$ lies completely in the single exciton subspace which allowed us to represent the time evolution of every non-zero element of $\hat{\rho}_{S}$ as the expectation value of a system operator, as shown in (76),(77). This difference led also to a change of the requirements for the form of the effective action of the system that is obtained after applying some stochastic unraveling approach to it. Instead of having to derive a time-local action at the cost of introducing integrals over new variables that can later be interpreted as Gaussian random variables, and trying to obtain a time-local stochastic equation of motion for $\hat{\rho}_{S}$ Stockburger and Grabert (2002); Imai et al. (2015); Stockburger and Mak (1999); Stockburger and Grabert (2001), we have derived an action that is completely real and linear in the quantum fields at the cost of introducing integrals over new variables that are weighted with pseudo-probability density functions. Even though in this case one can derive an exact equation of motion for the classical variables this method is numerically expensive. The rapid increase of the number of trajectories that are needed for the calculation of some observable origins from the new pseudo-probability density functions, which can modify the weight and the total sign of every contribution to the observable.

To overcome this problem we have derived an approximation of the equation of motion of the classical variables, where only some moments of the new variables weighted with the corresponding pseudo-probability density functions appear. We have shown that the quality of the approximation in the weak coupling limit depends on the form of the dissipation kernel of the Feynman-Vernon influence functional describing the effect of the environment on the system. If this kernel dissipates slowly in time without changing its sign, the corrections to the TWA give exact results even at large time scales ($t/\tau_{S}\gg 1$). This makes our approach applicable in situations, where the slow decay of the dissipation kernels does not allow the application of the Markovian approximations to the environment.

The method still remains simple enough to apply it to large systems. For a system composed of $\mathcal{N}$ sites, where each of them is linearly coupled to a separate environment, we can describe the time evolution of the system density matrix by a time local equation of an $\mathcal{N}\times\mathcal{N}$ matrix and $\mathcal{N}$ time local equations of vectors, where every vector describes the effect of one of the environments on the system. The dimension of each of those vectors is equal to the number of exponentially decaying functions, whose linear combination can describe the noise and dissipation kernel part of the Feynman-Vernon functional of the corresponding environment. In the limit of large networks the complexity of the approach depends mainly on the complexity of the time local equation of motion of the $\mathcal{N}\times\mathcal{N}$ matrix. It follows that the numerical effort to calculate a single trajectory of the observable is numerically as expensive as the calculation of some Markovian master equation in Lindblad form of the same system. So the increase of the numerical effort in comparison to the master equation depends on the number of trajectories of the observable which is needed to obtain a good approximation of its mean value. This number is of the order of $10^{4}-10^{5}$ and is required to sample correctly the initial state of the system which is described by a pseudo-probability density function. This is the price that we have to pay to take into account the non-Markovian effects of the environment on the system.

Acknowledgements.

We acknowledge support from the European Union (EU) through the Collaborative Project QuProCS (Grant Agreement No. 641277).

Appendix A Proof of the mapping of the action to a time local expression

We have to show that

[TABLE]

where $\rho$ and $\Sigma$ are defined in (40), (41) and the number of elements in $L\cup\tilde{L}$ is $\mathcal{M}$. To do this we integrate out all $\varphi$-fields by the use of (29). The equation of motion for the $\phi$-fields, defined by the set of $\delta$-functions, and its solution are given by:

[TABLE]

We replace $\bm{\phi}_{\tau}$ in the last term of the second line of (A) with the solution in (120). In the next step we integrate out the $x_{m,\tau}$ variables, which requires the replacement $\int d\bm{W}(\tau)=\sum_{\tau}\sqrt{{\scriptstyle\Delta}t}\bm{x}_{\tau}$. This gives the following identity

[TABLE]

where $m\in\{1,\ldots\mathcal{M}\}$ and the function $m(t,t^{\prime})$ is defined as:

[TABLE]

The second line of (A) transforms to:

[TABLE]

In the next step we have to decompose $m(t,t^{\prime})$ into a part that does and a part that does not depend on $(t-t^{\prime})$, where the first part has to be identified with the noise kernel of the action. To do this we choose the matrix $A$ such that the set of its eigenvalues coincides with $\{\lambda_{l}\}_{l\in L\cup\tilde{L}}$. We assume, that each of the eigenvalues of $A$ has an algebraric and geometric multiplicity of one. It follows, that there exists an invertible matrix $S$ and a diagonal matrix $D$ ($S,D\in\mathbb{R}^{\mathcal{M}\times\mathcal{M}}$) such that $S^{-1}DS=A$. By the use of the definitions

[TABLE]

we can show, that

[TABLE]

This allows us to decompose $m(t,t^{\prime})$ as a sum of the following functions:

[TABLE]

After integrating out the $\phi_{0}$-fields, the first line of (123) modifies to

[TABLE]

which cancels with the second term of (129). It follows, that (123) transforms to

[TABLE]

A direct comparison of the last equation with (A) and use of the fact, that the eigenvalues of $A$ are non-degenerate leads to the following equation for $\bm{v},\bm{\varepsilon},b_{1},\ldots b_{\mathcal{M}}$:

[TABLE]

Appendix B Proof of the linearization of the action in the quantum variables

To show that (60) is fulfilled, we have to integrate out the $x^{n}_{m,\tau}$ variables in the right-hand side of the equation. If we isolate only the $x$-terms and the terms that are non-linear in $\eta$ from this expression, we will obtain the following identity for every $n\in\{1,\ldots\mathcal{N}\}$:

[TABLE]

We will present a set of transformations, that have to be applied to every $(n,m,\tau)$ contribution to the last expression. For better legibility we will omit these indices. We introduce the initially unknown function $f(x,y,z)$, which has to fulfill the following equation:

[TABLE]

where the integrals over $x,y,z$ are from $-\infty$ to $\infty$. The function $f(x,y,z)$ can be obtained by an inverse Fourier transformation:

[TABLE]

We have to point out that $\alpha$ and $\beta$ are positive constants, which is a necessary condition to perform the following transformations. We can integrate out $\varphi$ and then apply the following variable transformations: $(\Re\eta,\Im\eta)\rightarrow(\Re\eta\beta^{-1/2}\alpha^{1/4},\Im\eta\beta^{-1/2}\alpha^{1/4})$, $(\Re\eta,\Im\eta)\rightarrow(\rho\cos(\varepsilon),\rho\sin(\varepsilon))$ with $\rho=|\eta|$. By the use of the relation $y\cos(\varepsilon)+z\sin(\varepsilon)=\sqrt{y^{2}+z^{2}}\sin(\varepsilon+\delta)$, where $\delta=\arcsin(y/\sqrt{y^{2}+z^{2}})$, we can integrate out the $\varepsilon$ variable. The function $f(x,y,z)$ is then equal to:

[TABLE]

We insert the result back in (135) and apply the following variable transformations: $(x,y,z)\rightarrow$ $(x\alpha^{1/2},y\beta^{1/2}\alpha^{-1/4},z\beta^{1/2}\alpha^{-1/4})$ and $(y,z)\rightarrow$ $(r\cos(\theta),r\sin(\theta))$. The right-hand side of (135) is then equal to:

[TABLE]

which completes the proof of (60). By comparing the left side of (135) and (138) we see, that a change of the sign of the $i\varphi x\alpha^{1/2}$ expression does not turn (138) into it its complex conjugate as it is the case for the left side of (135). This is related to the fact, that only specific combinations of random variables have a nonzero expectation value.

Appendix C Details of the application of the mapping procedure

To apply the mapping described in II.3, where the noise and dissipation kernels of $i\mathcal{S}^{n}_{B,SB}$ are given by (88),(89) and the sum over $l$ in (89),(93),(94) is neglected, we use the matrix ($n$-superscript will be omitted)

[TABLE]

which has the same eigenvalues $\lambda_{1,2}$ as in (99). The equations for the coefficients of the vectors $\bm{v},\bm{\varepsilon}\in\mathbb{R}^{2}$ and for the diagonal elements of the matrix $B={\rm diag}[b_{1},b_{2}]$ that were derived in Appendix A, are equivalent to the following set of equations:

[TABLE]

The first two equations can be derived from (42) by using

[TABLE]

The last two equations origin from the constraint that the part of $m(t,t^{\prime})$ in (43) depending on the difference of the two arguments, is equal to $\mathcal{F}(t-t^{\prime})$ which is given in (89). In the second line of (144), (145) we have used the assumption that we work at high temperatures and we have neglected all sums over $l$ in (89), (93), (94).

For the case $a^{\prime}_{a}=0=a_{1}$ we can solve the four equations in the high temperature limit by setting $\varepsilon_{1}=b_{2}=v_{2}=0$ and for the case $\Omega=0$ we can set $b_{2}=v_{2}=0$. A direct calculation shows that $\bm{\varepsilon}^{T}e^{At}Bd\bm{W}(t)=(a^{\prime})^{1/2}g(t)\mu dW_{1}(t)$ where $a^{\prime},g(t),\mu$ are defined in (111), (114), (117).

The initial distribution of the $\phi^{n}_{m,0}$-variables, described by $\rho^{n}_{\Phi}(\bm{\phi}^{n}_{0})$, is given in (40) where the covariance matrix $\Sigma^{n}$, defined in (41), takes the following form ($n$-index is omitted):

[TABLE]

Appendix D Alternative derivation of the CTWA and calculation of multitime correlation functions

For the following discussion it will be useful to include the terms

[TABLE]

in $i\mathcal{S}+\sum_{n}i\mathcal{S}^{n}_{B,SB}$ which will modify Eq. (72) by adding $-\bm{\nu}_{t}{\scriptstyle\Delta}t$ on its right side. We can use the formal solution of this equation

[TABLE]

and construct $\bm{\psi}(t)\bm{\psi}^{*T}(t)$ for the case $\bm{\nu}(t)=\bm{\nu}^{*T}(t)=0$. We can take the average over all $r^{n}_{m,\tau}$, $\theta^{n}_{m,\tau}$ variables ( denoted by $\langle\ldots\rangle_{r,\theta}$) and neglect the mixed terms $\bm{\psi}_{TWA}(t)\bm{\psi}^{*T}_{QC}(t)$, $\bm{\psi}_{QC}(t)\bm{\psi}^{*T}_{TWA}(t)$. Then we can apply the following self-consistent approximation:

[TABLE]

where we have used (79) to show the relation

[TABLE]

The $\tilde{h}$ matrix is defined in the same way as in (49) with the difference that the $\bm{\phi}$ terms are obtained from (51) by replacing $|\psi_{n,t}|^{2}$ term with the $nn$-component of $\langle\bm{\psi}(\tau)\bm{\psi}^{*T}(\tau)\rangle_{r,\theta}$ from (151). The differential of (151) is then equal to (82).

In the end we will briefly discuss the possibility to to use this approach to calculate multitime correlation functions. As example we will calculate the nonlinear response function $(t_{3}>t_{2}>t_{1})$

[TABLE]

where $\hat{a}^{\times}_{n_{j}}(t)\hskip 1.42262pt\bullet=[\hat{a}_{n_{j}}(t),\bullet]$ and \hat{a}^{{\color[rgb]{1,1,1}*}}_{n_{j}}(t) is an operator in the Heisenberg picture. To calculate (152) we can use the idea proposed in Polkovnikov (2010). We include again (149) in the definition of the action and choose the $\bm{\nu}(t)$, $\bm{\nu}^{*}(t)$ to be equal to

[TABLE]

and \eta^{*}_{n_{2},t_{2}}\eta^{*}_{n_{1},t_{1}}\eta^{{\color[rgb]{1,1,1}*}}_{n_{0},t_{0}} is replaced with

[TABLE]

where $F$ is defined such that:

[TABLE]

The $F$-function can again be represented as a product of a real probability density function, a normalization factor and a ${\rm Sign}$-function, and the integration over the \nu^{{\color[rgb]{1,1,1}*}}_{n_{j}}, $\nu^{*}_{n_{j}}$-variables can be carried in exactly the same way as the integration over over $\bm{\psi}_{0}$, $\bm{\psi}^{*}_{0}$. A direct expansion of $e^{i\mathcal{S}}$ in powers of $\nu_{j}$, $\nu^{*}_{j}$ shows that in the limit ${\scriptstyle\Delta}t\rightarrow 0$ the new and the old expression are equal. Our task now is to calculate the average of $\bm{\psi}_{t}$ over all $r^{n}_{m,\tau}$, $\theta^{n}_{m,\tau}$ variables. To do this we can use the formal solution (150) for $\bm{\psi}$ and neglect $\langle\bm{\psi}_{QC}(t)\rangle_{r,\theta}$ since $\langle d\chi^{n}_{m,\tau}\rangle_{r,\theta}=0$. The $|\psi_{n}(\tau)|^{2}$ term that appears in the diagonal elements of the $\tilde{h}(t)$-matrix can be replaced with the $nn$-term if $\langle\bm{\psi}(\tau)\bm{\psi}^{*T}(\tau)\rangle_{r,\theta}$ defined in (151). If we calculate the differential of this approximation we will obtain the following equation (the $\langle\ldots\rangle_{r,\theta}$ brackets are neglected):

[TABLE]

The only difference between the last equation and the corresponding TWA lies in the definition of the memory term in the diagonal elements of $\tilde{h}(t)$.