Signatures of spatially correlated noise and non-secular effects
in two-dimensional electronic spectroscopy
James Lim,1, a) David J. Ing,2, a) Joachim Rosskopf,1 Jan Jeske,2 Jared H. Cole,2 Susana F. Huelga,1, b) and Martin B. Plenio1, c) 1)Institut fur Theoretische Physik, Albert-Einstein-Allee 11, Universitat Ulm, D-89069 Ulm, Germany 2)Chemical and Quantum Physics, School of Applied Sciences, RMIT University, Melbourne, Victoria 3001, Australia
We investigate how correlated fluctuations affect oscillatory features in rephasing and non-rephasing twodimensional (2D) electronic spectra of a model dimer system. Based on a beating map analysis, we show that non-secular environmental couplings induced by uncorrelated fluctuations lead to oscillations centered at both cross- and diagonal-peaks in rephasing spectra as well as in non-rephasing spectra. Using an analytical approach, we provide a quantitative description of the non-secular effects in terms of the Feynman diagrams and show that the environment-induced mixing of different inter-excitonic coherences leads to oscillations in the rephasing diagonal-peaks and non-rephasing cross-peaks. We demonstrate that as correlations in the noise increase, the lifetime of oscillatory 2D signals is enhanced at rephasing cross-peaks and non-rephasing diagonal-peaks, while the other non-secular oscillatory signals are suppressed. We discuss that the asymmetry of 2D lineshapes in the beating map provides information on the degree of correlations in environmental fluctuations. Finally we investigate how the oscillatory features in 2D spectra are affected by inhomogeneous broadening.

arXiv:1701.00045v1 [quant-ph] 31 Dec 2016

I. INTRODUCTION
In the first step of the light-harvesting process1,2, the neutral electronic excitations (excitons) created by light absorption are transferred through the molecular system until free charge carriers are generated by exciton dissociation1,2. The interaction between electronic and vibrational degrees of freedom governs the exciton transfer dynamics, such as coherent and incoherent features in energy transport35.
Two-dimensional electronic spectroscopy (2DES) has been employed to study the exciton transfer dynamics in the light-harvesting systems on a sub-picosecond timescale6. For various natural712 and artificial1320 systems, oscillatory signals observed in 2D experiments have been interpreted as a signature of quantum coherences within the molecular system, coherence that is generated by laser pulses but sustained by the intrinsic dynamics. These coherences can originate in principle from both electronic and vibrational degrees of freedom. The electronic states of the light-harvesting systems are coupled to their vibrational environments with a moderate coupling strength, such that electronic coherences are not completely washed out by the environmentinduced noise. The intra-pigment vibrations of the lightharvesting systems exhibit underdamped vibrational motions on a picosecond timescale, which can leads to vibrational coherences in the electronic ground state manifold21,22. A vibronic (electronic-vibrational) coupling between electronic states and underdamped vibrational
a)These authors contributed equally to this work. b)Electronic mail: susana.huelga@uni-ulm.de c)Electronic mail: martin.plenio@uni-ulm.de

motions leads to a mixing of electronic and vibrational degrees of freedom, inducing vibronic coherences in the electronic excited state manifold2230. A model dimer system has been widely employed to investigate how electronic-vibrational interactions are reflected in spectroscopy2733.
For various light-harvesting systems, 2DES has demonstrated the presence of long-lived quantum coherences, which are sustained beyond the lifetime of optical coherences between electronic ground and excited states717. To identify the microscopic origin of the long-lived oscillatory signals observed in 2D experiments, several hypotheses have been formulated to explain how quantum coherences are sustained under a noisy environment at ambient conditions. Coherent vibronic coupling has been shown to induce long-lived vibronic coherences when the vibrational frequency is resonant with the energy-level difference between exciton states2228 and the vibronic coupling strength is moderate in magnitude17,27. In this case, oscillatory 2D signals originate from the combination of excited-state vibronic coherence and ground-state vibrational coherence, as they originate from the same mechanism, namely the vibronic Hamiltonian. Recent 2D experiments on J-aggregates of cyanine dyes17 confirmed the validity of this theoretical approach. Other experimental results on a synthetic dimer15,34 pointed towards the need of a threshold in the vibronic coupling strength to allow vibronic mixing to be relevant under environmental effects.
Correlated fluctuations in the transition energies of neighboring pigments have been suggested as an alternative mechanism where purely electronic coherences in the excited state manifold induce long-lived oscillatory 2D signals3,35,36. In the correlated fluctuation model, where underdamped vibrational motions are not consid-

2

ered, uncorrelated noise between electronic ground and excited states leads to the finite homogeneous broadening of 2D spectra, while highly correlated noise between different excited states induces long-lived inter-excitonic coherences. This purely electronic coherence model is similar in spirit to the decoherence-free subspaces in quantum information science37,38. Correlated fluctuations have been shown to enhance the lifetime of coherent features in electronic motions, such as population dynamics36,3943. Correlated fluctuations have also been considered in the simulations of 2D electronic spectroscopy44,45. Using the secular Redfield theory, it was found that uncorrelated and correlated fluctuation models give very similar absorption spectra and rephasing diagonal-peaks, but show large differences in rephasing cross-peak regimes44. The enhancement of the lifetime of 2D oscillations was observed for rephasing diagonal- and cross-peaks, but the diagonal-peak oscillations were attributed to the overlap of non-oscillatory diagonal-peaks with nearby oscillating cross-peaks44. It was found that correlations in the noise can alter the relative phase between diagonaland cross-peak oscillations45. Correlations in inhomogeneous broadening have also been considered in the context of 2DES where correlated and anti-correlated static disorder lead to cross-peaks strongly elongated along diagonal and anti-diagonal directions, respectively46. A similar feature was observed in the simulations of 2D infrared (IR) spectroscopy in the slow bath limit, while different correlation models lead to similar Lorentzian 2D lineshapes in the fast bath limit47. However, it is unclear what the microscopic origin of such correlations is. Quantum mechanics/molecular mechanics (QM/MM) simulations of photosynthetic systems, such as the FennaMatthews-Olson (FMO) complex48,49 and phycoerythrin 545 (PE545) complex50, for instance, have shown no evidence of spatially correlated fluctuations. This is contrary to the 2D experiments on the FMO complex9,10, colloidal semiconductor nanoplatelets18 and J-aggregates of porphyrins19, where the presence of correlated fluctuations was suggested. The discrepancy between theory and experiment shows the need for further investigations of how the degree of correlations in the noise is reflected in experimental observables, so that the presence of correlated fluctuations can be verified or ruled out based upon experimental results.
In this work, we employ the Bloch-Redfield equation36,51 where the degree of correlations in the noise is parameterized by a continuous variable to investigate how the correlations affect the oscillatory features in rephasing and non-rephasing 2D spectra of a model dimer system. We do not include underdamped vibrational modes in the model to focus on the influence of spatial noise correlations on optical responses without vibronic effects. The employed Bloch-Redfield equation includes non-secular environmental couplings between exciton populations and inter-excitonic coherences, and also the non-secular interaction between different interexcitonic coherences. These non-secular terms are disre-

garded in the secular approximation, where the dynamics of exciton populations are decoupled from those of interexcitonic coherences. The secular terms describe the relaxation between exciton populations and the decay of inter-excitonic coherences, independently.
We show that in the absence of noise correlations, non-secular environmental couplings induce oscillatory 2D signals centered at both cross- and diagonal-peaks in rephasing spectra as well as in non-rephasing spectra. In Ref. 52, it was found that rephasing diagonalpeaks show notable oscillatory features when non-secular couplings are taken into account in simulations. The origin of the oscillations was attributed to the nonsecular interaction between exciton populations and coherences52,53. It was also shown that the diagonal-peak oscillations are suppressed when the secular approximation is employed, even though both diagonal- and crosspeaks showed similar amplitudes of oscillations52. In this work, to clarify the contribution of non-secular effects to rephasing diagonal- and non-rephrasing crosspeak oscillations (non-secular oscillations) in the presence of spectral overlap of diagonal- and cross-peaks, we provide a quantitative beating map analysis in terms of the eigenstates of the Liouville space operator with associated Feynman diagrams. The beating map analysis helps to clarify whether the oscillation of a given 2D peak originates from itself or merely from the overlap with nearby oscillating peaks. For a homodimer, we show that the non-secular oscillations can be induced by the nonsecular coupling between different inter-excitonic coherences, even if their dynamics are decoupled from those of exciton populations. For a heterodimer, where the dynamics of exciton populations and coherences are all coupled to one another, we show that the non-secular oscillations are dominated by the mixing of different interexcitonic coherences, rather than population-coherence mixing suggested in Refs. 52 and 53. In addition, we show that the uncorrelated noise can induce asymmetric lineshapes of homogeneously broadened 2D peaks in the beating map, elongated along excitation or detection axis, when the broadening is dominated by relaxation, rather than pure dephasing noise. For the FMO complex, small pure dephasing rates have been identified as a condition for the oscillations of rephasing cross-peaks54. Finally we show that as the degree of correlations in the noise increases, the lineshape of absorption and 2D electronic spectra, including both diagonal- and cross-peaks, is significantly changed, as the homogeneous broadening starts to be dominated by pure dephasing, rather than relaxation. We show that the noise correlations induce long-lasting oscillations in the rephasing cross-peaks and non-rephasing diagonal-peaks with suppressed nonsecular oscillations, leading to symmetric 2D lineshapes in the beating map. We discuss that the asymmetry of 2D lineshapes can provide information on the degree of correlated fluctuations in J-aggregates17 and colloidal semiconductor nanoplatelets18, for which asymmetric lineshapes elongated along the excitation axis were observed

3

(a) 1b11b11 

1 b21 b21 

.....

.....

J

1

2





k b1k b1k

k b2k b2k

(b) 1b11b11

1 b21 b21







.....

.....

J

1

2





 k b1k b1k

k b2k b2k

FIG. 1. A schematic representation of (a) local phonon baths and (b) a shared phonon bath. In (a), a phonon mode b1k is locally coupled to site 1. In (b), a phonon mode b1k is coupled to both sites 1 and 2 with the relative electron-phonon coupling strengths quantified by  and , respectively.
in 2D experiments.

II. THE MODEL

To investigate how correlations in the noise are reflected in two-dimensional electronic spectroscopy, we consider a dimer system consisting of two sites coherently coupled by an electronic coupling J. The Hamiltonian of the dimer system is modeled by

He = 11+1- + 22+2- + J (1+2- + 1-2+), (1)

with |g1

k denoting e1|  112 and

the energy 2- = 111 

level of |g2 e2|

the site k, 1- = represent the an-

nihilation operators of the electronic excitation at sites

1 and 2, respectively, with |gk and |ek denoting the

ground and excited states of the site k, respectively, with

11k = |gk gk| + |ek ek|. The electronic eigenstates of

He are expressed as

|g = |g1, g2 ,

(2)

| 1 = - sin() |e1, g2 + cos() |g1, e2 ,

