# Exact Stochastic Unraveling of an Optical Coherence Dynamics by Cumulant Expansion

## Abstract

A numerically exact Monte Carlo scheme for calculation of open quantum system dynamics is proposed and implemented. The method consists of a Monte-Carlo summation of a perturbation expansion in terms of trajectories in Liouville phase-space with respect to the coupling between the excited states of the molecule. The trajectories are weighted by a complex decoherence factor based on the second-order cumulant expansion of the environmental evolution. The method can be used with an arbitrary environment characterized by a general correlation function and arbitrary coupling strength. It is formally exact for harmonic environments, and it can be used with arbitrary temperature. Time evolution of an optically excited Frenkel exciton dimer representing a molecular exciton interacting with a charge transfer state is calculated by the proposed method. We calculate the evolution of the optical coherence elements of the density matrix and linear absorption spectrum, and compare them with the predictions of standard simulation methods.

## 1Introduction

The theory of absorption spectra of molecular aggregates depends on a correct evaluation of the time evolution of the investigated system, either in a formalism of the wave function or the density operator. In Condon approximation, i. e. when the transition dipole moment of the investigated system can be assumed independent of the nuclear degrees of freedom (DOF), one can naturally separate the studied system into an electronic one, which is affected by light and treated explicitly, and a nuclear bath, which is treated by some effective theory and does not directly interact with light. Linear absorption spectrum requires an evaluation of the time-dependent polarization, and consequently of the optical coherence elements of the electronic density matrix. For a two level system linearly coupled to a bath of harmonic oscillators, this problem can be solved exactly for an arbitrary coupling strength [1]. The theory is thus very well established particularly in the limit of strong system-bath coupling (SBC), when the resonance coupling between individual transitions in an aggregate can be neglected. In the weak SBC limit, perturbative treatments in second order can be applied, leading to various forms of quantum master equations [2]. An interesting observation in this context is the convergence of the second order results and the exact solution for a two level system. The exact solution is provided by a second order equation [1].

Many practically important cases, especially those related to photosynthetic antennae, fall in between the two treatable limits. In the weak SBC limit one can successfully base the theory on the eigenstate basis of the electronic Hamiltonian, the so-called *excitonic basis*. In the strong SBC limit, on the other hand, one can rely on the so-called *local basis*, i. e. the basis of states local to the molecules involved in the aggregate. No such preferred basis is readily available for an intermediate SBC strength. The intermediate coupling results in a renormalization of the energy eigenstates of the Hamiltonian with respect to the excitonic basis, which leads e. g. to a presence of non-zero off-diagonal elements of the stationary long time reduced density matrix [5]. From the point of view of the reduced density matrix propagation, these off-diagonal elements are a result of a coupling between equations of motion for the state populations and the coherences between the basis states, and/or a coupling between different coherences. The weak coupling limit is often accompanied by the so-called *secular approximation*, which neglects such coupling. The effects of this coupling are correspondingly termed *non-secular effects*. In some particular cases, non-secular effects can be quite dramatic, leading to a noticeable localization of the eigenstate basis, exhibited e. g. by a temperature dependent absorption band shift from the values corresponding to a delocalized excitonic basis to those corresponding to a local basis [7]. Similar effects lead to significant changes of fluorescence depolarization dynamics in photosynthetic antenna as was demonstrated in Ref. [8]. A theory describing correctly the dynamics of the reduced system interacting with a bath under arbitrary SBC strength would have to be able to predict a smooth change of the basis in which the long time reduced density matrix is diagonal from the excitonic bases (weak coupling) to the local basis (strong coupling). We will denote the basis in which the reduced density matrix diagonalizes at long times as the *preferred basis* in analogy with the preferred basis problem studied in decoherence theory [9].

