Exciton coupling induces vibronic hyperchromism in light-harvesting complexes

In photosynthesis, light absorption and charge separation take place in specialised multi-chromophoric systems. While the reaction centre (RC) complexes where primary charge separation occurs are very similar in different photosynthetic systems, there are a wide variety of light-harvesting antennas [1]. Clearly, it is possible to optimise the absorption of light and make excitation transport to the RC efficient in many different ways. One of the strategies used in most photosynthetic systems is to organise the pigments energetically in an energy funnel with the RC in the bottom of the sink [2]. In addition, the pigment molecules are usually close enough to enable sufficient resonant interaction and thereby fast excitation transport. Such general optimisation strategies have been widely reported and are well understood in photosynthetic antenna systems.

One result of the relatively dense packing of the pigment molecules in the antenna complexes is the emergence of delocalised Frenkel exciton states [3–6]. The latter has led to the recognition of the role of memory effects and coherence in Frenkel exciton dynamics in photosynthetic light-harvesting. This was recently revitalised by novel coherent multidimensional spectroscopy measurements revealing long-lived oscillatory features in the so-called Fenna–Mathews–Olson (FMO) light-harvesting antenna complex [7]. The oscillations were interpreted as electronic coherences, and their long lifetime was taken as evidence for possible coherent excitation transport. It was argued that parallel coherent excitation transport pathways may enable interference-based quantum optimisation possibilities that make excitation transfer through the FMO complexes more efficient. Analogous observations were also reported for light-harvesting complexes in photosynthetic marine algae [8]. This initiated numerous new theoretical investigations of the role of coherence in excitation energy transport [9–11]. An important conclusion drawn from these efforts concerned the optimal parameter region of the transport—it was realized that the most efficient transport takes place in the intermediate regime between fully coherent and incoherent transport, where the inter-pigment resonance interaction and the system-bath interaction are not so different [12, 13]. Interestingly, the efficiency curve has a rather flat maximum, meaning that optimal transport can be achieved by quite a broad set of parameters.

In contrast to the large number of simulations following excitation dynamics [10, 14–17], few studies have made the effort to calculate the coherent multidimensional spectroscopy signals, which would directly correspond to what is observed experimentally [18–20]. It turned out that reproducing the long-lived oscillations in calculated 2D spectra as electronic coherences is not straightforward. Even studies based on hierarchy equations of motion (HEOM), where bath memory effects are included, could not reach agreement with the experimental observations [21, 22]. Additional assumptions had to be made. For example, assuming correlated static and dynamic disorder at different pigment sites led to long-lived electronic coherences [23, 24]. However, neither correlated nuclear motions nor inhomogeneous broadening has any independent support. Quite the opposite. Thorough molecular dynamics studies did not show any evidence for such correlations [25]. A recent study using the HEOM approach produced long-lasting electronic coherences [26]. The authors used a spectral density that fitted the fluorescence line-narrowing experiments well [27] for frequencies above $50\;\text{c}{{\text{m}}^{-1}}$. However, the low-frequency region, which is particularly important for electronic dephasing and is most sensitive to temperature, had an unrealistically small amplitude.

Initially, the possible vibrational origin of the oscillatory features in the two-dimensional (2D) signal in FMO was discarded since the nuclear modes in the frequency region of the observed oscillations ($\sim 150\;\text{c}{{\text{m}}^{-1}}$) are very weakly coupled to the electronic transition, leading to a negligible Huang–Rhys factor S. In the $180\;\text{c}{{\text{m}}^{-1}}$ region, there are a few modes with the total S reaching almost 0.05 [27, 28], but even this is too weak to produce significant oscillatory features in 2D signals. The transition energies of the bacteriochlorophyll (BChl) molecules in FMO are quite well established [29]. The 0–1 vibronic transition energy of the $180\;\text{c}{{\text{m}}^{-1}}$ mode region on the red-most BChl is almost in resonance with the 0–0 transition energy of one of its neighbours. The consequences of such resonance for natural light-harvesting have not been fully analysed before. A vibronic exciton model where one mode per pigment is treated explicitly [30] was applied to calculate the electronic 2D response functions of FMO [31]. It was shown that, owing to the resonance, the weak transition of mainly vibronic origin at the lowest energy BChl obtains a significant additional oscillator strength due to mixing with the strong purely electronic transition at the neighbouring molecule. Furthermore, the long-lived oscillations observed in the 2D experiments were explained by the fact that the energy fluctuations of the vibrational levels of a molecule are correlated, and consequently the corresponding coherences between vibronic levels have a long lifetime. The intensity borrowing is quite robust, and the moving to some extent out from the resonance due to the inhomogeneous broadening does not influence the results appreciably.

