Nonperturbative modeling of fifth-order coherent multidimensional spectroscopy in light harvesting antennas

Recent advances in coherent multidimensional spectroscopy have boosted interest in exciton coherences in light harvesting complexes. We present nonperturbative calculations of two-dimensional (2D) electronic spectroscopy from a fifth-order phase-matching direction. The calculations show clear patterns that correspond to the electronic structure of one- and two-exciton manifolds of a Fenna–Matthews–Olsson light harvesting complex. Such signals can provide new information about the coherent properties of antenna pigment protein complexes.


3
However, as pointed out by Brüggemann et al [35], in the case of multilevel systems, there are fifth-order Liouville pathways that give signal in the same phase-matching directions as the conventional 2D third-order photon-echo. In multichromophoric systems, those pathways pass double-excited manifold of states and some of the pathways can be related to excitation annihilation [36]. Of course, the fifth-order signals have their own phase-matching directions where no third-order signal is present. In this work, we take a closer look at one such fifth-order signal in the Fenna-Matthews-Olsson (FMO) antenna pigment protein complex. The complex contains seven bacteriochlorophyll a (BChl) molecules.

Setup and signal directions
We use the method developed by Brüggemann et al [35] where the calculations of the 2D signal follow the same logic as an experiment. In an experiment, three beams of fs laser pulses in triangular configuration are focused to the sample plane as shown in figure 1(A). We point out that in real experiments the angles between the beams are just a few degrees and most of the 2D experiments use a so-called boxcar geometry where beams form a square. The phase-matched signal direction −2k 1 + 2k 2 + k 3 , investigated in this work, is highlighted in red. The waveforms of the three pulses that have passed the sample are drawn in green color. The time positions of the pulse maxima and a certain time point t of the signal are shown by blue triangles. In a real experiment, a fourth pulse is used as a local oscillator for recovering the phase information of the field [37]. This pulse propagation coincides with the signal phase-matching direction and comes usually before the other pulses. The time in the figure runs from right to left. Delay between pulses 1 and 2 is called coherence time τ . τ > 0 for pulse ordering 1-2 and τ < 0 for 2-1. Time delay T is called population time and is defined as a delay between 2 or 1, whichever is last of the two, and 3. T > 0 in our calculations, which means that pulse 3 is always last. This, however, is not a principal restriction.
There are many possible phase-matching k-vector combinations that correspond to nonlinear signals. The even-numbered combinations usually do not give signal because of symmetry considerations, and if they do, then the signal is in a very different spectral range compared to the incoming laser beam. An example of such signal is second harmonic generation. The odd-numbered combinations, however, lead to photon-echo signals, which can be fairly strong. Figure 1(B) shows a photograph taken in the laboratory during a fourwave-mixing experiment on CS 2 . Three laser spots marked by numbers 1-3 are surrounded by photon-echo signals. The third-order signals are marked by yellow triangles. The signs correspond to two-pulse echo signals (one pulse interacts twice). s mark three-pulse photon-echo signals. The third-order signals are surrounded by the fifth-order signals. Green circles mark the fifth-order signals, where one pulse interacts three times. For example, the three upper right spots from up to down correspond to −2k 1 + 3k 3 , −k 1 − k 2 + 3k 3 and −2k 2 + 3k 3 , respectively. The signals where two different pulses interact twice are marked by pentagrams. The dominating red pentagram is the spot −2k 1 + 2k 2 + k 3 that is studied in this work.
In the calculations, we use a random distribution of FMO complexes in the x y-plane ( figure 1(A)). Generalization to 3D sample volume is possible and would add possible pulse propagation effects [38,39]. We assume that all pulses are vertically polarized. The electric field from the laser pulses at a certain FMO with a position x r is given by where A is the field amplitude and E p (t) is the normalized Gaussian field envelope function of the pulse p with the center position at t p and carrier frequency ω. k p x r gives the spatial phase of the pulse at the position of an FMO complex r .