Standard approaches to calculation of the dynamics, based on the second order master equations, often provide unreliable results for non-secular effects. Outside the secular approximation, one cannot, for example, guarantee positivity of the density matrix (see e. g. [6] and the references therein). There is a big variety of different methods to solve the dynamics of the molecular aggregates apart from the second order master equation theory. Among stochastic methods we find many path integral techniques [10], Monte Carlo methods [14], and also the stochastic wavefunction method [21], quantum state diffusion method [22], multi-configuration time-dependent Hartree method [23], mixed quantum-classical Liouville equation [24] and some other approaches [25]. In the formalism of master equations, the Nakajima-Zwanzig identity [26] and the Hashitsumae-Shibata-Takahashi (also known as convolutionless) identity [28] are often used as a starting point of higher order perturbation expansion theories [30]. Some exact results can be obtained for particular model systems, for example for the problem of a single molecule with one electronic transition. This problem is analytically solvable through second cumulant of Magnus expansion [31] exact for harmonic bath [34]. An exact solution for a system of multiple coupled transition can be obtained by the Hierarchical equations of motion (HEOM) [35]. Both, path integral methods and HEOM reach the preferred basis in the long-time evolution [38].

The HEOM method became very popular recently for calculation of excited state dynamics in photosynthetic systems as it combines feasibility with accuracy [39]. It has been implemented on modern parallel computers [44] and graphics-processing units (GPU-HEOM) [45]. One limitation of the method is that the calculations become more difficult with decreasing temperature [39].

In this paper, we propose another method which provides a formally exact solution to the reduced density matrix problem for harmonic baths. The method is based on a stochastic unraveling of the equation of motion for the reduced density matrix in the resonance coupling term. The leading idea is to cover the resonance coupling term in Hamiltonian by stochastic unraveling rather than doing it with the SBC as it is usual in ordinary stochastic methods. The evolution of the system’s state is modeled by an ensemble of trajectories in the space of the projectors on the states in the system’s Hilbert space. This projector space is known in the theory of non-linear spectroscopy as the *Liouville space* (see Ref. [46]). Each trajectory from the ensemble can be assigned a sequence of resonance coupling-free evolution operators that remains after the unraveling. The resulting expression is related to the high order non-linear response functions, and it can be evaluated analytically with the knowledge of the bath correlation function. The properly weighted sum over trajectories gives an exact result for the system’s reduced dynamics.

The proposed method can offer an advantage over the existing exact methods in systems with strong system-bath coupling and comparatively weak resonance coupling since the strong coupling to the bath does not increase the computational cost. It can be also very well used in systems with complicated spectral densities, which require more work in other methods like the HEOM.

We apply our method to the case of a heterodimer, motivated by the works on an interaction between charge transfer (CT) states and the excitons in photosynthetic reaction center [7]. In the absorption spectrum of this system, one can observe a large blue shift of the lowest energy band with increasing temperature [47]. Previous works used either the local basis [7] or the excitonic basis [5] as a starting point of their theory. It was concluded that the large reorganization energy of the CT state is the reason for the temperature dependent shift of the absorption band. It can be shown that within the weak SBC theory, the necessary condition for the band shift is the difference of the reorganization energies of the two involved types of states, the excitonic and the CT states. In second order theories, the non-secular terms, and correspondingly the shift, vanishes when the reorganization energies in the dimer are the same [48]. Because in both Refs. [7] the description of the shift is provided by theories which are effectively outside their respective range of validity, it is important to compare to the shifts predicted by exact theories. We therefore concentrate on a dimer in which the one dipole forbidden local state is characterized by a large reorganization energy and a zero transition dipole moment (playing thus a role of the CT state), and the second state is optically allowed, characterized by a moderate reorganization energy (playing a role of an excitonic state). Due to a resonant interaction between these two excited states, we observe two peaks in the absorption spectra, which shift as a function of the system parameters and various approximations discussed.

The paper is organized as follows: In Section ? we introduce the model dimer and its description. Section ? is devoted to the theory of absorption spectrum and in Section ? we introduce our stochastic unraveling of the equations of motion. The application of the unraveling is demonstrated and compared with other methods for calculation of the reduced dynamics in Section ?.

## 2Model System and Formalism

We illustrate the proposed method on a molecular dimer which represents the interacting CT state - exciton system. Considering just the linear absorption, this problem is equivalent to a standard problem of a molecular heterodimer, and we will therefore formulate it as such. The special characteristics of the CT state - exciton problem enter only through a specific set of parameters used for numerical demonstrations.