An analogous resonance between the 0–0 and 0–1 transitions of the neighbouring pigments was recently used in an excitonic dimer model, where the adiabatic potential energy surfaces were calculated [32, 33]. The authors argue that because of the near resonance between the levels, the adiabatic approximation breaks down and non-adiabatic coupling has to be taken into account. The non-adiabatic coupling causes mixing of the states, leading to intensity redistribution much like that in [31]. Excitation dynamics due to non-adiabatic coupling in the context of light-harvesting complexes was analysed in terms of curve crossing and surface hopping [34, 35]. Schröter and Kühn studied the interplay between non-adiabatic dynamics and Frenkel exciton transfer in aggregates where both the ${{\text{S}}_{1}}$ and ${{\text{S}}_{2}}$ transitions were considered [36]. Excitation annihilation has also been described as non-adiabatic coupling between one- and two-exciton manifolds [37, 38]. Tiwari et al [32] investigated a different situation where the potentials are nested and the crossing points are far from the equilibrium, resembling the situation in internal conversion. Using 2D signal calculations, the authors argue that the long-lived oscillating features of the 2D spectroscopy of light-harvesting systems could be explained in terms of the ground-state coherent nuclear motions excited via enhanced transitions with strong vibronic character.

Several other studies have recently addressed the issue of how to separate the electronic and vibrational quantum beats in electronic 2D spectroscopy. A comparison of the oscillatory patterns and oscillation phases in two representative systems—displaced oscillator and excitonic dimer—was reported [39]. The model was further developed for a dimer with the ground-state vibrational levels included [40]. The ground and excited state vibrational coherences were thoroughly compared in the vibronic exciton dimer model [41]. The oscillatory phonon features in colloidal quantum dots were modelled [42], and it was shown that following separately positive and negative frequency components of the population time Fourier transform gives additional control over the Liouville pathways, thereby enabling a more clear distinction between the electronic and vibrational beatings [43].

Here, we point out that the method used by Christensson et al [31], the so-called one particle approximation (OPA), only includes excited state vibrations. In order to account for the ground-state vibrational states, at least the two particle approximation (TPA) [44, 45] has to be employed. TPA is exact in the case of a dimer, and has been used for analyses of 2D spectra [32, 40]. Various analogous approaches have been applied for modelling a vibronic dimer [46, 47], and so an obvious question arises: what effects does one miss by using simple OPA calculations, as in [31], compared to the exact solution of the problem? Here, we address the issue by carrying out a comparison of OPA and TPA using a model dimer resembling two neighbouring low-energy BChl molecules in the FMO. The numerically exact reference is provided by the multi-configuration time-dependent Hartree (MCTDH) method [48, 49].

The article is organised as follows. We start by presenting the basics of the vibronic exciton model and the OPA as well as MCTDH approaches. The theory is followed by comprehensive model calculations of a vibronic heterodimer for various parameter sets. We show that despite the dominantly monomeric character of the heterodimer, the intensities of the vibrational features in the fluorescence spectrum are significantly affected. This has consequences in how to use fluorescence line-narrowing spectroscopy to experimentally determine the spectral densities [53]. In the final part, the results are discussed and conclusions formulated.

2.1. The Frenkel exciton Hamiltonian

In the following, we will use the Frenkel exciton Hamiltonian to describe the coupled electronic (zero, $\left| 0 \right\rangle$, and one-exciton, $|m\rangle$, space) and nuclear degrees of freedom, $\mathbf{Q}=\left\{ {{\mathbf{Q}}_{m}} \right\}$, [3]

