{'fraction_non_alphanumeric': 0.031023156511742485, 'fraction_numerical': 0.01601248152405978, 'mean_word_length': 4.760192980796519, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 9, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
The nature of the low energy band of the Fenna-Matthews-Olson complex: vibronic signatures

Based entirely upon actual experimental observations on electron-phonon coupling, we develop a theoretical framework to show that the lowest energy band of the Fenna- Matthews-Olson (FMO) complex exhibits observable features due to the quantum nature of the vibrational manifolds present in its chromophores. The study of linear spectra provides us with the basis to understand the dynamical features arising from the vibronic structure in non-linear spectra in a progressive fashion, starting from a microscopic model to finally performing an inhomogenous average. We show that the discreteness of the vibronic structure can be witnessed by probing the diagonal peaks of the non-linear spectra by means of a relative phase shift in the waiting time resolved signal. Moreover, we demonstrate the photon-echo and non-rephasing paths are sensitive to different harmonics in the vibrational manifold when static disorder is taken into account. Supported by analytical and numerical calculations, we show that nondiagonal resonances in the 2D spectra in the waiting time, further capture the discreteness of vibrations through a modulation of the amplitude without any effect in the signal intrinsic frequency. This fact generates a signal that is highly sensitive to correlations in the static disorder of the excitonic energy albeit protected against dephasing due to inhomogeneities of the vibrational ensemble.


I. INTRODUCTION
The dynamics of photosynthetic complexes, such as the Fenna-Matthews-Olson (FMO) complex, have received rapidly increasing attention as an example of how nature appears to favor an intermediate noise regime, as suggested by the discovery of experimental evidence for longlasting coherent oscillations present in signals arising from non-linear spectroscopy protocols 1,2 , and the theoretical discovery that environmental fluctuations play a crucial role in explaining the very high efficiencies for excitation energy transport in such complexes [3][4][5][6] . The FMO complex is a trimer, with individual subunits composed of eight strongly coupled chromophores that serve as a link between the Chlorosome antennae and the reaction center of green sulfur bacteria 7,8 . Chromophores have a spatial extension that supports the existence of intramolecular vibrations [9][10][11][12] . These vibrations play an important role when the chromophore structures are perturbed due to the electronic charge redistribution upon photon excitation.
Non-linear spectroscopy has proved to be useful in order to unveil the dynamics involved in excitonic transfer of light harvesting complexes, due to the fact that it is sensitive to excitonic quantum superpositions, i.e., excitonic coherences. Two dimensional (2D) non-linear spectroscopy can resolve the third-order polarization of the electronic system, from the heterodyne detected spectrally resolved signal S(t 1 , t 2 , ω 3 ), arising from photoexcitation of three consecutive pulses with wave vectors k 1 , k 2 and k 3 , separated by time intervals t 1 and t 2 13 . Experimentally, it is customary to Fourier transform the time t 1 dimension to yield S(ω 1 , t 2 , ω 3 ) in order to generate a 2D spectra parametrized by the waiting time t 2 .
This spectra allows inquiring into the dynamics of excitation transfer in photosynthetic light harvesting complexes.
The appreciation of the role of discrete vibrational modes on the excitation transfer process has evolved over time 14 , and additional signatures of their presence in the 2D signal will contribute to a thorough understanding. The main purpose of this paper is to show that the experimental evidence points towards the presence of a discrete vibronic structure in the lowest energy band of FMO, and that the dynamical contribution of these vibronic modes can be extracted from the information available through non-linear spectra protocols. Particularly, the multiplicity of states resolved in the low energy band in the FMO has been the subject of intense debate. It has not been possible to date, to uniquely assign its nature to electronic contributions. In this manuscript we propose a model for the lowest energy band in the FMO that includes explicitly the most significant resonances of its environment spectral density, which is consistent with the available experimental evidence achieved through hole burning 11,12,15 , fluorescence line-narrowing 12,16 , and accumulated photon echo 17 . We explore the possibility to detect the coupling of intra-molecular low energy vibrations of chromophores, to exciton dynamics, and shed light on the nature of the lowest energy band of the FMO (to be referred as B825 from its absorbance maximum). We show that the experimental evidence can be explained with the existence of a representative vibronic resonance with energy ≃ 36 cm −1 and, beyond current suppositions that neglect higher order contributions, its overtone. Interestingly, this resonance has also been specifically addressed in resonant Raman spectra of the special pair BChla in purple bacteria reaction centers 18 , inferred through the Stokes shift in LHC-II protein complexes Chla pigment of higher plants 19 , and exists within the energy range expected for low-frequency vibrations of helical structures in proteins 20,21 . This strongly suggests that its origin are either the absorbing chromophores or a ubiquitous protein feature such as the α-helices both of which, are therefore not specific to the FMO structure.
In order to introduce our model we comment on the outcome of different experimental examinations in the B825, and qualitatively reproduce their results with the assumption of the presence of a single excitonic state, a discrete vibronic dynamical contribution, and its overtone. The dynamical contribution of this vibronic resonance is studied through the phase difference among signals obtained in the neighborhood of an electronic diagonal peak in the 2D spectra, which allows to distinguish the fundamental from higher vibrational harmonics. By unraveling the underlying vibrational dynamics in the 2D spectroscopy time domain, we study the characteristics of the beating signal that arise due to this electronphonon interaction (the phonon term will be used for the intramolecular vibrations of our interest). Accordingly, we explore the dynamics of non-diagonal contributions in the 2D spectra, to show that a modulation of the amplitude in the excitonic coherence signal arises due to these low frequency vibrational modes. We obtain an analytical expression for the excitonic coherence that explains this beating pattern as a result of vibrational wave-packet motion. We highlight that the coherent oscillations are extremely sensitive to the degree of correlation in the static disorder of neighboring chromophores and propose a physical picture involving the spatially extended α-helices of the protein scaffold, that is able to explain the observed long lifetime for excitonic coherence. On the other hand, the modulation of the signal amplitude due to vibrational wave packet motion does not alter the intrinsic frequency of oscillations of the polarization signal in 2D spectra and therefore, produces no dephasing due to inhomogeneities of the vibrational manifold. The latter effect allows us to explain recent experimental results where, even though structural modification in the FMO chromophores were accomplished, the dephasing rate of the electronic coherence in such a complex remained unaltered 22 .
A. Brief history of the FMO low energy band B825 Due to the energetic difference between the first exciton and higher lying excitonic states, the lowest energy band in the FMO complex has become the optical transition in which to study and understand, the extremely long lifetime of excitonic coherences found in this photosynthetic complex 1 . Nevertheless, the results obtained with different experimental techniques have not reached a consensus on the nature of this energy band.
The first hole burning and fluorescence line-narrowing experiments performed in FMO 12 , showed two components in the B825, whose absorption differed by ∆ω ≈ 30 cm −1 . The possibility that these frequencies were zero-phonon and phonon side bands was considered but ruled out at the time, because the energy difference between emission and absorption maxima in the B825 was ∆E ≃40 cm −1 , and could not be reconciled with the calculated Huang-Rhys factor s = 0.3 from measured absorbance difference, that resulted in a Stokes shift 2sω = 2(0.3)30 ≈ 20 cm −1 = ∆E. Rather, these levels were regarded as excitonic states arising from delocalized excitonic components of the FMO trimer 12 , a nature corroborated by theoretical calculations of the whole FMO trimer 23 that suggested the existence of equivalent BChl's 7 (using the numbering scheme of Ref. 7 ) on different units as the origin of the three excitonic states on the B825. Even though linear dichroism confirmed delocalization on B825 due to the broad linewidth in tripletminus-singlet spectra 24 , the linear dichroic absorption with respect to a continuous microwave field showed that the angle among the triplet and Q y transition dipoles in the B825 is constant throughout the band 25 , and close to that observed in single BChla 26 . Since the triplet and singlet orientations are sensitive to excitonic delocalization, the aforementioned equivalence among the B825 and single BChla, suggested that such a band was composed of localized excitations within a single BChla. This last turn of events was even further corroborated when the simultaneous fitting of absorption spectra, triplet-minussinglet and linear dichroic absorbance, concluded that the B825 was composed of localized excitons on the BChl 3 27 . Shortly afterwards, Stark hole burning studies 15 compared the purple bacteria B800 band (composed of almost completely localized excitonic states) with the B825 to provide further support that in the B825, the inhomogeneous broadening among equivalent BChls over the trimer generates excitons localized on single monomers of the FMO. The low energy band was successfully fitted with three gaussians, that were thought to represent the contributions of the localized lowest exciton states of each subunit of the FMO trimer 28 . The view that the contributions in B825 would come from energetic inequivalence arising from structural heterogeneity was also supported by the absence of excitonic dipole moment redistribution, probed by polarized hole burning for this energy band 11 . In that work, two important results were exposed. First, a greater homogeneous linewidth was found with greater energies within the B825, and confirmed the trend found by accumulated photon echo, of higher dephasing rates as the involved pulses were blue-shifted 17 . Second, a very rich structure of satellite low energy holes regarded as intramolecular vibrations (with the biggest components having energies of 36 and 72 cm −1 ) in the B825 was first observed.
The trend of increased linewidth at higher energy in the B825 was attributed then, to cascade excitation Förster energy transfer among individual chromophores on different subunits. However, Monte Carlo simulations of the Hamiltonian under the hypothesis of these three different chromophores being BChl 3 29 , were not able to simultaneously fit the lineshape and the gradual variation of the homogeneous linewidth and dephasing times observed. Almost concurrently with hole-burning experiments, the technique of fluorescence line-narrowing 16 successfully determined the phonon spectral density and confirmed the hole-burning experiments 11 about the existence of discrete phonons with frequencies 36 and 70 cm −1 moderately interacting with excitons. The possibility that multiple dynamical components on this band could arise from these modes, was first suggested.
Difference fluorescence line narrowing 30 provided further evidence of three contributions on the B825, and suggested once again, the possibility of exciton delocalization in the lowest exciton state since new simulations of the B800 in purple bacteria harvesting structure, showed a delocalization length across 2-3 pigments 31 . Recently, molecular dynamics simulations 32 showed an appreciable narrowing of the of site energies in a trimer compared to a monomer. This narrowing renders infeasible that static structural inhomogeneities in the trimer are the cause of the B825 multiple dynamic contributions.
At present, the consensus view is that of a single excitonic state in the low energy band to describe the FMO optical spectra. This view was corroborated with more detailed calculations of the FMO electronic Hamiltonian 33 , that have shown increasing success in fitting optical spectra, when linewidth assignments of such transitions involved a fit of the spectral density coarse grained envelope with a global Huang-Rhys factor s ≃ 0.5. However, no model has been able explain simultaneously the facts that within the B825, 1) there is a single electronic excitation where at least two representative resonances appear, 2) a marked trend of greater homogeneous linewidths (or dephasing rates) while probing higher energies was found. Beyond these points (highlighted in a recent review 34 ), the fact that 3) a very rich intramolecular vibrational structure is observed but never explicitly accounted for. This suggests the lack of a complete understanding of the low energy band of FMO.

