Coherence oscillations measured in two-dimensional (2D) electronic spectra of pigment-protein complexes may have electronic, vibrational, or mixed-character vibronic origins, which depend on the degree of electronic-vibrational mixing. Oscillations from intrapigment vibrations can obscure the inter-site coherence lifetime of interest in elucidating the mechanisms of energy transfer in photosynthetic light-harvesting. Huang-Rhys factors (*S*) for low-frequency vibrations in Chlorophyll and Bacteriochlorophyll are quite small (*S* ≤ 0.05), so it is often assumed that these vibrations influence neither 2D spectra nor inter-site coherence dynamics. In this work, we explore the influence of *S* within this range on the oscillatory signatures in simulated 2D spectra of a pigment heterodimer. To visualize the inter-site coherence dynamics underlying the 2D spectra, we introduce a formalism which we call the “site-probe response.” By comparing the calculated 2D spectra with the site-probe response, we show that an on-resonance vibration with Huang-Rhys factor as small as *S* = 0.005 and the most strongly coupled off-resonance vibrations (*S* = 0.05) give rise to long-lived, purely vibrational coherences at 77 K. We moreover calculate the correlation between optical pump interactions and subsequent entanglement between sites, as measured by the concurrence. At 77 K, greater long-lived inter-site coherence and entanglement appear with increasing *S*. This dependence all but vanishes at physiological temperature, as environmentally induced fluctuations destroy the vibronic mixing.

## I. INTRODUCTION

In pigment-protein complexes (PPCs) devoted to light harvesting, light absorption generates molecular excitons,^{1} or singly excited states delocalized over the coupled pigment molecules (sites). Two-dimensional electronic spectroscopy (2D-ES)^{2} has become the method of choice for measuring exciton dynamics on ultrafast time scales.^{3} The technique provides information on electronic couplings, homogeneous and inhomogeneous lineshapes, and the kinetics of exciton relaxation and energy transfer from one state to another. In recent years, long-lived oscillatory signals have been measured in the 2D spectra of diverse photosynthetic complexes,^{4–10} including at physiological temperature.^{11,12} The oscillations^{13} are attributed to coherent superpositions of electronic, vibronic, and/or vibrational states. The evidence that electronic coherences may survive on sufficiently long time scales to affect the mechanism of electronic energy transfer (EET) in the photosynthetic apparatus has stimulated a large body of work at the intersection of biophysical chemistry and quantum information.^{14–17}

Understanding nature’s mechanisms for energy transfer over nanoscale distances is an important milestone toward replicating the robust and efficient solar energy conversions observed in photosynthesis. Long-lived coherence allows rapid energy displacement as a part of the unitary evolution of a superposition state. Further work pointed to coherence as a means for enhancing EET unidirectionality^{18,19} and robustness to energetic disorder.^{20} The environment, too, plays a role in efficient energy transfer. Both by dissipating energy and by destroying superpositions before they evolve back into their initially prepared state, environmentally induced fluctuations assist^{21–23} in delivering excitation to the lowest-energy pigment. Ishizaki and Fleming^{19} calculated that electronic coherence in the Fenna-Matthews-Olson (FMO) complex isolated from green sulfur bacteria could persist for 700 fs at 77 K and 300 fs at 300 K. At physiological temperature, the theory and experiment agree. However, oscillations measured at 77 K persist in the spectra for at least 1.5 ps.^{12} It has been suggested that the long-term oscillations are due to the influence of intra-molecular vibrations,^{24–29} and that in some systems, electronic-vibrational interactions enhance EET.^{30–34} Considerable interest arose in developing analyses that distinguish between vibrational and electronic components of oscillations in 2D-ES,^{35–46} as well as time resolved measurements selective to inter-site coherence and excitation delocalization. Polarization-dependent 2D-ES^{7,8,47} experiments can be site-selective, if transition dipoles on different pigments have sufficiently different orientations. Scholes and Smyth have suggested using excited-state absorption to measure electronic delocalization that is strongly dependent on a vibrational coordinate,^{48} and Tempelaar *et al.* have proposed time-resolved fluorescence experiments as a means to monitor the decay of the coherence length in molecular aggregates.^{49}

Mechanisms that balance inter-site electronic coherence with environmentally induced decoherence on the energy transfer time scale are especially relevant if they can survive and be measured in systems that mimic the conditions of natural light-harvesting. Most photosynthetic complexes fall in the so-called intermediate coupling regime, where neither the environment nor the inter-pigment couplings are strong enough to dominate EET. A feature of the Huang-Rhys factors for the photosynthetic pigments Chlorophyll (Chl) and Bacteriochlorophyll (BChl) is that they are remarkably small (*S* ≤ 0.05)^{50–54} compared to typical values in molecular dyes, where *S* ≥ 0.5 is common, and to model dimers where a clear distinction was made between vibronic and electronic coherences.^{44,45} Given the very small *S* values in natural systems, it is tempting to conclude that intramolecular vibrations are likely to play little role in either the spectroscopy or the quantum dynamics. However, as Jonas and co-workers pointed out,^{26} this may not be the case. Our goal in this paper is to clarify the role of vibrations in the small *S* regime in terms of their influence on two-dimensional spectra and on the underlying quantum dynamics. To this end, we carried out a systematic study of the influence of *S* ≤ 0.05 for a single mode in an electronically coupled dimer. We aim to establish the phenomenology with a model that is simple enough for exact numerical calculations, incorporating the vibrational degrees of freedom (DOFs) in the spectral density and not in the system Hamiltonian. This is important because the small Huang-Rhys factors make the vibronic resonance very fragile with respect to environmentally induced fluctuations, and so it is essential to treat this environmental influence to avoid overestimating the influence of near-resonant vibrations and to allow the comparison of dynamics at very different temperatures. Of course, in a real system, there may be many vibrational modes that play a role. Recently, Schulze and Kühn^{55} examined this issue using a zero temperature, numerically exact multilayer multiconfiguration time-dependent Hartree method on a dimer chosen from the FMO complex. Their calculation showed that a large number of weakly coupled modes participate in the dynamics, although the zero temperature assumption makes their role at physiological temperature difficult to predict.

