Signatures of three-exciton correlations in the coherent and incoherent nonlinear optical response of photosynthetic complexes

We present a simulation study of population and coherence relaxation dynamics in the third-order nonlinear responses of the Fenna–Matthews–Olson (FMO) light-harvesting photosynthetic complex. Three-exciton correlations have been incorporated and the relevant six-point relaxation superoperator is calculated to second order in coupling to a phonon bath. The dynamics of three-exciton correlations induce new peaks that are absent or much weaker when the three-exciton variables are factorized. Such features were observed in a recent experiment (Engel G S et al 2007 Nature 446 782–6).


Introduction
The extremely efficient solar-chemical energy conversion in photosynthetic systems [1]- [3] has triggered intensive activity aimed at new energy solutions [4]- [19]. Exciton couplings ensure the efficient transfer of energy from light-capturing antenna pigments to the reaction center where charge separation occurs [1]- [3], [14]. The exciton coupling pattern in these systems is critical to understanding the mechanism of solar to chemical energy conversion.
Third-order nonlinear optical techniques [6]- [8], [20] are most suitable for probing energy transfer [13,14]. However, although tremendous progress has been made in recent years, there is still a long way to go in order to achieve a clear physical picture of the energy and coherence transfers in the photosynthetic process. The nonlinear optical responses are usually modeled either by analytical methods based on response functions or by numerically solving equations of motion (also referred to as a quasi-particle approach) [26]. The response-function method is more suitable for cases where coherent dynamics are dominant. For example, in the photonecho technique, there are three response functions in the coherent limit with rotating wave approximation. The key advantage of the response-function-based approach is that it provides an intuitive picture of the third-order nonlinear response in the coherent limit.
When incoherent dynamics, such as coherence and energy transfer, are included, the nonlinear responses become more complicated and cannot be one-to-one mapped with the three response functions. In this situation, employing response functions to understand these responses is only qualitatively valid. The alternative type of quasi-particle approach is thus more desirable to handle this complicated situation. This approach involves solving the equations of motion numerically [27]- [29] with finite pulses, or quasi-analytically in the impulsive limit [26]. However, when solving the equations of motion, such as the nonlinear exciton equations (NEE) [21], with this approach, people usually have to make approximations to keep 3 calculations tractable. For example, in the NEE formalism, a local-basis three-exciton correlation variable is explicitly factorized into lower-order correlations. In other cases, the variable is formally included but the relevant relaxation operator is constructed from low-order relaxation matrices. It can be shown that three-exciton correlations are equally missing when either factorizing the three-exciton variable or using a simplified or constructed relaxation operator.
In this work, we employ a numerical approach to solving the NEE without factorizing the three-exciton variable. Furthermore, we follow [21] and calculate the three-exciton relaxation operators using Markovian second-order perturbation theory with respect to system-bath coupling, instead of using simplified or constructed relaxation operators. Our simulations show that the three-particle correlations play critical roles in the relaxation dynamics of Fenna-Mathews-Olson (FMO) [22]- [25]. Some new peaks arise from the inclusion of the threeexciton correlations and the relevant relaxations. These qualitatively reproduce some features in recent experiments [14]. Because three-exciton correlations are always present in thirdorder nonlinear response, we expect that the full incorporation of three-point correlations also significantly modifies the nonlinear responses of other condensed matter systems, particularly when system excitons are strongly coupled.