II. THE MODEL
Our immediate purpose is to explain the possible nature of the multiple resonances found with increasing linewidth for blue-shifted probing of the B825. At the same time, the model must be consistent with the approximate 30 cm −1 absorption energy difference among these transitions, with the observed Stokes shift and with the fluorescent line narrowing spectral density. The dipole moment of such transitions must not exhibit variations in polarization while scanning this energy band. Let us examine these issues in detail. The most recent Huang-Rhys factor s = 0.5 describing the global electron-vibration interaction, is better able to explain the Stoke's shift of 40 cm −1 observed, occurring mainly due to the 36 cm −1 mode, i.e., 2sω = 2(0.5)36 ≈ 40 cm −1 . The whole phonon wing Huang-Rhys factor is an upper bound for the interaction strength among this mode and electronic excitations. In particular, hole burning spectra 11 assigned Huang-Rhys factor s = 0.12 for the most representative resonances, namely, resonances having energies of 36 and 72 cm −1 . Fluorescence line narrowing 16 reported for both these resonances, a Huang-Rhys of s = 0.01, on a procedure that might overestimate the weight of the phonon wing. As the authors explicitly point out, it is heuristically assumed that individual modes just refine details over the phonon wing, and therefore minimally contribute to build up the maximum of the spectral density function. The result of either model differ beyond the purely practical aspect of a fitting procedure, since a discrete contribution may lead to coherent evolution that will result in crucial features to be studied below. Therefore, if the coupling of this mode is upper bounded by the global Huang-Rhys factor s = 0.5, and lower bounded by the reported hole-burning Huang-Rhys value s = 0.12, a perturbative analysis is dubious. Moreover, the vibronic contribution does not represent a redistribution of dipole moment, allowing a single polarization throughout the B825. Hence, we proceed with the explicit inclusion of this mode.
The Hamiltonian for electron-vibration interaction reads where σ i and b (+) k are usual spin 1/2 Pauli matrices, and annihilation (creation) operators of the kth boson field, the electron-vibration coupling strength √ s k ω k is given in terms of the phonon frequency ω k and Huang-Rhys factor s k . In the present section, our aim is to consistently describe the homogeneous broadening trend as the probing laser field is blue shifted when vibrations are explicitly included. To that end we consider a single vibration and the laser interaction with the charge distribution of the chromophore under the rotating wave approximation (RWA) in a frame rotating at the laser frequency, we obtain a hamiltonian: where ∆ = ω 0 − ω L and Ω are laser detuning from the energy of the zero phonon line (ZPL) electronic transition and Rabi frequency, respectively. Now, the broadening of the resonances is of utmost importance in the historical discussion and is required for the analysis of spectra, as it influences the interaction of the system with the environment. The homogeneous linewidth of the resonances will primarily depend on the different processes inducing decoherence beyond the ensemble average. Accumulated photon echo 17 fitted the dephasing rate in 1.8-50 K range with the expression: for the homogeneous dephasing time T 2 , and highlight the main processes involved: population relaxation, pure dephasing from spectral diffusion promoted by the surrounding protein, and a thermally induced activation process. Given that the FMO complexes used in optical experiments lack the RC, the first term in eq.(3) is lower bounded by the fluorescent relaxation. The second term, is attributed to dynamics of tunneling among energy minima of different structural configurations in the protein, and characterized by a rate that obeys a power law temperature dependence ∝ T 1.3 , typical of glassy hosts. The third terms involves thermally activated transitions not specified in 17 , that accounts here for the thermal equilibration of the vibrational manifold. Note that eq.(3) only describes thermal relaxation, and therefore is best suited to the red edge of the B825. A fit to the red edge photon echo signal was accomplished with free parameters {T 2 , γ z , γ, δE} = {420 ps,0.09 ns −1 /K 1.3 , 23.6 ns −1 , 15 cm −1 }. In our model the most prominent thermally activated process occurs among states whose energies involve boson quanta of energy 36 cm −1 which restricts δE = 36cm −1 . A good agreement among both sets of parameters is accomplished for a constant γ = 70.09ns −1 =2.36 cm −1 (results not shown). In a first model, these mechanisms can be described under the Born-Markov approximation with super- with L 1 = √ Γσ − for population relaxation (σ −(+) are lowering (rising) Pauli spin ladder operators) and L 2 = 2(n + 1)γb, L 3 = √ 2nγb + describing thermalization to temperature T = 1/βk of the mode having an equilibrium average quantan = (e βω − 1) −135 . The factor 2 in this latter pair is introduced since the dephasing produced is half the rate of relaxation/absorption processes. The protein configuration tunneling induced dephasing is accounted by the operator L 4 = γ z T 1.3 σ z . The magnitude of rates Γ, γ and γ z are hence provided from the dephasing rate obtained through accumulated photon echo technique 17 .