In a recent paper,^{56} Ishizaki and coworkers examined the effect of environmentally induced fluctuations on the mixing between electronic and vibrational DOFs in a heterodimer inspired by BChls 3 and 4 in FMO. They included the influence of an underdamped vibration at the acceptor pigment site, with *S* = 0.025, and numerically simulated the stimulated emission, ground state bleach, and excited state absorption contributions to the 2D spectra, showing that oscillatory signals originate from dynamics on both the electronic ground and excited state. They showed that if the vibrational frequency is close to resonant with the electronic energy gap, oscillations persist in the spectra at 77 K–despite the small Huang-Rhys factor chosen. These disappear at physiological temperature, as increased fluctuations destroy the electronic-vibrational mixing. Moreover, they observed that large amplitude beating of vibrational origin emerges in the 2D spectra, but with no accompanying effect on EET dynamics. This raises the important question of the extent to which, for a given value of *S*, the vibration-electronic mixing in PPCs contributes to coherence between sites (correlation between the excitation state of separate pigments) or else generates only local effects such as a vibrational wavepacket. The distinction is important because only signatures of inter-site coherence indicate that coherent EET mechanisms are possible.

In this paper, we explore the dependence of the 2D spectra on the magnitude of the Huang-Rhys factor while holding the Stokes shift (from the environmental reorganization) constant, and we consider vibrational modes both resonant and off-resonant with the exciton energy gap. Our results are intended to be useful for interpretation of experiments on photosynthetic systems, and so we confine our calculations to the range of small *S* values typical of Chls and BChls in pigment-protein complexes. As Ishizaki and coworkers demonstrated,^{56} the observation of distinct beats in a 2D spectrum guarantees neither that electronic coherence is involved in the energy transfer nor that intramolecular vibrations significantly enhance the energy transfer rate. In order to explore this point further, we introduce a formalism for calculating measures of inter-site coherence and entanglement dynamics. This enables a comparison of the electronic or vibronic coherence in the site basis, as a function of the Huang-Rhys factor, with the features observed in a 2D electronic spectrum.

## II. MODEL

We simulate EET dynamics in a model dimer with coupling to a harmonic vibrational mode in the ground and excited states, under the influence of environmentally induced fluctuations with Gaussian statistics. Figure 1 shows a schematic of the excitons that result from electronic and vibronic mixing in this dimer. We use the following notation: for pigment *m*, |*ϕ _{ma}*〉, with

*a*=

*g*or

*a*=

*e*describes the ground and first excited electronic states, respectively. $| \chi a \nu m \u3009$ describes the nuclear configuration of pigment

*m*with vibrational occupation number

*ν*on electronic state

*a*. We consider only the ground or first excited electronic states of each pigment, although higher excited states are sometimes of consequence in nonlinear spectroscopic signals.

^{57,58}As illustrated in Fig. 1, the interaction between $| \varphi 1 e \u3009| \chi e 0 1 \u3009$ (the vibrational ground state on the donor pigment) and $| \varphi 2 e \u3009| \chi e 1 2 \u3009$ (a vibrationally excited state on the acceptor) creates the mixed-character vibronic excitons, $| e 1 + \u3009$ and $| e 1 \u2212 \u3009$. The details of the theoretical method used are described in Ref. 56. The Hamiltonian for the dimer is given by

where $ H \u02c6 m a ( x m ) $ is the diabatic Hamiltonian for the environmental and nuclear DOFs, *x _{m}*, of pigment site

*m*in state

*a*. The Franck-Condon transition energies are then obtained by an average over ground state environmental and nuclear DOFs,

where 〈…〉_{mg} denotes the canonical average with respect to the environmental and nuclear DOFs associated with the |*ϕ _{mg}*〉 state.

In this work, we have set the Franck-Condon transition energies to be Ω_{1} = 600 cm^{−1} and Ω_{2} = 400 cm^{−1}, with electronic coupling *J*_{12} = − 50 cm^{−1}, such that $\Delta = ( \Omega 1 \u2212 \Omega 2 ) 2 + 4 J 12 2 =223.6 cm \u2212 1 $ is the gap between the two electronic energy eigenstates. The electronic coupling, approximate energy gap, and the environment parameters (described below) are taken from values for BChls 3 and 4 in the FMO complex, which were generated by simultaneous fits of measured linear absorption and 2D rephasing, nonrephasing, and polarization-dependent electronic spectra.^{47} We describe the influence of the environmental and vibrational DOFs on the electronic transition of pigment *m* using the relaxation function,^{59} Ψ_{m}(*t*), of the collective energy gap coordinate defined by $ u \u02c6 m = H \u02c6 m e ( x m ) \u2212 H \u02c6 m g ( x m ) \u2212\u0127 \Omega m $. It should be noted that the relaxation function is independent of temperature by definition, when Gaussian fluctuations in electronic transition energies and harmonic vibrations (i.e., the linear response regime) are considered.^{56,59,60} In this work, we assume that the relaxation functions for the two pigments are identical, i.e., Ψ_{1}(*t*) = Ψ_{2}(*t*) ≡ Ψ(*t*). We take the relaxation function to have two components: an exponential decay due to the collective influence of overdamped environmental DOFs, with reorganization energy *λ*_{env}, and the influence of a single underdamped vibration modeled as a Brownian oscillator^{61} with vibrational frequency *ω*_{vib} and relaxation rate *γ*_{vib}. The total relaxation function, including both the environmental and vibrational contributions is then