Each molecule of the dimer is considered to be a two-level system, either in the ground electronic state or the excited electronic state . Index numbers the molecules. We introduce the local basis of collective states , , , and we remove the double-excited state , since it is far off-resonant and does not contribute in calculation of linear absorption spectrum. The total Hamiltonian is written as sum of Hamiltonian of non-interacting monomers and the resonance coupling interaction term

where are energies of collective states , is the collective ground state energy and is the resonance coupling between electronic states. We denote , the potential electronic surfaces (PES) of the bath DOF around molecule in the ground and the excited states, respectively, and represents the kinetic term of nuclear DOF around molecule . It is more convenient to work with collective potential and kinetic terms defined as

The PES and depend on the generalized coordinates , and we define a collective coordinate

We split the Hamiltonian, Eq. (Equation 1), into the system, system-bath and bath parts in a standard manner [6]

The angular brackets denote the bath averaging which is defined as

where the trace is performed over the bath DOF, the density matrix of the bath is assumed to be of the

representing the canonical equilibrium, and the symbol denotes an arbitrary operator. In the definition of , Eq. ( ?), we used the so-called energy-gap operator defined as

The dynamics of the open quantum system is described by the reduced density matrix (RDM)

where is the density matrix of the total system.

Before we proceed with the discussion of the bath model, we introduce the so-called superoperator formalism, which is advantageous for the description of open quantum systems. We define Liouville superoperator (the Liouvillian) as

where index can hold values , , , and . We also define evolution superoperators

where is an ordinary evolution operator.

The bath is represented by harmonic oscillators in our model, further referred to as harmonic oscillator model (HOM). We pay special attention to three cases: Case (no bath is present) demonstrates the use of proposed method on calculation of pure quantum state dynamics. Case can serve as a good test of the method, because time dynamics of complete density matrix can be found explicitly, and we can compare between the explicit calculation and calculation of reduced dynamics by the proposed method. Case is chosen as a typical example of an open quantum system with irreversible dynamics. For , we simply put

For , we use definition

with the creation and annihilation operators

The potential term of the harmonic oscillator in the excited state is shifted by with respect to its potential in the ground state. In explicit calculation, we solve the Liouville – von Neumann equation to obtain dynamics of . The RDM, Eq. (Equation 5), is then obtained by tracing over DOF of the harmonic oscillator. In a calculation by the proposed Monte-Carlo method, however, we treat the system as open, i. e. the bath is described via the so-called *energy gap correlation functions* (EGCF) [46]

and their integrals – the lineshape functions

We use the energy gap correlation function [46]

for . This can be derived from the explicit form of bath operators, Eqs. (Equation 8- ?), if we assume , , and . The parameter

represents the reorganization energy of the system, and is the thermodynamic beta with temperature . The EGCF, Eq. (Equation 10), provides no relaxation or temperature; is the temperature of the initial condition from which the system gets excited. Note, that in our manuscript, we use the convention in which reorganization energy has dimension of energy and the EGCF has dimension of its square. In spectroscopy, it is also very common to use frequency that corresponds to the reorganization energy in place of and EGCF has dimension of frequency squared. Both conventions differ by the factor in appropriate power.

For the case , we use EGCF of an overdamped harmonic bath [46]

The parameter is given by the characteristic time of damping of the oscillators . The so-called Matsubara frequencies are defined as .

## 3Theory of Absorption Spectra

Linear optical properties of a given system are fully characterized by the Fourier transform of the linear response function [46]. In particular, the linear absorption coefficient can be evaluated as

where is frequency-dependent refractive index and is the speed of light. Let us describe the interaction of the dimer system with electric field by semi-classical Hamiltonian in dipole approximation

where we introduced the transition dipole moment operator . We rewrite using its absolute value and the polarization vector . Now, we can express the response function as

where we introduced the transition dipole moment operator projected on the polarization vector

and where is the initial density matrix of the system.

We can notice that all we need to calculate the absorption spectrum of the molecular dimer are the coherence elements of the evolution superoperator. The elements are often neglected which corresponds to the secular approximation. In Section ?, where we compare our method with standard methods, we plot directly these elements of evolution superoperator.

## 4Stochastic Unraveling of Coherent Dynamics with Pure Dephasing

### 4.1Basic Principle

