Modelling of Oscillations in Two-Dimensional Echo-Spectra of the Fenna-Matthews-Olson Complex

Recent experimental observations of time-dependent beatings in the two-dimensional echo-spectra of light-harvesting complexes at ambient temperatures have opened up the question whether coherence and wave-like behaviour plays a significant role in photosynthesis. We perform a numerical study of the absorption and echo-spectra of the Fenna-Matthews-Olson (FMO) complex in chlorobium tepidum and analyse the requirements in the theoretical model needed to reproduce beatings in the calculated spectra. The energy transfer in the FMO pigment-protein complex is theoretically described by an exciton Hamiltonian coupled to a phonon bath which account for the pigments electronic and vibrational excitations respectively. We use the hierarchical equations of motions method to treat the strong couplings in a non-perturbative way. We show that the oscillations in the two-dimensional echo-spectra persist in the presence of thermal noise and static disorder.


Introduction
Light harvesting complexes are the part of photosynthetic systems that channel energy from the antenna to the reaction centre. One of the most studied complexes is the Fenna-Matthews-Olson (FMO) complex [1,2,3] which is part of the photosynthetic apparatus of green sulphur bacteria. The observation of long lasting beatings, up to 1.8 ps at T = 77 K in time-resolved two-dimensional (2d) spectra [4,5,6] has generated enormous interest in theoretical modelling the process of energy transfer through lightharvesting complexes. Recently molecular-dynamics simulations started to develop a more microscopic picture of the dynamics [7,8], but an atomistic and fully quantummechanical model of the whole process reproducing the experimentally obtained spectra is not available. Therefore theoretical models treat the light-harvesting pigment-protein complex as an exciton system coupled to a bath which accounts for the vibrational modes of the pigments.
The choreography of the exciton-energy transfer can be elucidated with 2d echospectroscopy using three laser pulses which hit the sample within several femtoseconds time spacings. The experimentally measured 2d echo-spectra of the complex show a variety of beating patterns with oscillation periods roughly consistent with the differences in eigenenergies of the excited system obtained from fits to the absorption and 2d echo-spectra [5]. At long delay times between pump and probe pulses the system displays relaxation to the energetically lowest lying exciton eigenstates as expected from the interactions between the electronic degrees of freedom and the vibrational modes of the complex which drive the system into a state of thermal equilibrium [9].
One main goal of the research in this field is to clarify the role of a disordered and fluctuating environment in the energy-transfer efficiency. Theoretical quantumtransport studies show that the energy transfer in these systems results from a balance between coherent and dissipative dynamics [10,11]. Calculations at room temperature also show that the optimal energy-transport efficiencies appear as a compromise between coherent dynamics and thermal fluctuations [12].
Another open question in the theoretical study of these systems is the description of experimental observables, i.e. the absorption and the 2d echo-spectra. Most of the theoretical studies analyse the beatings in the population dynamics and compare them with the ones observed in the experimentally measured 2d echo-spectra. However, the 2d echo-spectrum reflects different information than the population dynamics and needs to be calculated and analysed separately in detail. In principle techniques such as quantum state and process tomography made out of a sequence of 2d echo-spectra can be used to map out the complete density matrix [13]. In contrast to energytransfer efficiency studies where an initial excitation enters the complex at specific sites close to the antenna, in 2d echo-spectra the whole complex is simultaneously excited. Also the two-exciton manifold yields prominent contributions to the signal, resulting in negative regions in the 2d echo-spectra. The non-pertubative calculation of 2d echo-spectra presents a considerable computational challenge due to the presence of two excitons giving raise to excited state absorption and the requirement to consider ensemble averages over differently orientated complexes with slightly varying energy levels. Previous calculations have employed Markovian approximations [14,15] or exclude the double exciton manifold [16]. In addition the systematic study of beatings in a series of 2d echo-spectra requires effective means to calculate a huge number of such spectra. So far no theoretical method has been able to describe the longlasting beatings in the time-resolved 2d spectra [17,14,16]. One possible explanation for the persistence of long coherence times has been the sluggish absorption of the reorganisation energy by the molecule, which requires theoretical descriptions that go beyond the Markovian approximation and the rotating wave approximation [18]. The hierarchical equations of motions (HEOM), first developed by Tanimura and Kubo [19] and subsequently refined in [20,21,22,23], show oscillations in the dynamics of the exciton populations that persist even at temperature T = 300 K [24,25]. The HEOM include the reorganisation process in a transparent way and are directly applicable for computations at physiological temperatures. A recent calculation at temperature 77 K of the 2d echo-spectra with the HEOM method has been performed by Chen et al [17], which does not display the strong beatings observed experimentally. Chen et al conclude that the agreement between theory and experiments needs to be improved. Here, we analyse under what conditions the theoretical calculations give better agreement with the experimental results. We calculate the absorption and 2d echo-spectra in a temperature range between 77 K to 277 K using the HEOM and investigate the role of pigment and protein vibrations in the calculated spectra. We systematically discuss the parameters needed to obtain results that agree better with the experimentally measured spectra. We discuss the role of static disorder for the shapes of the peaks and the appearance of beatings in the time-resolved 2d echo-spectrum and their robustness against temperature changes and static disorder.
The paper is organised as follows: the exciton model and the HEOM are introduced in sec. 2. The absorption spectra and the line shapes obtained from the HEOM model of the FMO complex are analysed in sec. 3 where we compare different methods for performing rotational averages. In sec. 4 we discuss the 2d spectrum and analyse how disorder affects it. Furthermore we analyse the persistence of beatings in the spectra as a function of the delay time for several temperatures and with and without static disorder.