where $ \omega \u0303 vib = ( \omega vib 2 \u2212 \gamma vib 2 ) 1 / 2 $, and *S* is the dimensionless Huang-Rhys factor for the underdamped vibration. Here, we set the environmental reorganization time scale to $ \gamma env \u2212 1 =50fs$ and the vibrational dephasing time to $ \gamma vib \u2212 1 =2ps$. This form of the relaxation function (Eq. (3)) yields a broad spectral density *J*(*ω*) with a narrow peak in the vicinity of *ω* = *ω*_{vib}, via a Kubo formula,^{59,62}

Here, we vary *S* and *ω*_{vib} but constrain *λ*_{tot} = *λ*_{env} + *Sω*_{vib} = 35 cm^{−1}, so that the Stokes shift Ψ(0) = 2*ħ λ*_{tot} is fixed in every scenario. We consider the effect of environment structure on the dynamics while holding approximately constant the dissipation associated with the response of all over- and underdamped bath modes, thus avoiding a systematic increase in the mean-squared fluctuations of the electronic transition energies. However, we note that in our approach, the effects of increasing *S* are not purely due to the added presence of the vibration, but also include the relative decrease in contributions from overdamped bath modes.

We compute the dynamics of observables via the reduced density matrix, for which the unobserved environmental DOFs have been traced out: *ρ*_{exc}(*t*) = Tr_{env}[*ρ*(*t*)]. For ease of notation, we will write |*g*〉 = |*ϕ*_{1g}〉|*ϕ*_{2g}〉 for the collective ground state, |*m*〉 = |*ϕ _{me}*〉|

*ϕ*〉 for the electronic state with pigment

_{ng}*m*in its first electronic state and pigment

*n*in its ground state, and |12〉 = |

*ϕ*

_{1e}〉|

*ϕ*

_{2e}〉 for the doubly excited state. In deriving the equation of motion, following Appendix A of Ref. 56, we begin from the response function and the symmetrized correlation function of the collective energy gap coordinate. The response function is simply expressed in terms of the relaxation function as Φ(

*t*) = − (

*d*/

*dt*)Ψ(

*t*). The symmetrized correlation function can be expressed in terms of the spectral density as

For the model described by Eq. (3), the integrand of Eq. (5) possesses four singularities $\omega =\xb1 \omega \u0303 vib \xb1i \gamma vib $, in addition to an infinite number of singularities *ω* = ± *ν _{k}* (

*k*= 1, 2, …), where

*ν*= 2

_{k}*πk*/

*βħ*are the Matsubara frequencies. Hence, Eq. (5) may be recast into the form of an infinite series. However, in order to save on the computational cost of calculating 2D electronic spectra for various waiting times, we approximate Eq. (5) as

where Ψ_{env}(*t*) denotes the overdamped environmental component of the relaxation function given in Eq. (3). In the first term of the right hand side of Eq. (6), we have assumed that the classical fluctuation-dissipation relation applies to the environmental dynamics by omitting the quantum components of the fluctuations, characterized by the Matsubara frequencies. The drawback of this assumption is that we deviate from the exact detailed balance condition and thus, a breakdown of the positivity condition for the population is possible. In this work, however, we focus on coherence dynamics under the influence of the environmentally induced fluctuations, so that this drawback is of no consequence. In the second term, on the other hand, the hyperbolic cotangent function indicates that the quantum zero-point energy of the nuclear vibration is properly taken into account. In deriving the expression, it has been assumed that the vibrational component of the spectral density does not resonate with the Matsubara frequencies. This assumption is reasonable as long as we consider slow vibrational dephasing times and correspondingly sharp peaks in the spectral density.^{56} Thus, both the response function and the symmetrized correlation function are expressed in terms of exponential functions, as required in Eq. (A8) of Ref. 56.

In the simplest picture, the electronic coherence oscillations appear on the cross-peak in the rephasing spectrum and on the diagonal peaks in the nonrephasing spectrum.^{3} Therefore, to better separate the oscillatory dynamics from population relaxation, we focus on the rephasing 2D spectra. Assuming infinitely broadband laser pulses, the rephasing 2D spectrum is given by

In Eq. (7), *R _{R}*(

*t*

_{3},

*t*

_{2},

*t*

_{1}) is the third-order responser function for the rephasing signal expressed as

^{63}

where *ρ*_{exc}(−∞) is the equilibrium reduced density operator for the electronic ground state. The operator $G ( t i ) $ describes the time evolution and relaxation during the intervals between pulses *t _{i}*. For simplicity, we assume parallel transition dipoles of unit magnitude for both sites. In this case, the operators that describe interaction with the identically polarized laser pulses are

*μ*

_{←}= |

*g*〉〈1| + |

*g*〉〈2| + |1〉〈12| + |2〉〈12| and

*μ*

