Extension of the Nakajima-Zwanzig approach to multitime correlation functions of open systems

TL;DR

This paper extends the Nakajima-Zwanzig approach to multitime correlation functions in open quantum systems, enabling combined short- and long-time simulations through memory kernels and stochastic unraveling.

Contribution

The authors develop a method to compute multitime correlation functions using memory kernels and stochastic unraveling, bridging short- and long-time dynamics in open quantum systems.

Findings

Validated approach with 2D-spectra simulations for donor-acceptor models

Compared results with reduced hierarchy equations of motion

Applicable to periodically driven two-level systems in equilibrium

Abstract

We extend the Nakajima-Zwanzig projection operator technique to the determination of multitime correlation functions of open quantum systems. The correlation functions are expressed in terms of certain multitime homogeneous and inhomogeneous memory kernels for which suitable equations of motion are derived. We show that under the condition of finite memory times these equations can be used to determine the memory kernels by employing an exact stochastic unraveling of the full system-environment dynamics. The approach thus allows to combine exact stochastic methods, feasible for short times, with long-time master equation simulations. The applicability of the method is demonstrated by numerical simulations of 2D-spectra for a donor-acceptor model, and by comparison of the results with those obtained from the reduced hierarchy equations of motion. We further show that the formalism is also applicable to the time evolution of a periodically driven two-level system initially in equilibrium with its environment.

Equations

$$\mathcal{L}={\mathcal{L}}_S+{\mathcal{L}}_B+{\mathcal{L}}_{SB}.$$

$$\displaystyle\hat{T}\Big{[}e^{\int^{t}_{t_{0}}ds(B(s)+C(s))}\Big{]}=\hskip 142.26378pt$$

$$\displaystyle=\hat{T}\Big{[}e^{\int^{t}_{t_{0}}dsB(s)}+\int^{t}_{t_{0}}d\tau e^{\int^{t}_{\tau}dsB(s)}C(\tau)e^{\int^{\tau}_{t_{0}}d\tilde{s}(B(\tilde{s})+C(\tilde{s}))}\Big{]}$$

$$\displaystyle=\hat{T}\Big{[}e^{\int^{t}_{t_{0}}dsB(s)}+\int^{t}_{t_{0}}d\tau e^{\int^{t}_{\tau}ds(B(s)+C(s))}C(\tau)e^{\int^{\tau}_{t_{0}}d\tilde{s}B(\tilde{s})}\Big{]},$$

$$\begin{array}[]{ll}Be^{\mathcal{L}Bt}=e^{B\mathcal{L}t}B,&B\in\{\mathcal{P},\mathcal{Q}\},\\ Be^{C\mathcal{L}t}=e^{\mathcal{L}Ct}B=B,&(B,C)\in\{(\mathcal{P},\mathcal{Q}),(\mathcal{Q},\mathcal{P})\}.\end{array}$$

$${\rho }_S(t)$$

$$U(t)$$

$$V(t)$$

$$\begin{array}[]{l}U(t)={\rm tr_{B}}\big{[}e^{\mathcal{L}\mathcal{P}t}R\big{]}\\ +\int^{t}_{0}d\tau\int^{\tau}_{0}d\tau^{\prime}{\rm tr_{B}}\big{[}e^{\mathcal{L}\mathcal{P}(t-\tau)}\mathcal{L}\mathcal{Q}e^{\mathcal{L}\mathcal{Q}(\tau-\tau^{\prime})}\mathcal{L}\mathcal{P}e^{\mathcal{L}\tau^{\prime}}R\big{]}.\end{array}$$

$$\begin{array}[]{rcl}{\rm tr_{B}}[\mathcal{L}\mathcal{P}X]&=&\mathcal{L}_{\bar{S}}{\rm tr_{B}}[X],\end{array}$$

$$\displaystyle{\rm tr_{B}}\big{[}e^{\mathcal{L}\mathcal{P}t}X\big{]}$$

$$u_{\bar{S}}(t)$$

$$U(t)=u_{\bar{S}}(t)+{\int }_0^{t}d\tau {\int }_0^{\tau }d{\tau }^{\prime }{u}_{\bar{S}}(t-\tau )K(\tau -{\tau }^{\prime })U({\tau }^{\prime }),$$

$$K(t)={\rm tr_{B}}\big{[}\mathcal{L}e^{Q\mathcal{L}t}Q\mathcal{L}R\big{]}={\rm tr_{B}}\big{[}\mathcal{L}_{SB}e^{Q\mathcal{L}t}Q\mathcal{L}_{SB}R\big{]}.$$

$$V(t)$$

$$I(t)$$

$${\rho }_S(t)$$

$$⟨A_N(t_{N\dots 1})\dots A_1(t_1)⟩$$

$$U_{(t_N,\dots ,t_1)}^{A_{N-1}\dots A_1}$$

$$V_{(t_N,\dots ,t_1)}^{A_{N-1}\dots A_1}$$

$$K_{(t_N,\dots ,t_1)}^{A_{N-1}\dots A_1}=$$

$$\displaystyle{\rm tr_{B}}\big{[}\mathcal{L}_{SB}e^{Q\mathcal{L}t_{N}}QA_{N-1}\ldots e^{Q\mathcal{L}t_{2}}QA_{1}e^{Q\mathcal{L}t_{1}}Q\mathcal{L}_{SB}R\big{]},$$

$$I_{(t_N,\dots ,t_1)}^{A_{N-1}\dots A_1}=$$