Hamiltonian
The Frenkel exciton Hamiltonian [21], which describes the third-order optical response of aggregates, is given by the sum of the electronic part H e , the phonon HamiltonianĤ ph , the exciton/phonon couplingĤ int and the system-field interaction: Here,B n (B † n ) are exciton annihilation (creation) operators on site n, q α and p α are respectively nuclear coordinates and momenta, h mn represents the one-exciton Hamiltonian, U mn,kl describes the exciton-exciton scattering potential, m α and ω α are respectively the reduced mass and frequency of the α phonon mode, andh mn,α andŪ mn,kl;α denote exciton/phonon interactions originating from the q α -dependence of h mn and U mn,kl , respectively. We shall represent the optical electric field in equation (1) as Here, E + α (E − α = (E + α ) * ) is the envelope of the positive (negative)-frequency component of the αth pulse (α = 1, 2, 3, 4) centered at time instant t α , with carrier frequency ω α , polarization unit vector e α and wavevector k α .
The commutation relations for the exciton operators are where P mk,nl is a tetradic matrix describing the deviation of the exciton commutation relations from those of bosons (P mk,nl = 0). It follows from the Jacobi identity for the commutators that P mn,kl = P nm,kl = P mn,lk .
For two-level chromophores (hard-core bosons), we have Equation (6) then becomes For bosons (e.g. vibrations) and for Wannier excitons in semiconductors P mn, pq = δ m 1 q 1 δ p 1 n 1 δ m 2 p 2 δ n 2 q 2 + δ m 2 q 2 δ p 2 n 2 δ m 1 p 1 δ n 1 q 1 , where m = (m 1 , m 2 ), n = (n 1 , n 2 ), p = ( p 1 , p 2 ) and q = (q 1 , q 2 ) denote electron-hole pairs. Finally, the polarization that is directly related to the signal in nonlinear response is given byP where µ (1) n is the transition dipole moment vector of the chromophore located at r n and µ (2) n,kl is the transition dipole moment correction for the creation of two excitons. For hard-core bosons, we usually accept µ (2) n,kl = 0 and thus the second term in equation (12) is irrelevant.

The nonlinear exciton equations (NEE)
Using the Hamiltonian equation (1) and the Heisenberg equation of motion, we obtain the NEE for the following exciton variables relevant to third-order optical response [21]: N mn ≡ B † nB m , Z mn, j ≡ B † jB mBn , where B n and Y mn describe respectively single-and two-exciton dynamics, N mn contains exciton populations and Z mn, j describes three-exciton correlations. The equations of motion for these variables are given by where B,Ẏ ,Ṅ andŻ in equations (14)-(17) denote relaxation terms induced by coupling with phonons. We assume that phonons induce simple time-local (Markovian) dephasing and transport. We then havė where mn , R Y mn,kl , R N i j,mn and R m n , j mn, j are the local-basis relaxation matrices for the NEE variables B m , Y mn , N mn and Z mn, j , respectively. Here, the most complicated and lengthy calculation arises from the relaxation matrices R m n , j mn, j . It contains 16 terms within our approximation. A sample calculation for one of the 16 terms is given in appendix B. All other relaxation matrices can be obtained in a similar way but with far fewer calculations. These can be found in [26].  [20,26] Computing the Z -variable dynamics is the most expensive task in solving the NEE. By factorizing the Z variable in different ways, we obtain different approximate NEE with different numbers of NEE variables. A major motive for various factorization schemes is to simplify the calculation at the price of losing some types of exciton correlations. We shall compare five levels of approximations. In the lowest level, mean-field factorization (i), the three-point correlation in equation (14) is factorized as (27) and the Y and N variables are similarly factorized into products of single-exciton quantities. All correlations among excitons are neglected in this level and the equations for the B variables, equation (14), are now closed.

Factorization schemes for NEE variables
In the second, exciton coherent dynamics limit (ii), the three-point correlation is factorized according to and thus only the B and Y variables are relevant. In this limit, the correlation among two excitons is taken into account but population transport dynamics is missed. The simplest factorization that incorporates population transport is (iii) where only B and N variables survive but two-exciton correlations, Y , are ignored. A more general factorization that includes mean-field interaction, two-exciton correlation and population dynamics is given by (iv) In this case, we retain only B, Y and N variables and eliminate only Z . The fifth-type NEE (v) fully incorporates the Z variable without additional factorizations. We had also considered an alternative approach that keeps the Z variable but uses a simplified relaxation matrix. However, it can be shown (see appendix A) that using a simplified relaxation matrix is equivalent to factorizing the Z variable. Thus exciton correlations are still neglected even if we retain the Z variables in NEE.