$$\begin{eqnarray}&&{{H}_{\text{agg}}}\left( \mathbf{Q} \right)={{H}^{(0)}}\left( \mathbf{Q} \right)+{{H}^{(1)}}\left( \mathbf{Q} \right),\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&{{H}^{(0)}}\left( \mathbf{Q} \right)=\underset{m}{\mathop{\mathop{\sum }^{}}}\,{{H}_{m,g}}\left( {{\mathbf{Q}}_{m}} \right)\left| 0 \right\rangle\left\langle 0 \right|{{\equiv\mathcal{E}}_{0}}\left| 0 \right\rangle\left\langle 0 \right|,\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&{{H}^{(1)}}\left( \mathbf{Q} \right)=\underset{mn}{\mathop{\mathop{\sum }^{}}}\,\left[ {{\delta }_{mn}}\left( {{\mathcal{E}}_{0}}+{{U}_{m,e}}\left( {{\mathbf{Q}}_{m}} \right) \right)+{{J}_{mn}} \right]\left| m \right\rangle\left\langle n \right|.\end{eqnarray} \tag{ 3 }$$

Denoting the monomeric electronic states by $\left| {{a}_{m}} \right\rangle$, where $a=\left( g,e \right)$ are the ground and excited states, one has $\left| 0 \right\rangle={{\mathop{\prod }^{}}_{m}}\left| {{g}_{m}} \right\rangle$ and $\left| m \right\rangle=\left| {{e}_{m}} \right\rangle{{\mathop{\prod }^{}}_{n\ne m}}\left| {{g}_{n}} \right\rangle$. The nuclear motion will be described in the displaced oscillator model, i.e. the ground-state Hamiltonian and excited state coupling read for each site (${{E}_{m}}$: electronic energy)

$$\begin{eqnarray}&&{{H}_{m,g}}\left( {{\mathbf{Q}}_{m}} \right)=\underset{\xi }{\mathop{\mathop{\sum }^{}}}\,\frac{\hbar {{\omega }_{\xi }}}{2}\left( -\frac{{{\partial }^{2}}}{\partial Q_{m,\xi }^{2}}+Q_{m,\xi }^{2} \right),\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&{{U}_{m,e}}\left( {{\mathbf{Q}}_{m}} \right)={{E}_{m}}+\underset{\xi }{\mathop{\mathop{\sum }^{}}}\,\hbar {{\omega }_{\xi }}{{g}_{\xi }}{{Q}_{m,\xi }}.\end{eqnarray} \tag{ 5 }$$

Here, ${{\omega }_{\xi }}(\xi =1,\ldots ,{{N}_{\text{vib}}})$ is the vibrational frequency of mode ${{Q}_{m,\xi }}$ (note the use of dimensionless units), which is assumed to be identical for the different monomers. The same approximation is made for the linear coupling constant, ${{g}_{\xi }}$, which relates to the Huang–Rhys factor as ${{S}_{\xi }}=g_{\xi }^{2}/2$.

The absorption and emission spectra will be calculated in Condon approximation for the monomeric transition dipole matrix elements ${{d}_{m}}$, summed to give the aggregate dipole operator according to

$$\begin{eqnarray}&&d=\underset{m}{\mathop{\mathop{\sum }^{}}}\,{{d}_{m}}\left| m \right\rangle\left\langle 0 \right|+\text{h}\text{.c}..\end{eqnarray} \tag{ 6 }$$

2.2. The n-particle approach

The problem of coupled exciton-vibrational dynamics can be approached by an n-particle approximation scheme [44]. To this end, we introduce vibrational states for the different potential energy surfaces according to $\left| {{M}_{m,a}} \right\rangle$, where it is understood that ${{M}_{m,a}}$ could be a multi-index in cases where several vibrational modes per monomer need to be taken into account. The eigenstates $\left| {{\alpha }^{\left( 1 \right)}} \right\rangle$ of ${{H}^{\left( 1 \right)}}\left( \mathbf{Q} \right)$ can be expanded as

$$\begin{eqnarray}&&\left| {{\alpha }^{\left( 1 \right)}} \right\rangle=\underset{m,\mu }{\mathop{\mathop{\sum }^{}}}\,{{C}_{m,\mu }}\left( {{\alpha }^{\left( 1 \right)}} \right)\left| m\mu \right\rangle+\underset{mn}{\mathop{\mathop{\sum }^{}}}\,\underset{\mu \nu }{\mathop{\mathop{\sum }^{}}}\,{{C}_{mn,\mu \nu }}\left( {{\alpha }^{\left( 1 \right)}} \right)\left| m\mu ,n\nu \right\rangle+\ldots \end{eqnarray} \tag{ 7 }$$

with the one-particle states