$$\displaystyle{\rm tr_{B}}\big{[}\mathcal{L}_{SB}e^{Q\mathcal{L}t_{N}}QA_{N-1}\ldots e^{Q\mathcal{L}t_{2}}QA_{1}e^{Q\mathcal{L}t_{1}}Q\rho_{0}\big{]}.$$

$$\mathcal{B}{e}^{\mathcal{L}t}\mathcal{Q}$$

$$\mathcal{Q}{e}^{\mathcal{L}t}\mathcal{P}$$

$$\mathcal{Q}{e}^{\mathcal{L}t}\mathcal{Q}$$

$$+{\int }_0^{t}d\tau {\int }_0^{\tau }d\tilde{\tau }\mathcal{Q}{e}^{\mathcal{L}\mathcal{Q}(t-\tau )}\mathcal{L}\mathcal{P}{e}^{\mathcal{L}(\tau -\tilde{\tau })}\mathcal{P}\mathcal{L}{e}^{\mathcal{Q}\mathcal{L}\tilde{\tau }}\mathcal{Q},$$

$$\mathcal{Q}{e}^{\mathcal{L}t}\mathcal{Q}{\rho }_0$$

$$+{\int }_0^{t}d\tau \mathcal{Q}{e}^{\mathcal{L}\mathcal{Q}(t-\tau )}\mathcal{L}\mathcal{P}{e}^{\mathcal{L}\tau }\mathcal{Q}{\rho }_0,$$

$$\displaystyle{\rm tr_{B}}\big{[}e^{\mathcal{L}t_{2}}A_{1}e^{\mathcal{L}t_{1}}\mathcal{B}\big{]},$$

$$\displaystyle{\rm tr_{B}}\big{[}e^{\mathcal{L}t_{3}}A_{2}e^{\mathcal{L}t_{2}}A_{1}e^{\mathcal{L}t_{1}}\mathcal{B}\big{]},$$

$$U_{(t_2,t_1)}^{A_1}$$

$$+U(t_2-{\tau }_2)K_{({\tau }_2,{\tau }_1)}^{A_1}U(t_1-{\tau }_1),$$

$$U_{(t_3,t_2,t_1)}^{A_2{A}_1}$$

$$+U(t_3)A_{2}U(t_2-{\tau }_2)K_{({\tau }_2,{\tau }_1)}^{A_1}U(t_1-{\tau }_1)$$

$$+U(t_3-{\tau }_3)K_{({\tau }_3,{\tau }_2)}^{A_2}U(t_2-{\tau }_2)A_{1}U(t_1)$$

$$+U(t_3-{\tau }_3)K_{({\tau }_3,{\tau }_2)}^{A_2}U(t_2-{\tau }_2-{\tau }_2^{\prime })K_{({\tau }_2^{\prime },{\tau }_1)}^{A_1}U(t_1-{\tau }_1)$$

$$+U(t_3-{\tau }_3)K_{({\tau }_3,t_2,{\tau }_1)}^{A_2{A}_1}U(t_1-{\tau }_1)$$

$$=U(t_3)A_2{U}_{(t_2,t_1)}^{A_1}$$

$$+U(t_3-{\tau }_3)K_{({\tau }_3,{\tau }_2)}^{A_2}{U}_{(t_2-{\tau }_2,t_1)}^{A_1}$$

$$+U(t_3-{\tau }_3)K_{({\tau }_3,t_2,{\tau }_1)}^{A_2{A}_1}U(t_1-{\tau }_1).$$

Full text

I Introduction

The simulation of the dynamics of an open quantum system coupled to an infinitely large environment TheWork is still a problem that attracts a lot of attention since there is a need for the development of reliable and fast numerical methods. Approaches based on Redfield- or Lindblad-like equations (see, e.g., Vulto_1999 ; Jang_2002 ; Jang_2004 ; Mohseni_2008 ; Castro_2008 ; Braig_2003 ; Mitra_2004 ; Koch_2005 ; Haenggi_2013 ), (self-consistent) perturbation expansions in some small parameter within the Keldysh formulation Aligia_2006 ; Monreal_2010 or projection operator techniques Nakajima ; Zwanzig ; Shibata do not cover the whole system parameter range of interest. These gaps can be filled by the use of numerically expensive approaches like time dependent density matrix numerical renormalization groups Daley_2004 ; White_2004 ; Schmitteckert_2004 , multilayer multi-configuration time-dependent Hartree method in second quantisation representation Wang_2003 ; Wang_2011 , real time quantum Monte Carlo methods Muehlbacher_2008 ; Muehlbacher_2004 and iterative path summation schemes Weiss_2008 ; Segal_2010 .

Although being exact, their computational cost increases exponentially in time, which often requires their combination with other methods in order to obtain the stationary state of the system. In Wilner_2013 ; Cohen_2011 the Nakajima-Zwanzig generalised quantum master equation is used to extract the specific memory kernels from the early time evolution of the system, which was initially obtained by the use of one of the exact approaches mentioned above. The memory kernels are then used for the calculation of the system dynamics for arbitrary long times.

The advances in nonlinear optical spectroscopy increased the need to develop efficient methods for calculating system multitime correlation functions.

One of the important examples is the two-dimensional (2D) spectroscopy of a photosynthetic pigment-protein complex known as Fenna-Mathew-Olsen (FMO) complex, which is obtained from the knowledge of two three-time correlation functions Engel_2007 .

In this work we extend the Nakajima-Zwanzig projection operator approach to the calculation of multitime correlation functions (MTCF), which requires the introduction of multitime homogeneous and inhomogeneous kernels. By having the information about the kernels in some finite time range we are able to calculate the MTCFs for an arbitrary set of times. We will see that the formalism can also be applied to problems which, at first glance, do not require the calculation of MTCFs, namely the time evolution of a periodically driven system being initially in equilibrium with its environment.