Solving the NEE with three-exciton correlations
To retain the three-point correlation, we next calculate the three-point relaxation operatoṙ Z microscopically. We calculate other relaxation operatorsḂ,Ẏ ,Ṅ in equations (23)- (26) according to the method described in [21] and summarized in [26] 7 The equation of motion for the phonon-assisted variables Z where s can be either bath coordinate q or momentum p.
where m α and ω α are the reduced masses and frequencies of nuclear modes. β is defined as 1/(k B T ). N (s) m j;α and Y (s) mn;α are respectively the phonon-assist variables for other NEE variables N nm and Y nm . Z (q) mn, j;α in equation (32) was calculated by employing a Green's function technique (see appendix B). Substituting Z (q) mn, j;α into equation (31) and employing the Markovian approximation, we recast equation (31) aṡ where the relaxation matrix R m n , j mn, j is a sum of 16 terms whenŪ mn,kl;α are neglected. The detailed procedure for calculating R m n , j mn, j is given in appendix B.

Two-dimensional (2D) spectroscopy
In this section we demonstrate the signatures of three-point correlations in 2D spectroscopy [26]. It is a coherent four-pulse phase-locked time-domain approach, where two of the three time delays t 1 , t 2 and t 3 between the three excitation pulses k 1 , k 2 , k 3 and the detection heterodyne pulse k 4 are varied to give 2D signals. Two techniques, S I and S II (detected respectively along −k 1 + k 2 + k 3 and k 1 − k 2 + k 3 directions), involve population dynamics and thus are suitable for probing population transfers. We focus on the S I (rephasing) signal represented by three Feynman diagrams, as shown in figure 1. These pathways are ground-state bleaching (GSB), excited-state emission (ESE) and excited-state absorption (ESA), where e i and f i (i = I, II, . . .) represent respectively single-and two-exciton manifolds and g represents the ground state. The vertical arrows indicate coherence transfer during t 1 or t 3 , and population/coherence transfer or coupling during t 2 . These diagrams reduce to the coherent-limit diagrams (i.e. no vertical arrows) when there is only free propagation during each time delay. We now examine the relation between the three pathways of S I to the NEE formalism. For hard-core bosons, we have µ (2) m,kl = 0 and therefore, according to equation (12), the time-domain S I signals can be written as where B (3) m is the B m variable calculated to third order in the external field. The signal S I (t 1 , t 2 , t 3 ) is dominated by the three-point correlation function Z . Equation (14) does not depend on the two-exciton variable Y . We have numerically shown that the dependence on the N variable is weak (it enters only via interaction with the field due to Pauli blocking). For soft-core bosons, the signal only depends on the Z variable. The dependence of the signal on Y and N variables is only through the Z variable, as in equation (17). The dynamics of the Z variable therefore contain all three pathways.
We shall show in the next sections that the three-point correlations strongly affect the relative weights of three pathways. Our simulations show that the conventional ESE-type population transfer pathway is overestimated when the three-point correlations are factorized. However, a new ESA-type population transfer pathway can be identified. Different treatments of three-point correlations assign different weights to the three pathways. These pathways interfere and completely or partially cancel each other. The absence of a certain type of population/coherence transfer signal (e.g. the ESE signal below the diagonal line) does not necessarily indicate their real absence. The physical processes are indeed there but the signals might have been canceled by signals with opposite signs from other pathways (e.g. the ESA pathway).