Multiexcitonic density matrix
Each BChl is modeled as a three-level system having the ground state S 0 , the first excited state S 1 and a higher excited state S n that is responsible for the excited state absorption (ESA) from the S 1 state at the same spectral region as the S 1 absorption. The diagonalization of the excitonic Hamiltonian gives rise to delocalized exciton states, which can be divided into manifolds according to the number of excitations in the system. Dipolar transitions are allowed only between the neighboring manifolds. We only consider one-and two-exciton manifolds here. For the one-exciton manifold, we obtain |α 1 = m C α 1 (m)|m , and for the two-exciton manifold, we obtain |α 2 = m,n C α 2 (m, n)|m, n using the respective exciton expansion coefficients C α 1 (m) and C α 2 (m, n). It should be noted that the two-exciton states include the monomeric higher excited state S n via |m, m .
The laser-pulse-induced time evolution of the system is calculated using a multiexciton density matrix theory. The details of the approach are presented in [40]- [42]. The electric field of the laser pulse is included explicitly in the time evolution scheme. Relaxation processes within the manifolds are described using the Redfield relaxation theory highlighted below.
The relaxation rate from the N -exciton state α N to the N -exciton state β N is given by where n( ) denotes the Bose-Einstein distribution, andh (α N , β N ) is the energy difference between the two states. The excitonic spectral density J (α N β N , β N α N ; α N ,β N ) describes how efficiently the energy generated by the interstate relaxation between |α N and |β N can be deposited into vibrational modes of both the chromophore and the protein matrix. It is given by for the one-exciton states, and by for the two-exciton states. We use a local uncorrelated spectral density j (ω), which is obtained from fluorescence line narrowing spectroscopy [43]. The exciton-exciton annihilation rate k α 2 →β 1 [44], i.e. the rate from the two-exciton state α 2 to the one-exciton state β 1 , is calculated as The rate is proportional to the square of the overlap of a one-exciton wavefunction at a molecule and the doubly excited molecular part S n of the two-exciton wavefunction at the same molecule. Contributions from all molecules are added and multiplied by the internal conversion rate k IC between S n and S 1 of the BChl.

Equation of motion and signal
The equation of motion of the multiexciton density matrixρ (r ) of the rth complex is [41,42] ∂ ∂tρ The superoperators for the free time evolution of the excitons and the dissipation are given by L (r ) 0 and D (r ) . The field part of the propagator is given by the commutator with the field times the dipole operator: [41,42]. The polarization of the r th complex is obtained as a sum over all states of the density matrix α N |ρ (r ) |β N ±1 = ρ (r ) times the respective elements of the dipole operator α N |μ (r ) |β N ±1 = µ (r ) α N β N ±1 (N , N ± 1 is 0 for ground state, 1 for one-and 2 for two-exciton states). A discrete Fourier transformation of the polarizations to the k-space gives the polarization in the k out -direction, The polarization depends on the time delays τ and T between the pulses in E r and the signal time t defined in figure 1(A). Finally, the 2D Fourier transform is calculated to obtain the 2D spectra in the k out -direction [45], In the following, only the real part of the signal is discussed. The signal can be further divided into rephasing (τ > 0) and nonrephasing (τ < 0) contributions.