(3)

| 2 = cos() |e1, g2 + sin() |g1, e2 ,

(4)

|f = |e1, e2 ,

(5)

with



=

1 2

tan-1(2J/(1

- 2))

for

1



2.

The

asso-

ciated eigenvalues are given by

g = 0,

(6)

1

=

1 2

1 + 2 -

(1 - 2)2 + 4J 2 ,

(7)

2

=

1 2

1 + 2 +

(1 - 2)2 + 4J 2 ,

(8)

f = 1 + 2,

(9)

where 1 < 2 and the bi-exciton binding energy is not considered for the sake of simplicity, leading to f = 1 + 2. Here |g represents a common ground state, | 1 and | 2 are low and high energy single exciton states, respectively, and |f is a bi-exciton state.
The phonon environment coupled to the electronic
states is modeled by independent harmonic oscillators

Hp = ( kb1kb1k + kb2kb2k),
k

(10)

where the dephasing interaction between electronic states and phonon environment is chosen to be of the form

He-p = 1+1- 

gk((b1k + b1k) + (b2k + b2k))

k

+ 2+2- 

gk((b2k + b2k) + (b1k + b1k)),

k
(11)

where the phonon mode b1k (or b2k) of frequency k is coupled to both sites 1 and 2 with relative coupling strengths quantified by dimensionless scaling fac-
tors 0    1 and   1 - 2 (or  and ), respectively. The electron-phonon couplings gk are modeled by a shifted Ohmic spectral density, as will be discussed below. When  = 1, leading to  = 0, the phonon environment is reduced to local phonon baths, where the phonon modes bjk are locally coupled to site j, inducing spatially uncorrelated noise. When  = 0 and  = 0, the phonon environment is reduced to a shared phonon bath, as the phonon modes are coupled to both sites 1 and 2, leading to spatially correlated noise. A schematic representation of the local and shared phonon baths is displayed in Fig. 1. The degree of correlations in the noise is quantified by a correlation length  defined by exp(-d/) = 2  [0, 1] with d denoting the distance between sites 1 and 2.
In this work, we employ the Bloch-Redfield formalism36,51 to describe the dynamics of the dimer system based on the Hamiltonian above. This formalism

4

is well suited for describing the effect of correlated fluctuations36, including the existence of partial cor-
relations, as summarized in Appendix A. The phonon
environment is modeled by a shifted Ohmic spectral density54

J ()

=

 

2

+

 ( -

s)2

+

2

+

 ( +

s)2

, (12)

where  denotes the reorganization energy, defined by

=

 0

dJ

(

)-1,

and



is

the

bath

relaxation

rate. For the FMO complex, small pure dephasing rates

have been identified as a condition for the oscillations of

rephasing cross-peaks54. In this work, we take the shift

s of the phonon spectral density to be resonant with the exciton splitting, i.e. s = | 2 - 1|, so that the homogeneous broadening is dominated by relaxation, rather than

pure dephasing, as discussed in Appendix A. In simula-

tions, we take a fast bath relaxation rate of  = (50 fs)-1,

corresponding to a broad spectral density, to avoid vi-

bronic effects induced by underdamped modes, such as

vibronic progressions or mixing in absorption55 and 2D

spectra. With the Bloch-Redfield equation, we simulate

2D electronic spectra in the impulsive limit with the as-

sumption that the transition dipoles of sites 1 and 2 are

mutually orthogonal, as discussed in Appendix B.

III. RESULTS
In this section, we provide a beating map analysis of rephasing and non-rephasing spectra of a model dimer system to demonstrate how non-secular environmental couplings and correlated fluctuations affect the oscillatory features in 2D spectra. We will consider two cases in the simulations: a homodimer, where both sites 1 and 2 have the same site energy 1,2 = 12500 cm-1, and a heterodimer, where the two sites have different site energies, 1 = 12600 cm-1 and 2 = 12400 cm-1. For both cases, the electronic coupling between sites is taken to be J = 100 cm-1, and the phonon bath is modeled by typical values encountered in natural photosynthetic systems3:  = 50 cm-1,  = (50 fs)-1 and T = 77 K. We will also show how the beating map is changed by static disorder. The numerical results, which will be provided in this section, are investigated analytically in Appendix D in full detail.
Before we provide a beating map analysis, we demonstrate in Fig. 2 that correlated fluctuations can modify absorption and 2D lineshapes when homogeneous broadening is dominated by relaxation, and induce long-lived 2D oscillations within the Bloch-Redfield formalism. For the homodimer, Fig. 2(a) shows the absorption spectrum when correlations in the noise are negligible (i.e. local phonon baths). The high-energy absorption peak centered at 2 = 1.26  104 cm-1 is broader than the lowenergy peak centered at 1 = 1.24  104 cm-1. This is due to the fact that the homogeneous broadening is dominated by relaxation in the model, where the broadening

of the low-energy peak is dominated by pure dephasing only, while the high-energy peak is broadened by both relaxation and pure dephasing. A larger broadening of the high-energy peak makes its amplitude smaller than the low-energy peak. As shown in Figs. 2(b) and (c), the amplitude and homogeneous broadening of two peaks become similar as correlations in the noise increase, where the broadening is dominated by pure dephasing with suppressed relaxation.
In Fig. 2(d), the real part of the rephasing spectra at waiting time t2 = 0 is displayed for the case that correlations in the noise are absent. The high-energy diagonal peak R22 (cf. four peaks marked by black circles in Fig. 2(d)) is hardly visible, as expected from the small amplitude of the high-energy absorption peak shown in Fig. 2(a). The amplitude of upper-diagonal cross-peak R12 is larger than lower-diagonal cross-peak R21. This is due to the excited state absorption (ESA) signals. In Fig. 2(j), the sum of ground state bleaching (GSB) and stimulated emission (SE) contributions to the rephasing spectra is displayed where the main peak is R11, as its amplitude and broadening along excitation and detection axes are governed by the lineshape of the low-energy absorption peak (cf. Feynman diagrams in Fig. 2(j)). In Fig. 2(k), the ESA contribution is displayed where the main peak is R12 centered at detection frequency of 3 = 2. Here the coherence |f 1| between bi-exciton and low-energy exciton state leads to the ESA peak centered at 3 = f - 1 = 2, as described in the Feynman diagram in Fig. 2(k), leading to a large amplitude and a narrow linewidth along 3-axis, similar to the low-energy absorption peak. Figs. 2(e) and (f) show that the rephasing lineshape becomes more symmetric and high-energy diagonal peak R22 starts to have a larger amplitude as correlations in the noise increase.
Figs. 2(g)-(i) show how the dynamics of lower-diagonal cross-peak R21 are affected by the degree of correlations in the noise. For uncorrelated noise, the cross-peak R21 shows oscillatory dynamics up to t2  300 fs, as shown in Fig. 2(g). As the correlation length  increases, the lifetime of the oscillations in R21 is increased as shown in Figs. 2(h) and (i), describing partially and fully correlated noise, respectively. These results demonstrate that correlations in the noise can enhance the lifetime of excited state coherences, as expected from Refs. 36, 39 44, leading to persistent oscillatory 2D signals when /d  . The advantage of the current formalism is that we can tune  to cover the two extreme cases of completely uncorrelated and perfectly correlated noise.
To investigate in detail how correlations in the noise affect the oscillatory 2D signals, we use a beating map analysis, which visualizes the lineshape of oscillatory 2D signals in the (1, 3) domain as a function of the beating frequency 2. To this end, we extract oscillatory components from the total 2D spectra that contain both damped oscillations and non-oscillatory components. The non-oscillatory components include exponential and static t2-transients (cf. Figs. 2(g) and (h)).

5

Uncorrelated fluctuations

Partially correlated fluctuations

Fully correlated fluctuations

(a)

(b)

(c)

(d)

(e)

(f)

R12

R22

R11

R21

(g)

(h)

(i)

R21

(j)

(k)

R11 (GSB) R11 (SE)

|g g| |g g|

GSB+SE

|1 g | |g g| | g 1| |g g|

|1 g | |1 1| | g 1| |g g|

R12 (ESA)

|1 1|

| f 1|

|1 1|

| g 1|

ESA

|g g|

FIG. 2. Absorption (Abs.), the real part of rephasing 2D spectra at waiting time t2 = 0 and the t2-transient of the lowerdiagonal cross-peak R21 centered at (1, 3) = ( 2, 1) with 1 = 1.24  104 cm-1 and 2 = 1.26  104 cm-1. In (a), (d), (g), (j), (k), we consider local phonon baths characterized by a short correlation length  = 10-3d, leading to e-d/  0. In (b),
(e), (h), we consider an intermediate case where  = 3d, leading to e-d/  0.7. In (c), (f), (i), we consider a shared phonon
bath characterized by a long correlation length  = 103d, leading to e-d/  1. In (j), the sum of GSB and SE contributions
to 2D spectra shown in (d) is displayed with the Feynman diagrams responsible for the main peak R11. In (k), the ESA
contribution to (d) is displayed with the Feynman diagram for the main peak R12. The ESA contribution makes R12 stronger than R21 in both (d) and (e). Here we employed 1 = 2 = 12500 cm-1, J = 100 cm-1,  = 50 cm-1,  = (50 fs)-1 (cf.   106 cm-1), s = 200 cm-1 and T = 77 K.

In 2D experiments, the oscillatory components are ex-
tracted from raw 2D spectra by fitting multi-exponentials
to the raw t2-transients for each (1, 3) value, or by fitting 2D decay-associated spectra (2DDAS) to the raw 2D data1113,17. In this work, we directly calculate the oscil-
latory components by removing time-evolution operator
components leading to non-oscillatory 2D signals, which
will be detailed in Appendix D. By avoiding the fitting

procedure in simulations, one can avoid potential artefacts and numerical errors in the beating map. Throughout this work, the response function that only contains the oscillatory components is denoted by S(1, t2, 3), while the total response function that contains both oscillatory and non-oscillatory signals is represented by S(1, t2, 3). The oscillatory component S(1, t2, 3) is generally expressed as a sum of complex-valued damped

6

(a) S(1, t2, 3), t2 = 0

(b)

R21

(c) (e) S(1, 2, 3), 2 = 200 cm-1
(d) R21
FT

R21 S = S(1, t2, 3)
S-S R21
S = S(1, t2, 3) R21

FIG. 3. A schematic representation of beating map calculation. In (a), 2D spectra S(1, t2, 3) at t2 = 0 are displayed (cf. Fig. 2(d)). In (b), t2-transient of the cross-peak R21 is shown (cf. Fig. 2(g)). The transient consists of (c) nonoscillatory component S - S, including exponential and static t2-transients, and (d) oscillatory component S = S(1, t2, 3) (see text). By extracting the oscillatory components S from the raw 2D spectra S for each (1, 3) value and Fourier transforming S with respect to t2, one can obtain the beating map in the (1, 2, 3) domain, as shown in (e), where the beating frequency 2 is taken to be the exciton splitting of | 2 - 1|. In (a)-(d), we display the real part of S for the sake of simplicity, but the imaginary part of S is also included in the computation of the beating map S(1, 2, 3) (see text). We note that in this work, S0.1 is displayed (cf. Fig. 4(b)), instead of S (cf. (e) and Fig. 5(a)), to make small amplitudes more visible.