Exciton model of the FMO complex
The FMO complex described by Fenna, Matthews and Olson [2,3] is a pigmentprotein complex that channels the energy in the photosynthetic apparatus in green sulphur bacteria. The FMO complex is arranged as a trimer with the different subunits interacting weakly with each other. We thus restrict our study to a single subunit which contains eight bacteriochlorophyll molecules (BChls) [26,27]. The BChls pigments, which are wrapped in a protein environment, can be electronically excited and are coupled due to dipole-dipole interactions. The BChls form a network that guides the energy from the antenna to the reaction centre. Since the eighth BChl is only loosely bound it usually detaches from the others when the system is isolated from its environment to perform experiments [26,28]. Therefore, we do not take the eighth BChl into account and describe the energy transfer using the seven site Frenkel exciton model [29,30,31]. The single exciton manifold is given by with N sites = 7. The state |k corresponds to an electronic excitation of BChl k with site energy ǫ k . Dipole-dipole couplings between different chromophores yield the inter-site couplings J kl , which lead to a delocalisation of the exciton over all BChls. The site energies of the BChls ǫ k = ǫ 0 k + λ k consist of 'zero-phonon energies' ǫ 0 k , and a term due to the reorganisation Hamiltonian H reorg = N sites k=1 λ k |k k|. The reorganisation energy λ k = ( i ω i q 2 i /2) k denotes the energy which is removed from the system to adjust the vibrational coordinates to the new equilibrium value [31,24,32]. Here, q i denotes a dimensionless displacement of the ith oscillator measured in units of /(m i ω i ), where m i and ω i are mass and frequency of the ith oscillator. The vibrational modes of the BChl are modelled as a set of harmonic oscillators The vibrational environment couples linearly to the exciton system H ex−phon = N sites k=1 |k k| and induces fluctuations in the site energies. We neglect correlations between the vibrations at different sites and model the coupling strength by a structureless spectral density where we assume identical couplings λ k = λ at each site. This choice of spectral density facilitates the use of the HEOM method since the time-dependent bath correlations then assume an exponential form. For the FMO complex a variety of spectral densities are discussed in the literature, where some are a parametrisation of experimental fluorescence spectra [33], while others are extracted from molecular dynamics simulations [34]. We use the model Hamiltonian for chlorobium tepidum given in Ref. [14]. The Hamiltonian originates from previous calculations in [35,36]. Cho et al. [14] changed the coupling constant between BChl 5 and BChl 6 and obtained new site energies from a fit to experimental results. The Hamiltonian of Ref. [14] has been used in previous studies of 2d echo-spectra [14,16,17] in units of cm −1 , where we have added an energy shift λ S = 12075 cm −1 to the diagonal elements of the exciton Hamiltonian H 1ex e to shift the positions of the eigenenergies to the frequency range of the figures in Ref. [14]. We set the parameters of the spectral density to λ = 35 cm −1 and γ −1 = 50 fs. These parameters are chosen to provide best agreement to the smooth part of the spectral density in Ref. [33], where the low frequency component of the spectral density was fitted to fluorescence line narrowing spectra [37]. Since a single Lorentz-Drude peak is not sufficient to accurately reproduce the smooth part of the spectral density of Ref. [33], we choose the parameters to obtain agreement within a factor of 1.2 to 1.5 in the energy range from 100 to 150 cm −1 , which corresponds to the typical difference between the exciton energies.
The HEOM approach rewrites time non-local effects into a hierarchy of coupled equations of motion for a set of auxiliary matrices σ n . The superscript n is a vector of N sites entries, where n k gives the order to which the bath of BChl k is taken into account. The reduced density matrix ρ(t) = tr phonon R(t) describing the exciton system is obtained after tracing out the vibrational degrees of freedom. The time evolution of the reduced density matrix is given by L e is the Liouville operator of the coherent exciton dynamics and we define e k as the kth unit vector with N sites = 7 entries. The operators Φ k and Θ k are defined by their action on a test operator A where The auxiliary matrices are initialised to σ n (0) = 0 for n = 0 and σ 0 (t) = ρ (t) [24]. The hierarchy is truncated for sufficiently large N max = N sites k=1 n k , where the diagonal coupling in (5) becomes dominant. For the parameters used in the calculation, reorganisation energy λ = 35 cm −1 and phonon relaxation time scale γ −1 = 50 fs, we verified that the relative differences in the interesting energy region of the absorption spectrum using N max = 4 versus N max = 12 are below 0.13 percent. In the 2d echo-spectrum we compared the rephasing stimulated emission diagram for N max = 4 and N max = 12 and found a difference of less than 0.07 percent of the amplitude. We therefore consider a truncation at N max = 4 sufficient to reach converged results. The derivation of the hierarchy stated in (4) and (5) relies on a high temperature approximation. For a temperature of 77 K the high temperature approximation is not valid and low temperature correction terms have to be included [25]. Within the low temperature correction we replace where ν 1 = 2π/β is the first bosonic Matsubara frequency following [25]. Agreement with [38] and proper thermalization to the Boltzmann distribution [12] confirms that the first order correction in (8) is indeed sufficient. A further improvement of the low temperature correction can be attained by including higher Matsubara terms. Recently, also a Padé spectrum decomposition has been proposed [39,40] which shows faster convergence than the Matsubara decomposition. The numerical evaluation of the HEOM requires to propagate a large number of auxiliary matrices. Although the number of auxiliary matrices can be reduced by a scaled HEOM approach [41,42] as well as by advanced filtering techniques [41], the HEOM are computationally very demanding, which sets limits to treatable system sizes. For 2d echo-spectroscopy one has to propagate the set of auxiliary matrices not only for a huge number of time steps (typically of the order of 10 5 ), but also the dimension of the auxiliary matrices gets enlarged since excited state absorption requires to extend the exciton Hamiltonian to the two exciton manifold. We overcome the computational limitation by using a GPU (Graphics Processing Unit) implementation where the computation time goes down with the number of cores on a single GPU [12]. For example on a single NVIDIA C2050 graphics board we obtain a 450 fold speedup compared to an equivalent single-core CPU implementation. With our efficient GPU algorithm the computation time for a 2d spectrum is reduced to 50 minutes on a single GPU.