We will start with the closed system, whose evolution superoperator is a solution of the Liouville equation

The reason for separation of into and is that for the case , we can solve the problem with the bath exactly via the cumulant expansion [32]. Provided there is initially no entanglement between the system and the bath, we can write as a time-ordered exponential using time-dependent perturbation theory

The interaction picture is taken with respect to the , i. e. . The assumption about the non-entangled initial state is quite natural for systems which have an optical energy gap and reside in the equilibrium state corresponding to electronic ground state before they are excited by a short laser pulse.

The proposed stochastic scheme is the following: We generate trajectories, where system exhibits random jumps between projectors on electronic states on Liouville space of electronic states. The jumps are generated in such a way that they reconstruct the action of electronic -coupling, i. e. of . Between the jumps, the system evolves according to . We introduce a discretization of the time axis into intervals . The model is exact in the limit . In every time step , there is a probability of the jump between states and the same probability of the jump . In addition to the time evolution according to , the trajectory weighting factor is multiplied by a complex number , which we will call “coherent factor” (CF) further on in the text. The coherent factor assures the correct stochastic unraveling. For each jump between bra-states, the trajectory gets a factor of , while for a jump between ket-states, it gets a factor of Hence

where and are the numbers of the jumps between bra-states and ket-states in the given trajectory, respectively.

If we introduce jump superoperators as

we can describe a trajectory with jumps in times , , by a sequence of jumps , , , , where every index should be replaced by “details” of the jump, i. e. it should specify if it is a jump in bra or ket vector, and it should state between which of the states the jump occurs. The total evolution superoperator can be then written as

The index numbers the trajectories, and is number of trajectories.

To show that Eq. (Equation 14) gives correct result for evolution superoperator, we will investigate the individual terms of the expansion, Eq. (Equation 12). We can see that the term “” is covered by trajectories with no jumps, which have the probability

where is the probability that no jump occurs in time interval . If the trajectory starts in a projector reads as

The term

of Eq. (Equation 12) is represented by trajectories with one jump at time , which constitute the individual terms of the sum in Eq. (Equation 17). A trajectory with one jump evolves according to the evolution superoperator , which corresponds to the time evolution in the projector for time , the action of which transfers the projector into the projector and a time evolution in this projector for time . The factor “” (complex unity) is included in the . The probability of such a trajectory with time of the jump is given by

It yields the correct ratio

The Liouvillian is hence not explicitly present in the sum over trajectories, but it is included by ratio of trajectories with particular number of jumps, as

There is also a trajectory with a jump between the projectors and for each trajectory with a jump from to . It comes from the second term of the commutator, Eq. (Equation 6). The trajectory gets additional minus sign, and the CF is therefore , see Eq. (Equation 13). One easily verifies that trajectories with multiple jumps reconstruct the higher order terms of the expansion, Eq. (Equation 12).

### 4.2Bath influence

For the case of a closed system, the superoperators in Eq. (Equation 14) are obtained explicitly. For open systems, we perform a trace over bath DOF and the cumulant expansion (CE), and we get complex factors in terms of the lineshape functions, Eq. (Equation 9). The reduced evolution superoperator can then be written as

The factor

can be evaluated analytically using second order cumulant expansion in a manner similar to the evaluation of non-linear response functions (see e. g. [46]).

### 4.3Some Numerical Considerations

We described the basic principle of the method in the previous section, and we showed that it is equivalent to the time evolution via the expansion, Eq. (Equation 12). The sum over trajectories gives the evolution superoperator in some fixed time . is, however, calculated up to a normalization constant, which depends linearly on the number of trajectories and according to Eq. (Equation 15) also on time. We would like to generate trajectories to a maximum time , and use them to evaluate for all times in such a way that the trajectories are not generated for each time independently. One possibility is simply to use scaling of the normalization, Eq. (Equation 15), for times of every trajectory. For technical reasons, we decided to use a different method. We included the trajectory in summation only at times for which , where is the time of the last jump in the trajectory, and we ignored the trajectory in evaluation of the times . This also leads to the correct result, because the ratio of the trajectories that have a jump in the interval to the total number of trajectories is proportional to the scaling of normalization factor Eq. (Equation 15) with time