A. Linear spectra
The fluorescence intensity in the line narrowing signal is proportional to the excited state population in the stationary state of the chromophore σ + σ − t→∞ , obtained from eq.(4), and presented in Fig.1(a). First, note that with both the upper and the lower bounds of the Huang-Rhys factor, three fluorescent transitions (highlighted by arrows) can be resolved: the zero-phonon line, the phonon side band, and the overtone, which coincides with the 70-72 cm −1 resonance reported in 11 . It should be noted that there is no need to include another vibration 70-72 cm −1 since a single mode and its harmonic lead to the most representative resonances at ≃ 36 and ≃ 72 cm −1 . Increasing the Huang-Rhys factor to s = 0.5 results in a higher overlap between levels having different boson quanta in the excited and ground electronic states, and therefore a redistribution of fluorescence intensity increases the contribution of higher harmonics. It should be noted that even a moderate Huang-Rhys factor s = 0.12 results in three dynamical contributions at T =77 K, as shown in the inset of Fig.1(a). At this temperature, the homogeneous phonon side band is commensurate with the zero-phonon line, and therefore, even in the presence of ensemble inhomogeneities, both transitions should be of importance in the description of the FMO in 2D spectroscopy (T =77 K, Ref. 1 ).
Having captured these three dynamical contributions, we proceed with the analysis of the individual linewidths. The hole burning experiments, sample the absorbance A(ω) proportional to the Fourier transform (FT) of a two-time correlation function, explicitly: which is presented in Fig.1(b) (calculated with the quantum regression theorem according to the Lindblad form master equation eq.(4)) for the relevant burning frequencies. In this figure, it is shown that higher energy on the burning beam generates greater linewidths of the holes in the absorption spectrum. With this model, the linewidths, say γ 0 , γ 1 , γ 2 , of the ZPL, the phonon sideband, and overtone, respectively, show a difference that persists for the whole range of temperatures measured as shown in the inset of this figure. This is in complete agreement with photon echo experiments (1.8-50 K, 17 ). The broadening can be understood from transitions involving n-th vibrational states, that become more populated with increasing temperature and whose thermalization rate is proportional to the excitation number L 2(3) |n n| ∝ n. The model also captures two important facts found in hole burning and photon echo experiments. Firstly, the greatest difference between linewidths occurs at low temperatures with a trend to converge as temperature is increased 11 . Secondly, before convergence is reached to equal linewidths, the three components are widening in the same fashion 17 . This behavior illustrates the prevalence of different mechanisms depending on the temperature. At low T , spontaneous emission and thermalization decay are the key mechanisms, the former being independent of temperature. At intermediate tem- malization as the dominating dephasing mechanism induces equal increment in all linewidths with temperature.
In practice, the reddest photon echo signal according to eq.(3) from accumulated photon data 17 , extrapolates to a dephasing rate of 3.3 cm −1 . Based in our theoretical results, this rate is ≈ γ 0 /2 = 6.2/2 cm −1 (see the inset Fig.1(b) at T=77 K). The difference of the data with our ZPL linewidth can be traced back to mixing of echoes originating from the ZPL and the side-band within the ensemble.
At this stage, we conclude that by inclusion of a vibrational mode with frequency of 36 cm −1 and an Huang-Rhys factor s ≥ 0.12, we are able to capture three dynamical contributions with equal polarization, whose position in absorption spectra and linewidths trend follow the experimental observation in the lowest energy band of the FMO. Without inclusion of a such mode, the absorptive response showing three dynamical contributions below T=10 K 11 cannot be explained relying solely on an inhomogeneously broadened single excitation electronic contribution.
Additional information is available in the 2D nonlinear spectra that might better characterize this low energy vibration, and further support (or reject) its role in the additional dynamical component on the lowest energy band of the FMO.