Absorption spectra of the FMO complex
Spectroscopy provides a tool that gives direct insight into the energy states of a quantum system. Absorption spectroscopy is used to identify the energy eigenstates and their interaction strengths with the radiation field, expressed in terms of dipole moments. The energy eigenstates and the dipole moments are commonly identified with the peak position and the peak areas respectively. The spectral peaks are broadened due to the interaction of the energy eigenstates with the remaining degrees of freedom of the system. For an ensemble of light-harvesting complexes the coupling of the electronic dynamics to the vibrational modes and the presence of fluctuations in the energy levels, due to e.g. the motion of the surrounding proteins, broadens the absorption peaks.
To calculate absorption-spectra we have to expand the system and add the ground state |0 to our basis. The Hamiltonian is thus extended to an 8 × 8 block matrix where we have chosen the zero exciton energy H gs e = 0 and H 1ex e is the 7 × 7 matrix defined in (3).
Using the same block structure as in (9) the dipole matrix that governs the photonexciton conversion is constructed as where T and the zero to single excitation dipole block is defined using the seven dipole moments µ k that give the excitation of BChl k by the laser pulse.
The absorption spectrum is calculated as where ρ 0 = |0 0| is the density operator for no excitation and we assume δ-shaped laser pulses. The different methods for performing the rotational average denoted by · rot are detailed below. The time evolution of the dipole operators in (12) is given by In the FMO complex the seven BChl molecules have a fixed orientation with respect to each other. Within this configuration each BChl has a dipole vector d k that characterises its interaction with light. We consider the dipole vectors of the BChls to be along the connection line of two nitrogen atoms N B −N D [43] whose positions we obtain from the protein data bank [26,44]. The norm of d k accounts for the strength of the dipole and we assume d k = d since all BChl molecules are identical. For the simulations we choose d = 1. Given a polarisation direction of the laser l = (sin β cos γ, sin β sin γ, cos β), where γ ∈ [0, 2π) and β ∈ [0, π], the dipole moments in (11) read Typical experiments are performed using spectroscopy of a sample solution and many FMO complexes are illuminated at the same time. Each of them is randomly orientated with respect to the electric field in the laser pulse. If we assume a coherent dynamics expression (12) is readily evaluated in eigenbasis as where E j is the jth eigenenergy and v j k = E j |k denotes the overlap of the kth state in site basis with the jth eigenstate. The integrationally averaged dipole moments in eigenbasis areμ Thus, we obtain the integrationally averaged squares of dipole moments in the eigenbasis  [14].
However, we can also read the result asμ 2 where the three polarisation vectors l i are orthogonal to each other. This now implies that the rotational average can be performed using three orthogonal directions for l instead of sampling over the whole unit sphere. Note that this calculation is valid only if the interaction of system and bath does not interfere with the rotational average.
In the following we discuss how the line-shapes of the peaks depend on the theoretical method and approximations commonly used for performing the timedependent propagation in (12) before we return to the different methods of performing the rotational average.