In order to calculate the required multitime kernels we construct a set of equations. If the problem can be solved efficiently by the hierarchy equations of motion (HEOM) method Tanimura_06 ; Tanimura_09 , then the input information for the equations can be obtained by slight modification of the method. Otherwise one can use a stochastic unravelling approach which is well suited for this task as long as the kernels decay to zero sufficiently fast.

The main advantage of this two-step approach is that it gives us the possibility to calculate MTCFs for problems that can not be described by HEOM. Even if the HEOM method is applicable it can be still more efficient to calculate the multitime kernels via the HEOM method and then the multitime propagators, than the direct calculation of the MTCFs.

The paper is structured as follows. In Sec. II.1 we derive briefly the solution of the Nakajima-Zwanzig equation. The generalisation of the problem to multitime correlation functions is presented in Sec. II.2, and in Sec. II.3 we derive the rules for constructing equations for the multitime kernels. In Sec. II.4 we apply the formalism to a periodically driven system being initially in equilibrium with its environment. The input information needed for the solution of the equations for the multitime kernels is calculated by use of a stochastic unravelling approach, which is presented in Sec. II.5. The reliability of the method is tested in Sec. III. In Sec.  III.1 the numerical results for the problem proposed in Sec.  II.4 are presented, and in Sec. III.2 the 2D-spectra of a donor-acceptor model is calculated. In both cases the results are compared with those obtained from the HEOM approach. Finally, conclusions about the advantages and drawbacks of the method are given in Sec. IV.

II Theory

II.1 Nakajima-Zwanzig equation

We consider an open system $S$ coupled to some bath $B$. The total Hilbert space is $\mathcal{H}=\mathcal{H}_{S}{}\otimes\mathcal{H}_{B}$. The Liouvillian operators, that describe the system, the bath and the system-bath coupling, are denoted by $\mathcal{L}_{S}$, $\mathcal{L}_{B}$ and $\mathcal{L}_{SB}$, respectively. The Liouvillian for the total system is thus given by

[TABLE]

For simplicity we assume that all operators do not depend on time, but the results can be extended for time dependent Liouvillians as well. We also define the projection operators $\mathcal{P}={\rm tr_{B}}[\ldots]\otimes R$ and $\mathcal{Q}=1-\mathcal{P}$, where ${\rm tr_{B}}[\ldots]$ denotes the trace over the bath degrees of freedom and $R$ is an arbitrary density operator for the bath with the properties ${\rm tr_{B}}[R]=1$ and $\mathcal{L}_{B}R=0$.

In the following we will often use the identity

[TABLE]

where $B(t)$ and $C(t)$ are any time dependent superoperators and $\hat{T}$ is the time ordering operator. The pair $(B,C)$ will be replaced by $(\mathcal{P}\mathcal{L},\mathcal{Q}\mathcal{L})$, $(\mathcal{L}\mathcal{P},\mathcal{L}\mathcal{Q})$ or $(\mathcal{Q}\mathcal{L},\mathcal{P}\mathcal{L})$, $(\mathcal{L}\mathcal{Q},\mathcal{L}\mathcal{P})$. The last set of identities that we will need is

[TABLE]

The density matrix of the open system $\rho_{S}(t)$ is given by

[TABLE]

where the initial state of the total system is denoted by $\rho_{0}$ and $U$, $V$ represent the homogeneous and the inhomogeneous propagators, respectively:

[TABLE]

By applying Eq. (2a) for $(B,C)=(\mathcal{L}\mathcal{P},\mathcal{L}\mathcal{Q})$ to U(t)={\rm tr_{B}}\big{[}e^{\mathcal{L}t}\mathcal{P}R\big{]}, substituting $\mathcal{P}RX=RX$ for every operator $X$ acting on $\mathcal{H}_{S}$ and then using Eq. (16b) we obtain

[TABLE]

Since for every operator $X$ acting on $\mathcal{H}$ we have

[TABLE]

where $\mathcal{L}_{\bar{S}}\equiv\mathcal{L}_{S}+\langle\mathcal{L}_{SB}\rangle$ with $\langle\mathcal{L}_{SB}\rangle\equiv{\rm tr_{B}}[\mathcal{L}_{SB}R]$, it follows that

[TABLE]

Equation (6) can then be rewritten as

[TABLE]

where the memory kernel is given by

[TABLE]

The last expression of the previous equation is obtained by use of $\mathcal{L}_{B}R=0$ and of the relation ${\rm tr_{B}}[\mathcal{Q}X]=0$ which holds for all operators $X$ acting on $\mathcal{H}$.

By the use of Eq. (2b) for $(B,C)=(\mathcal{Q}\mathcal{L},\mathcal{P}\mathcal{L})$ and ${\rm tr_{B}}[e^{Q\mathcal{L}t}\mathcal{Q}X]=0$ the equation for the inhomogeneous propagator [Eq. (5b)] becomes

[TABLE]

The kernel $I(t)$ is also known as inhomogeneity. Combining Eqs. (4) and (11a) and then using Eq. (9) we obtain the solution of the Nakajima-Zwanzig equation:

[TABLE]

At first glance it seems unnecessary to express $\rho_{S}{\scriptstyle}$ in terms of $U$ and $V$ as it is done in Eq. (4), because for the time evolution of the system density matrix it is sufficient to solve only Eq. (II.1). The alternative form in Eq. (4) is preferred if we are interested in the calculation of multitime correlation functions.