Comparison with experiment
We have calculated the 2D signals by solving the NEE perturbatively in the optical field. Alternative numerical nonperturbative schemes are possible [27]- [29]. In this work, we calculate the 2D signals by solving the equations of motion for four NEE variables and by directly selecting the relevant spatial Fourier components [30]. The method employed is applicable to any pulse sequence, whether well separated or overlapping, and includes the pulse shapes. An approximate numerical orientational averaging scheme was employed to significantly reduce the calculation effort. Inhomogeneous broadening is added by introducing disorder to the diagonal elements of the Hamiltonian. The Hamiltonian parameters were taken from previous publications [31,32]. The Gaussian optical pulses are centered at 12 362 cm −1 with a 16 fs FWHM. -Ω 1 /1000 (1/cm) Comparison among experiment (top row), full theory (without factorization of three-exciton correlation, 2nd row) and two different lowerorder theories (3rd and 4th rows). The spectra in the 3rd row are calculated by factorizing the Z -variable Green's function into single-exciton and two-exciton Green's functions [26]. The spectra in the 4th row are calculated according to the factorization scheme in equation (30). Only the full theory (2nd row) reproduces the strong ESA-type blue peaks appearing in the experiments.
For comparison, the 2D signals [14] (S I + S II ) measured in recent experiments are shown in panels (a)-(d) of figure 2. The calculated 2D spectra are shown in panels (e)-(h). Each panel is first normalized to the strongest peak and displayed from −1 (blue) to 1 (red), and then is shown in an arcsinh scale [14]. We are most interested in the cross peaks in each 2D spectrum. One set of cross peaks is above the diagonal line (mainly negative) and the other set is below the diagonal line (mainly positive). The striking feature is the appearance of strong negative cross peaks (A)-(C) in panels (e)-(h), which qualitatively agree with the experiments in panels (a)-(d). These negative peaks are actually the dominant cross peaks in almost all the experimental results presented. The strong inhomogeneous feature in the experiments is not reproduced. We attribute this to the simple model of diagonal disorder used here.
These strong negative peaks only appear when the three-exciton correlation is fully included in the simulations and vanish or become much weaker otherwise. Panels (i)-(l) are calculated by employing the nonequilibrium Green's function where the Z -variable Green's function is factorized into single-exciton and two-exciton Green's functions [26]. -Ω 1 /1000 (1/cm) In figure 3, we compare the simulated S I , S II signals and their sum, S I + S II . Compared to the total signals shown in panels (a)-(d), the S I signals shown in panels (e)-(h) still have dominant ESA peaks arising from the three-exciton correlation. Some negative ESA peaks are even stronger than those in the total S I + S II signals. For S II signals shown in panels (i)-(l), we find that, in S II signals, the positive cross peaks below the diagonal line are stronger than their counterparts in S I signals and thus can clearly show the energy transfer.
To summarize, we showed that the three-point correlations result in strong negative cross peaks. The positive peaks (ESE pathway) below the diagonal line are commonly used to reveal the relaxation dynamics of FMO. The negative ESA peaks can provide further complementary information to the positive ESE peaks. They are related to population transfer and population/coherence coupling during t 2 , and to coherence transfer during t 3 .