$$\begin{eqnarray}&&\left| m\mu \right\rangle=\left| {{e}_{m}} \right\rangle\left| \mu ={{M}_{m,e}} \right\rangle \underset{n\ne m}{\mathop{\mathop{\prod }^{}}}\,\left| {{g}_{n}} \right\rangle|{{0}_{n,g}}\rangle \end{eqnarray} \tag{ 8 }$$

and the two-particle states

$$\begin{eqnarray}&&\left| mn,\mu \nu \right\rangle=\left| {{e}_{m}} \right\rangle\left| \mu ={{M}_{m,e}} \right\rangle\left| {{g}_{n}} \right\rangle\left| \nu ={{N}_{n,g}} \right\rangle \underset{k\ne m,n}{\mathop{\mathop{\prod }^{}}}\,|{{g}_{k}}\rangle |{{0}_{k,g}}\rangle \end{eqnarray} \tag{ 9 }$$

In principle, equation (7) will contain further terms, but for the present case of a molecular heterodimer, the two-particle ansatz is already exact. The restriction to the first term is the so-called OPA.

Using the eigenstates, equation (7), the absorption spectrum will be calculated in the zero temperature limit as (${{I}_{0}}$ normalisation constant)

$$\begin{eqnarray}&&I\left( \omega \right)={{I}_{0}}\omega \underset{{{\alpha }^{\left( 1 \right)}}}{\mathop{\mathop{\sum }^{}}}\,\frac{\Gamma {{\left| \left\langle {{\alpha }^{\left( 1 \right)}} \right|d\left| {{\alpha }^{\left( 0 \right)}}=0 \right\rangle \right|}^{2}}}{{{\left( \omega -{{\omega }_{{{\alpha }^{\left( 1 \right)}}0}} \right)}^{2}}+{{\Gamma }^{2}}},\end{eqnarray} \tag{ 10 }$$

where $\Gamma$ is a parameter mimicking the finite line width (dephasing time) of the real system, $\left| {{\alpha }^{\left( 0 \right)}} \right\rangle$ denotes the eigenstates of ${{H}^{\left( 0 \right)}}$, with $\left| {{\alpha }^{\left( 0 \right)}}=0 \right\rangle$ being the overall ground state (see below), and ${{\omega }_{{{\alpha }^{\left( 1 \right)}}{{\alpha }^{\left( 0 \right)}}}}$ the transition frequency. The emission spectrum is calculated as

$$\begin{eqnarray}&&F\left( \omega \right)={{F}_{0}}\underset{{{\alpha }^{\left( 1 \right)}},{{\alpha }^{\left( 0 \right)}}}{\mathop{\mathop{\sum }^{}}}\,P\left( {{E}_{{{\alpha }^{\left( 1 \right)}}}} \right)\frac{\Gamma {{\left| \left\langle {{\alpha }^{\left( 1 \right)}} \right|d\left| {{\alpha }^{\left( 0 \right)}} \right\rangle \right|}^{2}}}{{{\left( \omega -{{\omega }_{{{\alpha }^{\left( 1 \right)}}{{\alpha }^{\left( 0 \right)}}}} \right)}^{2}}+{{\Gamma }^{2}}},\end{eqnarray} \tag{ 11 }$$

where $P\left( {{E}_{{{\alpha }^{\left( 1 \right)}}}} \right)$ is the Boltzmann population of a one-exciton vibrational state with energy ${{E}_{{{\alpha }^{\left( 1 \right)}}}}$ and ${{F}_{0}}$ is a normalisation constant.

2.3. The MCTDH approach

In principle, the n-particle approach provides a systematic route to the exact eigenstates of the one-exciton Hamiltonian. However, with increasing aggregate size and number of vibrational coordinates, one will soon face a dimensionality bottleneck and the problem will become numerically intractable. An efficient alternative approach is the MCTDH method, which rests on the expansion of the time-dependent state vector into a basis of time-dependent Hartree products that are composed of single particle functions (SPFs) [48, 49]. Applications to exciton dynamics and spectroscopy have been given in [50, 51].

First, the state vector is expanded in terms of the one-exciton basis, i.e

$$\begin{eqnarray}&&\left| \Psi \left( \mathbf{Q};t \right) \right\rangle=\underset{m}{\mathop{\mathop{\sum }^{}}}\,{{\chi }_{m}}\left( \mathbf{Q};t \right)\;\left| m \right\rangle.\end{eqnarray} \tag{ 12 }$$