Background
In 2D spectroscopy, the heterodyne detected signal S(ω 1 , t 2 , ω 3 ) allows straightforward comparison with the FT of the theoretical analogue for a two level electronical system with energy difference ω eg among excited and ground states |e , |g respectively. In this case, the response functions R 2 and R 3 (R 1 and R 4 ) usually termed as rephasing (non-rephasing) paths, allow calculation of signal in the direction of observation k s = −k 1 + k 2 + k 3 (k s = k 1 − k 2 + k 3 ) surviving the rotating wave approximation. In order to present the data in a single quadrant, the FT for the non-rephasing and rephasing signals have conjugate variables {−ω 1 , ω 3 } and {ω 1 , ω 3 }, respectively. In detail, where upper (lower) signs correspond to rephasing (nonrephasing) terms, respectively. In the above expression, µ is the electronic transition dipole moment and the function f n (t 3 , t 2 , t 1 ) depends on the lineshape function g(t) evaluated in combinations of these three time intervals 13 . The lineshape function of a Brownian oscillator describes the contribution to the response function from a discrete vibration with frequency ω 13 . In this particular case, the function g(t) generates through f n (t) 36 , an increase of the excitonic frequency by the reorganization energy sω, ω eg → ω eg + sω, produces a global reduction of the signal by the Debye-Waller factor squared, e −2s coth(βω/2) , and has additional time dependent terms of the form e s(i sin(ωt)−coth(βω/2) cos(ωt)) . The Fourier transforms required for the 2D plots are performed here analytically from a Taylor series expansion e Os ≈ 1+Os+ (Os) 2 +· · · . The exponential function has a rapid convergence on this expansion for s coth(βω/2) ≪ 1, partially fulfilled for the Huang-Rhys factors of interest at the actual temperature of experiments (0.12 ≤ s ≤ 0.5, ω = 36 cm −1 , T=77 K; s coth(βω/2) = 0.37 − 1.53). A first order expansion of this kind was developed in Ref. 36 . In the following we further develop the expansion to second order, to study the possibility of resolving the higher harmonics (conspicuous in the absorption and fluorescence intensity) and develop a set of features able to better characterize vibronic resonances in the 2D spectra. The decay of the time resolved heterodyne signal in the impulsive limit, is equivalent to the linewidths of spectrally resolved observation under continuous-wave illumination. Henceforth, the effects of the bath in the 2D spectra will be introduced in the aforementioned expansion for terms proportional to e ±iωt1 , e ±iωt3 , e ±2iωt1 , e ±2iωt3 , with the Fourier transform argument ω 1(3) + iγ 1 on the former two, and ω 1(3) + iγ 2 on the latter pair, to include the homogeneous broadening for side-bands and overtone. The terms lacking of any of these factors, are Fourier transformed with an argument ω 1(3) + iγ 0 modeling the width of the zero phonon line. The addition of these imaginary contributions in the Fourier transform argument is consistent with a decay to the stationary state in the waiting time domain due to the substitution e iωt2 → e iωt2 e −γ1t2 and e 2iωt2 → e 2iωt2 e −γ2t2 .