oscillations, i.e. S(1, t2, 3) = k Ak(1, 3)e(ivk-k)t2 with frequencies vk and associated damping rates k. We evaluate the beating map by Fourier transforming
S(1, t2, 3) with respect to the waiting time t2



S(1, 2, 3) =

dt2S(1, t2, 3) exp(-i2t2) ,

0

(13)

where 2 is the beating frequency. Here we consider a complex-valued response function S(1, t2, 3), rather than only its real or imaginary part, so that we retain the

full information of the oscillatory signals. In this way,

we can distinguish positive and negative frequency com-

ponents that oscillate in the form of exp(i |v| t2 - t2) and exp(-i |v| t2 - t2), respectively, with an overall decay rate of . The positive and negative frequency

components are reflected in the beating map as the

Lorentzian functions centered at 2 = |v| and 2 = - |v|, respectively, with a width of  along 2-axis. A schematic representation of the beating map evaluation

is shown in Fig. 3. A separate analysis of positive and negative frequency components has been employed to distinguish electronic and vibrational coherences for a model quantum dot system56 and experimentally estimate the Hamiltonian and decoherence rates of an atomic vapour57.
In Fig. 4, we show the resulting beating map of the rephasing spectra of a homodimer with the parameters used in Fig. 2. Figs. 4(a)-(c) show the case of local phonon baths considered in Figs. 2(d) and (g). In Fig. 4(a), the beating map at a negative frequency of 2 = - | 2 - 1| is displayed, which is dominated by the upper-diagonal cross-peak R12. This is due to the interexcitonic coherence in the form of | 2 1|, where | 1 and | 2 denote lower and higher energy single exciton states, respectively. The inter-excitonic coherence leads to the negative frequency component, as 2 > 1. It is notable that there are weak diagonal-peaks centered at R11 and R22, which are not artefacts of the beating map calculations. Here the maximum value of S  S(1, 2, 3) is normalized to 1, i.e. 0  S  1, and S0.1 is displayed instead of S, so that the small amplitudes are more visible in the beating map. In Fig. 4(b), the beating map at a positive frequency of 2 = | 2 - 1| is shown, where the lower-diagonal cross-peak R21 is induced by the inter-excitonic coherence in the form of | 1 2|. Interestingly, all the other peaks R11, R12 and R22 are visible in Fig. 4(b) and the amplitude of R11 is comparable to that of R21.
To understand the lineshape of oscillatory signals in more detail, Fig. 4(c) displays the amplitudes of the cross- and diagonal-peaks as a function of the beating frequency 2. The cross-peaks R12 and R21 are centered at the negative and positive beating frequencies, respectively. The large amplitude and broadening of the R12 peak explains the reason why R12 is visible in both Figs. 4(a) and (b). Interestingly, the diagonal-peak R11 has comparable amplitudes at both positive and negative frequencies, contrary to the cross-peaks R12 and R21. We note that the diagonal-peak R11 in the beating map does not originate from the overlap of the homogeneously broadened cross-peaks R12 and R21. As shown in Fig. 5(a), where S is displayed instead of S0.1, the homogeneous broadening of the cross-peaks is not large enough to dominate the amplitude of the diagonalpeak R11. More specifically, the distance between R11 and R21 is the same to that between R21 and point A, marked by a purple circle. Since the regime around point A has no overlap with other peaks, the amplitude of the homogeneously broadened cross-peak R21 at the position of R11 can be approximately estimated by the value of S at point A (cf. Fig. 5(b)). Similarly, the amplitude of the homogeneously broadened cross-peak R12 at the position of R11 can be approximately estimated by the value of S at point B, marked by a light blue circle (cf. Fig. 5(c)). The contribution of the cross-peaks R12 and R21 to the amplitude of R11 is more than two times smaller than the amplitude of R11, implying that the diagonal-peak does

7

Uncorrelated fluctuations

Partially correlated fluctuations

Fully correlated fluctuations

(a)

neg.

(d)

neg.

(g)

neg.

R12

R11

R21

S0.1

(b)

(e)

(h)

pos.

pos.

pos.

(c)

S

R12

(f) R12

(i)

R12

R21

R21

S R11

R21 R22

FIG. 4. The beating map of complex-valued rephasing spectra that visualizes the lineshape of oscillatory 2D signals at a
frequency of 2. In (a)-(c), we consider local phonon baths with the parameters used in Fig. 2(a). In (a), the lineshape of oscillatory signals in the form of exp(i2t2 - t2) is displayed with a negative frequency of 2 = - | 2 - 1| = -200 cm-1 and an overall decay rate of . Here the maximum value of S  S(1, 2, 3) is normalized to 1, i.e. 0  S  1, and S0.1 is displayed instead of S, so that the small amplitudes are more visible in the beating map. In (b), the lineshape of oscillatory signals with a positive frequency of 2 = | 2 - 1| = 200 cm-1 is displayed. It is notable that the lineshape of R21 is asymmetric with a larger homogeneous broadening along 1-axis when compared to the broadening along 3-axis. In both (a) and (b), there are oscillations centered at the lower diagonal-peak R11. In (c), the amplitudes of the cross-peaks R12 and
R21 and diagonal-peaks R11 and R22 in the beating map are displayed as a function of the beating frequency 2. Here S is displayed instead of S0.1. The cross-peaks R12 and R21 are centered at negative and positive beating frequencies, respectively.
On the other hand, the diagonal-peak R11 has comparable amplitudes at both positive and negative frequencies. In (d)-(f), we
consider an intermediate case with the parameters used in Fig. 2(b). In (d) and (e), the lineshape of the cross-peaks R12 and
R21 is more symmetric and the diagonal-peak R11 is less visible when compared to the case of the uncorrelated noise shown in
(a)-(c). In (g)-(i), we consider a shared phonon bath with the parameters used in Fig. 2(c). In (g) and (h), the diagonal-peaks
R11 and R22 are not visible, and the lineshape of the cross-peaks R12 and R21 is symmetric along the 1- and 3-axes. In (i), all the peaks have very narrow linewidths along 2-axis, as the overall decay rate  of the oscillatory 2D signals is very low due to the highly correlated noise (cf. Fig. 2(i)). In (i), the diagonal-peaks R11 and R22 have very small amplitudes, but this
is due to the homogeneous broadening of the cross-peaks R12 and R21 along 1- and 3-axes, as the diagonal-peaks are not visible in (g) and (h).

8

(a)

B

pos.

0.12

R12 R21 A
R11
S

0.1 0.08 0.06 0.04 0.02 0

(b)

S

R11

R21 0.033
A

(c)

S

R11

0.106

R12
0.011 B

FIG. 5. The cross sections of the rephasing beating map shown in Fig. 4(b). In (a), S is displayed instead of S0.1 (cf. Fig. 4(b)). To demonstrate that the oscillations in the diagonal-peak R11 do not originate from the overlap of the homogeneous broadening of the cross-peaks R12 and R21, in (b), S is shown as a function of 1 for 3 = 1 = 1.24  104 cm-1, while in (c), S is displayed as a function of 3 for 1 = 1 = 1.24  104 cm-1. The sum of the values of S at points A and B is  0.044, which is more than two times smaller than the value of R11 ( 0.106). This implies that if one assumes that the cross-peaks R12 and R21 have the Lorentzian lineshapes, the amplitude of the diagonal-peak R11 cannot be explained by the overlap of the homogeneously broadened cross-peaks R12 and R21.
not originate solely from the overlap of the cross-peaks.
So far we have analyzed the beating map for the case that correlations in the noise are absent (cf. Figs. 4(a)(c)). We now show how correlations in the noise change features of the beating map (cf. Figs. 4(d)-(i)). Figs. 4(d)-(f) show the beating map in the presence of partially correlated noise (cf. Figs. 2(e) and (h)), while Figs. 4(g)-(i) display the case of fully correlated noise (cf. Figs. 2(f) and (i)). It is notable that the overall 2D lineshapes become more symmetric as correlations in the noise increase. For instance, the asymmetric lineshape of the cross-peak R21 elongated along 1-axis becomes more symmetric as the correlation length  increases, as shown in Figs. 4(b), (e) and (h). Note also that the amplitude of the diagonal peak R11 is suppressed as the correlation length  increases. Indeed, R11 is not visible at all for fully correlated noise, as shown in Figs. 4(g) and (h). These results demonstrate that uncorrelated noise

can induce oscillations in the rephasing diagonal-peaks and make the lineshapes of the rephasing cross-peaks asymmetric in the beating map. Conversely, correlations in the noise suppress these features, leading to symmetric lineshapes of the rephasing cross-peaks in the beating map with suppressed diagonal oscillations.
In Fig. 6, we now show that the characteristics of the rephasing beating map of a heterodimer is similar to that of the homodimer, shown in Fig. 4. We also show that the qualitative features of non-rephasing beating maps are significantly affected by correlations in the noise, as is the case for rephasing beating maps. Fig. 6(a) shows rephasing beating maps of the heterodimer in the absence of correlations in the noise at negative and positive 1.22 1.24 1.26 1.28 frequencies 2 =  | 2 - 1|. Here the oscillations occur at both cross-peaks R12 and R21 as well as at the diagonal-peaks R11 and R22. Note that the amplitude of the upper diagonal-peak R22 is more visible when compared to the case of the homodimer shown in Fig. 4(b). Fig. 6(b) shows the rephasing beating maps in the presence of highly correlated noise, where the oscillatory 2D signals occur only at the cross-peaks, as in the case of the homodimer. Figs. 6(c) and (d) show the nonrephasing beating maps of the homodimer considered in Fig. 4. In the absence of correlations in the noise, oscillatory non-rephasing signals occur at both diagonal-peaks N11 and N22 as well as at the cross-peaks N12 and N21, as shown in Fig. 6(c). In the presence of highly correlated noise, the oscillations occur only at the diagonalpeaks N11 and N22, as shown in Fig. 6(d). Note that the asymmetric lineshape of the diagonal-peak N22 becomes more symmetric as correlations in the noise increase. The non-rephasing beating map of a heterodimer shows similar features, as demonstrated in Figs. 6(e) and (f). In Fig. 6(e), the diagonal-peak N22 shows a seemingly discontinuous lineshape due to the interference of the oscillatory signals from the SE and ESA contributions, where each contribution leads to continuous 2D lineshapes in the beating map (not shown here).
These results demonstrate that in the absence of correlations in the noise, the oscillations in rephasing and non-rephasing spectra can appear at both cross- and diagonal-peaks with asymmetric lineshapes in the beating map. These features are suppressed as correlations in the noise increase, leading to oscillatory signals centered only at rephasing cross-peaks and non-rephasing diagonal-peaks with symmetric lineshapes. This is contrary to the simulated 2D spectra based on the Redfield equation within the secular approximation, where the oscillations occur only at rephasing cross-peaks and nonrephasing diagonal-peaks, even in the absence of correlations in the noise21. This suggests that the non-secular environmental couplings in the Bloch-Redfield equation, which couple the dynamics of an inter-excitonic coherence to that of the other inter-excitonic coherences and exciton populations36, may be responsible for our observations, as suggested in Refs. 52 and 53 for rephasing spectra without a quantitative description. We note