In the next step, the nuclear wave function is written in MCTDH form as follows

$$\begin{eqnarray}&&{{\chi }_{m}}\left( \mathbf{Q},t \right)=\underset{{{j}_{1}}\ldots {{j}_{D}}}{\overset{{{n}_{{{j}_{1}}}}\ldots {{n}_{{{j}_{D}}}}}{\mathop{\mathop{\sum }^{}}}}\,C_{{{j}_{1}},\ldots ,{{j}_{D}}}^{\left( m \right)}\left( t \right)\phi _{{{j}_{1}}}^{\left( m \right)}\left( {{Q}_{1}};t \right)\ldots \;\phi _{{{j}_{D}}}^{\left( m \right)}\left( {{Q}_{D}};t \right).\end{eqnarray} \tag{ 13 }$$

Here, $C_{{{j}_{1}},\ldots ,{{j}_{D}}}^{\left( m \right)}\left( t \right)$ are the time-dependent expansion coefficients weighting the contributions of the different Hartree products, which are composed of SPFs, $\phi _{{{j}_{k}}}^{\left( m \right)}\left( {{Q}_{k}};t \right)$, for the kth nuclear degree of freedom in state m. Overall, there are $D={{N}_{\text{agg}}}\times {{N}_{\text{vib}}}$ nuclear degrees of freedom.

Although it is in principle possible to obtain eigenstates with the MCTDH method, being a time-dependent approach it is more suited to solving the time-dependent Schrödinger equation. Therefore, we will use this to calculate the absorption spectrum according to the time-dependent reformulation of equation (10) [3]

$$\begin{eqnarray}&&I\left( \omega \right)={{I}_{0}}\omega \;\operatorname{Re}\mathop{\int }^{}_{0}^{\infty }dt\;{{e}^{i\omega t-\Gamma t}}\left\langle {{\Psi }_{0}} \right|d{{e}^{-i{{H}_{\text{agg}}}\left( \mathbf{Q} \right)t/\hbar }}d\left| {{\Psi }_{0}} \right\rangle,\end{eqnarray} \tag{ 14 }$$

where d is the dipole operator, equation (6), and $\left| {{\Psi }_{0}} \right\rangle$ is the ground-state wave function. Since there is no coupling in ${{H}^{\left( 0 \right)}}$, the ground state is given as a Hartree product, i.e.

$$\begin{eqnarray}&&\left| {{\Psi }_{0}} \right\rangle=\phi _{0}^{\left( 0 \right)}\left( {{Q}_{1}} \right)\ldots \phi _{0}^{\left( 0 \right)}\left( {{Q}_{D}} \right)\left| 0 \right\rangle\end{eqnarray} \tag{ 15 }$$

with $\phi _{0}^{\left( 0 \right)}\left( {{Q}_{j}} \right)$ being the respective ground-state wave function for the jth mode. All wave packet propagation simulations have been performed using the Heidelberg program package [52]. The MCTDH dimer setup, which has been applied here, includes three electronic states (one ground state and two singly excited states). Since the ground state has been used for the preparation of the initial wave packet only, one SPF per mode was sufficient to describe it properly. The actual propagation of the wave packet (propagation time 6 ps) takes place on the two excited states; here, four SPFs were necessary to obtain converged results. For the primitive basis we used harmonic oscillator (20 points) discrete variable representation within the interval [−3.5:3.5] for all modes. For all calculations, the multi-set method was used.