II.2 Multitime correlation functions

The multitime correlation function of an arbitrary set of system operators $\{A_{j}\}_{j\in\mathbb{N}}$ applied at times $t_{j\ldots 1}\equiv t_{j}+\ldots+t_{1}$ is given by

[TABLE]

where we have introduced the $N$-time homogeneous and inhomogeneous propagators:

[TABLE]

Every $N$-time propagator can be expressed as a function of $n$-time kernels with $n\leq N$ and of the propagators $U$ and $V$. The multitime homogeneous $K^{A_{N-1}\ldots A_{1}}$ and inhomogeneous $I^{A_{N-1}\ldots A_{1}}$ kernels are defined as

[TABLE]

The procedure to obtain the desired expressions can be entirely summarised in applying the following reduction rules:

[TABLE]

where $\mathcal{B}\in\{\mathcal{P},{\rm tr_{B}}\}$. In addition, one always has to decompose the system operators as $A=\mathcal{P}A\mathcal{P}+\mathcal{Q}A\mathcal{Q}$.

We will apply these rules for the two- and three-time propagators. Starting from

[TABLE]

and applying the reduction rules for $\mathcal{B}=\mathcal{P}R$ we obtain the following equations:

[TABLE]

Here and in all other equations we integrate over all $\tau_{j}$ variables from [math] to $t_{j}$ and over all $\tau^{\prime}_{j}$ variables from [math] to $t_{j}-\tau_{j}$ if the integration range is not shown explicitly.

The desired equations for the two- and three-time inhomogeneous propagators are obtained by applying the reduction rules to Eq. (17) for $\mathcal{B}=\mathcal{Q}\rho_{0}$:

[TABLE]

A diagrammatic representation of the terms contributing to Eqs. (18) and (19) is shown in Figs. 1 and 2. From Fig. 1 we can see that $U^{A_{1}}$ ($U^{A_{2}A_{1}}$) are obtained by constructing all possible combinations between $U$, $A_{1}$, $K^{A_{1}}$ (and $A_{2}$, $K^{A_{2}}$, $K^{A_{2}A_{1}}$), such that $A_{1}$ (and $A_{2}$) appear only once in every combination. In addition, the sum of the first two diagrams and the next two diagrams in Fig. 1b give the first and the second term of Eq. (18b). Also the diagrams contributing to $V^{A_{1}},V^{A_{2}A_{1}}$ can be obtained by taking all diagrams from Fig. 1a, 1b respectively and replacing the last homogeneous propagator by an inhomogeneous one or by replacing the homogeneous kernels containing $A_{1}$-superscript and the homogeneous propagator on their right side with the corresponding inhomogeneous kernels. The sum of the first three diagrams and the next three diagrams in Fig. 2b give the first and the second term of Eq. (19b).

II.3 Equations for the multitime kernels

We will consider the case of having time-independent Liouvillians but the results can be easily generalised to the time-dependent case. Looking at the rules for constructing diagrams we can conclude that every $N$-time homogeneous propagator ($N>1$) contains an $N$-time homogeneous kernel, the first and last arguments of which are convoluted with $U$. We can always derive an equation for $K^{A_{N-1}\ldots A_{1}}$ by applying a time derivative operator to the first and the last argument of $U^{A_{N-1}\ldots A_{1}}$. For the $N=1$ case we just have to take the second time derivative of $U$. A closer look at the time derivatives of a multitime propagators will allow us to cancel a significant amount of terms, which will make the resulting equation numerically more stable. First, we define the following system operators:

[TABLE]

where $\mathcal{L}_{\bar{S}\bar{B}}\equiv\mathcal{L}_{SB}-\langle\mathcal{L}_{SB}\rangle$ and $\mathcal{L}_{\bar{S}}\equiv\mathcal{L}_{S}+\langle\mathcal{L}_{SB}\rangle$. The $i,f$ subscript shows that we have applied $\mathcal{L}_{\bar{S}\bar{B}}$ at the beginning/end of the expression before taking the trace.

By use of the fact that

[TABLE]

we can see that after applying the time derivative w.r.t. $t_{N}$ at both sides of the equation for $U^{A_{N-1}\ldots A_{1}}$, all terms proportional to $\mathcal{L}_{\bar{S}}$ cancel out, such that only $U^{A_{N-1}\ldots A_{1}}_{f}$ remains on the left-hand-side of the equation. The same argument is valid also for the case of applying $\partial_{t_{1}}$ on both sides of the new equation. For the two-time homogeneous propagator $U^{A_{1}}$ we obtain for example the following set of equations:

[TABLE]

For the $N=1$ case we derive an equation similar to the Volterra equation of the second kind for $K$, that was defined in Zhang . The derivation of an equation for the mutitime inhomogeneous kernels $I^{A_{N-1}\ldots A_{1}}$ can be carried out similarly, the only difference being that only $\partial_{t_{N}}$ has to be applied on both sides of the equation for $V^{A_{N-1}\ldots A_{1}}$ since it contains always an $N$-time inhomogeneous kernel convoluted on the left side with $U$.

II.4 Periodically driven systems

Consider a system which initially is in a unique steady state with its environment. At $t=0$ a periodic force with period $t_{P}$ is turned on, that is applied at the system. We can simulate numerically this problem by starting from an arbitrary product state $\tilde{\rho}_{S}\otimes R$, letting it evolve in time until it relaxes to its unique steady state and then turning on the periodic force. We denote the Liouvillians describing the total system before and after turning on the periodic force by $\mathcal{L}^{\prime}$ and $\mathcal{L}(t)$, respectively. The state of the system is then given by:

[TABLE]

The time the system, initially prepared in the state $\tilde{\rho}=\tilde{\rho}_{S}\otimes R$, needs to relax to equilibrium will be denoted by $t_{R}$. If we express the unit operator $\mathbb{1}$ as a sum of $\mathcal{P}$ and $\mathcal{Q}$, and apply to the last equation the reduction rules given in Eq. (16), then we obtain the following result:

[TABLE]

where we have defined

[TABLE]

with $t_{32}\equiv t_{3}+t_{2}$. Equation (26) is just the extension of Eq. (5a) to time dependent Liouvillians. In order to calculate $\tilde{U}$ we also have to use a similar extension of Eq. (10):

[TABLE]

The propagators $U$, $\tilde{U}$ describe the time evolution of a system described by $\mathcal{L}^{\prime},\mathcal{L}(t)$, which is initially in the product state $\rho_{S}\otimes R$ with arbitrary $\rho_{S}$. The second argument of $\tilde{U}(t_{2},t_{1})$ gives the initial phase of the function of the periodic force, while the difference $t_{2}-t_{1}$ is the actual evolution time. It follows that we need to know $\tilde{U}(t_{2},t_{1})$ only in the range $t_{1}\in[0,t_{P})$. This property, also valid for $\tilde{K}(t_{2},t_{1})$, can be seen directly from Eqs. (26) and (29) and formally reads

[TABLE]

where $m\in\mathbb{Z}$ and $\mathcal{B}\in\{\tilde{U},\tilde{K}\}$. This fact is important since it allows us to calculate $\tilde{U}(t_{2},t_{1})$ by knowing $\tilde{K}(t_{2},t_{1})$ only in the range $t_{1}\in[0,t_{P})$, $t_{2}\in[t_{1},f(t_{1})]$, where $f(t_{1})$ is chosen such that $\tilde{K}(t_{2},t_{1})=0$ for $t_{2}>f(t_{1})$.

In order to better understand how we have to choose the relaxation time $t_{R}$, and to explain why this is not the time the system density matrix needs to reach its steady state, we split $\rho_{S}(t)$ in terms of homogeneous and inhomogeneous propagators as it is done in Eq. (4), and solve Eq. (24) for the inhomogeneous kernel $\tilde{I}$, which is defined as

[TABLE]

This definition is just the extension of Eq.(11b) to time dependent Liouvillians. The result is

[TABLE]

The constraint that the system was initially in its steady state means that $\tilde{I}(t,0)$ is independent of $t_{R}$ for all $t\geq 0$. If we define by $t^{\prime}$ the time that the propagator $U$ needs to reach its steady state and by $t^{\prime\prime}$ the time, where $\mathcal{K}(t,0,\tau_{1})=0$ $\forall\tau_{1}>t^{\prime\prime}$, then any $t_{R}>t^{\prime}+t^{\prime\prime}$ gives the same result in Eq. (32). We are free to set $t_{R}\rightarrow\infty$ and replace $U(t_{R}-\tau_{1})\tilde{\rho}_{S}$ by the steady state of the system $\rho_{S}(0)$, which is not driven by a periodic force. Thus, we obtain the following result:

[TABLE]

Finally, we mention that the argumentation of the previous subsection can also be applied to Eq. (25) in order to obtain an equation for $\mathcal{K}(t_{3},0,t_{1})$.

II.5 Stochastic unravelling method

In the following we always assume that the system is coupled to Gaussian environments such that the coupling is linear in the environmental fields. This allows us to integrate out analytically the reservoir degrees of freedom. All (multitime) propagators will then contain the same time nonlocal contribution $u_{NL}$ of the form:

[TABLE]

where $\Upsilon$ refers to the system part of the system-bath coupling operators and the index $j$ denotes the different environments to which the system is coupled. We have introduced the superoperator notation $f^{\times}A\equiv[f,A]$ and $f^{o}A\equiv\{f,A\}$. The dissipation and noise kernels $g_{D,j}$ and $g_{N,j}$ contributing to $u_{NL}$ are defined as:

[TABLE]

where for the spectral densities $J_{j}(\omega)$ we have to use the definition given in Eq. (47) or Eq. (57) depending on the problem we are interested in. In the following we will apply the stochastic unravelling method to the case of having a single element in the sum over $j$ and the index $j$ will be omitted.

Our goal is to make the action local in time. The first step to achieve this is to eliminate the $\theta$-function in Eq. (34). Since $g_{N}(t)$ is symmetric in $t$, we can replace $(\theta g_{N})(t)$ with $\frac{1}{2}g_{N}(t)$. The elimination of $\theta$ from $g_{D}$ requires the introduction of the following Fourier transformation:

[TABLE]

where the improper integral over $\varepsilon$ is calculated by the Cauchy principal value method. If $J(\omega)$ and $f(\omega)$ go to zero for large values of $\omega$ we can discretize the integrals over $\omega$ in Eqs. (36) and (37) by a finite sum of terms. Then we can make the $u_{NL}$ local in time at the cost of introducing a finite number of Gaussian integrals and Eq. (34) transforms to:

[TABLE]

The functions $\xi,\chi,\tilde{\chi}$ are given by:

[TABLE]