Uncorrelated fluctuations

9

Rephasing / Heterodimer

Nonrephasing / Homodimer

Nonrephasing / Heterodimer

(a)

(c)

(e)

neg.

neg.

neg.

R12

R22

N12

N22

R11

R21

S0.1

N11

N21

pos.

pos.

pos.

(b)

(d)

(f)

neg.

neg.

neg.

pos.

pos.

pos.

Fully correlated fluctuations

FIG. 6. The rephasing beating map of a heterodimer and the non-rephasing beating map of homo- and heterodimers. In (a) and (b), we consider a heterodimer modeled by 1 = 12600 cm-1, 2 = 12400 cm-1, J = 100 cm-1,  = 50 cm-1,  = (50 fs)-1, s = | 2 - 1|  283 cm-1 and T = 77 K. In (a) and (b), the correlation length is taken to be  = 10-3d (local phonon baths) and  = 103d (a shared phonon bath), respectively, for which the rephasing beating maps at negative and positive frequencies 2 =  | 2 - 1| are displayed. In (c) and (d), where  = 10-3d and  = 103d, respectively, the nonrephasing beating maps of a homodimer are displayed with the model parameters used in Fig. 2. In (e) and (f), where  = 10-3d and  = 103d, respectively, the non-rephasing beating maps of a heterodimer are displayed with the model parameters used in
(a) and (b).

that our observations are not sensitive to model parameters, as shown in Appendix C and Fig. 7, where uncorrelated static disorder is taken into account in 2D simulations. For small static disorder with a full width at half maximum (FWHM) of 50 cm-1, asymmetric 2D

lineshapes, rephasing diagonal-peak and non-rephasing cross-peak oscillations are visible in simulations. For larger static disorder with a FWHM of 100 cm-1, 2D lineshapes start to be elongated along the diagonal due to inhomogeneous broadening, but rephasing diagonal-peak

10

and non-rephasing cross-peak oscillations are still visible. In Appendix D, we investigate this issue analytically to clarify how the non-secular terms and correlations in the noise affect oscillatory features in 2D spectra and how the non-secular effects can be described quantitatively with Feynman diagrams. For the homodimer, we show that non-secular oscillations are induced by non-secular interaction between different inter-excitonic coherences, as the dynamics of exciton populations are decoupled from those of coherences. For the heterodimer, where the populations and coherences are all coupled to one another, we show that non-secular oscillations are mainly induced by the mixing of different coherences, rather than population-coherence mixing, contrary to the suggestions in Refs. 52 and 53.
IV. DISCUSSION
In Ref. 17, three of the present authors demonstrated that the experimentally observed asymmetric lineshape in the rephasing beating map of J-aggregates cannot be explained by a correlated fluctuation model within the secular approximation. The asymmetric lineshape of the rephasing cross-peak of J-aggregates was found to originate from the fast population relaxation from higher to lower energy excitons17. The present results based on the Bloch-Redfield equation beyond the secular approximation further support the claim that when the oscillatory 2D signals have a long lifetime, and the lineshapes in the beating map are sufficiently asymmetric, long-lived beating signals are not dominated by correlated fluctuations. In Ref. 17, it was found that homogeneous broadening dominates the 2D lineshapes of J-aggregates and the exciton splitting is of the order of  700 cm-1. For such a large exciton splitting, non-secular effects are unlikely to induce notable signatures in oscillatory 2D signals, which are in line with the results shown in Ref. 58 based on quantum process tomography.
In Ref. 18, the experimentally measured 2D spectra of colloidal semiconductor nanoplatelets were reported, where heavy- and light-hole excitons exhibit lower and higher energy peaks in 2D spectra. It was found that the 2D lineshapes of the semiconductor system are dominated by homogeneous broadening, and the broadening of the higher energy exciton is approximately three times larger than that of the lower energy exciton, leading to asymmetric 2D lineshapes in the (1, t2, 3) domain: a beating map analysis in the (1, 2, 3) domain was not provided in Ref. 18. Given that highly asymmetric lineshapes were observed in experiments, our theoretical study predicts that highly correlated noise is unlikely to be present in the semiconductor system, and purely electronic coherences are unlikely to induce long-lived 2D oscillations. This is in line with the experimental observations where the lifetime of oscillatory 2D signals is similar to that of the optical coherences of the heavy- and lighthole excitons18. The authors of Ref. 18 concluded that

partially correlated noise is present in their system. We note that the exciton splitting of the semiconductor system was found to be in the range of 1200  1600 cm-118, depending on the sample preparation. For such a large exciton splitting, our results predict that non-secular effects are unlikely to induce notable features in oscillatory 2D signals. This is in line with the experimental 2D spectra where oscillatory features are present only at rephasing cross-peaks and non-rephasing diagonal-peaks18.
In Refs. 52 and 53, it was suggested that the nonsecular interaction between exciton populations and inter-excitonic coherences may induce oscillatory features in the rephasing diagonal-peaks of the photosystem II reaction center and the FMO complex. Our results support the claim that electronic coherences can induce such diagonal oscillations in the rephasing spectra, mediated by non-secular couplings, as the exciton splittings of the photosynthetic systems are relatively small, typically in the range of 100  200 cm-152,53. However, our quantitative analysis demonstrates that non-secular effects may be dominated by the interaction between inter-excitonic coherences, rather than the mixing of exciton populations and inter-excitonic coherences, depending on the model parameters. A detailed quantitative analysis with a beating map may be helpful for the identification of the microscopic origin of the oscillatory 2D signals of photosynthetic complexes, at least for simulated 2D spectra. We note that in Refs. 52 and 53, it was suggested that non-secular terms may be related to the functional relevance of the inter-excitonic coherences in exciton transport, as these non-secular terms couple the dynamics of the inter-excitonic coherences to that of exciton populations. Our results demonstrate that non-secular effects are suppressed as correlations in the noise increase, which suggests the possibility that there could be a trade-off between non-secular effects and the lifetime of purely electronic coherences. A further theoretical investigation based on non-Markovian quantum master equations59,60, and numerically exact methods such as TEDOPA23,26 and hierarchical equations of motion (HEOM)6163, could be helpful for the identification of the trade-off relation.
In our simulations, each time interval, i.e. t1, t2, t3, of the response function was described independently using the Bloch-Redfield equation, with the Born-Markov approximation where the phonon bath is in its equilibrium state in the electronic ground state manifold for all times. We note that the non-equilibrium dynamics of the phonon bath can induce the correlations between different time intervals and lead to much richer spectral lineshapes beyond the Lorentzians64,65. In the spatially correlated noise model, slow bath relaxation or a strong coupling of electronic states to a phonon environment can induce notable temporal correlations, which can be studied using numerically exact methods, such as HEOM54,61,62. Such a study of temporal and spatial correlations is beyond the scope of this work. Temporal correlations also play an important role in a vibronic model where underdamped vibrational modes modulate

11

2D lineshapes66. We also note that absorption lineshapes computed by using exact methods can show quantitative differences from those computed by using approximate methods55,67,68. Exact simulations of absorption and 2D spectra will be helpful for fully characterizing the signatures of correlated fluctuations. Such studies, which go well beyond the scope of the present work, will be pursued in a forthcoming work.
V. SUMMARY AND CONCLUSIONS
In this work, we investigated the influence of nonsecular couplings and spatial noise correlations on oscillatory 2D signals in rephasing and non-rephasing spectra. We employed the Bloch-Redfield formalism where we can tune the degree of correlations in the noise, such that we can cover the two extreme cases of uncorrelated and fully correlated noise. We performed a beating map analysis to identify the signatures of non-secular effects and noise correlations in oscillatory 2D spectra.
For uncorrelated noise, we found that non-secular couplings induce the mixing of exciton populations and inter-excitonic coherences, which lead to oscillations centered at rephasing diagonal-peaks and non-rephasing cross-peaks. With a developed quantitative method, we showed that the mixing of different inter-excitonic coherences is mainly responsible for the 2D oscillations induced by non-secular couplings. We also showed that the uncorrelated noise can induce asymmetric lineshapes of 2D peaks elongated along the excitation or detection axis.
For correlated noise, we showed that the non-secular effects are suppressed by correlations in the noise. This spatially correlated noise can induce long-lasting 2D oscillations centered at rephasing cross- and non-rephasing diagonal-peaks, but with suppressed oscillatory features in rephasing diagonal- and non-rephasing cross-peaks. We also showed that correlations in the noise enforce symmetry onto 2D lineshapes, hinting that the degree of asymmetry in 2D lineshapes could be used to estimate to what degree the noise is spatially correlated. Our results demonstrate that a detailed analysis of the oscillatory features in 2D electronic spectra may provide information on the structure of vibrational environments, such as correlations in the noise and the strength of non-secular environmental couplings.
ACKNOWLEDGEMENTS
This work was supported by the EU STREP project PAPETS and QUCHIP, the ERC Synergy grant BioQ, the Deutsche Forschungsgemeinschaft (DFG) within the SFB/TRR21 and an Alexander von Humboldt Professorship, and the state of Baden-Wurttemberg through bwHPC. This research was undertaken with the assistance of resources from the National Computational Infrastructure (NCI), which is supported by the Australian

Government.

Appendix A: Bloch-Redfield equation

The Bloch-Redfield equation is expressed as

d dt

=

- i [He, ] +

2

-sj V qjkV 

j,k=1

(A1)

+V qjkV sj - V q^jkV sj + sj V q^jkV  ,

where  denotes the density matrix of the dimer sys-

tem, V =

4 n=1

|n

an| is a unitary operator with

{|n } representing the electronic eigenstates of the sys-

tem Hamiltonian He, defined by He |n = n |n , in

an arbitrary basis {|an }, and the other terms are given

by

s1 = 1+1-, s2 = 2+2-, an| qjk |am = an| V skV |am
an| q^jk |am = an| V skV |am

1 2

Cj k (m

-

n),

1 2

Ckj (n

-

m),

(A2) (A3) (A4)
(A5)

where the spectral functions Cjk() are defined by

Cjk() =

1
2


d ei
-

eiHp / Bj e-iHp / Bk , (A6)

B1 =
k
B2 =
k

gk((b1k + b1k) + (b2k + b2k)), gk((b2k + b2k) + (b1k + b1k)).

(A7) (A8)

More specifically, we consider the site basis, given by
|a1 = |g1, g2 , |a2 = |e1, g2 , |a3 = |g1, e2 , |a4 = |e1, e2 . The spectral functions Cjk() are reduced to