_{→}= [

*μ*

_{←}]

^{†}. We have used the superoperator notation $ \mu \xd7 O \u02c6 = [ \mu , O \u02c6 ] $ to simplify the appearance of commutators in Eq. (8).

### A. Site-probe response

The excitation population on site *m* is given by

where the reduced density operator $ \rho exc ( 1 ) ( 0 ) $ describes the chosen initial condition for the electronic DOFs in the single-excitation manifold. For example, if the initial population is localized on the donor pigment, then $ \rho exc ( 1 ) ( 0 ) =|1\u3009\u30081|$. Ishizaki and coworkers demonstrated that for a dimer under the influence of environmentally induced fluctuations, a vibronic resonance can greatly impact the 2D-ES measurement without strongly affecting the relaxation to equilibrium of an excitation initially located on the donor pigment.^{56} They called attention to the high sensitivity of nonlinear measurements to coherent dynamics that may only play a small role in the overall energy transfer process.

In this work, we examine the correspondence between the spectroscopic measurement and the coherent aspects of dynamics in the site basis. To this end, we introduce the theoretical formalism for a site-specific nonlinear measurement, which we call the “site-probe response,”

where *m* = 1, 2 indicate projection on the donor and acceptor pigments, respectively. This quantity is analogous to the optical response in Eq. (8). Whereas site population dynamics represent energy transfer as a given site population relaxes to equilibrium, the site-probe response represents dynamics of correlation between pump-induced optical polarizations and the site population, just as the third-order optical response represents dynamics of correlation between pump and probe optical polarizations. In terms of the “doorway-window” picture of Fried and Mukamel,^{61,64} the site-probe response employs the same two-pulse “doorway function” as is used to calculate the 2D spectrum, $D= \mu \u2192 \xd7 G ( t 1 ) \mu \u2190 \xd7 $. After time propagation with $G ( t 2 ) $, the “window function” (or measurement) is given by $ W m =|m\u3009\u3008m|$, the projection onto the excited state localized at site *m*. This site projection replaces the third pulse interaction and resulting signal field ($W=\mu G ( t 3 ) \mu \u2190 \xd7 $ in Eq. (8)) and provides a site-specific measure of the system dynamics after the pump excitation. With a Fourier transform over the “doorway” coherence time *t*_{1}, we obtain the site-probe response as a function of absorption frequency *ω*_{1}.

Just as oscillations in the 2D spectrum show coherent energy flow between excited states, oscillations in the site-probe response show coherent flow between pigment sites. The site-probe oscillations depend on the presence of coherences between excitons, as they are a direct consequence of the unitary dynamics of coherent exciton superpositions. However, not all exciton coherences produce site-probe oscillations. For instance, in a coupled dimer with a small mixing angle, the excitons are localized on the sites. Hence, [|*m*〉〈*m*|, *H*] ≈ 0, and thus, the site projection measurement is time-independent except for dissipation, even if the excitons were prepared in superposition. In this scenario, the site-probe response is just proportional to the excited state population, and hence to the contribution to the spectrum from either stimulated emission or excited state absorption. However, as *J*_{12} and the mixing angle increase, dynamics in the site-probe response increasingly represent a mixture of exciton population and coherence dynamics. Due to this mixing, there is no simple correspondence between the phase of the site response oscillations and oscillations at any single point in a 2D spectrum; nonetheless, the presence of site-probe oscillations reports when coherent energy transfer occurs. The site projection measurement is blind to the vibrational DOFs, so that while a localized vibrational wavepacket may give rise to long-lived oscillations in the 2D-ES, the site-probe response is insensitive to purely vibrational effects.

## III. RESULTS AND DISCUSSION

### A. Simulated 2D-ES

We simulate the rephasing 2D electronic spectra, Eq. (8) and site-probe response, Eq. (10) across waiting time *t*_{2}, for an electronic heterodimer with coupling to a vibrational mode in the ground and excited electronic manifolds. Our simulations are split into two cases. In the first case, the vibrational frequency matches the pure electronic energy gap (*ω*_{vib} = Δ = 223.6 cm^{−1}). This leads to two closely spaced vibronic excitons $| e 1 + \u3009$ and $| e 1 \u2212 \u3009$ in addition to |*e*_{0}〉 (Fig. 1). In the second case, the vibration is off-resonant, *ω*_{vib} = Δ + 100 cm^{−1}. This study is limited to small Huang-Rhys factors, 0 ≤ *S* ≤ 0.05, which correspond to the range of values measured for individual vibrational modes in BChl and Chl in low-temperature glasses.^{53} In the FMO complex, several low-frequency vibrations^{51,52} have been identified with Huang-Rhys factors 0.008 ≤ *S* ≤ 0.015; the strongest low-frequency mode measured has *S* = 0.05.

Figure 2(a) shows the simulated linear absorption spectrum of the dimer for a range of *S*-values and with the constraint *λ*_{tot} = *λ*_{env} + *Sω*_{vib} = 35 cm^{−1}. Maintaining a constant Stokes shift (2*ħ λ*_{tot}) prevents a strong *S*-dependence in the linear absorption peak positions and widths and is a key to making a useful comparison of spectra across a range of vibronic couplings.