Vibronic features in 2D spectra: Diagonal peaks
Genuine vibronic sidebands in 2D spectra are challenging to discriminate from the electronic contributions. Fortunately the exciton 1, |ψ 1 , in the FMO, has an energy separation > 150cm −1 from all other excitonic states which allows its spectral resolution from the other electronic contributions 1,37 . The 2D spectra can be used in order to characterize the most salient features of the vibronic spectra in this lowest energy band, such as the resonances found at 36 and 72 cm −1 .
The result of the second order expansion of the homogeneously broadened 2D spectra for exciton 1 interacting with a vibrational mode at different waiting times t 2 , is presented in Fig.2(a)-(b) for the rephasing and non-rephasing contributions, respectively. A hyperbolic sine scale is used to highlight small features. First, note that in these figures, an asymmetry is prominent in the absorption dimension. This is highlighted in the zero-phonon emission profile (continuous line), where the difference in height of the peaks at ω 3 = ω eg and ω 1 =ω eg ± ω accounts for Stokes and anti-Stokes contributions due to differences among populations of the initial canonical state in the vibrational manifolds, ρ = n e −βnω |n n| /Tr{·}.
Second, notice that the rephasing and non-rephasing contributions to the response function are sensitive to different processes on the vibrational manifold at a given waiting time. In detail, at zero waiting time the vibrational structure around the ZPL peak is prominent in the rephasing spectra, where even the overtone is conspicuous at ω 1 =ω eg + 2ω, while such structures are almost inexistent in the non-rephasing signal at the same waiting time. For a waiting time t 2 ∼ = π/ω, corresponding to half a vibrational period, the behaviors are exchanged leading to a non-rephasing contribution with very well resolved side-bands 36 .
A better contrast of the side-bands is possible in homogeneously broadened spectra, when instead of following the established procedure eq.(6) of adding the rephasing and non-rephasing responses (Fig.2(c)), we proceed with their subtraction (Fig.2(d)).The relative phases of rephasing and non-rephasing signals induce a stronger ZPL when added, but strongly suppress the weight of the ZPL to make the side-band structure gain relative weight when subtracted. This example allows to point out that the relative phases of the components that lead to the total signal are sensitive to algebraic operations over the individual contributions. However, this scenario where the average among the ensemble is not performed, neglects the static inhomogeneity of importance in the ensemble required to individually address rephasing and non-rephasing contributions in different spatial directions.
In the following, the exciton energy inhomogeneities are studied, considering the second order expansion of the line width function. The localized nature of exciton 1 (at pigment 3) has been confirmed with numerical simulations 39 , using a site energy distribution width of 110 cm −1 which due small exchange narrowing leads to a value of σ E = 102 cm −1 for the total width of the B825 band. Recent numerical simulations of the FMO 2D spectra used a comparable width for the exciton 1 energy distribution (100 cm −1 ) 40 . The diagonal 1-1 peak for rephasing and non-rephasing signals is presented in Fig.3(a)-(b), respectively, where the inhomogeneous broadening is introduced according to a gaussian distribution of the excitonic energy with standard deviation σ E = 102 cm −1 . Interestingly, a very dramatic change occurs when moving from the homogeneous to the inhomogeneously broadened ensemble, where each path (rephasing/non-rephasing) is sensitive to addition of signals arising from different processes. On the one hand, Fig.3(a) shows the zero-phonon emission line (continuous) of the rephasing component, in order to highlight side-bands located at ω 1 − ω 3 = ±ω. On the other hand, the non-rephasing part, Fig.3(b), is sensitive to the overtone, highlighted by the continuous lines over the surface at the ZPL and overtones at ω 1 − ω 3 = ±2ω. A better contrast of the vibronic structure is accomplished in the non-rephasing 2D spectra due to the drastic reduction of the ZPL. Hence, the inclusion of inhomogeneous broadening does not simply regards to a general widening of the homogeneous resonances breath, since it affects in a different fashion to each of the resonances. In short, a discussion of this reduction will allow us to understand the qualitative differences among homogeneous and inhomogeneously broadened spectra.
The inhomogeneous broadening in the excitonic frequencies corresponds to widening along the diagonal in the 2D spectra, with a ZPL maximum at ω 1 = ω 3 . We have chosen the waiting time long enough in Fig.3(a)-(b), to avoid the t 2 dependence on the 2D spectra. These results remain valid for different waiting times as will be demonstrated below. The expression for the zeroth order (waiting time independent) contribution of the ZPL reads: where the upper and lower signs correspond to rephasing and non-rephasing expressions, respectively. The real part of the rephasing contribution is always positive for any value of the exciton 1 energy E ψ1 , while that of the non-rephasing part, is positive for E ψ1 − ω 1 < γ 0 and negative otherwise. The imaginary part of the rephasing signal is zero at the diagonal, while that of the nonrephasing, is positive for E ψ1 < ω 1 and negative otherwise. This situation leads to have at a single point in the diagonal of the 2D spectra, only positive contributions for the rephasing part that will add up to build the ZPL peak. On the other hand, in the nonrephasing path the signals have varying signs depending on the exciton energy, which when added up in the inhomogeneous sample, will lead to a decreased ZPL. This fact is quantified in Fig.3(c)-(e), where the weighted contribution e −(E ψ 1 −ωeg ) 2 /2σ 2 E S(ω 1 , ω 3 , E ψ1 ) of the complete second order expansion as a function of the difference E ψ1 −ω eg are presented. The Fig.3(c) shows that the signal at ω 1 =ω 3 =ω eg has alternating (equal) signs for real and imaginary parts of the non-rephasing (rephasing) signal. The trend continues (Fig.3(d)) for points regarding the phonon side-band (ω 1 =ω eg + ω, ω 3 =ω eg ) where the real (imaginary) part of the rephasing signal presents over a great range of the whole inhomogeneous distribution, positive (negative) values that constructively add to  3. In (a) and (b), rephasing and non-rephasing 2D spectra of exciton 1 diagonal peak of FMO complex, broadened by excitonic energy inhomogeneities. In (c), (d) and (e) are presented the real (continuous, blue online) and imaginary (dashed, red online) contributions to the inhomogeneous signal at points (ω1 = ω3 =ωeg), (ω1 =ωeg + ω, ω3 =ωeg) and (ω1 =ωeg + 2ω, ω3 =ωeg) respectively, for complexes whose energy differ E ψ 1 −ωeg from the ensemble average as explained in the text. In (f) and (g), rephasing and non-rephasing 2D spectra of exciton 1 diagonal peak of FMO complex, broadened by both excitonic and vibrational mode energy inhomogeneities. σE = 102 cm −1 , σω = 10cm −1 , t2 = 2 ps.
build up a appreciable side band in the rephasing path. However, rephasing and non-rephasing parts have positive and negative contributions along the inhomogeneous distribution from the signals arising at the overtone resonance (ω 1 =ω eg +2ω, ω 3 =ω eg ) which lead to a reduction commensurate in both. The concomitant reduction of the side band and drastic reduction of the ZPL for the nonrephasing direction, explains the surprising enhancement in the contrast for the signal arising from overtone, due to the (usually troublesome) sample inhomogeneities.
The inhomogeneities are not only ascribed to the exciton energy, but also the vibrational manifold. For instance, hole burning can select a subensemble having an optical transition frequency, but the inhomogeneities in the satellite structure due to the vibrational manifold will not be different from those of any other subensem-ble. Therefore, the width of the satellite structure in these experiments reflects that of the ensemble. The low energy satellite structure in the B825 was observed in hole burning 11 , with a width of 10 cm −1 for the ω = 36 cm −1 resonance. Hence, we proceed to perform an average of the vibrational frequency, according to a gaussian with σ ω = 10 cm −1 , and present the results for rephasing and non-rephasing paths in Fig.3(f)-(g), respectively. In the rephasing contribution, this average smears out the vibronic structure, while some of it remains in the non-rephasing path. Since the experiment involve macroscopic samples, the vibrational manifold ensemble inhomogeneities will lead to a reduced structure in the 2D spectra, and a simple inspection will not permit (in general) a straightforward inference of the presence of discrete vibronic contributions. Furthermore, the sampling     of the signal in t 1 is done during a finite time that imposes a minimum frequency, hence a maximum resolution among frequencies. Currently, studies of 2D spectra tracking the vibrational wave-packet motion, used a scanning range for t 1 from -100 to 100 fs, which allows a resolution in ∆ω 1 ≃ 52 cm −1 for rephasing/non-rephasing signals, or half of it, ∆ω 1 ≃ 26cm −1 for the total signal 38 . Accordingly, the discovery of vibronic structure in the ensemble becomes difficult relying on the inspection of a single 2D spectrum alone. However, the vibronic contributions can be detected from the waiting time t 2 , which induces features that allow us to infer the presence of discrete phonon modes. The additional structure whose peaks lie at ω 1 − ω 3 = ±nω, imply rotations occurring in the vibrational manifold that, as will be seen, induce a phase shift in the signal on waiting time dimension. To illustrate our point, we start with the second order line-width expansion at the homogeneous level and focus on the real part of the signal concerning the product of absorptive (non-dispersive) spectra in ω 1 and ω 3 .
For instance, at ω 1 =ω eg + ω, ω 3 =ω eg the signal is constructed from the Liouville paths contribution shown in the Feynman diagrams of Fig.4(a), where time evolves from bottom to top over the vertical lines that represent the ket and bra evolution. The wiggly arrows denote interactions with electromagnetic pulses which excite (deexcite) the system when they are directed into (out of) the ket-bra vertical lines. The additional kets or bras, are used to clarify the transitions involved after a given excitation. The diagrams R 1 and R 2 (R 3 and R 4 ) involve rotations in the excited (ground) state manifold for the waiting time domain. Observation of the point ω 1 =ω eg + ω is generated when the bra (ket) e n+1 ( e n+1 ) is excited with the first pulse. On the other hand, at ω 3 =ω eg the evolution involves the same boson quanta in both bra and ket. Therefore a flip of the boson state n + 1 → n is required in order to end up in a population ket-bra combination before the trace operation is performed. For diagrams R 3 and R 4 this rotation can be accomplished when the second pulse de-excites the system to g n+1 and g n+1 , respectively. However, for R 1 and R 2 the rotation must happen in the waiting time, and therefore, the signal is sensitive to the the required flip in the vibrational manifold. This rotation can be explained along the mentioned mathematical procedure, but can also be understood by considering the Frank-Condon principle.
After pulsed excitation, the system is left in a superposition of eigenstates of the Hamiltonian (1), which also represent a vibronic progression |ñ whose energy differ by ω among consecutive levels in a given electronic manifold. As the thermal population involves the lowest energy levels, the contribution to the signal will primarily arise from the ZPL and the phonon side-band, approximately described by |φ ≈ α 0 sin(ωt) 0 + α 1 cos(ωt) 1 . Simply stated, the signal at ω 1(3) =ω eg ± nω is mainly contributed by the |ñ bosons state population, as a witness of excitation and dynamics of the vibrational wavepacket motion. The phase of the signal will vary as the 2D spectrum is swept reflecting the phase difference among populations of the levels being scanned. In the waiting time, Fig.4(b), the signal presents the discussed phase shift ≈ π among rephasing and non-rephasing contributions (remind Fig.3(a)-(b)) . However, this panel also shows a phase shift among ZPL and side band signals of π, to resemble the just mentioned population inversion. A population picture can be misleading for the phase difference among the side band and overtone, since this latter resonance develops a signal with strength pro-portional to s 2 , small enough to be contaminated by the tails of the other resonances. The right panel of Fig.4(b), shows that if the breath of the homogeneous broadened resonances is amplified in accordance with loosing the mode discreteness, the phase shift among different points in the 2D spectrum is lost when scanning the neighborhood of a diagonal peak. Hence, a smeared out interference of several indistinguishable contributions in a given area of the 2D spectrum reduces the phase shift differences.
In Fig.4(c) left panel, is shown the result of averaging static inhomogeneity to the signals resolved in the t 2 dimension. Here it can be first noticed the reduction of the ZPL and side band of the non-rephasing compared to rephasing path. Secondly, this procedure only slightly reduces the phase difference among rephasing and non-rephasing components. Third, the phase differences among resonances on a given path is persistent among ZPL, side band and overtone. A variation in the phase leads also to a variation of the slope ∂ t2 Re[S](ω 1 , ω 3 ) at a fixed t 2 , before the settlement of the stationary state. This variation is presented in the right panel of Fig.4(c), and specially highlights (by filling the area underneath) the change of sign of the signal slope for a fixed waiting time as a function of the coordinate ω 1 . Even though the actual values of this phase are affected by the presence of the ensemble inhomogeneities, its variation remains, and since its origin is traced back to the discreteness of the vibrational manifold, it can be used as a tool to diagnose the presence of vibronic discrete structure beyond the unavoidable presence of static disorder.
The 2D spectra consist of diagonal and non-diagonal peaks. We have covered the behavior of signatures in both homogeneous and inhomogeneously broadened samples for the diagonal peaks. Next, we consider a prominent signature that arises in the waiting time evolution of non-diagonal excitonic peaks. The single excitation peak on the model just presented is unable to capture the behavior of excitonic coherence. Moreover, it has been previously stated 1,37 , that extracting the dynamics of non-diagonal peaks in nonlinear spectra, that is the excitonic coherences, has the advantage of less congested spectra. Therefore, it is an important contribution to understand the features that will arise in the cross peak dynamics on the FMO light harvesting complex.
In the following, we investigate the excitonic coherence of a dimer interacting with localized vibrations in order to address the FMO exciton 1-2 cross peak, when the coupling to the ω = 36 cm −1 mode is accounted for. The sites with the greatest contribution on excitons 1 and 2, are chromophores historically 7 named as sites 3 and 4 7 (states |3 and |4 ). The Hamiltonian restricted to a single excitation for these two sites, each interacting with intramolecular vibrations is: where the last two terms correspond to the electronic Hamiltonian, and are annihilation (creation) operators of the ith vibration. Diagonalization of the electronic Hamiltonian leads to excitonic states describes the limiting cases of localized (|V | ≪ δ, a ≈ 0) and fully delocalized (|V | ≫ δ, a ≈ ±1) excitons. In order to determine the effect of the vibrational modes in the exciton dynamics in the particular case of the FMO excitons 1 and 2, we proceed with a numerical solution constrained to the experimentally determined parameters. The dipole-dipole couplings in agreement with linear absorption, circular dichroism and linear dichroism 33 , and the experimentally determined exciton transition energies 37 of the full FMO, were used in an iterative procedure as described in Ref. 37 to determine the site energies of 12142 cm −1 and 12315 cm −1 for sites 3 and 4, respectively. Then, the value a = 0.28 is determined from the dipole-dipole coupling (-53 cm −1 , 33 ) between sites 3 and 4. We constrain the reduced dimer model to match the delocalization degree a = 0.28 and the experimentally determined exciton 1 and 2 energy difference δ ψ1ψ2 = 158 cm −137 , which results in renormalized site energies δ = 135 cm −1 and coupling strength V = −41cm −1 . Given that our interest is the 1-2 cross peak of FMO whose beating signal is a manifestation of excitonic coherences, Fig.5(a) shows the result of the numerical calculation of excitonic coherence. In the absence of electron-vibration coupling, the amplitude of the coherence should be monotonically decreasing, since a dimer model will only have a single excitonic frequency in the one excitation manifold. The most prominent manifestation of the presence of the vibronic structure is observed here, as a modulation of the excitonic coherence amplitude. In order to provide a figure of merit beyond a qualitative estimation for this modulation, an analytical   study is much preferable.
In the excitonic basis, the projector |3 3| = 1 1+a 2 (|ψ 1 ψ 1 | + a 2 |ψ 2 ψ 2 | + a |ψ 1 ψ 2 | + a |ψ 2 ψ 1 |) with a similar expression for |4 4| by the exchanges ψ 1 → −ψ 2 and ψ 2 → ψ 1 , shows that in situations where δ is much greater than the Coulomb interaction, the Hamiltonian eq.(10) is diagonal in the excitonic basis. With this approximation the excitonic coherence Tr{|ψ 1 ψ 2 | ρ(t)} can be calculated using the total density matrix evolution ρ(t) = e iHt ρ(0)e −iHt . Given that our interest lies in the dynamics on the waiting time domain, we suppose an initial state ρ(0) which is a product of the ith oscillator thermal state ρ bi and an electronic superposition proportional to the site transition dipole moments, ρ e ≃ i,j µ i µ j |i j| = i µ ψi µ ψj |ψ i ψ j | to best describe the situation after the arrival of two broad band low intensity pulses within the single excitation subspace. After tracing out the exciton manifold and then noting that the trace over the oscillators can be arranged as a product of traces due to commutation among operators describing vibrations in different chromophores, we make use of the displacement operator This procedure leads to the interaction picture displacement operator U i (t) = e iH bi t U i e −iH bi t correlation function Tr bi U −1 i (−t)U i ρ bi , which for harmonic oscillators, is a standard result 13 . Hence follows: where g i (t) are the linewidth functions (remind eq.(8)) addressing the properties of vibrations in the ith site. As shown in the second line of the above equation, if oscillators are identical, the addition of the linewidth functions exp[g 3 (t) + g * 4 (t)] results in a real argument exp[−2s coth(β ω/2)(1 − cos(ωt))] presented (dashed) in Fig.5(a), that quantitatively reflect the amplitude modulation of the excitonic coherence. Also the complete solution, eq.(11) in dotted lines shows excellent agreement with full numerical simulations. This underline the fact that disregarding excitonic population transfer terms mediated by vibrational modes ∝ a 1+a 2 |ψ i ψ j | seems correct in the B825, for the description of the vibronic wave packet motion. Moreover, this agreement in the current situation where V ≃ δ, highlights that in a more general context, the periodic modulation of the excitonic coherence, i.e., of the signal amplitude, is a signature of coherent dynamics occurring in the quantized vibrational manifold. The analytical result shows that the modulation of amplitude in the excitonic coherence will present a frequency equal to that of the prominent vibrational mode, namely ω, and will have a depth proportional to the coupling of such vibration with the electronic transitions, via the Huang-Rhys factor s. Therefore, for practical purposes, the assumption that the Hamiltonian eq.(10) is diagonal in the exciton basis when accounting on vibronic effects on the 1 and 2 excitonic coherence in the FMO, is well supported. We proceed with this assumption to investigate whether the modulation of excitonic coherence amplitude can be readily observed in the 1-2 FMO cross-peak from 2D spectra. The diagrams in Fig.5(b) are the contributions for the excitonic cross peak that present excited state evolution in the waiting time at a frequency equal to the energy differenceω ψ2 −ω ψ1 = δ ψ1ψ2 , with corrections arising from the line-width functions. These oscillating terms carry dynamical information that provide access to the interference of the most salient dynamical contributions, and read: where µ ψi,f is the dipole moment of the transition |ψ i → |f , and f 1(2) (t 1 , t 2 , t 3 ) correspond to functions associated with the cumulant expansion of the linewidth functiond of the paths shown in Fig.5(b) under the assumption of a diagonal exciton-phonon interaction 41 : These expressions involve photon-echo rephasing contributions with linewidth functions of the singly excited exciton states |ψ 1 , |ψ 2 and doubly excited states |e 3 , e 4 = |f . The transformation of the coordinates X i in the exciton basis yields to , which allows calculation of the linewidth function g ij (t) in terms of the site linewidth function g(t). For identical vibrations in each site, it can be written g ij = C ij g(t), When a → 0, the above expressions fulfill f 1 (t 1 , t 2 t 3 ) = f 2 (t 1 , t 2 t 3 ), while µ ψ1,f = µ ψ2 = µ 4 and µ ψ2,f = µ ψ1 = µ 3 to produce the expected vanishment of non-diagonal peaks in absence of electronic coupling among sites. Following the same procedure used for the diagonal peak in the homogeneously broadened FT signal, we perform a first order Taylor expansion of e fi(t) and present in Figs.5(c)-(e) the result of the 2D non-diagonal peak following eq.(12). The 2D spectra at early waiting times in Fig.5(c), shows appreciable contribution from the side-band structure. At the coordinate {ω 1 , ω 3 } = {ω ψ1g ,ω ψ2g } (Fig.5(d)) the modulation of the signal envelope is absent and supports that its principal contribution will mainly arise from ZPL transitions. On the other hand, the signal arising from the satellite structure {ω 1 , ω 3 } = {ω ψ1 + ω,ω ψ2 } shows an appreciable modulation of its amplitude. This modulation has a frequency ≃ ω, that can be accounted for in virtue of eq.(11) to the vibrational wave packet motion.
The modulation in the signal of the cross peaks makes it a suitable witness for discrete vibronic structure in the environment. We follow with the detailed study of the progression from homogeneous to inhomogeneous broad-ening of the signal, which leads to interesting conclusions.
C. Ensemble inhomogeneities: On the static nature of the difference among site and excitonic coherence decay Two different scenarios are shown in Fig.6(a) for the introduction of static disorder. In the first case (dashed line), the signal is averaged assuming uncorrelated static fluctuations, distributed according to the product of probability density functions from individual energies P (ω 1 )P (ω 2 ), using comparable inhomogeneous widths on each. In such a case, the damping of the oscillations occurs in a time-scale ≈ 100 − 200 fs, which does not represent the experimentally measured result 22,42 . By virtue of the microscopic model proposed here, the homogeneous decay presented in Fig.5(d)-(e) is similar, while the decay due to uncorrelated variation in the inhomogeneous average in Fig.6(a) dashed line, is far too fast, compared to the appreciable visibility experimentally found for oscillations well beyond t 2 = 1 ps 22,42 . Hence it follows that there must exist a degree of correlation in the static energy variations in order to reproduce the experimental evidence 22,42 . Such correlation has been recently outlined as a process of paramount importance for understanding the signal arising from macroscopic ensembles 43 . Studies from first principles with molecular dynamics 44 , suggest that the energy shift in individual BChls in the FMO are due to two contributions: 1) local, arising from different geometries of each BChl deformed within their fitting into the protein scaffold, 2) environmental, induced by the charge density of the protein due to its crystal structure. The latter process is generally accounted for by the charge density of α-helices which extend across spatially extended areas involving in general, neighboring BChls. For example, BChla 3, 4 and 5, are situated close to the α-helix 5 44 . Therefore, the electric fields generated by the protein are created in spatially extended regions including several chromophores within a single monomer, and is expected that a set of neighboring chromophores will experience Stark shifts from conformational changes of a single α helix, i.e., a correlation on energy shifts. If the shifts share sign and approximate magnitude, this situation will lead within the ensemble to small deviation from the average in the differenceω ψ2 −ω ψ1 .
The fact that the main contributions in excitonic states |ψ 1 and |ψ 2 come from the neighboring sites 3, and 4 and 5 respectively, will imply in this scenario a degree of correlation in the energy of these excitonic states. Hence, FMO monomers with exciton energiesω ψ1 +ǫ andω ψ2 +ǫ ′ fulfilling ǫ ′ 2 ≃ ǫ 2 ≈ σ 2 E , will result in ǫ ≃ ǫ ′ and (ǫ − ǫ ′ ) 2 ≪ σ E . With this assumption, is straightforward to average the signal according to the product of gaussian distributions P (ǫ), P (ǫ − ǫ ′ ) which describe the correspondent variations of exciton 1 energy and of the relative difference among excitons 1 and 2 energies, having standard deviation σ ǫ = σ E and σ ǫ−ǫ ′ ≪ σ ǫ . The result of the correlated averaging procedure is presented in Fig.6(a) (continuous line) which shows an appreciable amplitude of oscillations in the signal at t 2 = 2 ps, in good agreement with experimental evidence. It is important to highlight that this model reproduces the inhomogeneous width recorded in absorption spectra or in the ω 1 and ω 3 dimension in 2D spectra, while concomitantly, it is able to reproduce the enhanced lifetime of oscillation in the waiting time domain. This last conclusion requires a more thorough analysis. Remember the Feynman diagrams Fig.5(b), where the signal in the coherence time t 1 evolves at a frequency cycling among ground and excited electronic states, while that in the waiting time promotes oscillations among different excitonic ket-bra combinations. The local field will therefore lead in t 1 domain to a short dephasing decay time (as happens with the dashed line in Fig.6(a), ≈ 100−200 fs) proportional to the magnitude of the Stark shift due to the static fluctuations. The correlation in the contributing chromophores of the excitons allows (as presented in Fig.6(a)) the existence of one order of magnitude difference in the decay of the signal when compared to non-correlated inhomogeneity. Hence, regardless of the magnitude of the Stark shift but sensitive to the gradient of the field producing it, the dephasing rate in the waiting time t 2 will decay slower than in t 1 when correlated fluctuations in the Stark shifts are present. The extended local field induced Stark shift picture is further supported by the recent experimental finding 42 that cross-peaks associated with excitonic states whose main contributions involve neighboring chromophores have signals (in the waiting time domain) with smaller dephasing rates. It cannot be explained based on heuristic arguments that shifts should be correlated by having equal signs at a given time. A more detailed analysis based on the actual charge distribution over the α helices may ad-dress microscopically this possibility, but will not to be treated here.
The oscillation of the signal when correlated static inhomogeneities are present in Fig.6(a), resemble extremely well the behavior of the signal in the FMO 1-2 cross peak reported in 22,42 . This has been achieved by doubling the Huang-Rhys factor of this representative vibration well within the bounds 0.12 ≤ s < 0.5, s = 0.24, and decreasing the homogeneous width of the mode γ → γ/8 which is acceptable from the fact that fluorescence line narrowing or spectral hole burning already give an inhomogeneously broadened ensemble average of the vibrational manifolds. The lower bound s = 0.12 has been used to infer the minimal effect of the mode of interest and keep consistent with experimental evidence. Moreover, the discussion concerning the diagonal peaks enhancement due to inhomogeneity and phase shift remain valid for s = 0.24. Quantitatively, the enhanced Huang-Rhys factor results in an increased contrast of the vibronic bands compared to the ZPL. Experimentally, the modulation of electronic coherences has been ascribed solely to overlap from peaks 1-2 and 1-3, which would produce a beating according to the recently estimated Hamiltonian 37 , at a frequency of 38 cm −1 . However, this result highlights that such beating could arise solely from the discrete contribution of the most prominent vibronic mode. Hence, it must be underlined that many processes are observed at once in the 2D spectra and that beating of the amplitude in signals can arise from excitonic interference and electron-vibrational coupling. Recently, resolution of the waiting time dynamics at different cross peaks has been achieved 42 , and therefore, the data is available to explore whether an amplitude modulation having a frequency ≈ 36 cm −1 is a general feature.
In order to clarify whether vibrational coherence is responsible for the long lasting oscillations present in the non-linear response of the FMO complex, structural modifications were implemented 22 , whose effect was intended to change the vibrational manifold frequencies.
Interestingly, none of the modifications produced appreciable changes in the electronic coherence dephasing rate, nor the oscillation frequency of the signal. Note however that our result implies a change in the period of the modulating envelope, and not that of the underlying oscillations when structural modifications are generated. Moreover, the usual mechanism by which ensemble oscillations are damped from averaging dephased contributions, does not induce dephasing among signals regarding modes of different frequency. Fig.6(b) shows that increasing the vibrational mode frequency inhomogeneities, will only decrease the beating depth without enhancing the dephasing rate of the signal. We model the signal averaging according to a gaussian probability density function, with σ ω of 10 cm −1 in agreement with 11 . This figure shows that using a physically feasible Huang-Rhys factor is enough to present an appreciable modulation of the signal amplitude whose depth decreases as the vibrational mode inhomogeneities are larger. Accordingly, static vibrational mode heterogeneities will lead to a decreased modulation of the signal amplitude but will not influence the overall decay of the signal. In this way the robustness to structural modifications of the signal can be understood due to this novel non-dephasing effect over the ensemble inhomogenities.