Absorption spectra obtained by different propagation methods
We calculate the absorption spectra of the FMO complex at the temperature T = 77 K. We propagate the exciton-bath system during 2048 fs with a numerical time step of 2 fs using the HEOM with reorganisation energy λ = 35 cm −1 and a time scale for the bath correlations given by γ −1 = 50 fs. Each HEOM propagation takes only 1 second on a GPU.
In figure 1 (a), we compare the results obtained for the absorption spectrum calculated using the HEOM and the secular Born-Markov approximation (SBM) following Refs. [45,12]. The line shapes for each eigenenergy can be calculated in the eigenbasis I j = ∞ 0 e iωt tr µ E j (t) µ E j (0) ρ 0 rot . In the SBM simulation the complete spectrum is identical to the sum of the seven positive absorption peaks (not shown). Due to the rotating wave approximation the terms tr µ E j (t) µ E l (0) ρ 0 for j = l are zero. For the HEOM, these terms have non-vanishing contributions and the sum of the line shapes I k do no longer yield the spectrum. Thus we have to calculate the line shape of the jth exciton peak as . Including these terms leads to more complicated line shapes with shoulders and negative contributions as shown in figure 1 (b) for the peak of the second eigenstate. We find that the more complicated line shapes arise from the non-secular approximations. Comparisons of different shapes of the spectrum obtained by different propagation methods and different spectral densities have been presented in [46,47] for the B850 ring of the light-harvesting complex LH2.
We observe in figure 1 (a) that the energies for the HEOM peaks are shifted compared to the peaks obtained in the SBM approximation. Since we do not take into account the time-dependent Lamb shift in the SBM approximation (see eq. (8) in Ref. [45]) the corresponding SBM peaks are located at the exciton energies of Hamiltonian (3), marked by the vertical lines. Note that the eigenenergies of Hamiltonian (3) include the reorganisation energy λ added to the "bare excitonenergies". Due to the dissipation of the reorganisation energy within the vibrational environment, taken into account by the HEOM [24], the peaks of the absorption spectrum shift to lower energies. For a fast phonon relaxation γ −1 → 0 (Markov limit) the reorganisation energy dissipates instantaneously and the peaks are shifted to lower energies by approximately λ and then coincide with the bare exciton energies. Figure  1  longer phonon relaxation times, the reorganisation energy dissipates slowly during the exciton dynamics and the peaks stay closer to the eigenenergies of (3). While the spectra shown above are calculated using fitted dipole moments from [14], we will now analyse the changes due to rotational averaging on the absorption spectrum. In figure 2 (a) we show the results for the spectrum using the HEOM for λ = 35 cm −1 and γ −1 = 166 fs. We compare different methods of rotational average discussed in the previous section, that is, a single propagation from integrationally averaged dipole moments usingμ eq. (16), a Monte Carlo (MC) set of 300 orientations and µ FIT . The results using three propagations with perpendicular laser directions xyz are not shown since the curve is almost identical to the MC spectrum for 300 realisations. The MC method matches the experimental conditions closest, since the experiment shows the averaged spectrum of many randomly orientated complexes. To avoid propagating density matrices for many random orientations, the xyz approach can be chosen instead as both spectra show excellent agreement. Neither the integrated dipole moments, nor the fitted ones yield comparable results to the MC method and only in the Markovian limit it is possible to run the propagation with pre-averaged dipole moments. Here, we chose a longer bath correlation time γ −1 = 166 fs for the comparison of different rotational averages as the differences vanish for decreasing γ −1 .
In the following we always use several orientations for rotational averages and do not resort to pre-averaged dipole moments.