Ignoring trajectories at times thus provides a correct normalization.

### 4.4Connection to the Feynman-Vernon Influence Functional

The connection can be drawn between the described method and the well-known Feynman-Vernon influence functional [49]. The path integral for a time evolution of a reduced density matrix can be written in a time-discretized form on a Hilbert space of a finite dimension, similarly to Eq. (5) in [51]:

where , is the initial system reduced density matrix and , , , number states from the system’s Hilbert space. The sums run through the whole Hilbert space of the system. The time is discretized into steps of size . The influence functional has a form

where .

Our method, to which we will refer to as Stochastic Unraveling of the Resonance Coupling (SURC) in the rest of this work, relies on an approximation of the expression in the Eq. (Equation 21) using the Trotter expansion . This expression after further approximation yields

If the path integral was to be performed without any importance sampling at this point, at each time step, we would have trajectories that jump between the states and and gain factor and those that stay in the same state and gain factor “1”. Trajectories with too many jumps, however, tend to cancel each other. The SURC can thus be viewed as an importance sampling in which we prefer the trajectories with less jumps by factor , and we increase correspondingly the phase change if the jump occurs. This allows us to trivially perform most of the summations in Eq. (Equation 21), because terms for most of the cases. Luckily, our influence functional has a particularly simple form: If there is a sequence of consecutive states , ,, , for which , all factors with indices between and can be expressed as one factor in terms of the lineshape functions and the states can be excluded from the expression for the influence functional.

## 5Numerical Demonstrations

In this section, we study the dynamics of optical coherences in a molecular dimer in order to test the precision and numerical stability of our method. We will demonstrate how the SURC works for a simple system with no bath (case of HOM) and for an exactly solvable problem with simple bath model (case of HOM). We also calculate absorption spectra of a model dimer with full harmonic bath (case of HOM) and compare the results obtained by our method with those obtained by the time-dependent Redfield tensor (described in Ref. [6]), its secularized version and the HEOM.

### 5.1Hierarchical Equations of Motion

The HEOM as described in Ref. [38] is an exact method for systems with overdamped bath, Eq. (Equation 11). To achieve its convergence at low temperatures becomes increasingly difficult, since the computational cost increases exponentially with both the system size and the number of included Matsubara frequencies. We denote the calculations with HEOM method using a particular number of Matsubara frequencies by an abbreviation HEOM . For practical HEOM calculations on systems larger than a dimer, it is vital to employ the high-temperature approximation (HTA) to correlation function of Eq. (Equation 11). This yields [39]

and it requires . We denote the calculations with this method by abbreviation HEOM HTA. Compared to HEOM 0, i. e. to neglecting the Matsubara part of the EGCF completely, the HTA extends the range of validity of HEOM towards lower temperatures while at the same time it preserves the crucial exponential time-decay of the correlation function. This is essential for constructing the reduced hierarchy in a computationally efficient manner. Arbitrary spectral densities and lower temperatures are readily implemented in SURC using the exact EGCF, whereas the treatment by HEOM requires a decomposition of the spectral density into shifted peaks [52].

### 5.2Stochastic Unraveling of the Resonance Coupling

In order to test the ideas of Section ?, we first demonstrate that the method correctly reproduces coherent quantum dynamics. We set the same excitation energies of both molecules and a non-zero resonance coupling .

The bath is not present, and the dynamics can therefore be solved exactly. We calculate the element of the evolution superoperator, and we compare the exact dynamics with the dynamics obtained by SURC. The results are presented in Fig. ? We can see the results of three runs of SURC, each of them performed with trajectories. The time step is chosen so that it smoothly discretizes the line shape functions, the oscillations given by the difference between electronic transition energies of the monomers and the oscillations due to the coupling . For trajectories, numerical noise leads to accuracy breakdown for times longer than approximately , and this value is independent of when it is chosen to give a smooth discretization. We chose to discretize the interval into 2048 steps which gives the time step . In calculations where the time step would be larger than 0.25 fs, we chose instead. As the Fig. ? demonstrates, SURC dynamics is very close to the exact one before the dynamics is overwhelmed by noise.