Signatures of three-point correlations
To reveal the origin of the negative cross peaks, we now calculate S I signals by incorporating different types of elements of the exciton-basis relaxation matrices, R N e 4 e 3 ,e 2 e 1 and R e 3 e 2 ,e 1 e 6 e 5 ,e 4 , -Ω 1 /1000 (1/cm)  respectively, for the NEE variables N and Z . Figure 4 shows the calculated 2D signals at the same delay time t 2 = 600 fs. Panel (a) shows the spectrum with full R N e 4 e 3 ,e 2 e 1 and R e 3 e 2 ,e 1 e 6 e 5 ,e 4 . Panel (b) shows the spectrum where we eliminate only the population transfer elements in R N e 4 e 3 ,e 2 e 1 . Comparing panels (a) and (b), we find that the strengths of both peaks (B) and (A) are considerably changed. This indicates that the ESA peaks are related to the population transfer during t 2 . This feature can be understood from the Feynman diagrams in figure 1. Only the ESEtype peaks, such as peak (D) in all panels, were previously investigated for population transfer. The origin of this type of positive peak can be understood by the population transfer during t 2 in the ESE diagram in figure 1. However, the same population transfer also occurs during t 2 in the ESA diagram of figure 1. There is no particular mechanism to block the appearance of While the ESE-type peak (C) and ESA-type peak (D) are well resolved in panels (a) and (b) due to the artificial introduction of a big U mn,kl , the corresponding peaks (C) and (D) in panels (c) and (d) largely cancel each other due to a zero U mn,kl . Therefore, the strong negative ESA peaks may cancel the positive ESEtype population transfer peaks, which leads to an underestimation of population transfer. Panels (e)-(h) repeat the calculations in (a)-(d) except with a lowerorder theory, where the Z -variable Green's function is factorized into the product of single-and two-exciton correlations [26]. Here we do not observe the same cancelation as described for panels (a)-(d).
such ESA-type peaks (negative) representing the population transfer. Therefore, even from a simple diagrammatic analysis, we expect that these peaks above the diagonal lines are relevant to population transfers. However, the neglect of three-point correlations significantly reduces their strength.
Panel (c) is obtained by neglecting the coherence transfer elements in R e 3 e 2 ,e 1 e 6 e 5 ,e 4 . The spectrum is again strongly modified as compared to the full calculation in panel (a). This suggests that the coherence transfer in three-point correlation relaxation also plays an important role in the dominant negative peaks above the diagonal line. Comparing with the effect of coherence transfer during t 3 , we further find that the largest modification to the full calculation in panel (a) comes from the secular approximation in R N e 4 e 3 ,e 2 e 1 , where peak (B) becomes significantly stronger, as shown in panel (f).This suggests that the population-coherence couplings neglected in secular approximation also play an important role in determining the strength of ESA peaks.
Within the secular approximation, if we further eliminate the coherence transfer elements in R N e 4 e 3 ,e 2 e 1 , we obtain panel (e), where the strong negative peaks do not change a lot but the ESE-type population transfer peaks (D)-(F) become clearly resolved. This implies that the coherence transfer during t 2 , as included in panel (f), plays a role in reducing the ESE population transfer peaks. If we only eliminate the population transfer within the secular approximation for R N e 4 e 3 ,e 2 e 1 , then we obtain panel (d), where the conventional population transfer peaks, such as peak (D), are almost gone. This is different from panel (b), where we still have a weak peak (D) due to the inclusion of population and coherence couplings beyond secular approximation. However, in both panels (b) and (d), the negative ESA peaks exist even when we neglect the population transfer in R N e 4 e 3 ,e 2 e 1 . This indicates that the dominant ESA peaks depend not only on population transfer during t 2 , but also on other factors, such as coherences during t 2 and coherence transfers during t 3 . Therefore, like the conventional ESE-type peaks where population transfer and coherence are always coupled during t 2 , ESA peaks also have this property. The dependence of ESA peaks on coherence during t 2 can also be understood intuitively from the three Feynman diagrams in figure 1. From the diagrams, we know that the third pulse generates ge i (GSB and ESE) and f i e i (ESA) coherences. The generation of these coherences will rely on initial conditions, i.e. the coherence among different states and the population distribution at the end of t 2 . From this perspective, the t 2 -dependence of the ESA cross peaks can also be understood, although these cross peaks are mainly introduced by the three-exciton correlation that is finite only when t 3 > 0.
Here we have studied the dependence of ESA negative peaks on the various mechanisms given by the relaxation matrices of NEE variables N and Z . The strength of these peaks strongly depends on the population/coherence transfers and population-coherence couplings during t 2 , and on the coherence transfer during t 3 . We emphasize that the underlying reason for the existence of these dominant ESA cross peaks is the three-exciton correlation. If we employ a simplified relaxation matrix for the Z variable, we can still have coherence transfer during t 3 . But we will lose the dominant ESA peaks even with R N e 4 e 3 ,e 2 e 1 untouched. In our simulations, we also find that, with simplified R e 3 e 2 ,e 1 e 6 e 5 ,e 4 , the 2D signals will not be so sensitive to the different treatments of R N e 4 e 3 ,e 2 e 1 .