C11() = C22() = C(), C12() = C21() = 2C(),
C() = 2gk2 [(n(k) + 1)( - k)
k
+n(k)( + k)] ,

(A9) (A10)
(A11)

where n(k) = (exp( k/kBT ) - 1)-1 is the mean phonon number of a phonon mode with a frequency of k at temperature T , while (x) denotes the Dirac delta function. Based on the fact that 0  2 =
2 1 - 2  1, we introduce a correlation length  to quantify the degree of spatial correlations in the noise, defined by exp(-d/) = 2, where d denotes the spatial distance between sites 1 and 2. When  d, C11() = C22() = C() and C12() = C21()  0, leading to local (or spatially uncorrelated) noise, while

12

when  d, Cjk()  C() for all j and k, leading to fully correlated noise. An intermediate case of   d
leads to partially correlated noise. The correlated noise
is known to enhance the lifetime of electronic coherences in the single excitation subspace36,3944. This is contrary
to the anti-correlated noise defined by  = - 1 - 2
in Eq. (11), which is known to suppress the lifetime of excited state coherences40,43. In this work, we do not
consider the anti-correlated noise, as we are interested in
the scenario that correlations in the noise enhance the
lifetime of excited state coherences, leading to long-lived
oscillatory 2D signals. Therefore, the spectral functions
Cjk() in the presence of spatially correlated noise can be summarized as

C11() = C22() = C(), C12() = C21() = e-d/C().

(A12) (A13)

In the continuous limit of phonon modes, leading to a phonon bath, the spectral function C() is reduced to

2J ()(n() + 1)  C() = 2J (||)n(||)
lim0 2J ()n()

 > 0,  < 0,  = 0,

(A14)

where J () is the phonon spectral density that describes
the phonon mode density weighted by the electronphonon coupling strength gk and satisfies J (0) = 0. In this work, J () is modeled by a shifted Ohmic spec-
tral density described in Eq. (12). The shift s of the Ohmic spectral density can make C(0) C (| 2 - 1|) (cf. Eq. (A14)), for instance, when s  | 2 - 1|, such that the pure dephasing rates proportional to C(0) are
smaller than the relaxation rate between single exciton states | 1 and | 2 . We note that the pure dephasing and relaxation rates are not only determined by the spec-
tral function C(), but also by the system parameters
of the electronic Hamiltonian, described by sj and V in Eqs. (A1), (A4) and (A5).

Appendix B: 2D electronic spectroscopy
In 2D experiments, three excitation pulses interact with a molecular system and the resultant third-order optical response of the system is measured as a function of the time delays between the first and second, the second and third, and the third excitation pulse and the emitted signal from the molecular system. These time delays are called coherence time t1, waiting time t2 and rephasing time t3, respectively. The Fourier transformation of the response function with respect to t1 and t3 leads to 2D spectra as a function of excitation frequency 1 and detection frequency 3. When the pulse duration of the excitation pulses is short enough, the excitation fields can be approximately described by the Dirac delta function in the time domain, for which the thirdorder optical response function can be described within

the rotating wave approximation in the impulsive limit. This is equivalent to the assumption that the laser spectrum is broad enough to cover the electronic states of the dimer system in the frequency domain. When the laser spectrum is not broad enough for a given system, one needs to take into account the pulse duration explicitly in 2D simulations69. For the waiting times longer than the pulse duration, it was found that the finite pulse duration mainly acts as a frequency filter70.
Within the rotating wave approximation in the impulsive limit, rephasing 2D spectra are formally expressed as





SR(1, t2, 3) =

dt1

dt3 e-i(1 t1 -3 t3 )

0

0

 [RGSB + RSE - RESA],

(B1)

where RGSB, RSE and RESA denote the ground state bleaching (GSB), stimulated emission (SE) and excited state absorption (ESA) contributions to the rephasing spectra, respectively:
RGSB(t1, t2, t3) = tr[-u(t3)[+u(t2)[u(t1)[eq-]+]]], (B2)
RSE(t1, t2, t3) = tr[-u(t3)[u(t2)[+u(t1)[eq-]]+]], (B3)
RESA(t1, t2, t3) = tr[-u(t3)[+u(t2)[+u(t1)[eq-]]]], (B4)

with eq representing the equilibrium state in the electronic ground state manifold, u(t) is a formal representation of the propagator, determined by the Bloch-Redfield equation in this work. Here  denote the transition dipole operators of the molecular system, describing the optical transition between ground and excited states by the excitation pulses

+ = (e^  d1)1+ + (e^  d2)2+, - = (e^  d1)1- + (e^  d2)2-,

(B5) (B6)

where e^ denotes the polarization direction of the excitation pulses, which are all assumed to be parallel in this
work, while dk represents the transition dipole moment of site k. In 2D simulations, we take into account the
rotational averaging of the dipole moments dk with respect to the polarization direction e^, as we are considering 2D measurements of an ensemble of dimers. We assume that the sites 1 and 2 have mutually orthogonal transi-
tion dipoles with the same magnitude, i.e. d1  d2 = 0 and d1  d1 = d2  d2. In the GSB pathway, the system is in the ground state during waiting time t2, while in the SE and ESA pathways, the system is in the single excitation subspace during t2. Within our model, the oscillatory 2D signals originate only from the SE and ESA contributions, as we are not considering ground state vibrational coherences induced by underdamped vibrational motions.

13

Similarly, non-rephasing 2D spectra can be formally expressed as





SN (1, t2, 3) =

dt1

dt3 ei(1 t1 +3 t3 )

0

0

 [NGSB + NSE - NESA],

(B7)

where the GSB, SE and ESA contributions are expressed as

NGSB(t1, t2, t3) = tr[-u(t3)[+u(t2)[-u(t1)[+eq]]]], (B8)
NSE(t1, t2, t3) = tr[-u(t3)[u(t2)[u(t1)[+eq]-]+]], (B9)
NESA(t1, t2, t3) = tr[-u(t3)[+u(t2)[u(t1)[+eq]-]]]. (B10)

Appendix C: Inhomogeneous broadening
Here we demonstrate how oscillatory features in the beating map are affected by inhomogeneous broadening. In Fig. 7, we consider uncorrelated disorder, where the site energies 1 and 2 of a dimer are described by two independent Gaussian distributions centered at the average values of 1 and 2 , respectively. Here we consider the same full width at half maximum (FWHM) for both Gaussian distributions, and employ the model parameters of the heterodimer used in Fig. 6, where
1 = 12600 cm-1 and 2 = 12400 cm-1. Figs. 7(a) and (b) show the rephasing and non-rephasing beating maps, respectively, for the case of uncorrelated fluctuations (i.e.  = 10-3d) with the inhomogeneous broadening modeled by a FWHM of 50 cm-1. Compared to Figs. 6(a) and (e), where the inhomogeneous broadening is not considered, the overall 2D lineshapes in Figs. 7(a) and (b) become broader due to the inhomogeneous broadening. However, the oscillatory features in the rephasing diagonal-peaks and non-rephasing cross-peaks (i.e. nonsecular effects) and the asymmetric 2D lineshapes elongated along 1-axis are still visible. As the FWHM increases further, the 2D lineshapes are elongated along the diagonal (1 = 3), but non-secular effects are still visible for a FWHM of 100 cm-1, as shown in Figs. 7(c) and (d). This is somewhat relevant, as the static disorder of the FMO complex has been modeled by a FWHM of  100 cm-1 in other works54,71. On the other hand, for the case of correlated fluctuations, the rephasing beating map shows strong elongation of a cross-peak along diagonal, as shown in Fig. 7(e), while the non-rephasing beating map shows a relatively symmetric 2D lineshape of a diagonal peak, as shown in Fig. 7(f). This is due to the difference in the phase distributions of the rephasing and non-rephasing spectra in the (1, 3) domain46,72.

Appendix D: Diagonalization of the Liouville space operator

Within the Bloch-Redfield formalism, the oscillatory signals in rephasing spectra are induced by excited state coherences described by the stimulated emission (SE) and excited state absorption (ESA) contributions in the theory of 2D spectroscopy6 (cf. Eqs. (B3) and (B4)). Here we show how the SE contribution to the oscillatory rephasing signals can be described quantitatively to identify the role of non-secular couplings and spatial noise correlations in the beating map. The analytical approach, which will be presented below, can be generalized to the ESA contribution to the rephasing spectra as well as to the SE and ESA contributions to the non-rephasing spectra. The analysis is based on the diagonalization of the Liouville space operator. This approach can be generalized to the other quantum master equations beyond the Bloch-Redfield equation employed in this work.
The lineshape of 2D spectra along excitation axis 1 is determined by the dynamics of optical coherences between ground state and singly excited states during the coherence time t1. In the exciton basis, the optical coherences are expressed as g1(t1) |g 1| + g2(t1) |g 2| where the time evolution of g1(t1) and g2(t1) is governed by

d dt1

g1(t1) g2(t1)

=

X11 X12 X21 X22

g1(t1) g2(t1)

,

(D1)

where the super-operator X describes both the Hamiltonian dynamics and decoherence. For the Bloch-Redfield equation summarized in Appendix A, the elements of X are given by

X11

=

-

1 4

C (0)(2

-

(1

-

e-d/ )s2 (2))

-

1 4

C (-

)(1

-

e-d/ )s2 (2)

+

i

1,

X22

=

-

1 4

C (0)(2

-

(1

-

e-d/ )s2 (2))

-

1 4

C (

)(1

-

e-d/ )s2 (2)

+

i

2,

X12

=

1 8

(C (0)

-

C (

))(1

- e-d/)s(4),

X21

=

-

1 8

(C (0)

-

C (-

))(1

-

e-d/ )s(4),

(D2)
(D3)
(D4) (D5)

with s()  sin(),  quantifies the delocalization of excitons in the site basis, and  = | 2 - 1| denotes the exciton splitting between | 1 and | 2 (see Eqs. (3) and (4)). Here C() represents the spectral function determined
by the phonon spectral density, as shown in Eq. (A14).
The lineshape of 2D spectra along the excitation axis
1 can be represented analytically by using the eigenstates of the super-operator X, defined by Xxk = kxk. In the exciton basis, the eigenvalue equation is given by

X11 X12 X21 X22

x(gk1) x(gk2)

= k

x(gk1) x(gk2)

,

(D6)

Rephasing (pos.)

14

Uncorrelated noise / Small disorder Uncorrelated noise / Large disorder

(a)

(c)

(e)

pos.

pos.

Correlated noise / Large disorder
pos.

(b)

neg.

(d)

neg.

(f)

neg.

Nonrephasing (neg.)