Broadening of the spectral lines
If we compare the only rotationally averaged absorption spectra in figure 1 (a) with the experimentally measured ones [14], we find that the rotationally averaged peaks are too narrow. Furthermore six peaks are well pronounced while in experiments several peaks overlap and the third exciton eigenstate is not individually resolved. This result is expected, since there are several additional effects besides rotational averaging, which broaden the peaks and bring the simulation results closer to the experimental data. If we increase either λ (not shown) or the relaxation time γ −1 , the peaks in the spectrum become broader, see figure 1 (b). Both parameters change the spectral density and hence the coupling to the vibrational modes. The third mechanism is static disorder caused by the slowly fluctuating protein environment. We model the static disorder by adding a Gaussian distributed noise of a given standard deviation (SD) to each diagonal term in (3). The resulting spectra are shown in figure 2 (b) where we obtain broader peaks and show with increasing disorder the diminishing of the third peak in the spectrum. In addition we observe that the broadening changes the relative peak heights of the second and fourth peak yielding a better agreement with the experimental results.
For the absorption spectrum changing the spectral density and inhomogeneous broadening yield similar results and it is not possible to distinguish which mechanism is relevant for the FMO complex. It is in the 2d echo-spectra where the difference becomes visible since the disorder results in a broadening along the diagonal frequencies only, while changing the exciton-phonon coupling broadens the peaks in all directions [48].