To show that the SURC method is exact also for an open system with harmonic bath, we calculate element of the evolution superoperator for a molecular dimer with simple bath that allows exact solution of the model. The bath is represented by a single undamped harmonic oscillator (case of HOM) for each molecule. In the SURC calculation, harmonic oscillators are treated implicitly by EGCF, Eq. (Equation 10). This calculation is compared in Fig. ? with the exact explicit solution. The calculated system is a homodimer with parameters and the temperature of initial condition, Eq. (Equation 3), equals to . We performed four calculations for values of , , and . We can see that the correspondence of the exact solution with the SURC solution is very good. The time up to which we can calculate is again inversely proportional to .

Finally, let us compare the absorption spectra calculated for a bath model with by SURC and other frequently used theories. All the following SURC calculations are performed with the use of trajectories. This prolongs the time before the accuracy breakdown to approximately , and it required approximately 2.5 CPU-years of Intel Xeon CPU E5620 @ 2.40GHz per calculation. We have chosen the same number of time steps on the extended interval yielding . We use a dimer model with parameters . The big difference in reorganization energies is typical in situations where the molecules have very different surroundings (e. g. protein envelope and water) or for systems with both excitonic and CT states, as discussed in the Introduction. Only the site 1 has non-zero transition dipole moment, since the site 2 represents a CT state. The length of the dipole moment is arbitrary since it only changes normalization of the spectra.

Figs. ? and ? present absorption spectra of the studied dimer system. The Fig. ? shows absorption spectra in three cases: , and . For , this is a problem of non-interacting monomers which is exactly solvable [1], and all theories give the same result. With a gradual increase of , we see an increase of the excitonic splitting, noticeable as a shift of the higher of the peaks (lower transition frequency) to the red. The results of non-secular Redfield theory and SURC calculations are similar, while the secular approximation gives an exaggerated peak splitting. The HEOM, which is considered an exact theory in this regime of parameters, is in very good agreement with the SURC method. The Fig. ?A plots the positions of the lower frequency peak for the four theories and quantifies the difference in its shift with . In order to estimate the error of the SURC method, we divided the total 10 trajectories of every calculation into five independent simulation runs. For every run, we independently calculated the absorption spectrum and extracted the peak position. The error bars in Fig. ? are calculated as standard mean deviations of the peak position. Already here, we can conclude that the non-secular Redfield theory captures the decrease of the excitonic splitting with respect to the excitonic basis surprisingly well. This can serve as a support for the conclusions of Refs. [5] and [8] where non-secular effects in form of temperature dependent band shifts and changes in fluorescence depolarization dynamics were treated by Redfield theory.

In Fig. ? we plot the temperature dependence of the absorption spectra for the same model dimer. Here, we can notice that the SURC is in nearly perfect match with the HEOM HTA, and that the Redfield theory also captures the right tendency. We can notice a shift of the lower frequency peak to higher frequencies with increasing temperature. This demonstrates an increasing dynamic localization by the bath with increasing temperature. The secular Redfield gives opposite peak shift with respect to the other theories. This shift is caused by the slight temperature dependent changes in the lines shape, because the transition frequencies remain the same for all temperatures. For lower temperatures, we can notice a difference between the SURC and Redfield theory in the width of the wider peak. With increasing temperature, these two theories also increasingly differ in the lower frequency peak position, as is shown on Fig. ?B. As expected, the nearly perfect match of SURC and HEOM HTA gets worse at lower temperatures. Calculation by HEOM with increasing number of Matsubara frequencies shows a increasing degree of convergence towards the SURC result. There is a small discrepancy of between the HEOM HTA and SURC in high temperatures, where the theories should match. As the peak position is a very sensitive test, this is probably caused by small numerical differences between the program codes of HEOM and HEOM HTA. As can be seen on Fig. ?A, a difference on the same order between HEOM 2 and SURC at visible on Fig. ?C is almost not recognizable in the dynamics of the envelope of the optical coherence.