FIG. 7. The rephasing and non-rephasing beating maps of a heterodimer in the presence of inhomogeneous broadening. Here
we employed the model parameters used in Figs. 6(a), (b), (e) and (f): 1 = 12600 cm-1, 2 = 12400 cm-1 (the average site energies), J = 100 cm-1,  = 50 cm-1,  = (50 fs)-1, s  283 cm-1 (the exciton splitting for the average site energies) and T = 77 K. In (a) and (b), the rephasing beating map at a positive frequency of 2 = 283 cm-1 and the non-rephasing beating map at a negative frequency of 2 = -283 cm-1 are displayed, respectively, for the case that  = 10-3d (uncorrelated fluctuations) and the inhomogeneous broadening is modeled by Gaussian distributions with a FWHM of 50 cm-1. In (c) and
(d), the rephasing and non-rephasing beating maps are displayed, respectively, for the case that  = 10-3d (uncorrelated
fluctuations) and the inhomogeneous broadening is modeled by a larger FWHM of 100 cm-1. In (e) and (f), the rephasing and
non-rephasing beating maps are displayed, respectively, for the case that  = 103d (correlated fluctuations) and the FWHM is
taken to be 100 cm-1.

where the eigenvector xk in the Liouville space corresponds to an optical coherence x^k in the Hilbert space

x^k = x(gk1) |g 1| + x(gk2) |g 2| ,

(D7)

satisfying

d dt

x^k

=

k x^k

with

an

associated

eigenvalue

of

k. This implies that non-secular couplings X12 and X21

induce a mixing of two optical coherences |g 1| and

|g 2| in the exciton basis. The dynamics of the mixed

coherence is formally described by u(t1)[x^k] = ekt1 x^k.

The optical coherence created by the first excitation

pulse, i.e. |g g| -, can be represented as a superpo-

sition of x^k

2
|g g| - = |g
j=1

2
j |gj = kx^k,
k=1

(D8)

where gj represents the transition dipole strength between ground state |g and the j-th exciton | j for a given realization of the transition dipole moments of sites
1 and 2 (cf. Eqs. (B5) and (B6)). The coefficient k in Eq. (D8) describes the effective transition dipole strength
between ground state |g g| and mixed coherence x^k,

given by

1 2

=

x(g11) x(g21) -1 x(g12) x(g22)

g1 g2

.

(D9)

The dynamics of |g g| - during time t1 is then expressed as

2
u(t1)[|g g| -] = kekt1 x^k,
k=1

(D10)

and the Fourier transformation of ekt1 determines the lineshape of the homogeneously broadened 2D spectra along the excitation axis 1


dt1e-i1t1 u(t1)[|g
0

g|

-]

=

2 k=1

k

k - i1

x^k .

(D11)

Here the real and imaginary parts of the eigenvalue k

of x^k, denoted by Re[k] and Im[k], respectively, de-

termine the homogeneous broadening and peak location,

respectively, of the k-th Lorentzian peak along the exci-

tation axis. The non-secular couplings X12 and X21 in

15

Eq. (D1) can make the imaginary part of k deviate from the eigenvalue k of the system Hamiltonian (cf. Eqs. (7) and (8)), implying that 2D peak locations can be shifted by non-secular effects.
These results imply that the dynamics of the eigenstates x^1 and x^2 of the super-operator X lead to the lower and higher energy peaks, respectively, along the excitation axis (cf. Figs. 2-6). When the off-diagonal components X12 and X21 are comparable or larger in magnitude than the difference in the diagonal components X11 and X22, e.g. |X12| |X11 - X22|, the eigenstates x^k become a superposition of the optical coherences |g 1| and |g 2| in the exciton basis. We note that the difference in X11 and X22 is larger in magnitude than the exciton splitting, i.e. |X11 - X22|  |Im[X11 - X22]| = | 2 - 1|, indicating that, as expected, the non-secular effects are suppressed as the exciton splitting increases. The mixing is also suppressed by correlations in the noise, i.e. X11, X22 = 0 and X12, X21  0 as 1 - e-d/  0 in

Eqs. (D2)-(D5). So far we have analyzed the dynamics of the optical
coherences created by the first excitation pulse. We now consider the populations and coherences within the single excitation subspace created by the second excitation pulse. The population or coherence in the excited state manifold is expressed in the exciton basis as

2
ij (t2) | i j |,
i,j=1

(D12)

whose dynamics are governed by the super-operator Y ,

11(t2) Y11,11 Y11,22 Y11,12 Y11,21 11(t2)

d dt2

22(t2

) 

12(t2)

=

Y22,11 Y12,11

Y22,22 Y12,22

Y22,12 Y12,12

Y22,21 Y12,21

 

22(t2 12(t2

)  )

.

21(t2)

Y21,11 Y21,22 Y21,12 Y21,21

21(t2)

(D13)

For the Bloch-Redfield equation, the elements Yjk,lm of Y are given by

Y11,11 Y11,22 Y11,12 Y11,21 0 0

0

Y22,11 Y12,11

Y22,22 Y12,22

Y22,12 Y12,12

Y22,21 Y12,21

=

0  0

0 0

0 -i( 1 -

2)

0

 -2C(- )s2(2) 2C( )s2(2)

0 0

 

+

1

-

e-d/ 4

 

2C(- )s2(2) C(- )s(4)

-2C( )s2(2) -C( )s(4)

Y21,11 Y21,22 Y21,12 Y21,21

00

0

-i( 2 - 1)

C(- )s(4) -C( )s(4)

C (0)s(4)

C (0)s(4)



-C (0)s(4) -(C( ) + C(- ))s2(2) - 4C(0)c2(2)

-C (0)s(4) (C( ) + C(- ))s2(2)

, 

(C( ) + C(- ))s2(2)

-(C( ) + C(- ))s2(2) - 4C(0)c2(2)

(D14)

with c()  cos(). The first term on the right hand side in Eq. (D14) shows the Hamiltonian contribution to the system dynamics, which is proportional to the exciton splitting of 2 - 1. This implies that as the exciton splitting increases in magnitude, the difference in diagonal components of Y increases, and as a result the nonsecular interactions between exciton populations ii(t) and inter-excitonic coherences ij(t) with i = j, and those between different inter-excitonic coherences 12(t) and 21(t) are suppressed. The second term on the right hand side in Eq. (D14) describes decoherence within the single excitation subspace. The factor (1 - e-d/) in Eq. (D14) shows that in the long correlation length limit, i.e.  d, there is no decoherence in the single excitation subspace, and the dynamics of the singly excited states are governed by the Hamiltonian only. The dynamics of these singly excited states can be described by the eigen-

states y^l of the super-operator Y , satisfying

Y11,11 Y22,11 Y12,11
Y21,11

Y11,22 Y22,22 Y12,22 Y21,22

Y11,12 Y22,12 Y12,12 Y21,12

Y11,21 Y22,21 Y12,21 Y21,21

y1(l1) y2(l2) y1(l2)
y2(l1)

=

l

y1(l1) y2(l2) y1(l2)
y2(l1)

,

(D15)

which is expressed in the exciton basis as

2

y^l =

yi(jl) | i j |.

i,j=1

(D16)

The time evolution of the eigenstate y^l is formally expressed as u(t2)[y^l] = elt2 y^l where the real and imaginary parts of the eigenvalue l describe the decay and phase evolution, respectively, of the eigenstate y^l. The phase evolution leads to oscillatory 2D signals. For
the dimer system considered in simulations, we found
that two of the eigenvalues l have imaginary parts, which are approximately given by Im[1]  | 2 - 1| and Im[2]  - | 1 - 2|. The associated two eigenstates y^1 and y^2 are responsible for the oscillatory 2D signals with

16

positive and negative frequencies, respectively. The other where the two-dimensional amplitude Aklm(1, 3) de-

eigenstates y^3 and y^4 have negligible imaginary parts, implying that they are responsible for non-oscillatory 2D
signals, such as exponential and static t2-transients. The

scribes a Lorentzian peak centered at (1, 3) = (Im[k], -Im[m]), weighted by the effective transition dipole strength

time evolution of the SE contribution during the waiting time t2 can be expressed as
2

Aklm(1, 3)

=

k kl lm tr[- x^m ] (k - i1)(m + i3

)

,

(D23)

u(t2)[+u(t1)[|g

g| -]] = kekt1 u(t2)[+x^k]
k=1
(D17)

2

4

=

k ek t1

klelt2 y^l,

where the homogeneous broadenings along 1- and 3axes are determined by the real part of the eigenvalues k and m, respectively. Here kkllmtr[-x^m] denotes the rotational average (ensemble) of the effective tran-
sition dipole strength (cf. Appendix B). These results

k=1

l=1

(D18)

where +x^k =

4 l=1

kly^l

with

kl

representing

the

effec-

tive transition dipole strength between eigenstates x^k and

y^l. In 2D simulations, one can calculate the beating map

directly by removing the non-oscillatory components y^3

and y^4 from Eq. (D18), then Fourier transforming elt2

show that y^l=1,2 can induce oscillatory rephasing signals centered at (1, 3) = (Im[k], -Im[m]) when the associated transition dipole strength kkllm tr [-x^m] is not zero.
When the optical coherences x^k=1,2 are approximately represented by x^k=1,2  |g k|, the eigenstates y^l=1,2 can induce 2D oscillations centered at rephasing lower diagonal-peak R11 when the Feynman path-

which leads to the l-th Lorentzian peak along the 2-axis (cf. l = 1, 2).

way of |g g|  x^1  y^l  x^1  |g g| has a nonzero transition dipole strength. Here the transition from

Finally we consider the dynamics of the optical coherences 1g(t3) | 1 g| + 2g(t3) | 2 g| created by the third excitation pulse, whose dynamics during rephasing time
t3 are described by

d dt3

1g (t3 ) 2g (t3 )

=

X11 X12 X21 X22

1g (t3 ) 2g (t3 )

.

(D19)

x^1  |g 1| to y^l is allowed when y^l | 1 = 0, and the transition from y^l to x^1  | 1 g| is allowed when
1| y^l = 0. These conditions are not satisfied within the secular approximation36 where y^1 = | 1 2| and y^2 = | 2 1| and the super-operator Y is approximated by

The eigenstates of the super-operator X are given by

x^k =

x(gk1)


|1

g| +

x(gk2)


|2

g| ,

(D20)

11(t2) Y11,11 Y11,22 0

0  11(t2)

d dt2

22(t2

) 

12(t2)

=

Y22,11 0

Y22,22 0

0 Y12,12

0 0

 

22(t2 12(t2

)  )

.

21(t2)

0

0

0 Y21,21 21(t2)

(D24)

with the associated eigenvalues k. Here k (or x^k) is the complex conjugate (or adjoint) of k (or x^k). The

On the other hand, in the presence of non-secular couplings, the eigenstates y^l=1,2 become a mixture of differ-

SE contribution to the rephasing spectra (cf. Eq. (B3)) ent inter-excitonic coherences | 1 2| and | 2 1| and

is then expressed as

exciton populations | 1 1| and | 2 2|, which are for-

RSE(t1,

t2,

t3)

=

2 k,m=1

4 l=1

kekt1 klelt2 lmemt3 tr[-x^m], (D21)