Figures 2(b) and 2(c) show the real part of the simulated rephasing 2D spectra at 77 K and *t*_{2} = 0 fs for *S* = 0 and *S* = 0.05. Measured over *t*_{2}, the diagonal peaks DP1 and DP2 show the population dynamics in states |*e*_{0}〉 and $| e 1 \xb1 \u3009$, respectively (see Fig. 1). Oscillations arising from electronic coherence are expected to appear in the lower cross-peak, CP21,^{3} and the growth of this cross-peak indicates the energy relaxation from upper to lower electronic excitons. While the shape of DP2 and CP21 changes with *S* in Figs. 2(b) and 2(c), the peak positions remain the same. Vibrational dynamics are expected to contribute to both the diagonal and the cross-peaks. Due to the rapid environmental reorganization (∼50 fs) and weak electronic coupling (*J*_{12} = 50 cm^{−1}), we expect minimal electronic coherence oscillations in the absence of electronic-vibrational coupling.

We calculated the waiting time dynamics of 2D spectra for resonant and nonresonant vibration scenarios at 77 K and 300 K. Figure 3 displays the amplitude of CP21 as a function of *t*_{2}, with *ω*_{1} and *ω*_{2} taken from the maxima along a diagonal cut through DP2 and DP1 (respectively) at *t*_{2} = 0 fs. Oscillations appear in addition to an exponential rise. For *ω*_{vib} = Δ at 77 K (Fig. 3(a)) and *S* = 0.01 (black line), the cross-peak oscillates by about 2% of the *t*_{2} = 0 fs spectrum maximum even after 1 ps; oscillations of about 1% persist for *t*_{2} > 1 ps when the Huang-Rhys factor is as small as *S* = 0.005 (green line). For comparison, in the case of *S* = 0.05 (red line), the oscillation remaining after 1 ps has an amplitude of about 5% of the *t*_{2} = 0 fs spectrum maximum.

Ishizaki and coworkers showed that oscillations such as these originate from dynamics in both the ground and excited electronic manifolds.^{56} They suggested that the short-term (*t*_{2} < 300 fs) oscillations are due to the electronic coherence, while the late-time oscillation is due to vibrational wavepackets, created via the vibronic states $| e 1 \xb1 \u3009$. It is the resonance between $| \varphi 1 g \u3009| \chi g 0 1 \u3009\u2194| \varphi 1 e \u3009| \chi e 0 1 \u3009$ and $| \varphi 2 g \u3009| \chi g 0 2 \u3009\u2194| \varphi 2 e \u3009| \chi e 1 2 \u3009$ that creates the $| e 1 \xb1 \u3009$ states. The transition dipoles for $|g\u3009\u2194| e 1 \xb1 \u3009$ are stronger than for $| \varphi 2 g \u3009| \chi g 0 2 \u3009\u2192| \varphi 2 e \u3009| \chi e 1 2 \u3009$, due to mixing with the pure electronic transition.^{65} Thus, the amplitudes of oscillations in Fig. 3 should depend on this resonance. Figure 3(b) shows the cross-peak dynamics when *ω*_{vib} = Δ + 100 cm^{−1}. The oscillatory components of the signals are indeed significantly weaker than in the on-resonance case shown in Fig. 3(a), with amplitudes only nearing 1% of the spectrum maximum for *S* = 0.05 at 77 K.

The simulations at 77 K indicate that depending on the experimental signal-to-noise ratio, we expect that the low-frequency intra-pigment vibrations may easily impact waiting time oscillations, despite very small Huang-Rhys factors. Moreover, the very strongest intramolecular vibrations in BChl (*S* = 0.05) may influence the 77 K 2D spectrum even if the frequency is relatively far from resonant with the PPC electronic energy gap. However, as was found for the case of *S* = 0.025 examined in Ref. 56, the contribution of electronic-vibrational mixing to waiting time oscillations in 2D-ES is extremely sensitive to temperature. Figures 3(c) and 3(d) show cross-peak dynamics simulated at 300 K. Noting the changes in the vertical scale, oscillations at 300 K remain long-lived but have much smaller amplitudes. In the case of *ω*_{vib} = Δ and *S* = 0.05, oscillations of about 0.5% of the *t*_{2} = 0 fs spectrum maximum persist after 1 ps at 300 K. This calculation supports the claim that owing to the small Huang-Rhys factors and the influence of the environment, oscillations due to intramolecular vibrational coherence are small and unlikely to be detected in weakly coupled Chl and BChl complexes at physiological temperature.