IV. CONCLUSIONS AND PERSPECTIVES
We proposed a theoretical model for the low energy band of the FMO complex comprising a vibronic discrete structure and study the possibility to detect such vibronic dynamics with non-linear 2D optical spectroscopy. For experimentally determined parameters that characterize the electron-vibration coupling in FMO chromophores, the model is consistent with linear and non-linear optical spectra for the low energy band of the FMO complex. The main features of the vibrational side bands are variations of the signal phase in the neighborhood of diagonal peaks and a modulation of the amplitude in the waiting time resolved signal at cross peaks, with a frequency corresponding to that of the energy of the vibrational discrete mode. Complementarily, the robustness of excitonic coherence to dephasing can be further enhanced by correlations in the static energy disorder within a single monomer, presumably related with spatially extended electric field inhomogeneities that lead to correlated Stark shifts involving nearest neighbor chromophores.
We show that correlated static disorder can lead to sample inhomogeneities that induce about an order of magnitude enhancement of excitonic coherence lifetime as compared to the t 1 signal decay. When correlated static disorder is present, the discrete vibronic structure is not sensitive to ensemble inhomogeneities. Only under this conditions, the behavior of coherences at the single molecule level will represent the ensemble averaged signal. Lastly, the resilience of the waiting time domain os-cillations to structural modifications is understood in the 1-2 FMO cross peak as vibrations generate a modulation in the signal amplitude, and not a variation in the signal intrinsic frequency. The signal presents a dichotomic sensitivity to excitonic and vibrational features, where the former are extremely sensitive to inhomogeneities while the latter are almost unaltered. Noticeable effects in the waiting time 2D signal are enhanced by both discreteness and coupling strength of the prominent vibrations.
The fact that the mode here studied might resonantly couple other electronic transitions in the FMO, and therefore, be useful to enhance delocalization among excitons, is a possibility that remains to be explored. Our study suggests that a detailed understanding of the vibrational environment, even at low frequencies, is crucial to achieve a full understanding of the electronic dynamics of pigment-protein complexes.