Results and discussion
We apply the above method to the FMO complex of the green sulfur bacterium Chlorobium tepidum [46]. The parameters are the same as in [35]. Random orientations of the complexes are used as well as Gaussian-distributed static site energy disorder of 80 cm −1 (FWHM). A diagonalization with respect to the dipole-dipole coupling gives rise to delocalized exciton states. As in [33,35], we use a temperature of 77 K and wavelength 805 nm for the three pulses. To get a broad spectrum, we used a pulse length of t = 13 fs, corresponding to a bandwidth of ω = 2 ln 2/(π c t) = 1100 cm −1 . The bandwidth is somewhat broader than what was used in 2D experiments of FMO [33]. Figures 2-4 show the calculated 2D spectra of FMO in the fifth-order phase-matched direction −2k 1 + 2k 2 + k 3 . In figure 2, the total signal at T = 0.2 ps and T = 1 ps are presented, whereas in figures 3 and 4, the nonrephasing (τ < 0) and rephasing (τ > 0) signals are shown separately for the two T points. Vertical and horizontal lines in the plots indicate the energies of the 28 two-exciton states (vertical) and 7 one-exciton states (horizontal). For nonoverlapping pulses, during the time τ , all Liouville pathways are in ground-state to two-exciton state coherence. Therefore the two-exciton states become visible in ω τ . In ω t also the one-exciton states become clearly visible because of the Liouville pathways, which are in ground-state to one-exciton coherence during the t evolution. The lowest-energy one-excitons are strongest in the linear absorption spectrum. The α 1 = 1, 2, 3 are also particularly clearly visible in 2D spectra as the prominent features coinciding with the three lowest frequency horizontal lines.  The relaxation from the α 1 = 2, 3 to the α 1 = 1 level is also obvious on comparing the spectra at T = 0.2 ps and T = 1 ps-the weight of the α 1 = 1 contribution is rising, whereas the upper two become weaker. Comparing the nonrephasing and rephasing contributions to the total signal, one can clearly note the usual antidiagonal elongation of the nonrephasing signal, whereas the rephasing signal shows diagonal elongation. The effect is a direct analogue to a similar behavior in the third-order 2D spectra of excitonically coupled dimers [47].
We do not use response functions and Feynman diagrams in our calculations. Still, the diagrams provide a convenient tool for discussing the qualitative features and behavior of the results. In figure 5, we have presented all Feynman diagrams that contribute to the signal in fifth order. The corresponding rephasing and nonrephasing diagrams are paired up. All diagrams are in 'population' during the time period T . For diagram 1, this is a ground-state population. For the other diagrams, this can be a true population ρ α N α N or coherence between exciton states of the same manifold ρ α N β N , each corresponding to separate diagrams. All of these have their own dynamics during T , which is not explicitly depicted in the diagrams. For example, the quantum beats observed in the 2D spectra of FMO [7] are due to the evolution of coherences in one-exciton manifold. Population relaxation during the T leads to still another set of diagrams. In figure 5, we have included the possibility that diagram 3 can relax from two-exciton to oneexciton manifold during T , leading to diagrams 5 and 6. We point out that since the number of      The 2D signals in figures 2-4 are negative over the whole plot. Each additional order contributes to the response function a multiplier i [48]. This means that the sign rules in fifth order are opposite to the sign rules in third order. In fifth order, an even number of interactions from the right lead to the negative sign of the Feynman diagram. The negative diagrams 1 and 2 are in the one-exciton to ground-state coherence during the last time period t. Contrary to that, the positive diagrams 3 and 4 are during that time period in the two-exciton to one-exciton coherence. Obviously, the latter contribution is weaker. In a similar way, the third-order 2D spectra of the FMO are mostly positive [33,35]-clearly even there the negative ESA signal, due to the pathways where the last evolution is in the two-exciton to one-exciton coherence, is weak compared with the positive bleaching and stimulated emission pathways. Since we use the three-level system to describe BChls, the total transition dipole moment is not reduced by having one excitation present in the system. We conclude that the reason for the negative 2D signal in fifth order is in the specific level structure and transition dipole moments of the FMO. This explains also the change in the total signal intensity during the T evolution. During that time, diagram 3 can relax via excitation annihilation to diagrams 5 and 6. Diagram 6 is negative but gives weak signal because of the two-exciton to one-exciton coherence during the last time period t (see the above argument). However, the positive diagram 5 gives strong signal-it is in one-exciton to ground-state coherence. The general features of the 2D signal from diagrams 5, 1 and 2 are the same. Consequently, the growing positive diagram 5 reduces the total negative signal amplitude during the time evolution T .
The fifth-order 2D spectra in figures 2-4 resemble the expected third-order double quantum coherence 2D (2Q2D) spectra of such a system [49]. Also, the 2Q2D spectra would probe the two-exciton manifold structure in one frequency axis and one-exciton transitions in the other frequency axis. Since the 2Q2D is inherently a nonrephasing technique, we expect that those spectra would be elongated along the antidiagonal as the nonrephasing signals here. However, since the 2Q2D is third order, the signal will be positive as the conventional third-order 2D spectrum of FMO [33].
In our nonperturbative approach, all orders of the signal are inherently present in the same way as in experiment. This means that the seventh-order signals like −2k 1 + 2k 2 + 2k 3 − k 3 will contribute, if present. Analogously, in [35] it was found that at higher intensities, conventional 2D spectra calculated from the direction −k 1 + k 2 + k 3 can have contributions from the fifthorder signal. In relative terms, the difference between the seventh-and fifth-order contributions is smaller than the difference between the fifth-and third-order contributions. To judge the amount of possible higher order contributions, one should vary the field amplitude A in equation (1) as it was done in [35]. In fifth-order experiment, one needs also to be aware of the possible contributions from cascading third-order processes [50], which are not present in our calculations. The cascading-induced signals can be also included if the calculations are carried out for the 3D sample and the field generated by the earlier excited complexes is included in the incoming field equation (1).
We point out that very recently the experimental 3D fifth-order signal of a dye molecule was measured [51]. The phase-matching schema in that work is the same as that in our calculations if the pulse number 3 comes first. In this case, both pulse delays correspond to the evolution of inter-manifold coherences and one can carry out Fourier transformation along both times leading to 3D representation. To calculate such signals in multichromophoric systems like here, one should include possible transitions from two-exciton manifold to three-exciton manifold. In fact, also for the case studied here, there would be a pair of diagrams that involve three-exciton manifold. If we take diagram 3 (see figure 5) and let the third pulse interact from the left, we would obtain two diagrams with negative sign that generate signal from threeexciton to two-exciton coherence. Such calculations are beyond the scope of the present work and will be a topic of future studies.

Summary
We use nonperturbative computation schema based on numeric integration of the equation of the motion of multiexciton density matrix of the FMO pigment protein complex. A clear signal can be extracted in phase-matching directions owing to the interference of the polarizations generated by a large number of randomly placed unoriented energetically inhomogeneous complexes that interact with the laser field. In this paper, the signal in a fifth-order phasematching direction −2k 1 + 2k 2 + k 3 is analyzed. As in experiment, the signal field is used to construct electronic 2D spectra showing clear features that correspond to two-exciton and oneexciton energies. The results are discussed in terms of Feynman diagrams including exciton annihilation and will be useful in disentangling fifth-order experiments, once they are performed in such systems.