Two-dimensional echo-spectroscopy
In addition to the information provided by absorption spectroscopy, two-dimensional echo-spectroscopy is used to study time-resolved processes such as energy transfer or vibrational decay as well as to measure intermolecular coupling strengths [49,50,51]. Two-dimensional echo-spectroscopy consists of a sequence of three ultra-short pulses and a resulting signal pulse with fixed time-delay between the second and third pulse [52,53]. The theory of 2d echo-spectroscopy is explained in detail in [54,55,56]. The 2d echo-spectra are measured in two configurations which correspond to the rephasing (RP) and non-rephasing (NR) contributions of the third order response function . (17) Using the impulsive approximation which assumes δ-functions for the temporal envelopes of the pulses, the two components of the spectrum corresponding to  Figure 3. The six Feynman diagrams included in the 2d echo-spectrum calculation. Vertical axis represent the time propagation and the action of the dipole operators are shown by diagonal arrows. For a complete theoretical description see [55].
Sorting this expression for the rephasing and the non-rephasing parts yields six different pathways, which represent the stimulated emission, the ground-state bleach and the excited state absorption [54,55,56,57], shown schematically in figure 3. By taking the time between second and third pulse as a control parameter as done in (18) and (19) one can study the exciton dynamics which is probed by recording the emitted signal pulse. Due to the excited state absorption the dynamics is more involved compared to the one dimensional case. Including a double excited manifold, we construct an enlarged exciton Hamiltonian H 2d e which is block diagonal and of the form where H gs e = 0 represents the ground state, H 1ex e is the single-exciton Hamiltonian with N 2 sites elements as defined in (3), and the energy levels of the two excitons are described by H 2ex e . The two-exciton Hamiltonian H 2d e has [N sites (N sites − 1)/2] 2 entries which are constructed from the matrix elements of H 1ex e following [54,14]. For the FMO complex H 2d e is a 29 × 29 matrix. In site basis, we denote the two exciton states by |ij , where 0 < i < j ≤ N sites , and the diagonal entries of H 2ex e are given by The off-diagonal elements are given by Physically, this means that the |ij and |kl two-exciton states are coupled if they have one exciton in common while a second exciton is shifted from one site to another. The coupling constant is the same as for the transition of the corresponding exciton shift in the one exciton dynamics. Note also that the linear coupling to the independent baths, each of them associated to one BChl, leads to couplings between the one exciton state |i and the two exciton states |ij which do both influence the bath at site i and are in return influenced by the bath correlations. Following [31], we construct the dipole operator associated with the extended Hamiltonian (20) The zero to single excitation operators µ ± 1ex are defined as in (10) and (11). The single to double excitation operators are defined according to In order to calculate 2d spectra, we have to evolve the 29 × 29 density matrix to obtain the third order response function S (3) . The hierarchical method incorporates the effect of the baths through the auxiliary system of the σ-matrices. When calculating the 2d spectra, we propagate according to the Feynman diagrams shown in figure 3 and act on the density matrix and the auxiliary matrices with dipole operator to simulate the interaction with the laser field.
We observe in (17) that the 2d spectra depend on the fourth order dipole moments. To take into account the random orientation of the samples we need to compute rotational averages µ 2 i µ 2 j rot which are performed by sampling 10 laser polarisationvectors aligned to the vertices of a dodecahedron in one half of the coordinate space. We take the vertices (±1, ±1, 1), 0, ± 1 φ , φ , ± 1 φ , φ, 0 and ±φ, 0, 1 φ , where φ = 1 2 1 + √ 5 denotes the golden ratio. We have compared the dodecahedral summation with the full spherical averaging using a 10000-shoot MC and obtain differences of less than 1 percent of the MC results.