In its present form, the SURC is limited to short time evolution due to a breakdown of accuracy after the propagation time corresponding roughly to 1.5 times the period of the oscillations caused by the resonance coupling . This is in turn caused by the dynamic sign problem [53] originating from the phase factor assigned at each jump in a given trajectory. For the cases discussed here, the GPU-HEOM implementation [54] exceeds the SURC implementation in speed. However, the SURC is also massively parallelizable. Moreover, it has only minimal memory requirements independent of the system size and its computational cost does not increase with the complexity of the EGCF. For example, the treatment of spectral densities with shifted peaks corresponding to underdamped vibrational modes by HEOM can be a challenging task in some cases since the HEOM has computational difficulty growing as for EGCF of Eq. (Equation 11) and for spectral density with shifted peaks. Computational costs of HEOM HTA grow as both in memory and CPU time. Here is the number of sites and is the Hierarchy depth required for convergence [52]. In our calculations, typical hierarchy depth is . The SURC method scales in CPU and memory proportionally to the dimension of the reduced density matrix, which corresponds to for calculation of the absorption spectra and to in the block of single exciton states. Clearly, the initial cost of the SURC method is high in current implementation. When a complex EGCF is required for the description of a given bath, the cost of SURC calculation does not increase. The SURC may therefore offer advantage over the HEOM for larger systems with more involved spectral densities, or serve as a good check of the HEOM convergence.

## 6Conclusion

In this paper we proposed a Monte-Carlo method of evaluation of optical coherence dynamics based on stochastic unraveling of resonance coupling via cumulant expansion, which is exact for a harmonic bath. The method can be easily used with a bath specified by an arbitrary energy gap correlation function. We verified that the method correctly reproduces the coherent dynamics, and that it gives correct dynamics for exactly solvable model of molecular dimer with one harmonic oscillator coupled to each molecule representing the bath. We calculated absorption spectra of a model dimer for different values of a resonance coupling and demonstrated that the proposed scheme gives results very similar to the Hierarchical equations of motion. We investigated the temperature dependence of the spectrum for the same dimer, finding again good agreement between the Hierarchical equations of motion and the method introduced in this paper. Comparison with non-secular Redfield theory leads us to a conclusion that despite its well-known problems, such as the positivity breakdown under some parameters, the Redfield theory describes the non-secular effects in optical coherences rather well for the typical parameters used in our work.

### References