In the following application, the case of two sites (vibronic heterodimer, $m=\left( D,A \right)$) will be considered. The parameters are chosen to mimic the situation in the FMO complex. A number of studies proposed an electronic Hamiltonian of the complex, based on the modelling of a set of spectroscopic observables. For a review, see [54]. There is a general agreement that BChl 3 and 4 have the lowest site energies. The transition energy difference is proposed to be from about 110 to $180\;\text{c}{{\text{m}}^{-1}}$. Excitonic coupling between these two molecules is negative, with values from about −50 to $-75\;\text{c}{{\text{m}}^{-1}}$ [29]. Fluorescence line-narrowing spectroscopy has revealed that at around $180\;\text{c}{{\text{m}}^{-1}}$, there are a few vibrational modes in the FMO's lowest-energy BChl with a total Huang–Rhys factor of 0.03 [27]. We point out that the analogous experiments with BChl in triethylamine give a somewhat higher value of 0.04 [28]. Having these parameters in mind, we formulate our vibronic dimer as follows. The properties of the eigenstates and the spectra of the model will be scrutinised depending on the Coulomb coupling $J\equiv{{J}_{\text{DA}}}$ and the heterogeneity, $\Delta E\equiv{{E}_{\text{D}}}-{{E}_{\text{A}}}$. For reference, we use $\Delta E=180\;\text{c}{{\text{m}}^{-1}}$ and ${{J}_{\text{DA}}}/hc=-90\;\text{c}{{\text{m}}^{-1}}$. First, each site is coupled to one vibrational mode with parameters $\omega /2\pi c=180\;\text{c}{{\text{m}}^{-1}}$ and S = 0.05. In the following, we will start with a discussion of the differences between exact and OPA in the absorption spectrum for transitions to the one-exciton manifold. The focus will be on the dependence on detuning between the local excitation energies and the Coulomb coupling. In order to study the multi-mode effects, a second model including additional vibrational modes with varying frequencies and Huang–Rhys factors will be considered in an up to three-mode model ($\omega /2\pi c=365\;\text{c}{{\text{m}}^{-1}}$, S = 0.019 and $\omega /2\pi c=260\;\text{c}{{\text{m}}^{-1}}$, S = 0.023). The actual frequencies of the modes were chosen as approximately 1.5 and 2 times the original frequency $\omega /2\pi c=180\;\text{c}{{\text{m}}^{-1}}$. The corresponding Huang–Rhys factors were chosen to represent the total S of the nearby modes reported in [28]. Finally, the temperature dependence is addressed in terms of the emission spectrum.

3.1. The validity of OPA

The dependence of the exact and OPA absorption spectra on the Coulomb coupling ${{J}_{\text{DA}}}$ ($\Delta E/\hbar \omega =1$) and the detuning $\Delta E$ (${{J}_{\text{DA}}}/\hbar \omega =-0.5$) are compared in figure 1. This figure also contains the respective dependencies of the bare energy levels.

Zoom In Zoom Out Reset image size

First we notice that for ${{J}_{\text{DA}}}=0$, the two cases show different degeneracies of energy levels due to the restricted state space of the OPA. Looking at the dependence on ${{J}_{\text{DA}}}$ upon increasing its value, the energy levels are shifted and the degeneracies are partly lifted. Around ${{J}_{\text{DA}}}/\hbar \omega \approx 0.85$ and $\approx 1.35$, the energy levels approach each other, showing partially avoided crossings in the exact case. This pattern is not at all reproduced by the OPA. Overall, only the lowest state shows an agreement between the exact and OPA calculations. The same conclusion can be drawn from the $\Delta E$ dependence shown in figures 1(d), (e).

Comparing the quite different behaviour of the spectrum of ${{H}^{\left( 1 \right)}}$ for the exact and OPA models, the question of how this will reflect in the absorption spectra arises. Here, we notice from figure 1 that the differences between the exact and OPA are particularly large for the ${{J}_{\text{DA}}}$-dependence. Of course, whether or not deviations are visible depends on the distribution of the oscillator strength. Therefore, the positive coupling part of the ${{J}_{\text{DA}}}$-dependence is more visibly influenced. As a consequence, OPA gives reasonable results as far as the $\Delta E$-dependence is concerned since ${{J}_{\text{DA}}}<0$ for the present model. Generally, there is a trend that despite the radically different energy level structure in the two models, when it comes to the intensities of the spectral features, the models give surprisingly similar results. If one uses broader line shapes, the resulting spectra will not be very different for most of the used parameters. Interestingly, it is not true that OPA necessarily yields a simpler spectrum, as one would expect on the basis of the energy level structure of the one-exciton Hamiltonian. Indeed, one can state that the exact case shows a more pronounced collective effect, i.e. oscillator strength is essentially located in a single transition.

3.2. Multi-mode effects