Numerical results for the 2d echo-spectrum of the FMO complex
For the simulation of 2d spectra, we calculate the third order response function S (3) (t 1 , t 2 , t 3 ) in (17). We obtain 128 2 data points with a spacing of 16 fs in t 1 and t 3 direction using a numerical time step of 2 fs. This means that we have to perform over 131 000 propagation steps for each of the six Feynman diagrams. Depending on the delay time T delay = t 2 our simulations take 48 to 53 minutes on a single NVIDIA Fermi C2050 GPU which provides a 38 fold speed up compared to the CPU algorithm with filtering methods that reduce the hierarchy [17]. We always show the sum of the real parts of the rephasing and the non-rephasing spectrum. Figure 4 (a) shows a simulated 2d spectrum for a rotationally averaged sample with no disorder at T = 77 K. The figure is encoded using a linear colour scale. We divided the raw data by a constant to obtain an amplitude on a −10-40 range. The thin white lines have been added to indicate the expected peak positions. The peaks are denoted by their excitonic energy-level indices (i, j), with the first number referring to the ω 3 -axis and the second one to the ω 1 -axis. The spectrum shows similarities to the experimentally measured one in [14] insofar as it reproduces the large peak on the diagonal at the second eigenstate and a second near the fourth and fifth. Furthermore we find negative amplitudes above the diagonal and positive peaks below. As expected from the discussion of the absorption spectrum in section 3.2, the peak shapes differ from the experimental ones due to the lack of static disorder. Figure 4 (b) shows the peak oscillations for the diagonal peaks (1, 1) and (5, 5) and the off-diagonal peak (1, 5) below the diagonal.We renormalised the peak heights at T delay = 0 fs to unity such that the time resolved amplitude of the three peaks are conveniently put into the same panel. We choose these peaks because figure 4 shows peak (1,5) as the highest off-diagonal peak in the line where ω 3 equals the lowest lying eigenstate. We see that diagonal and off-diagonal peaks oscillate in time. This indicates that our model imitates the coherent energy transport within the FMO complex that has been observed experimentally [5,6]. Theoretically, an off-diagonal peak (i, j) is expected to undergo oscillations with a frequency corresponding to the difference of its excitonic eigenenergies (E j − E i ). For the (1, 5) peak we expect a period of 99 fs which is in good agreement to our numerical result yielding a best-fit period of 97 fs. When we calculate the peak heights, we integrate over squares of width ∆ω = 16 cm −1 centred around the white crosses. We follow the changes of the peak heights with increasing delay time by plotting the absolute value of the spectrum normalised with respect to the value at T delay = 0 fs. In order to analyse the robustness of the peak oscillations when the thermal fluctuations increase we plot in figure 5 (a) the oscillations of the (1, 5) cross peak for several temperatures. We compare with the experimental data in [6] and observe a qualitative similar behaviour with weaker oscillations and faster decay when temperature is increased. We do not obtain the picoseconds lasting beatings that have been observed experimentally. One possible factor is the shape of the spectral density that is used for the HEOM. We fit the oscillations in figure 5 (a) by an exponentially decaying sine. The fit is done separately for the stimulated emission and excited state absorption pathway, which couple to one or two baths respectively. Based on the damping of the oscillations on peak (1, 5) the temperature dependence of the dephasing rate of our simulation is given by γ ESA deph (T )/T = (0.66 ± 0.05) cm −1 /K and γ SE deph (T )/T = (0.50 ± 0.03) cm −1 /K. This is in good agreement to what we expect for the chosen spectral density (2k B λ/(γT ) = 0.46 cm −1 /K) and the value 0.52 cm −1 /K determined from the experimentally measured data [6].
In the experiments an anti-correlation in peak shapes between the diagonal peak width and the height [5], or alternatively the ratio of the diagonal peak-width over the anti-diagonal height [4] has been found. The width was measured by fitting a Gaussian to a diagonal cut through the spectrum and taking its standard deviation. Applying the same procedure to the simulation data yields also anti-correlations shown in figure 5 (b).
Next, we analyse the role of static disorder in the results for the 2d echo-spectra. When including static disorder, the peaks get elongated along the diagonal as observed when comparing the spectra at zero delay time shown in figures 4 (a) and 6 (a) and (b) for increasing values of disorder. The plots with disorder show better agreement with the experiment [14]. Figure 6 (c) shows the 2d spectra for SD= 25 cm −1 for a long delay time T delay = 4000 fs for which the system approaches its thermal stationary state. We observe that the excitation moves from the diagonal peaks into the regions with lower frequency ω 3 . This occurs as the system relaxes to lower lying eigenstates during the delay time. In experiments by Brixner et al [9] the decay of the population into the ground state of the system is also visible. The absence of the experimentally seen oscillation in the HEOM calculation in Ref. [17] has led Chen et al to conclude that the importance of static disorder requires careful reconsideration. First-principles simulations of static disorder in the FMO complex are computationally extremely time-consuming, since several hundred of realisations are required to reach converged results. In order to study the role of static disorder in the oscillations of the peaks of the spectrum and to ensure numerical convergence, we construct a dimer to model the (1, 5) peak of the FMO spectrum. For the model dimer we chose the Hamiltonian such that the first and fifth FMO complex eigenvalues are reproduced. The coupling between both sites is set to 50 cm −1 . The simulations with and without disorder show oscillations on the diagonal and off-diagonal peaks. The most striking effect of disorder is a decreased oscillatory amplitude. While the cross peak in the simulation without disorder has lost 45 percent of its amplitude in the first minimum, in the disorder case the loss is only 35 percent. From a fit, we observe that the periods of the disorder curves are slightly shorter T = 92 fs and T = 98 fs for the (1, 1) and (2, 1) curves while in the no-disorder plot T = 95 fs and T = 98 fs has been obtained, respectively. For the FMO complex we expect that oscillations are reduced in amplitude by a similar magnitude and therefore still persist within the model.