mally expressed as y^l = y1(l2) | 1 2| + y2(l1) | 2 1| + y1(l1) | 1 1| + y2(l2) | 2 2|. As shown in Appendix E, the conditions of y^l | 1 = 0 and 1| y^l = 0 can be satisfied for the homo- and heterodimers considered in our
simulations. For the homodimer, non-secular interac-

where y^l+ =

2 m=1

lmx^m

and

lm

denotes

the

effec-

tive transition dipole strength between eigenstates y^l and

x^m. The summations over k, l, m, where k  {1, 2},

l  {1, 2, 3, 4} and m  {1, 2}, lead to 16 different Feyn-

man pathways for the SE contribution to the rephasing

spectra. Since only y^1 and y^2 are responsible for oscilla-

tory 2D signals, there are only eight Feynman pathways

with l  {1, 2} contributing to the beating map. Thus,

the SE contribution to the oscillatory rephasing signals

in the (1, 3) domain can be expressed as

tion between populations and coherences is absent, as shown in Eq. (E1), where the eigenstates y^l=1,2 are the mixtures of different inter-excitonic coherences only, i.e. y^l=1,2 = y1(l2) | 1 2| + y2(l1) | 2 1|, as shown in Eqs. (E2) and (E3). In this case, the dynamics of exciton populations are decoupled from those of inter-excitonic coherences, but the non-secular interaction between coherences can induce rephasing diagonal-peak oscillations. For the heterodimer, all the exciton populations and inter-excitonic coherences are coupled to one another and induce a population-coherence mixing in the eigenstates

2

RSE(1, t2, 3) =

Aklm(1, 3)elt2 ,

k,l,m=1

(D22)

y^l=1,2, which also lead to oscillations centered at rephasing diagonal peaks. We found that for the model param-
eters used in our simulations, the mixing is dominated by

17

inter-excitonic coherences, and the contribution of exci- ences:

ton populations to y^l=1,2 is relatively small, as shown in Eqs. (E9) and (E10). A detailed quantitative description

y^1  0.991 | 1 2| - 0.133 i | 2 1| ,

(E2)

of non-secular effects with associated Feynman diagrams is provided in Appendix E.
Finally, we note that the asymmetric lineshape in the beating map originates from the fact that the lower and

y^2  0.133 i | 1 2| + 0.991 | 2 1| , y^3  0.717 | 1 1| - 0.717 | 2 2| , y^4  0.999 | 1 1| + 0.024 | 2 2| ,

(E3) (E4) (E5)

higher energy peaks have different homogeneous broadenings. When the spectral function satisfies C ( ) > C (0), as is the case of the model parameters used in this work, relaxation dominates the homogeneous broadening, and the higher energy peak shows a larger broadening than the lower energy peak, described by |Re[1]| < |Re[2]|. This leads to asymmetric lineshapes in the beating map, as shown in Figs. 4(a) and (b). As correlations in the noise increase, the super-operator X is governed by the pure dephasing noise described by C (0), as shown in Eqs. (D2)-(D5), where the relaxation described by C ( ) does not contribute to the homogeneous broadening, leading to |Re[1]|  |Re[2]| and symmetric lineshapes in the beating map, as shown in Figs. 4(g) and (h).

with the associated eigenvalues given by 1 = (-53 + 193 i) cm-1, 2 = (-53-193 i) cm-1, 3 = -105 cm-1 and 4 = 0. The first eigenstate y^1 is a superposition of the inter-excitonic coherences | 1 2| and | 2 1| due to the non-secular couplings Y12,21 and Y21,12. The imaginary part of the associated eigenvalue Im[1]  193 cm-1 shows that y^1 leads to a positive frequency component in the beating map with a beating frequency of 2  193 cm-1. Due to the non-secular effects, the beating
frequency is slightly different from the exciton splitting of 2 - 1 = 200 cm-1. Similarly, the second eigenstate y^2 is a superposition of the inter-excitonic coherences, but it has a larger amplitude in | 2 1| than in | 1 2|, contrary to y^1. This results in the imaginary part of the associated eigenvalue having the opposite

sign, Im[2]  -193 cm-1, implying that y^2 leads to a

negative frequency component in the beating map with

Appendix E: Non-secular effects

2  -193 cm-1. The eigenvalues of the other eigenstates y^3 and y^4 do not contain imaginary parts, implying

Here we apply the quantitative method developed in Appendix D to the model parameters of homo- and heterodimers considered in our simulations. We show that the mixing of inter-excitonic coherences is mainly responsible for the oscillations centered at rephasing diagonal peaks.

that they are associated with non-oscillatory 2D signals: 3 < 0 and 4 = 0 indicate that y^3 describes the relaxation of exciton populations, while y^4 is an equilibrium state within the excited state manifold.
Fig. 8 shows the Feynman diagrams of the SE contribution to the oscillatory rephasing signals for the homodimer. Figs. 8(a) and (b) show the Feynman diagrams

For the model parameters of the homodimer with responsible for the oscillations in the rephasing cross = 10-3d (i.e. local phonon baths), the off-diagonal peaks R21 and R12, respectively. Here the eigenstates

terms X12 and X21 in Eqs. (D4) and (D5) are zero, as sin(4) = 0 with  = /4 (cf. Eqs. (3) and (4)). This im-

y^1  | 1 2| - i | 2 1| and y^2  | 2 1| + i | 1 2| are approximately represented in terms of a small am-

plies that the eigenstates of the super-operator X (or X) are given by x^1 = |g 1| and x^2 = |g 2| (or x^1 = | 1 g| and x^2 = | 2 g|). In this case, the superoperator Y is reduced to

plitude 0 <  < 1. Note that the optical transitions x^2  y^1  x^1 in Fig. 8(a) and x^1  y^2  x^2 in Fig. 8(b) are allowed even in the absence of the small amplitude . On the other hand, Fig. 8(c) shows the Feynman diagram

responsible for the positive frequency component in the

11(t2) Y11,11 Y11,22 0

0  11(t2) rephasing diagonal-peak R11, where the optical transi-

d dt2

22(t2) 12(t2)

=

Y22,11 0

Y22,22 0

0 Y12,12

0 
Y12,21

22(t2) 12(t2)

,ntioonn-zx^e1ro, i.y^e1.

is x^1

allowed only = |g 1| 

if the small amplitude  is -i | 2 1| + | 1 2|  y^1,

21(t2)

0

0 Y21,12 Y21,21 21(t2) as a direct transition from |g 1| to | 1 2| is forbidden.

(E1) This implies that the mixing  of different inter-excitonic

as sin(4) = 0 (cf. Eq. (D14)), which shows that the dy-
namics of exciton populations ii(t2) are decoupled from those of inter-excitonic coherences ij(t2) with i = j.

coherences induced by the non-secular couplings Y12,21 and Y21,12 allows the transition x^1  y^1  x^1 in Fig. 8(c) to occur, leading to oscillations centered at the rephas-

However, the non-zero off-diagonal components Y12,21 ing diagonal-peak R11. Similarly, Fig. 8(d) shows the

and Y21,12 induce non-secular interactions between differ- Feynman diagram that induces the negative frequency

ent inter-excitonic coherences 12(t2) and 21(t2). Even though the off-diagonal components Y12,21 = Y21,12  53 cm-1 are an order of magnitude smaller than the

component in the rephasing diagonal-peak R11. We now consider the model parameters of the het-
erodimer with  = 10-3d (i.e. local phonon baths). In

difference in diagonal components, |Y12,12 - Y21,21|  400 cm-1, the eigenstates y^l=1,2 of the super-operator Y

this case, the off-diagonal components X12 and X21 are non-zero, which makes the eigenstates x^1 and x^2 of the

show a notable mixing of different inter-excitonic coher- super-operator X be in a superposition of the optical co-

18

(a) R21 (pos.)

(b) R12 (neg.)

eigenstates of the super-operator Y are given by

|g g| x^1 = |1 g | y^1  |1 2| - i |2 1| x^2 = | g 2|
|g g|
(c) R11 (pos.)
|g g| x^1 = |1 g | y^1  |1 2| - i |2 1|
x^1 = | g 1| |g g|

|g g| x^2 = |2 g | y^2  |2 1| + i |1 2| x^1 = | g 1|
|g g|
(d) R11 (neg.)
|g g| x^1 = |1 g | y^2  |2 1| + i |1 2| x^1 = | g 1| |g g|

FIG. 8. The Feynman diagrams of the SE contribution to the oscillatory signals in the rephasing spectra of a homodimer. In (a) and (b), the Feynman diagrams responsible for the oscillatory signals in the rephasing cross-peaks R21 and R12 are displayed, respectively. In (c) and (d), the Feynman diagrams responsible for the oscillatory signals in the rephasing diagonal peak R11 are shown, which lead to positive and negative frequency components, respectively. Here a small amplitude , satisfying 0 <  < 1, is employed to approximately represent the eigenstates y^1 and y^2 (see Eqs. (E2) and (E3)).

(a) R22 (pos.)

(b) R22 (neg.)

y^1  0.997 | 1 2| + 0.061e-1.59 i | 2 1|

(E9)

+ 0.011e-1.53 i | 1 1| + 0.011e1.61 i | 2 2| ,

y^2  0.061e1.59 i | 1 2| + 0.997 | 2 1|

(E10)

+ 0.011e1.53 i | 1 1| + 0.011e-1.61 i | 2 2| ,

y^3  0.169e1.55 i | 1 2| + 0.169e-1.55 i | 2 1| (E11)

+ 0.686 | 1 1| - 0.686 | 2 2| ,

y^4  0.999 | 1 1| + 0.005 | 2 2| ,

(E12)

with the associated eigenvalues given by 1  (-41 + 280 i) cm-1, 2  (-41 - 280 i) cm-1, 3  -70 cm-1 and 4 = 0. As is the case of the homodimer, the first two eigenstates y^1 and y^2 are responsible for oscillatory 2D signals, while the other eigenstates y^3 and y^4 induce exponential and static t2-transients, respectively. Note that the eigenstates are mixtures of exciton populations
and inter-excitonic coherences due to non-secular effects.
Fig. 9 shows the Feynman diagrams of the SE contri-
bution to the oscillations in the rephasing diagonal-peak
R22 of the heterodimer. Contrary to the case of the ho-
modimer, the mixing of the optical coherences |g 1| and |g 2| during coherence and rephasing times, described by a small amplitude |1| < 1, enhances the diagonal oscillations in the rephasing spectra.

|g g| x^2  1 |1 g | + |2 g | y^1  |1 2| + 2 |2 1|
x^2  | g 2| + 1 | g 1|
|g g|

|g g| x^2  |2 g | + 1 |1 g | y^2  |2 1| + 2 |1 2| x^2  1 | g 1| + | g 2|
|g g|

1H. van Amerongen, L. Valkunas and R. van Grondelle, Photo-
synthetic Excitons (World Scientific, 2000). 2R. E. Blankenship, Molecular Mechanisms of Photosynthesis
(Blackwell Science, 2002). 3A. Ishizaki, T. R. Calhoun, G. S. Schlau-Cohen, and G. R. Flem-
ing, Phys. Chem. Chem. Phys. 12, 7319 (2010).