with $\omega_{\nu}=2\pi\nu\hskip 0.28453pt\Delta\nu$, $a_{\nu}=2\pi\hskip 0.28453pt\Delta\nu$ for $\nu>0$, $a_{0}=\pi\hskip 0.28453pt\Delta\nu$ and $\nu_{c},\tilde{\nu}_{c}$ are properly chosen cutoffs of $J(\omega)$ and $f(\omega)$, respectively. By the use of Monte-Carlo integration techniques for the calculation of the Gaussian integrals we can interpret $x_{\nu},\tilde{x}_{\nu},y_{\nu},\tilde{y}_{\nu}$ as normally distributed random variables and $\xi(t),\chi(t),\tilde{\chi}(t)$ as Gaussian random variables with zero mean and the following correlation functions:

[TABLE]

All other correlations are equal to zero. This stochastic unravelling of $u_{NL}$ is just a specific realisation of the general idea explained in Stockburger_02 .

The calculation of an arbitrary multitime propagator $U^{A_{N-1}\ldots A_{1}}$ is carried out by replacing $u_{NL}$ by the last two lines of Eq. (II.5) in the definition of $U^{A_{N-1}\ldots A_{1}}$ and averaging the result over a large enough number of realisations of the normally distributed random variables $\{x_{\nu},\tilde{x}_{\nu},y_{\nu},\tilde{y}_{\nu}\}$. $U^{A_{N-1}\ldots A_{1}}_{m}$ ($m\in\{f,i,fi\}$) is calculated by multiplying every realisation of $U^{A_{N-1}\ldots A_{1}}$ on the left or/and on the right by $F(t_{N\ldots 1})$ and $F(0)$, respectively, where $F$ is given by

[TABLE]

III Results

III.1 Periodically driven system initially in equilibrium with its environment

We apply the formalism derived in Sec. II to the problem of a classically driven two-level system coupled to a bosonic environment Leggett1987 . The system is described by the following Hamiltonian:

[TABLE]

where $b^{\dagger}_{k},b_{k}$ are the bosonic creation and annihilation operators for the modes $k$ with frequency $\omega_{k}$, and $\lambda_{k}$ describes the strength of the interaction of the two-level system with its environment. The Pauli spin-$\frac{1}{2}$ operators are denoted by $\sigma_{j}$ $(j\in\{x,y,z\})$, the energy distance between the two levels of the spin is $\Delta$ and the classical driving force applied to the system is given by $\varepsilon(t)=V_{0}\sin(\Omega t)$. We set $R=\rho^{eq}_{B}=\frac{1}{N}\prod_{k}\exp[-\beta\omega_{k}b^{\dagger}_{k}b_{k}]$ with $\beta=1/T$ ($k_{B}=1=\hbar$), which allows us to describe the effect of the environment on the two-level system completely by the spectral density:

[TABLE]

This form of the spectral density together with the replacement of $\coth(x)$ by $\frac{1}{x}+\frac{2x}{x^{2}+\pi^{2}}+\frac{2x}{x^{2}+4\pi^{2}}$ allow us to compare our results with the HEOM method since the dissipation and noise kernels $g_{D}$ and $g_{N}$, defined in Eqs. (35) and (36), can be expressed by a finite sum of exponentially decaying functions.

We choose the following parameters (measured in units of $\Delta$): $J=0.7$, $V_{0}=0.5$, $\lambda=0.05$, $\gamma=1.3$, $\omega_{u}=0.7$, $\Omega=\pi/2$ and $T=0.35$. By use of the stochastic unravelling method we simulate $U_{l}(t)$ ($l\in\{i,fi\}$) in the range $t\in[0,12]$, $\tilde{U}_{l}(t_{2},t_{1})$ ($l\in\{i,f,fi\}$) in the range $t_{1}\in[0,t_{P})$, $t_{2}\in[t_{1},t_{1}+12]$ and $W_{fi}(t_{3},0,-t_{1})$ for all $(t_{3},t_{1})$ which fulfil the constraint $t_{1}+t_{3}\leq 12$. From the equations for the kernels

[TABLE]

we obtain $K(t),\tilde{K}(t_{2},t_{1}),\mathcal{K}(t_{3},0,t_{1})$ for the same range of times as $U_{l}(t),\tilde{U}_{l}(t_{2},t_{1}),W_{fi}(t_{3},0,-t_{1})$ respectively ($l\in\{i,f,fi\}$). Since $K(t),\mathcal{K}(t_{3},0,t_{1})$ are equal to zero outside this range and $\tilde{K}(t_{2},t_{1})=0$ for $t_{1}\in[0,t_{P}),t_{2}>t_{1}+12$, and the periodicity condition in Eq. (30), we are able to calculate $U(t),\tilde{U}(t_{2},t_{1}),W(t_{3},0,-t_{1})$ for arbitrary $t$, $(t_{2},t_{1})$ and $(t_{3},t_{1})$. The equations for $U,\tilde{U},W$ are obtained by the use of the reduction rules (16). Those for $U$ and $W$ are given explicitly in (9) and (25) and the equation for $\tilde{U}$ is given by

[TABLE]

In the following we denote the up- and down-state of the two-level system by $|1\rangle$ and $|0\rangle$. This means that $\sigma_{z}=|1\rangle\langle 1|-|0\rangle\langle 0|$, $\sigma_{x}=|1\rangle\langle 0|+|0\rangle\langle 1|$. An element of some system superoperator $M$ will be denoted by

[TABLE]

If we plot all elements of $U\scriptsize{(t)}$ we will see that they become constant for $t>200$. In addition, the elements of $U$ in its steady state obey the following relations:

[TABLE]