Kreisbeck *et al.*^{39} used a short-time Fourier transform to analyze 2D spectra simulated with a spectral density extracted from experiments studying FMO.^{52} They found that the oscillation frequency shifted from the electronic energy gap to vibrational frequencies over the waiting time and used this to distinguish between the different origins of oscillation. To further clarify our analysis, we also examine the differences between oscillations appearing at early (*t*_{2} = 0-400 fs) and late (*t*_{2} = 400-800 fs) waiting times. Figure 4 shows time-windowed Fourier transforms of cross-peak traces in Figs. 3(b) and 3(c), where *ω*_{vib} = Δ + 100 cm^{−1}. We have used a moving-average filter to partially remove the slow exponential rise and then transformed 400 fs segments of the waiting time oscillations. The resulting frequency-domain peak is especially broad due to the short *t*_{2} window, and the pure electronic oscillations (*S* = 0, dark blue line) are broad due to the rapid environmental reorganization ($ \gamma env \u2212 1 =50fs$). Nonetheless, the 100 cm^{−1} detuning of the vibration from the electronic energy gap is enough to distinguish the two frequencies in the *t*_{2} oscillations (red dashed lines). Late-time oscillations appear at the vibrational frequency for any *S* > 0 and at both 77 K (Fig. 4(b)) and 300 K (Fig. 4(d)); this supports ascribing them to vibrational wavepackets. The 300 K late-time oscillations are ∼70% weaker than at 77 K, which again demonstrates the reduced electronic-vibrational mixing in the presence of strong fluctuations. Early time oscillations at 300 K (Fig. 4(c)) also match the vibrational frequency. These physiological-temperature vibrational oscillations are likely too weak to be observed experimentally. Nonetheless, we can observe from the simulations that in this weakly coupled dimer (*J*_{12} = − 50 cm^{−1}), with rapid environmental reorganization, even the weakest of vibrations have more influence than electronic coherence at 300 K. It is only at 77 K (Fig. 4(a)) and *S* = 0 (dark blue line) that the frequency of the early time oscillation matches the electronic energy gap Δ, and the frequency tends strongly toward *ω*_{vib} as *S* increases. Given that in the FMO complex, many low-frequency vibrations have been measured with *S* ∼ 0.01; we take note that in Fig. 4(a), the simulation with *S* = 0.01 (black line) yielded a waiting time oscillation frequency midway between Δ and *ω*_{vib}. This suggests that for BChls 3 and 4 in FMO, the influence of the numerous low-frequency vibrations may be of a similar magnitude to the purely electronic oscillations observable at 77 K, even at early waiting times.

The purely vibrational effects in complexes with weak electronic coupling and small Huang-Rhys factors are extremely difficult to identify using 2D spectra alone. Although 2D Fourier maps of the waiting time oscillations have recently been used to assign a mixed character to coherences in the PSII reaction center^{9} (where the electronic coupling is large), in the case of weak electronic coupling and small Huang-Rhys factor, vibrational and electronic Fourier map patterns were not distinguishable.^{35} Polarization-controlled 2D-ES is a promising experimental solution and has already been applied to distinguish vibrational and electronic coherences in the bacterial reaction center.^{8} However, the polarization-control method is limited by its critical dependence on relative electronic transition dipole orientations.

### B. Comparison to site-probe response

When the vibrational and electronic character coherences are not distinguishable, the correspondence between oscillations in 2D spectra and inter-site coherence becomes unclear. A better measure of inter-site coherence would probe the correlation between absorptions and later excitation population at a specific pigment, as in Eq. (10). This “site-probe response” is intended as a useful benchmark for evaluating the correspondence between 2D-ES oscillations and underlying inter-site coherence in different model dimers. The site-probe response is insensitive to purely vibrational dynamics, because the nuclear DOFs are traced over at the end of the waiting time. The only DOF measured is the single excitation localized on pigment *m*, ensuring that the site-probe response manifests oscillations due to excited-state coherence between pigment sites. By comparing the oscillations in the simulated 2D spectra (Eq. (8)) with the site-probe response measurement, we can definitively distinguish between localized vibrational wavepackets and inter-site electronic or vibronic coherences in our simulations. In Fig. 5, we display the results of the site-probe response calculation for *ω*_{vib} = Δ and *ω*_{vib} = Δ + 100 cm^{−1} at 77 K and also at 300 K. The optical pump frequency *ω*_{1} was selected to match the DP2 feature in the 2D-ES, and the probe is at the donor site (*m* = 1 in Eq. (10)). Initially, the donor site-probe response becomes increasingly negative with *t*_{2}, because pumping $| e 1 \xb1 \u3009$ is positively correlated with population transfer from the donor to the acceptor (note that traces in each panel have been vertically offset for clarity). As was observed in the 2D-ES simulations, the amplitude of oscillations in the site-probe response increases with *S* and is also enhanced by the resonance, *ω*_{vib} = Δ. Comparing Figs. 3(a) and 5(a), we observe that for *S* ≤ 0.02, oscillation amplitudes are significantly more robust in the 2D spectrum than in the site-probe response. This same observation can be made from Figs. 3(b) and 5(b): while in the 2D-ES measurement, the amplitude of late-time oscillations increases steadily with increasing *S*, site-probe response oscillations are relatively insensitive to the vibration for *S* ≤ 0.02. Oscillations that, with increased *S*, appear in the 2D spectrum but not in the site-probe response must arise from vibrational wavepackets (both on the ground electronic state and localized on a single pigment in the excited state). The comparison between Figs. 3 and 5 confirms that in our simulations, the 2D-ES late-time oscillations at 77 K do not correspond to coherence between sites. Similarly for *ω*_{vib} = Δ + 100 cm^{−1}, Fig. 3(b) shows late-time oscillations that are nearly absent in Fig. 5(b), for all *S* ≤ 0.05. Detuning the vibrational frequency prevents sufficient mixing for long-lived vibronic delocalization, so that when *ω*_{vib} = Δ + 100 cm^{−1}, the late-time oscillations in the 2D spectra are vibrational in nature.

Examining the results for *S* > 0.02 in Fig. 5(a), we observe that increasing the vibronic coupling prolongs the lifetime of site-probe response oscillations, which implies a delay in irreversible decoherence and excitation localization. Since the Franck-Condon factor that weights the coupling between $| \varphi 1 e \u3009| \chi e 0 1 \u3009$ and $| \varphi 2 e \u3009| \chi e 1 2 \u3009$ scales as $ S $ for *S* ≪ 1, a larger *S* produces stronger vibronic mixing and more delocalized $| e 1 \xb1 \u3009$ states. The vibrational relaxation in this model is slow, so that coherences with more mixed character vibronic coherence persist for longer. Our simulation indicates that the Franck-Condon factor of 0.05 (*S* ≈ 0.056) measured for a 192 cm^{−1} mode in BChl^{53} may be large enough to induce long-lived vibronic delocalization effects in PPCs with a near-resonant electronic energy gap, albeit only at low temperature (see Fig. 5(c)).