Discussion
We have shown that different approximations to the three-exciton correlations can affect the relative weights of three pathways of the SI technique. In FMO, the three-point correlations give the dominant ESA-type cross peaks. However, the fact that one type of cross peak can dominate has more profound consequences. For coherent dynamics, the ESE and ESA pathways introduce respectively positive and negative cross peaks both above and below the diagonal lines. These cross peaks with different signs can thus cancel. Therefore, weak cross peaks do not necessarily 14 indicate weak coherence in the system. For incoherent dynamics, the population transfer from the ESE pathway introduces positive cross peaks only below the diagonal line due to the general transfer direction from higher excitonic states to the lower ones. The population transfer from the ESA pathway mainly introduces negative peaks above the diagonal line. However, it is also possible that ESA population transfer peaks occur below the diagonal line due to the many two-exciton states available. From Redfield theory we know that the ESE coherence is coupled with the ESE population transfer, making it more complicated to probe population transfer in photosynthetic systems. However, when the negative ESA-type (both coherence and population transfer) peaks become dominant, the positive population transfer signals from ESE would also be canceled or partially canceled by the ESA peaks. This will further complicate the probing of energy transfer. This also indicates that the absence of population transfer peaks for small t 2 does not necessarily mean the absence of population transfer. Thus, if the coherence time is long, then the population transfer peaks below the diagonal line are more likely to be weak. We shall demonstrate this with a simple dimer model.
Most simulations of FMO usually only consider the J -coupling (h mn ) and neglect the exciton-exciton coupling U mn,kl . In the following dimer model, we introduce a finite U mn,kl in order to separate the ESE and ESA pathways. Panels (a) and (b) in figure 5 are calculated with U mn,kl = 100 cm −1 respectively at t 2 = 0 fs and t 2 = 300 fs. The U mn,kl renormalizes the two-exciton state and thus the ESA contribution to cross peaks is shifted down, as given by peaks (A) and (C). ESE contributions are still at the regular cross peak positions (B) and (D) with positive signs (note that the GSB contribution overlaps with the ESE contribution but GSB peaks are not t 2 dependent). By setting U mn,kl = 0, we obtain panels (c) and (d). In (c), peaks (C) and (D) partially cancel. We expect this cancelation to become more evident if inhomogeneous broadening is added. However, after the cancelation, the resulting weak cross peak does not indicate the weak coherence in the system. The coherence is strong but the coherence signals from two pathways partially cancel each other and thus look weak. At finite t 2 , population transfer already occurs, as in panel (d). However, the similar cancelation will let us underestimate the population transfer. We repeat the above calculations with lower-order theory, as shown in panels (e)-(h). We find that in the lower-order calculation the ESA pathway is no longer dominant, and the cancelation effect is not obvious. If the dominant ESA peaks are more realistic in accounting for the dimer dynamics, then the lower-order calculation will overestimate the population transfer.
In summary, we demonstrated that three-point correlation effects can affect the relative strength of individual pathways of nonlinear responses of the photosynthetic FMO complex. Our simulations demonstrate general excitonic properties as long as exciton correlations are important. The importance of three-exciton correlation also arises from the type of experiments considered (third-order, weak-field limit), where three-exciton correlation is the highest correlation involved, and is the only correlation that is directly related to the signal. Our simulations show that the probing of relaxation dynamics is actually more complicated. More advanced techniques, such as those based on manipulating polarization configuration and pulseshaping techniques, may resolve these issues. We further showed that the negative peaks above the diagonal line contain more information regarding coherence, energy/coherence transfer, than the positive peaks below the diagonal line. This is because the latter energy transfer peaks may have been partially canceled by the ESA peaks. whereV mn,kl;α ≡ 2Ū mn,kl;α − 2 p P mn, pkh pl,α − 2 pq P mn, pqŪ pq,kl;α . (B. 3) The equation of motion for the phonon-assist variables Z where s can be either q or p, and we have used the notation By defining Z (q) kl, j;α = A (q) α , equation (B.4) can be written out in a unified way: where W (s) α is the rhs of the corresponding equation and h are linear operators representing the free evolution of the Z variable in the absence of relaxation. We further introduce the set of the Green functions of the Z operator, which satisfy the equations The solution of equation (B.8) can be represented in the form Using the general formal solution equation (B.10), we have We also know, by definition, that  where We kl, jα (τ ) 20 We next employ a simplified bath model where a single-exciton eigenstate is coupled to its own statistically independent bath with a continuous spectral density. Fluctuations of different chromophores are uncorrelated and they all have the same spectral density so that C mnkl (ω) = δ mn δ nk δ klC (ω), (B.17) whereC (ω) = π αd 2 mm,α [δ(ω − ω α ) − δ(ω + ω α )] is now independent of the system indices. The spectral density has the symmetryC (ω) = −C (−ω). Specific models for the spectral density (e.g. Ohmic, white noise, Lorentzian, etc) are often used for calculating the nonlinear responses of molecular systems. We assume that the bath has several collective modes α and recast the spectral density in the form where C α (ω) is a dimensionless spectral density of mode α and λ α is the coupling strength between. We define the auxiliary functions  In order to calculate this relaxation matrix, we need first to transform it into the exciton basis, This is the first term in the relaxation matrix. The other 15 terms are calculated in a similar way. The total relaxation matrix for Z is finally given by