and all other elements of $U$ are equal to zero. This assures that for every initial system density matrix $\tilde{\rho}_{S}$ the final steady state $U(t)\tilde{\rho}_{S}$ $(t>200)$ is the same. Taking into account that $\mathcal{K}(t_{3},0,t_{1})=0$ for $t_{1}>12$ we set $t_{R}=240$. In Fig. 3 we have plotted the time evolution of $U_{00,00}(t)$, which represents the occupation of the lower energy site given that the system was initially in $\rho_{S}(0)=|0\rangle\langle 0|\otimes\rho^{eq}_{B}$. The difference between $U_{00,00}(t)$ and the exact solution (obtained by the use of the HEOM approach) origins mainly from the large time step $(0.04)$ used in the calculation of $K$ and $U$.

In Fig. 4 we can see the time evolution of $\tilde{U}_{00,00}(t+t_{1},t_{1})$ for $t_{1}=0$ and $t_{1}=2$ which corresponds to the case of having a driving force of the form $\varepsilon(t)=\sin(\Omega t)$ and $\varepsilon(t)=\sin(\Omega t+\pi)$ respectively. For large enough $t$ $\tilde{U}(t+t_{1},t_{1})$ satisfies Eq. (53),(54) and also has the property that $\tilde{U}(t_{2},t_{1})=\tilde{U}(t_{2},0)$. This means that the steady state of a system being initially in $\rho_{S}(0)\otimes\rho^{eq}_{B}$ (and being evolved with $H$ given in Eq. (46) with $V_{0}=0.5$) does not depend on $\rho_{S}(0)$ but only on the initial phase of the driving force.

The time evolution of $W_{00,00}(t,0,-t_{R})$ is given in Fig. 5. From the inset in the figure we can see that the error increases by an order of magnitude in comparison to the previous two cases. The growth of the error origins from the second line of Eq. (25), where $\mathcal{K}$ is convoluted with the functions $U,\tilde{U}$, which , as shown in Fig. 3,4, deviate from the exact result. Even in this case the relative error remains below $2\%$ at short time scales and below $1\%$ at long time scales.

III.2 2D-spectra of a donor-acceptor model

As a second example we calculate the 2D-spectra of a system composed of a single donor and acceptor, each of them coupled to a different phononic bath. We reduce the description to the zero- and single exciton manifold. The total Hamiltonian is given by:

[TABLE]

where $\Delta_{1},\Delta_{2}$ are the energy levels of the states where only the donor $|1\rangle$ or the acceptor $|2\rangle$ are excited. The energy level of the ground state $|0\rangle$ is set to zero. The reorganisation energies $\Delta_{1,r},\Delta_{2,r}$ are defined as

[TABLE]

The second line of Eq. (55) describes the reservoir and the system-reservoir interaction in a similar way as it is done in Eq. (46). The spectral densities describing the effect of both environments on the system are given by:

[TABLE]

This means that we have already assumed that the initial state of the system is of the form $\rho_{S}(0)\otimes\rho^{eq}_{B}$, where $\rho^{eq}_{B}=\frac{1}{\mathcal{N}}\prod_{j}\prod_{k}\exp[-\beta\omega_{jk}b^{\dagger}_{jk}b_{jk}]$ with normalization constant $\mathcal{N}$.

The 2D-spectra $I(\Omega_{3},t_{2},\Omega_{1})$ is defined as a double Fourier transform of the rephasing $R_{rp}$ and nonrephasing $R_{nr}$ contributions to the third-order optical response function Mukamel_1995 :

[TABLE]

The operators $\mu_{L},\mu_{R}$ are contributions to the total dipole operator $\mu=\mu_{L}+\mu_{R}$, where $\mu_{L}=\mu_{1}|0\rangle\langle 1|+\mu_{2}|0\rangle\langle 2|$ and $\mu_{R}=\mu_{1}|1\rangle\langle 0|+\mu_{2}|2\rangle\langle 0|$. The system is initially in its ground state $\rho_{S}{\scriptsize(0)}=|0\rangle\langle 0|$.

In order to calculate $U^{\mu^{\times}_{R}\mu^{\times}_{R}}$ and $U^{\mu^{\times}_{R}\mu^{\times}_{L}}$ we use Eq. (18b). The terms on the right-hand side contain at most two time integrals, which substantially simplifies the numerical simulations. From Eqs. (9), (18a), (18b) we can derive equations for $K$, ($K^{\mu^{\times}_{L}}$, $K^{\mu^{\times}_{R}}$) and ($K^{\mu^{\times}_{R}\mu^{\times}_{R}}$, $K^{\mu^{\times}_{R}\mu^{\times}_{L}}$) respectively. They also contain only terms with at most two time integrals. The equations are solved for a finite time interval defined by the parameters $\tilde{t},\tilde{t}_{\mu},\tilde{t}_{\mu\mu}$ as follows: $K(t)$ is calculated for $t<\tilde{t}$, $K^{\mu^{\times}_{L}}_{(t_{2},t_{1})}$ and $K^{\mu^{\times}_{R}}_{(t_{2},t_{1})}$ - for $(t_{2},t_{1})$, which fulfil the condition $t_{2}+t_{1}<\tilde{t}_{\mu}$, $K^{\mu^{\times}_{R}\mu^{\times}_{R}}_{(t_{3},t_{2},t_{1})}$, $K^{\mu^{\times}_{R}\mu^{\times}_{L}}_{(t_{3},t_{2},t_{1})}$ - for $(t_{3},t_{2},t_{1})$ which fulfil the condition $t_{1}+t_{2}+t_{3}<\tilde{t}_{\mu\mu}$. The parameters $\tilde{t},\tilde{t}_{\mu},\tilde{t}_{\mu\mu}$ are chosen such that the kernels are zero outside the corresponding range. After solving the equations for the kernels we calculate $U$, which is used in the calculation of $U^{\mu^{\times}_{R}}$, $U^{\mu^{\times}_{L}}$. Finally, the three propagators are used in the calculation of $U^{\mu^{\times}_{R}\mu^{\times}_{R}}$, $U^{\mu^{\times}_{R}\mu^{\times}_{L}}$.