For *S* = 0.05 (red line) in Fig. 5(a), the site-probe oscillation shows clear early- and late-time components, with close to a *π* difference in phase. We find that early- and late-time components do not fit well to a sum of exponential sinusoids. Instead, we used the form

where *σ*(*t*) = [1 + *e*^{−k(t−t0)}]^{−1} is a logistic function that determines the amplitude and arrival timing for the late-time component. A least-squares fit of the *S* = 0.05 site-probe response trace in Fig. 5(a) was obtained with coefficient of determination,^{66} *R*^{2} = 0.9988. The early oscillatory component with *ω* = 223.5 cm^{−1} decays with *τ*_{3} = 145 fs. The late-time oscillatory component turns on after *t*_{0} = 261 fs, with *ω*′ = 222.3 cm^{−1} and decay time *τ*_{4} = 665 fs. The phase difference *ϕ* − *ϕ*′ between the components is 190°. This form suggests that the site-probe response dynamics result from two distinct processes. The two oscillation decay time scales both fall between the environmental reorganization time ($ \gamma env \u2212 1 =50fs$) and the vibrational dephasing time ($ \gamma vib \u2212 1 =2ps$), which implies a mixed-character vibronic origin for both components.

If it seems in Fig. 5(b) that increasing *S* prolongs the lifetime of oscillations, when *ω*_{vib} = Δ + 100 cm^{−1}, we must interpret the simulation carefully. These longer lifetimes are due in part to the adjustment in the overdamped environmental reorganization energy, *λ* = 35 cm^{−1} − *Sω*_{vib}. While this adjustment preserves the total amount of dissipation despite an increase in *S*, in order to do so, it must decrease the spectral density at *ω*≠*ω*_{vib} (see Eqs. (3) and (4)). This will cause an increase with *S* in the lifetime of oscillations with frequencies far from *ω*_{vib}. Nonetheless, the difference between the results in Figs. 5(a) and 5(b) indicates that this trivial dependence cannot be the only influence of the Huang-Rhys factor on the site-probe response and dynamical localization.

Figure 6 shows the time-windowed Fourier transforms (compared to Fig. 4) of the site-probe response traces. As in Fig. 4, the left hand panels in Fig. 6 show the oscillation frequencies for the waiting time period *t*_{2} < 400 fs, while the right hand panels show frequencies for 400 fs < *t*_{2} < 800 fs. In Figs. 6(a) and 6(b), the site-probe response oscillates at the energy gap Δ, not at *ω*_{vib} as seen in the 2D spectrum (Figs. 4(a) and 4(b)). We attribute this site-probe oscillation to the evolution of the $| e 1 \xb1 \u3009\u3008 e 0 |$ electronic coherence. At 300 K (Figs. 6(b) and 6(c)), very small amplitude site-probe oscillations instead appear near *ω*_{vib}. These are explained by the minimal amount of inter-site coherence oscillation associated with the vibrational wavepacket, when there is a small amount of residual vibronic mixing. By comparing the frequency content of the site-probe response and the 2D-ES cross-peak oscillations, it is clear that at 77 K, the 2D-ES oscillations (for the relatively weak couplings considered here) dominantly originate from electronic coherence only for systems with *very* weakly coupled vibrations, *S* < 0.005.

In all of the simulations described here, the majority of the energy transfer is accomplished within the environmental reorganization time scale, $ \gamma env \u2212 1 =50fs$. In the range of small Huang-Rhys factors that we consider, the rate of EET as measured in the calculated responses does not change significantly with increased *S*. We see this both in the site-probe measurement (Fig. 5) and in the initial decays of DP2 in the rephasing spectrum (not shown). This observation agrees with the findings of Ishizaki and coworkers:^{56} when both electronic coupling and electron-vibrational coupling are small, although multidimensional spectroscopies remain very sensitive to coherences of all types (vibrational, vibronic, and electronic), these coherences have little influence on the overall population transfer dynamics from one site to another via relaxation and energy transfer.

Our strategy of calculating time-dependent correlations between optical pump excitations and subsequent density matrix elements is not limited to probing site populations. To more clearly visualize the effect of vibrations on site-site entanglement and exciton delocalization, we calculate a quantity which we call the “concurrence-probe response.” Concurrence (*C*) is a widely used measure of entanglement between a pair of qubits.^{67,68} For a coupled dimer in which excitation is completely delocalized over the two sites, *C* = 1; for a localized excitation, *C* = 0. After two interactions with the laser field, the reduced density operator *ρ*_{exc}(*t*) has an “X-structure,”^{69} so that the concurrence in the single excitation manifold conveniently reduces to^{70}

Hence, similarly to Eq. (10), we can construct

Here, the doorway function $D= \mu \u2192 \xd7 G ( t 1 ) \mu \u2190 \xd7 $ and time-evolution $G ( t 2 ) $ are identical to those in the 2D spectrum (Eq. (8)) and site-probe response (Eq. (10)) calculations, but the window function is now taken from Eq. (12). Figure 7 shows waiting time traces of *C*_{12}(*ω*_{1}, *t*_{2}), for *ω*_{vib} = Δ and *ω*_{vib} = Δ + 100 cm^{−1}, at 77 K and 300 K. As with the site-probe response traces in Fig. 5, the optical pump frequency *ω*_{1} was chosen from the DP2 feature in the 2D-ES. This concurrence-probe response represents the time-dependent correlation between optical pump interactions at the upper exciton frequency and the subsequent amount of entanglement between the two pigment sites.