- R. Doll, D. Zueco, M. Wubs, S. Kohler, and P. Hanggi, Chem. Phys.
**347**, 243 (2008). - V. May and O. Kühn,
*Charge and Energy Transfer Dynamics in Molecular Systems*(Wiley-VCH, Berlin, 2001). - W. M. Zhang, T. Meier, V. Chernyak, and S. Mukamel, J. Chem. Phys.
**108**, 7763 (1998). - M. Yang and G. R. Fleming, Chem. Phys.
**275**, 355 (2002). - T. Mančal, L. Valkunas, and G. R. Fleming, Chem. Phys. Lett.
**432**, 301 (2006). - J. Olšina and T. Mančal, J. Mol. Model.
**16**, 1765 (2010). - T. Renger, Phys. Rev. Lett.
**93**, 188101 (2004). - O. Kühn, V. Sundström, and T. Pullerits, Chem. Phys.
**275**, 15 (2002). - M. Schlosshauer,
*Decoherence and the Quantum-to-Classical Transition*(Springer, Berlin, 2007). - J. T. Stockburger and H. Grabert, Phys. Rev. Lett.
**88**, 170407 (2002). - Q. Shi and E. Geva, J. Chem. Phys.
**119**, 12063 (2003). - R.-X. Xu, P. Cui, X.-Q. Li, Y. Mo, and Y. Yan, J. Chem. Phys.
**122**, 041103 (2005). - P. Huo, S. Bonella, L. Chen, and D. F. Coker, Chem. Phys.
**370**, 87 (2010). - J. Hubbard, Phys. Rev. Lett.
**3**, 77 (1959). - W. J. D. Beenken, M. Dahlbom, P. Kjellberg, and T. Pullerits, J. Chem. Phys.
**117**, 5810 (2002). - J. Shao, J. Chem. Phys.
**120**, 5053 (2004). - Y. Zhou, Y. Yan, and J. Shao, Europhys. Lett.
**72**, 334 (2005). - Y. Mo, R.-X. Xu, P. Cui, and Y. Yan, J. Chem. Phys.
**122**, 084115 (2005). - D. Lacroix, Phys. Rev. A
**72**, 013805 (2005). - J. Shao, Chem. Phys.
**322**, 187 (2006). - J. Dalibard, Y. Castin, and K. Molmer, Phys. Rev. Lett.
**68**, 580 (1992). - Y. Aharonov, L. Davidovich, and N. Zagury, Phys. Rev. A
**48**, 1687 (1993). - J. Schulze, M. Torbjörnsson, O. Kühn, and T. Pullerits, New J. Phys.
**16**, 045010 (2014). - S. Bai, W. Xie, L. Zhu, and Q. Shi, J. Chem. Phys.
**140**, 084105 (2014). - M. Mohseni, P. Rebentrost, S. Lloyd, and A. Aspuru-Guzik, J. Chem. Phys.
**129**, 174106 (2008). - S. Nakajima, Progress Of Theoretical Physics
**20**, 948 (1958). - R. Zwanzig, Physica
**30**, 1109 (1964). - N. Hashitsume, F. Shibata, and M. Shingu, J. Stat. Phys.
**17**, 155 (1977). - F. Shibata, Y. Takahashi, and N. Hashitsume, J. Stat. Phys.
**17**, 4 (1977). - M. Schröder, M. Schreiber, and U. Kleinekathöfer, J. Chem. Phys.
**126**, 114102 (2007). - W. Magnus, Commun. Pure Appl. Math.
**7**, 649 (1954). - R. Kubo, J. Phys. Soc. Jpn.
**17**, 1100 (1962). - S. Mukamel, Phys. Rev. A
**28**, 3480 (1983). - T. Mančal and F. Šanda, ArXiv 1011.3803v1 (2011).
- Y. Tanimura and R. Kubo, J. Phys. Soc. Jpn.
**58**, 101 (1989). - A. Ishizaki and Y. Tanimura, J. Phys. Soc. Jpn.
**74**, 3131 (2005). - M. Tanaka and Y. Tanimura, J. Phys. Soc. Jpn.
**78**, 073802 (2009). - A. G. Dijkstra and Y. Tanimura, Phys. Rev. Lett.
**104**, 250401 (2010). - A. Ishizaki and G. R. Fleming, Proc. Natl. Acad. Sci. U. S. A.
**106**, 17255 (2009). - A. Ishizaki and G. R. Fleming, J. Chem. Phys.
**130**, 234110 (2009). - A. Ishizaki and G. R. Fleming, J. Chem. Phys.
**130**, 234111 (2009). - B. Hein, C. Kreisbeck, T. Kramer, and M. Rodríguez, New J. Phys.
**14**, 023018 (2012). - C. Kreisbeck and T. Kramer, J. Phys. Chem. Lett.
**3**, 2828 (2012). - J. Strümpfer and K. Schulten, Journal of Chemical Theory and Computation
**8**, 2808 (2012). - C. Kreisbeck, T. Kramer, M. Rodriguez, and B. Hein, J. Chem. Theory Comput.
**7**, 2166 (2011). - S. Mukamel,
*Principles of nonlinear spectroscopy*(Oxford University Press, Oxford, 1995). - H. Huber, M. Meyer, H. Scheer, W. Zinth, and J.Wachtveitl, Photosynthesis Research
**55**, 153 (1998). - T. Mančal, L. Valkunas, E. L. Read, G. S. Engel, T. R. Calhoun, and G. R. Fleming, Spectroscopy
**22**, 199 (2008). - R. P. Feynman and F. L. Vernon, Annals of Physics
**24**, 118 (1963). - H. Grabert, P. Schramm, and G. L. Ingold, Phys. Rep.-Rev. Sec. Phys. Lett.
**168**, 115 (1988). - N. Makri and D. E. Makarov, J. Chem. Phys.
**102**, 4600 (1995). - Y. Tanimura and S. Mukamel, J. Phys. Soc. Jpn.
**63**, 66 (1994). - N. Makri, J. Math. Phys.
**36**, 2430 (1995). - C. Kreisbeck and T. Kramer, Exciton Dynamics Lab for Light-Harvesting Complexes (GPU-HEOM), 2013, electronic tool at nanohub.org, DOI: 10.4231/D3RB6W248.