In the following, we discuss the effect of the multiple modes coupling to the electronic transitions. The results of the MCTDH simulations of the absorption spectra are shown in figure 2. Overall, one can state that since for large negative ${{J}_{\text{DA}}}$ the oscillator strength is concentrated in the lowest state that has a dominantly electronic character, a pronounced effect of the further modes is seen for positive couplings only. For the case of FMO, this leads to quite similar spectra, as a function of the detuning $\Delta E$, as can be seen from the right column of figure 2. In the case of positive coupling, the details of the spectral change depend, of course, on the mode parameters, but the net effect is a broadening due to the more complex structure of the exciton-vibronic states.

Zoom In Zoom Out Reset image size

3.3. The temperature dependence of the emission spectra

The emission spectrum results for the one-mode model as a function of Coulomb coupling and detuning are shown in figure 3 for two different temperatures. Owing to the small Huang–Rhys factor, the spectra at 4 K are essentially showing the 0–0 transition and a small shoulder due to the 0–1 transition in figure 3(a). The intensity of these transitions gradually drops upon increasing ${{J}_{\text{DA}}}$ in the positive domain due to H-aggregate formation. This does not occur upon increasing the detuning in panel (c), where the spectrum becomes close to that of a monomer for large $\Delta E$. For $T=300\;\text{K}$, the spectrum is much more structured, as can be seen from figures 3(b), (d). The level structure with a complex intensity pattern makes it impossible to draw a priori conclusions on the thermal occupation. Notice that due to the Coulomb coupling, the spacing between the peaks above the 0–0 transition is not equal to a vibrational quantum. Since the lowest energy fluorescence feature is due to the ground-state vibrational structure, it is still shifted from 0–0 by the vibrational frequency.

Zoom In Zoom Out Reset image size

3.4. Excitonic distortions of fluorescence site-selection spectra

The ratio between the 0–0 and 0–1 emission intensities is investigated in more detail in figure 4. For a monomer, this ratio would be equal to the Huang–Rhys factor used in the calculations, S = 0.05 in the present case. This value is observed only for ${{J}_{\text{DA}}}=0$ (marked as a crossing of red arrows in figure 4) and for a finite ${{J}_{\text{DA}}}$ in the limit of large detunings (not shown). At this point, we should recall that Huang–Rhys factors and mode frequencies are usually obtained from low-temperature site-selected fluorescence spectra. Clearly, the outcome of such an experiment would be influenced by the excitonic coupling. The case of FMO would approximately correspond to the blue curve of figure 4 at the relative detuning of −1, marked as a crossing of blue arrows. Our calculations for that point give the 'effective observable' of about S = 0.034. We point out that the Huang–Rhys factor obtained from the FMO experiment is smaller than the S reported for BChl in solution. Our calculations give a straightforward explanation of this discrepancy. With the help of our results illustrated in figure 4, the Huang–Rhys factors obtained from the fluorescence site-selection spectroscopy of light-harvesting complexes and other molecular aggregates can be corrected. Obviously, the correction can easily be as large as 100%. In the case of the FMO in the $180\;\text{c}{{\text{m}}^{-1}}$ region, the observed S should be multiplied by a factor of 1.5. The correction factor becomes smaller for the higher frequency modes. Theoretical studies where the spectral density is extracted from atomic simulations and compared to the fluorescence spectra for benchmarking [25, 55], need to consider these correction factors.

Zoom In Zoom Out Reset image size

Finally, we note that the measured S values will be reduced in the case of negative ${{J}_{\text{DA}}}$ values (J-aggregate) and increased for the opposite sign (H-aggregate).

We carried out a thorough comparison of OPA and exact calculations of vibronic excitons in a dimer model. We found that the general behaviour of the calculated absorption spectra is surprisingly similar, despite the radically different energy-level structure in the two models. Multi-mode effects were discussed, depending on the mode parameters. In the case of negative coupling, which is relevant for the FMO complex, most of the oscillator strength is concentrated in the transition to a state of essentially electronic character, such that no pronounced effect on the spectrum arises. Calculations of the fluorescence spectra show that the Huang–Rhys factors obtained from fluorescence spectroscopy of light-harvesting complexes with significant excitonic couplings need to be corrected. Depending on the excitonic coupling strength and detuning, the correction can be as large as 100%. At the same time, the observed frequencies are not affected if measured at low temperature, where the transitions from higher energy levels do not provide any significant contribution.

The authors acknowledge support from KAW, the Swedish Energy Agency, the Swedish Research Council, the Deutsche Forschungsgemeinschaft (Sfb652), and the EU Erasmus program.