The low-temperature concurrence-probe response in Fig. 6 exhibits dramatic extinctions and revivals of the concurrence. Similarly, “entanglement sudden death” and “rebirth” have been observed in the dynamical behavior of open quantum systems;^{71,72} this effect can arise as a consequence of non-Markovianity in the environment^{73,74} and was found to be characteristic of the intermediate-coupling regime.^{75} With increasing *S* in our simulations, the contribution to the environment from the non-Markovian vibration increases relative to that of the overdamped environmental modes; the amplitude of revivals in Fig. 7 correspondingly increases. The effect of the vibration on the concurrence-probe response is somewhat more robust when *ω*_{vib} = Δ (Fig. 7(a)) than in the off-resonant case (Fig. 7(b)). At 300 K (see Figs. 7(c) and 7(d)), the entanglement revivals are small and much less dependent on *S*. We find these results in agreement with our analysis of the site-probe response: at 77 K, the *S*-weighted coupling and near-resonance between $| \varphi 1 e \u3009| \chi e 0 1 \u3009$ and $| \varphi 2 e \u3009| \chi e 1 2 \u3009$ create states that are more resistant to environment—induced localization. When electronic excitation energy localizes, information about the wavefunction at the other site has been irreversibly lost. The concurrence-probe oscillations demonstrate that when the vibrational and electronic DOFs are correlated through mixing, some of that information passes reversibly to the underdamped vibration and then returns to the electronic system. At 300 K, the environmentally induced fluctuations are powerful enough to destroy the mixing and therefore prevent strong revivals from appearing in the concurrence-probe response.

## IV. CONCLUSION

In order to better understand the effect of very weakly coupled vibrations on 2D electronic spectra and the coherence dynamics they measure, we have simulated 2D spectra of a model heterodimer. Our model incorporates environmental effects and an intra-pigment vibrational mode at each site via a phenomenological form of the relaxation function for the nonequilibrium excitation energy. Crucially, this approach correctly describes the destruction of vibrational-electronic mixing by environmentally induced fluctuations. We have systematically studied the impact of vibrations with very small Huang-Rhys factors, which are typical in photosynthetic pigment-protein complexes. We introduced the “site-probe response,” a formalism for effectively comparing simulated 2D-ES waiting time dynamics with the underlying dynamics at specific sites. We show that even the very weak electronic-vibrational interactions common in pigment-protein complexes can significantly affect the oscillations in a multidimensional spectrum so much that localized vibrational wavepackets obscure the underlying inter-site coherence dynamics at 77 K. Our simulations also demonstrate that at 77 K and for a Huang-Rhys factor comparable to the strongest low-frequency vibration in BChl, a vibration 100 cm^{−1} away from resonance with the electronic energy gap produces only localized vibrational coherence oscillations in the 2D spectra at late waiting time but that a near-resonant vibration significantly increases the inter-site coherence and delocalization remaining after 600 fs. However, by comparing simulated 2D spectra with the site-specific response calculation, we find that evidence of this vibronic resonance effect is obscured by the purely vibrational oscillations and is not discernible from the 2D spectrum alone. We confirm that, due to the fragility of vibronic mixing to the environmentally induced fluctuations, vibronic resonance effects vanish at 300 K.

The contrast between the 2D-ES simulations and the fictional site-probe response measurement simulated here highlights the potential impact of future experiments that access site-specific information, as opposed to the response from delocalized excitonic states, on our understanding of photosynthetic light-harvesting. Like the site-probe response discussed here, such experiments might measure the time-dependent correlation between excitons generated via Franck-Condon transitions and the energy subsequently available at a specific pigment. A spatially inhomogeneous observable, such as localized higher electronic states, transition dipole orientations, or intrapigment vibrations, plays the role of the site basis projection. The 2D electronic-vibrational spectroscopy^{76–78} recently developed in our laboratory, for example, is able to measure the site-probe response — so long as the probed vibration is confined on a single pigment. Such measurements, sensitive to changes in the location of excitation energy, enable a view of energy transfer not just as equilibration between eigenstates but as a flow across nanoscale distance. Besides insensitivity to vibrational wavepackets, such measurements might provide empirical evidence of how local variations in the spectral density^{79} affect EET mechanisms. While 2D electronic spectra provide an important measure of energy dynamics in pigment-protein complexes and other molecular aggregates, it should not be neglected that very weakly coupled low-frequency vibrations can easily mitigate the correspondence between long-lived oscillations in 2D spectra and inter-site coherence.

## Acknowledgments

This work was supported by the Director, Office of Science, Office of Basic Energy Sciences, of the U.S. Department of Energy, under Contract No. DE-AC02-05CH11231; the Division of Chemical Sciences, Geo-sciences, and Biosciences Division, Office of Basic Energy Sciences, through Grant No. DE-AC03-76SF000098 (at Lawrence Berkeley National Laboratory and University of California, Berkeley); and Grants-in-Aid for Scientific Research (Grant No. 25708003) from the Japan Society for the Promotion of Science. This research used resources of the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231. D.M.M. received a National Science Foundation Graduate Research Fellowship under Grant No. DGE-1106400.