Conclusions
We have presented calculations of the absorption spectra and the 2d echo-spectra for the FMO complex at different temperatures using the HEOM for solving the exciton and the vibrational bath dynamics and adding static disorder to account for the fluctuating environment. Important ingredients necessary to find good agreement with experiments for the absorption spectra using the HEOM propagation method is to perform the rotational average and to include static disorder. We have also discussed how other simplified propagation methods affect the line shapes and shift the peaks in the spectra. Our rotationally averaged 2d echo-spectra reflect the main features of the experimentally measured ones. At delay time T delay = 0 we find negative contributions above the diagonal and elongated peaks in the diagonal which are transferred to the lowest exciton states after long times on the picosecond scale. The HEOM method yields an oscillatory behaviour of peak amplitudes over a wide temperature range. At T = 77 K we obtain oscillations of cross peak amplitudes up to 300 fs with the expected frequency that corresponds to the difference in excitonic eigenenergies. In addition oscillations on the diagonal peaks are present. The long duration of oscillations in the experimental data (up to 1.8 ps in Ref. [6]) suggests that further refinements of the theoretical models are necessary, for example by adjusting the form of the spectral density which affects the population dynamics [38].
Remarkably, we find excellent agreement with the measured damping rate of the oscillations due to thermal fluctuations. This implies a good choice of the parameters of the model and might explain why we obtain different results compared to previous calculations, which were done for a higher value of γ −1 = 100 fs [17], leading to an increased dephasing by a factor of two. At ambient temperatures the HEOM have the advantage to require less terms to converge while approximate methods such as the secular and full Redfield approximations do not yield the correct time-evolution of the density matrix and overestimate the thermalisation rate.
In agreement with the experimental measurements we observe the anti-correlated oscillations of the peak heights and peak shapes for the diagonal peaks. The analysis of a model dimer shows that the cross peak oscillations persist even in the presence of static disorder which mainly results in a reduction of oscillation amplitude by about 30 percent.