FIG. 9. The Feynman diagrams of the SE contribution to the oscillatory signals in the rephasing spectra of a heterodimer. In (a) and (b), the Feynman diagrams responsible for the oscillations in the rephasing diagonal-peak R22 are displayed, which induce positive and negative frequency components, respectively. Here |1| < 1 and |2| < 1 are small amplitudes induced by non-secular effects.

4S. F. Huelga and M. B. Plenio, Contemp. Phys. 54, 181 (2013). 5A. Chenu and G. D. Scholes, Annu. Rev. Phys. Chem. 66, 69
(2015). 6D. M. Jonas, Annu. Rev. Phys. Chem. 54, 425 (2003). 7G. S. Engel, T. R. Calhoun, E. L. Read, T.-K. Ahn, T. Mancal,
Y.-C. Cheng, R. E. Blankenship, and G. R. Fleming, Nature 446,
782 (2007). 8E. Collini, C. Y. Wong, K. E. Wilk, P. M. G. Curmi, P. Brumer,
and G. D. Scholes, Nature 463, 644 (2010).

9J. R. Caram, N. H. C. Lewis, A. F. Fidler, and G. S. Engel, J.

Chem. Phys. 136, 104505 (2012).

herences |g 1| and |g 2| in the exciton basis, given by

10A. F. Fidler, E. Harel, P. D. Long, and G. S. Engel, J. Phys.
Chem. A 116, 282 (2012). 11E. Romero, R. Augulis, V. I. Novoderezhkin, M. Ferretti, J.

x^1  0.999 |g x^2  0.993 |g

1| + 0.005e1.69 i |g 2| , 2| + 0.115e-1.45 i |g 1| .

(E6) (E7)

Thieme, D. Zigmantas, and R. van Grondelle, Nat. Phys. 10, 676 (2014). 12F. D. Fuller, J. Pan, A. Gelzinis, V. Butkus, S. S. Senlik, D. E. Wilcox, C. F. Yocum, L. Valkunas, D. Abramavicius, and J. P.

The super-operator Y is reduced to

Ogilvie, Nat. Chem. 6, 706 (2014). 13F. Milota, V. I. Prokhorenko, T. Mancal, H. von Berlepsch, O.
Bixner, H. F. Kauffmann, and J. Hauer, J. Phys. Chem. A 117,

11(t2) Y11,11

d dt2

22(t2) 12(t2)

=

Y22,11 Y12,11

Y11,22 Y22,22 Y12,22

Y11,12 Y22,12 Y12,12

Y11,21 Y22,21 Y12,21


 

11 22 12

(t2 (t2 (t2

) )  )

6007 (2013). 14D. Hayes, G. B. Griffin, G. S. Engel, Science 340, 1431 (2013).
,15A. Halpin, P. J. M. Johnson, R. Tempelaar, R. S. Murphy, J.
Knoester, T. L. C. Jansen, and R. J. D. Miller, Nat. Chem. 6,

21(t2)

Y21,11 Y21,22 Y21,12 Y21,21

21(t2)

196 (2014).

(E8) 16Y. Song, S. N. Clafton, R. D. Pensack, T. W. Kee, and G. D.

which includes additional non-zero off-diagonal elements, such as Y12,22, inducing the non-secular coupling between

Scholes, Nat. Comm. 5, 4933 (2014). 17J. Lim, D. Palecek, F. Caycedo-Soler, C. N. Lincoln, J. Prior, H.
von Berlepsch, S. F. Huelga, M. B. Plenio, D. Zigmantas and J.

exciton populations and inter-excitonic coherences. The

Hauer, Nat. Comm. 6, 7755 (2015).

19

18E. Cassette, R. D. Pensack, B. Mahler, and G. D. Scholes, Nat.
Comm. 6, 6086 (2015). 19L. Bolzonello, F. Fassioli, and E. Collini, J. Phys. Chem. Lett. 7,
4996 (2016). 20A. De Sio, F. Troiani, M. Maiuri, J. Rehault, E. Sommer, J. Lim,
S. F. Huelga, M. B. Plenio, C. A. Rozzi, G. Cerullo, E. Molinari
and C. Lienau, Nat. Comm. 7, 13742 (2016). 21V. Butkus, D. Zigmantas, L. Valkunas, and D. Abramavicius,
Chem. Phys. Lett. 545, 40 (2012). 22V. Tiwari, W. K. Peters and D. M. Jonas, Proc. Natl. Acad. Sci.
USA 110, 1203 (2013). 23J. Prior, A. W. Chin, S. F. Huelga, and M. B. Plenio, Phys. Rev.
Lett. 105, 050404 (2010). 24A. W. Chin, A. Datta, F. Caruso, S. F. Huelga, and M. B. Plenio,
New J. Phys. 12, 065002 (2010). 25N. Christensson, H. F. Kauffmann, T. Pullerits, and T. Mancal,
J. Phys. Chem. B 116, 7449 (2012). 26A. W. Chin, J. Prior, R. Rosenbach, F. Caycedo-Soler, S. F.
Huelga and M. B. Plenio, Nat. Phys. 9, 113 (2013). 27M. B. Plenio, J. Almeida and S. F. Huelga, J. Chem. Phys. 139,
235102 (2013). 28A. Chenu, N. Christensson, H. F. Kauffmann, and T. Mancal,
Scientific Reports 3, 2029 (2013). 29V. Butkus, L. Valkunas and D. Abramavicius, J. Chem. Phys.
140, 034306 (2014). 30E. Basinskaite, V. Butkus, L. Valkunas and D. Abramavicius,
Photosynth. Res. 121, 95 (2014). 31M. Yang, J. Mol. Spec. 239, 108 (2006). 32S. Polyutov, O. Kuhn and T. Pullerits, Chem. Phys. 394, 21
(2012). 33M. Schroter, S. D. Ivanov, J. Schulze, S. P. Polyutov, Y. Yan, T.
Pullerits, O. Kuhn, Phys. Rep. 567, 1 (2015). 34H.-G. Duan, P. Nalbach, V.I. Prokhorenko, S. Mukamel, and M.
Thorwart, New J. Phys. 17, 072002 (2015). 35H. Lee, Y.-C. Cheng and G. R. Fleming, Science 316, 1462
(2007). 36J. Jeske, D. J. Ing, M. B. Plenio, S. F. Huelga and J. H. Cole, J.
Chem. Phys. 142, 064104 (2015). 37D. A. Lidar, I. L. Chuang and K. B. Whaley, Phys. Rev. Lett.
81, 2594 (1998). 38J. Jeske, N. Vogt, and J. H. Cole, Phys. Rev. A 88, 062333
(2013). 39P. Nalbach, J. Eckel, and M. Thorwart, New J. Phys. 12, 065043
(2010). 40A. Ishizaki and G. R. Fleming, New J. Phys. 12, 055004 (2010). 41X. Chen and R. J. Silbey, J. Chem. Phys. 132, 204503 (2010). 42D. P. S. McCutcheon and A. Nazir, Phys. Rev. B 83, 165101
(2011). 43J. Lim, M. Tame, K. H. Yee, J.-S. Lee and J. Lee, New J. Phys.
16, 018001 (2014). 44D. Abramavicius and S. Mukamel, J. Chem. Phys. 134, 174504
(2011).

45J. Seibt and T. Pullerits, J. Chem. Phys. 141, 114106 (2014). 46O. Rancova, R. Jankowiak, and D. Abramavicius, J. Chem. Phys.
142, 212428 (2015). 47R. Venkatramani and S. Mukamel, J. Chem. Phys. 117, 11089
(2002). 48C. Olbrich, T. Jansen, J. Liebers, M. Aghtar, J. Strumpfer, K.
Schulten, J. Knoester, and U. Kleinekathofer, J. Phys. Chem. B
115, 8609 (2011). 49S. Shim, P. Rebentrost, S. Valleau, and A. Aspuru-Guzik, Bio-
phys. J. 102, 649 (2012). 50L. Viani, C. Curutchet, and B. Mennucci, J. Phys. Chem. Lett.
4, 372 (2013). 51J. Jeske and J. H. Cole, Phys. Rev. A 87, 052138 (2013). 52D. Abramavicius and S. Mukamel, J. Chem. Phys. 133, 064510
(2010). 53G. Panitchayangkoon, D. V. Voronine, D. Abramavicius, J. R.
Caram, N. H. C. Lewis, S. Mukamel, and G. S. Engel, Proc. Natl.
Acad. Sci. USA 108, 20908 (2011).
54C. Kreisbeck and T. Kramer, J. Phys. Chem. Lett. 3, 2828
(2012). 55A. Gelzinis, D. Abramavicius, and L. Valkunas, J. Chem. Phys.
142, 154107 (2015). 56J. Seibt and T. Pullerits, J. Phys. Chem. C 117, 18728 (2013). 57H. Li, A. D. Bristow, M. E. Siemens, G. Moody and S. T. Cundiff,
Nat. Comm. 4, 1390 (2013). 58J. Yuen-Zhou, D. H. Arias, D. M. Eisele, C. P. Steiner, J. J. Krich,
M. G. Bawendi, K. A. Nelson, and A. Aspuru-Guzik, ACS Nano
8, 5527 (2014). 59J. Roden, W. T. Strunz and A. Eisfeld, J. Chem. Phys. 134,
034902 (2011). 60J. Iles-Smith, A. G. Dijkstra, N. Lambert and A. Nazir, J. Chem.
Phys, 144, 044110 (2016). 61Y. Tanimura, J. Phys. Soc. Jpn. 75, 082001 (2006). 62J. Strumpfer and K. Schulten, J. Chem. Phys. 134, 095102
(2011). 63F. Mascherpa, A. Smirne, S. F. Huelga, M. B. Plenio,
arXiv:1611.03377 (2016). 64T. Mancal and F. Sanda, Chem. Phys. Lett. 530, 140 (2012). 65J. Olsina and T. Mancal, Chem. Phys. 404, 103 (2012). 66A. Nemeth, F. Milota, T. Mancal, V. Lukes, H. F. Kauffmann
and J. Sperling, Chem. Phys. Lett. 459, 94 (2008). 67T.-C. Dinh and T. Renger, J. Chem. Phys. 142, 034104 (2015). 68J. Ma and J. Cao, J. Chem. Phys. 142, 094106 (2015). 69T. Brixner, T. Mancal, I. V. Stiopkin, and G. R. Fleming, J.
Chem. Phys. 121, 4221 (2004). 70P. Kjellberg, B. Bruggemann, and T. Pullerits, Phys. Rev. B 74,
024303 (2006). 71J. Adolphs and T. Renger, Biophys. J. 91, 2778 (2006). 72P. Hamm and M. T. Zanni, Concepts and Methods of 2D Infrared
Spectroscopy (Cambridge University Press, 2011).