We work with the following parameters (measured in units of $\Delta\equiv\Delta_{2}+\Delta_{2,r}$): $\Delta_{1}=1.9,\Delta_{2}=0.9,J=0.1,\lambda_{1}=\lambda_{2}=0.1,\omega_{1}=\omega_{2}=1.4,\gamma_{1}=\gamma_{2}=2.6$. The reorganisation energies are $\Delta_{1,r}=\Delta_{2,r}=0.1$ and the memory kernels are calculated in the range defined by $\tilde{t}=5$, $\tilde{t}_{\mu}=\tilde{t}_{\mu\mu}=4$. The 2D-spectra of the system are plotted for four different waiting times $t_{2}$ in Fig. 6. Since the plots obtained by a direct use of the HEOM look the same as those of Fig. 6, we compare the time evolution of the off-diagonal peaks in Fig. 7 in order to obtain information about the quantitative accuracy of the method. The relative error is below $3\%$, which in this case is accurate enough to reproduce all main features of the 2D-spectra.

In general, the error depends on the number of time arguments of the multitime propagator, from which the observable is derived. The higher this number is, the more integro-differential equations for the kernels have to be solved. Since the calculation of an $n$-time propagator requires not only the knowledge of the kernels but also of a set of $m$-time propagators ($m<n$), the error accumulates.

IV Concluding remarks

In this paper we have extended the Nakajima-Zwanzig projection operator technique to the calculation of multitime correlation functions, which required the introduction of multitime kernels. The applicability of the theory was demonstrated by simulating the time evolution of a driven two level system being initially in equilibrium with its environment, and by determining the 2D-spectra of a donor-acceptor model. It is important to mention that we have considered systems with environmental spectral densities of the form of Eq. (57) because the HEOM approach is well suited for such problems. If we work with an environment whose dissipation and noise kernels can not be approximated by a finite number of exponentially decaying functions, then the combination of stochastic unravelling and the equations for multitime kernels should become the preferable approach since its complexity, in contrast to the HEOM approach, will not increase as long as the kernels decay sufficiently fast.

In the examples of Sec.  III we have always assumed that all kernels are nonzero only for a finite range of times. However, this assumption is of course not fulfilled for all models of interest. As a trivial counterexample we can consider a spin-boson model with Hamiltonian

[TABLE]

The spectral density of the environment contains a $\delta$-peak which results in non-decaying multitime kernels. In general, we expect that the multitime kernels will decay sufficiently fast to zero if the spectral density of the environment is smooth enough. In cases where the stochastic unravelling method fails and the system-bath coupling $\lambda{}_{SB}$ is weak enough, we can still try to approximate all multitime kernels by expanding them in powers of $\lambda{}_{SB}$ and taking only the first few terms into account.

Besides the slow decay of the memory kernel, another problem for our approach could be the size of the system. For a system Hilbert space of dimension $\mathcal{N}$ we have to work with kernels and propagators which are represented by $\mathcal{N}^{2}\times\mathcal{N}^{2}$ matrices. This has to be compared with the $\mathcal{N}^{2}$-dependence of the HEOM method on the system size.

Additional problems arise from the fact that the calculation of an $m$-time propagator in its full time domain requires the knowledge of all $m^{\prime}$-time propagators and $m^{\prime\prime}$-time kernels $(m^{\prime}<m,m^{\prime\prime}\leq m)$ in their full time domain, which leads to accumulation of the numerical error by an increase of $m$.

If we want to calculate an $(m+n)$-time propagator in the time domain, where $n$ of its arguments are fixed, we can not guarantee that we will have to know only a set of kernels/propagators, whose time domain is at most $m$-dimensional. In the first example that we have considered in Sec. III.1 we have fixed the last two arguments of $W(t_{3},t_{2},t_{1})$ to $(t_{2},t_{1})=(0,-t_{R})$. But for the calculation of $W$ in its one-dimensional time domain we needed the pairs $\mathcal{K},U$ and $\tilde{K},\tilde{U}$, whose time domains were one- and two-dimensional, respectively. On the other hand, in the second example in Sec. III.2 we have fixed the second argument of $I$ and $U^{\mu^{\times}_{R}\mu^{\times}_{R}},U^{\mu^{\times}_{R}\mu^{\times}_{L}}$. For the calculation of $U^{\mu^{\times}_{R}\mu^{\times}_{R}},U^{\mu^{\times}_{R}\mu^{\times}_{L}}$ we used kernels/propagators whose time domains were at most two-dimensional.

In summary, we have presented a method for the calculation of MTCFs of systems which span finite Hilbert spaces. In the first step we calculate the kernels via a set of equations. The input information can be obtained via a modification of the HEOM method or via a stochastic unravelling method. In the second step the kernels are used for the calculation of MTCFs by use of equations which can be derived by a few simple reduction rules. Thus, the main advantage of the present method is that it can be applied to problems, where HEOM does not perform well, as long as the system is sufficiently small and the memory kernels decay sufficiently fast.

Acknowledgements.

HPB acknowledges support from the EU Collaborative Project QuProCS (Grant Agreement 641277).