Theory of 2D electronic spectroscopy of water soluble chlorophyll-binding protein (WSCP): Signatures of Chl b derivate

Riedl, Michael; Renger, Thomas; Seibt, Joachim

doi:10.1063/5.0200876

Author Notes

We investigate how electronic excitations and subsequent dissipative dynamics in the water soluble chlorophyll-binding protein (WSCP) are connected to features in two-dimensional (2D) electronic spectra, thereby comparing results from our theoretical approach with experimental data from the literature. Our calculations rely on third-order response functions, which we derived from a second-order cumulant expansion of the dissipative dynamics involving the partial ordering prescription, assuming a fast vibrational relaxation in the potential energy surfaces of excitons. Depending on whether the WSCP complex containing a tetrameric arrangement of pigments composed of two dimers with weak excitonic coupling between them binds the chlorophyll variant Chl a or Chl b, the resulting linear absorption and circular dichroism spectra and particularly the 2D spectra exhibit substantial differences in line shapes. These differences between Chl a WSCP and Chl b WSCP cannot be explained by the slightly modified excitonic couplings within the two variants. In the case of Chl a WSCP, the assumption of equivalent dimer subunits facilitates a reproduction of substantial features from the experiment by the calculations. In contrast, for Chl b WSCP, we have to assume that the sample, in addition to Chl b dimers, contains a small but distinct fraction of chemically modified Chl b pigments. The existence of such Chl b derivates has been proposed by Pieper et al. [J. Phys. Chem. B 115, 4042 (2011)] based on low-temperature absorption and hole-burning spectroscopy. Here, we provide independent evidence.

I. INTRODUCTION

In biomolecular complexes of bacteria and higher plants, which are often designated for specific functions, such as light harvesting, functional molecules are usually embedded in a protein scaffold. The interaction of such molecules, which, in the case of light-absorbing properties, are called pigments, with the protein environment can tune their spectroscopic properties, such as electronic excitation energies and excitonic couplings.^1–14 The water-soluble chlorophyll binding protein (WSCP), as a specific example of such a complex, is rather associated with photoprotection than with light harvesting.^15–17 It contains four chlorophyll units in an arrangement with quasi-D2 symmetry,^12,15 where two pairs of them form dimers with strong intra- and weak inter-dimer coupling. The symmetry of the arrangement of the chlorophylls in WSCP with a very similar protein environment of each pigment leads to comparable electronic properties of the individual pigments, particularly in view of their local transition energies (site energies). Furthermore, due to the arrangement of the chlorophylls within the protein scaffold and the resulting limitation of wave function overlap, short-range contributions to the site energies and excitonic couplings become negligible. Because of these aspects, WSCP represents a relatively simple model system to study pigment–pigment and pigment–protein interactions.^12,18–30 While, in such a context, findings from linear spectroscopy can be revealing already, nonlinear spectroscopic techniques, such as two-dimensional (2D) electronic spectroscopy, can yield more comprehensive insight into energetic and structural properties and into exciton-vibrational dynamics. Therefore, among other applications, it has often been used to study energy transfer and dissipation in molecular aggregates.^1,24,31–34 In a recent work by Fresch et al.,²⁹ unexpected differences in 2D spectra of WSCP containing the chlorophyll variants Chl a and Chl b have been reported. It was suggested that the differences could be caused by the following effect: the formation of a hydrogen (H) bond leads to an increased transition dipole moment and, as a consequence of the dipole–dipole interaction, to an increased excitonic coupling. Non-trivial differences between optical line shapes of Chl a and Chl b WSCP have been reported earlier using linear absorption and hole-burning spectroscopy.²² These differences have been interpreted by assuming the presence of a small amount of chemically modified Chl b, termed Chl b derivate, which exhibits a redshifted site energy compared to Chl b. To investigate these hypotheses quantitatively, we developed a theoretical description of 2D electronic spectroscopy and apply it to the linear and nonlinear experiments on Chl a and Chl b WSCP.²⁹ Our theory combines the concept of response functions with a time-convolutionless treatment of the system–bath (pigment–protein) interaction. The chosen partial ordering prescription (POP) relies on a specific ordering of the system–bath interaction in the framework of a cumulant expansion approach.³⁵ Even though this treatment neglects correlation effects between the vibrational dynamics during different time intervals of the excitation sequence resulting in the detected nonlinear signal,³⁶ it turned out to capture the relevant signatures of electronic excitations subject to our investigations. It furthermore facilitates a more efficient calculation of 2D spectra of larger systems because of the separability of the evolution during the different time intervals. In this context, implicit assumptions about thermal equilibration at different stages of an excitation pathway require closer inspection, as will be discussed in detail later.

This article is organized as follows: in Sec. II, we describe in detail how the response functions contributing to the third-order polarization are derived, thereby referring to the concept of POP and related approximations, e.g., concerning Markovianity of certain parts of the system–bath coupling, which have turned out to be appropriate for a comprehensive yet compact description.³⁵ In Sec. III, we discuss the features in calculated 2D spectra of Chl a WSCP and Chl b WSCP and their compatibility with the findings from the experiment and figure out in which respect 2D spectroscopy yields more comprehensive information about the investigated systems than linear absorption and circular dichroism (CD) spectroscopy. In particular, we investigate which modifications of the model assumptions in the description of Chl b WSCP are appropriate in view of reproducing findings from the experimental results. In this context, we also investigate whether the suggested mechanism of an enhanced inter-dimer coupling,²⁹ which leads to exciton delocalization within the tetrameric structure rather than in the separate dimers or the alternative derivate hypothesis,²² is appropriate for explaining the measured spectra.²⁹ In Sec. IV, we summarize our findings and discuss open questions. The experimental data of Ref. 29 were provided by Fresch et al.

II. THEORETICAL BACKGROUND

A. Components of Hamiltonian

We are aiming at the description of an aggregate of molecules with intermolecular excitonic coupling and interactions with an environment causing fluctuations of molecular excitation energies. For the treatment of such an open quantum system, the total Hamiltonian is separated into contributions of system, bath, and system–bath and system–field interactions,

\hat{H} = {\hat{H}}_{s} + {\hat{H}}_{b} + {\hat{H}}_{sb} + {\hat{H}}_{s - field}

⁠, which, in the given notation, are assumed to be represented in the exciton basis. In the basis of localized excited states (site basis), involving ground state |g⟩, single excitations |m⟩ localized at the mth monomer unit and double excitations |mn⟩ localized at monomer units m and n, the system Hamiltonian contains the site energies E_m, defined as the energy difference between |m⟩ and |g⟩ at the equilibrium positions of nuclei in |g⟩ as diagonal elements. E_m is composed of purely electronic contributions E_el,m and reorganization energies λ_m = ∑_iλ_m,i from vibrational modes i. The excitonic couplings J_mn determine the off-diagonal elements of the system Hamiltonian, which is given as

\begin{align} {\hat{H}}_{s, loc} & = \sum_{m} \sum_{n} (δ_{m n} (E_{el, m} + λ_{m}) + J_{m n}) | m 〉 〈 n | \\ + \sum_{k} \sum_{l} \sum_{m} \sum_{m} (δ_{k m} δ_{l n} (E_{el, m} + E_{el, n} + λ_{m} + λ_{n}) \\ + J_{k m} δ_{l n} + J_{k n} δ_{l m} + J_{l m} δ_{k n} + J_{ln} δ_{k m}) | k l 〉 〈 m n | . \end{align}

(1)

The bath Hamiltonian is separable from the basis of the system and depends on the frequencies ω_i and the position coordinates q_i of the involved bath oscillators in terms of

{\hat{H}}_{b} = \sum_{i} \frac{ω_{vib, i}^{2}}{2} (- \frac{\partial^{2}}{\partial q_{i}^{2}} + {\hat{q}}_{i}^{2}),

(2)

where ℏ has been set to 1. For excitation of monomer m, the harmonic potential of the bath modes exhibits a displacement

d_{i} = \frac{\sqrt{2 λ_{i}}}{ω_{vib, i}}

and we assume the same vibrational frequencies ω_vib,i as in the electronic ground state. The expansion of

\frac{ω_{vib, i}^{2}}{2} {({\hat{q}}_{i} + d_{i})}^{2}

leads to quadratic terms attributed to the bath Hamiltonian and the contribution of the reorganization energy to the system Hamiltonian, respectively. The mixed term enters the system–bath coupling Hamiltonian and accounts for the coupling of each localized state to the bath modes.^37–40 By introducing the gap coordinate operator

{\hat{u}}_{m} = \sum_{i} ω_{vib, i} \sqrt{2 λ_{m, i}} {\hat{q}}_{i},

(3)

the system–bath coupling Hamiltonian can be written as

{\hat{H}}_{sb,loc} = \sum_{m} {\hat{u}}_{m} | m 〉 〈 m | + \sum_{k} \sum_{l} ({\hat{u}}_{k} + {\hat{u}}_{l}) | k l 〉 〈 k l | .

(4)

A difficulty for theory is the equal magnitude of excitonic and exciton-vibrational couplings, where a measure for the latter is the local reorganization energy λ_m introduced above. Hence, there is no obvious small parameter that can be used for perturbation theory. However, such a parameter can be found by working in the exciton representation, in which the system Hamiltonian becomes diagonal upon transformation

{\hat{H}}_{s} = {\hat{T}}^{T} {\hat{H}}_{s,loc} \hat{T} = \sum_{α} ε_{α} | α 〉 〈 α |

⁠, with transformation matrix

\hat{T}

⁠. The system–bath coupling Hamiltonian in the exciton basis is obtained as

{\hat{H}}_{sb} = {\hat{T}}^{T} {\hat{H}}_{sb,loc} \hat{T} = \sum_{α, β} 〈 α | {\hat{T}}^{T} {\hat{H}}_{sb,loc} \hat{T} | β 〉 | α 〉 〈 β |

⁠. Often, the diagonal elements (α = β) of

{\hat{H}}_{sb}

are larger than the off-diagonal elements^33,40 due to partial localization of the excited states. In the present treatment, the diagonal elements will be taken into account exactly and a Markov and secular approximation will be used for the off-diagonal elements.³⁵

B. Time evolution and derivation of quantum master equation (QME)

For the derivation of the quantum master equation, we define the reduced statistical operator

\hat{ρ} = t r_{b} (\hat{W})

⁠, where the full density operator

\hat{W}

includes the probabilities of different components of a mixed state in the product basis of system and bath. For the statistical operator, we obtain the following closed differential equation:

\partial_{t} \hat{ρ} = - i t r_{b} ([\hat{H} (t), \hat{W}]) .

(5)

In the time evolution of the system, we take the interaction of the system with an external field into account, as described in the dipole approximation by the Hamiltonian

{\hat{H}}_{s - field} (t) = - \hat{\vec{μ}} \cdot \vec{E} (t)

involving the transformed transition dipole operator

\hat{\vec{μ}} = {\hat{T}}^{T} {\hat{\vec{μ}}}_{loc} \hat{T}

⁠, which contains transitions between ground state and singly excited states and between singly and doubly excited states with a representation in the localized basis as

\begin{align} {\hat{\vec{μ}}}_{loc} & = \sum_{m} ({\vec{μ}}_{m g} | m 〉 〈 g | + {\vec{μ}}_{g m} | g 〉 〈 m |) \\ + \sum_{m} \sum_{n} ({\vec{μ}}_{n g} | m n 〉 〈 m | + {\vec{μ}}_{m g} | m n 〉 〈 n | \\ + {\vec{μ}}_{g n} | m 〉 〈 m n | + {\vec{μ}}_{g m} | n 〉 〈 m n |) . \end{align}

(6)

Furthermore, the system–bath coupling Hamiltonian is represented in the interaction picture, with respect to both the system and the bath. Altogether, the time evolution of the components of the Hamiltonian is governed by

\begin{gathered} \partial_{t} {\hat{U}}_{s} = - i ({\hat{H}}_{s} + {\hat{H}}_{s - field} (t)) {\hat{U}}_{s}, \\ \partial_{t} {\hat{U}}_{b} = - i {\hat{H}}_{b} {\hat{U}}_{b}, \\ \partial_{t} {\hat{U}}_{s - b} = - i {\hat{U}}_{b}^{- 1} {\hat{U}}_{s}^{- 1} {\hat{H}}_{sb} {\hat{U}}_{s} {\hat{U}}_{b} {\hat{U}}_{s - b} . \end{gathered}

(7)

By representing also the full density operator in the interaction picture as

{\hat{W}}^{sb} = {\hat{U}}_{b}^{- 1} {\hat{U}}_{s}^{- 1} \hat{W} {\hat{U}}_{s} {\hat{U}}_{b}

⁠, one obtains

\hat{ρ} (t) = t r_{b} ({\hat{U}}_{s} {\hat{U}}_{b} {\hat{W}}^{sb} {\hat{U}}_{b}^{- 1} {\hat{U}}_{s}^{- 1}) = {\hat{U}}_{s} t r_{b} ({\hat{W}}^{sb}) {\hat{U}}_{s}^{- 1},

(8)

where the tracing leads to the reduced density matrix of the system

{\hat{ρ}}^{s} = t r_{b} ({\hat{W}}^{sb})

⁠. Time evolution with the full Hamiltonian

\hat{H} (t) = {\hat{H}}_{s - field} (t) + {\hat{H}}_{s} + {\hat{H}}_{b} + {\hat{H}}_{sb}

⁠, involving products of the time evolution operators specified in Eq. , leads to

\begin{gathered} i \partial_{t} {\hat{ρ}}^{s} (t) = t r_{b} ({\hat{U}}_{b}^{- 1} {\hat{U}}_{s}^{- 1} [{\hat{H}}_{sb}, \hat{W}] {\hat{U}}_{s} {\hat{U}}_{b}) \\ \partial_{t} {\hat{ρ}}^{s} (t) = - i t r_{b} ([{\hat{H}}_{sb}^{s - b} (t), {\hat{W}}^{sb} (t)]) . \end{gathered}

(9)

The quantum master equation (QME) is derived as described in Appendix A. By introducing the system part of the system–bath interaction of an excited monomer m in the exciton basis as

{\hat{K}}_{m} = {\hat{T}}^{T} | m 〉 〈 m | \hat{T}

and the correlation function of the bath components of the system–bath interaction of monomer units m and n as

C_{m n} (t) = t r_{b} ({\hat{u}}_{m} (t) {\hat{u}}_{n})

with the energy gap coordinate operator from Eq. , a reformulation of Eq. leads to an analogous expression as given in Ref. 41 for the QME in the interaction picture. A QME for the reduced density operator is obtained by applying an inverse transformation compared to the one which was specified in Eq. , resulting in

\begin{align} \partial_{t} \hat{ρ} (t) & = \partial_{t} ({\hat{U}}_{s}) {\hat{ρ}}^{s} (t) {\hat{U}}_{s}^{- 1} + {\hat{U}}_{s} {\hat{ρ}}^{s} (t) \partial_{t} ({\hat{U}}_{s}^{- 1}) + {\hat{U}}_{s} \partial_{t} ({\hat{ρ}}^{s} (t)) {\hat{U}}_{s}^{- 1} \\ = - i [{\hat{H}}_{s}, \hat{ρ} (t)] - i [{\hat{H}}_{s - field} (t), \hat{ρ} (t)] - \sum_{m} \sum_{n} \int_{t_{0}}^{t} d τ C_{n m} (t - τ) \\ \times ({\hat{K}}_{n} {\hat{U}}_{s} (t - τ) {\hat{K}}_{m} \hat{ρ} (τ) {\hat{U}}_{s}^{- 1} (t - τ) - {\hat{U}}_{s} (t - τ) {\hat{K}}_{m} \hat{ρ} (τ) \\ \times {\hat{U}}_{s}^{- 1} (t - τ) {\hat{K}}_{n}) + \sum_{m} \sum_{n} \int_{t_{0}}^{t} d τ C_{m n} (τ - t) \\ \times ({\hat{K}}_{n} {\hat{U}}_{s} (t - τ) \hat{ρ} (τ) {\hat{K}}_{m} {\hat{U}}_{s}^{- 1} (t - τ) - {\hat{U}}_{s} (t - τ) \hat{ρ} (τ) \\ \times {\hat{K}}_{m} {\hat{U}}_{s}^{- 1} (t - τ) {\hat{K}}_{n}), \end{align}

(10)

where

t r_{b} ({\hat{u}}_{m}) = 0

was used. The assumption

\hat{ρ} (τ) = {\hat{U}}_{s}^{- 1} (t - τ) \hat{ρ} (t) {\hat{U}}_{s} (t - τ)

implies that we neglect any dynamics beyond that of the isolated system between times t and τ. This assumption can be justified if the correlation function of the system–bath coupling, C_nm(t − τ), decays sufficiently fast, as assumed in the Markov approximation, leading to a standard Redfield description.⁴¹ Moreover, this approximation leads to an equivalence of our QME to the one which appears in the framework of the “partial ordering prescription” (POP) for the formulation of a second-order expression via cumulant expansion.^35,42 Note that the latter equivalence is rather fortuitous since POP provides an exact description of the diagonal part of the exciton-vibrational coupling in the exciton basis.³⁵ We furthermore expand the reduced density matrix in orders of the system–field interaction as

\hat{ρ} (t) = \sum_{n = 0}^{\infty} {\hat{ρ}}^{(n)} (t)

⁠, where

{\hat{ρ}}^{(n)} (t)

is in the nth “order” of

{\hat{H}}_{s - field}

⁠, thereby neglecting an interplay between dissipation and the external field. A separation of the different orders of the system–field interaction in the formulation of the QME leads to the requirement that the component multiplied with

{\hat{H}}_{s - field} (t)

is one order lower than the other components, which are all of the same order. We formulate the QME component-wise for each matrix element, thereby denoting the row and column indices as superscript and subscript indices, respectively. Using the notation

c_{m}^{α} = 〈 α | \hat{T} | m 〉

for the transformation coefficients and ω_αβ = ɛ_α − ɛ_β for the differences between eigenenergies of the exciton states, we obtain

\begin{align} \partial_{t} ρ_{α β}^{(n)} (t) & = - i ω_{α β} ρ_{α β}^{(n)} (t) - i {[{\hat{H}}_{s - field} (t), {\hat{ρ}}^{(n - 1)} (t)]}_{α β} \\ + \int_{0}^{t - t_{0}} A_{α β γ δ} (τ) d τ ρ_{γ δ}^{(n)} (t) . \end{align}

(11)

By introducing correlation functions in the eigenbasis of the system

C_{γ δ κ λ} (τ) = \sum_{m} \sum_{k} c_{k}^{γ} c_{k}^{δ} c_{m}^{κ} c_{m}^{λ} C_{m k} (τ)

and by taking into account the identities C_γδκλ(τ) = C_δγκλ(τ) = C_γδλκ(τ) = C_δγλκ(τ) for real-valued eigenenergies, the above A_αβγδ(τ) becomes

\begin{align} A_{α β γ δ} (τ) & = - δ_{δ β} \sum_{κ}^{κ \neq γ} C_{α κ κ γ} (τ) e^{- i ω_{κ γ} τ} - δ_{δ β} C_{α γ γ γ} (τ) \\ + C_{δ β α γ} (τ) e^{- i ω_{α γ} τ} + C_{δ β α γ} (- τ) e^{- i ω_{δ β} τ} \\ - δ_{α γ} \sum_{κ}^{κ \neq δ} C_{δ κ κ β} (- τ) e^{- i ω_{δ κ} τ} - δ_{α γ} C_{δ δ δ β} (- τ) . \end{align}

(12)

C. Formulation of population and coherence dynamics

For the element-wise formulation of the QME given in Eq. , an analytical solution can be obtained by applying the “variation of constants” approach, where the expressions containing the nth contribution from the series expansion of the reduced density matrix are attributed to the homogeneous differential equation, whereas the terms containing the contribution of the order n − 1 are treated as an inhomogeneity. In the following, we will apply the secular approximation,⁴¹ which consists in a separation of coherence contributions with α ≠ β, α = γ, and β = δ and population contributions with α = β and γ = δ. Depending on whether a coherence or a population is considered, different solutions are obtained. In the case of a coherence, taking into account the initial condition

ρ_{α β}^{(n)} (t_{0}) = 0

⁠,

ρ_{α β}^{(n)} (t) = - i \int_{t_{0}}^{t} e^{- i Ω_{α β} (t, T)} {[{\hat{H}}_{s - field} (T), {\hat{ρ}}^{(n - 1)} (T)]}_{α β} d T

(13)

is obtained, where we introduced

- i Ω_{α β} (t, T) = - i ω_{α β} (t - T) + \int_{T}^{t} \int_{0}^{t^{'} - t_{0}} A_{α β α β} (τ) d τ d t^{'}

⁠.

In the case of populations, we apply a Markov approximation in the sense of setting the upper integration border in Eq. to infinity,³⁵ which leads to

\begin{align} K_{α β} & = - δ_{α β} \sum_{γ}^{γ \neq β} \int_{0}^{\infty} d τ C_{β γ γ β} (τ) e^{- i ω_{γ β} τ} \\ - δ_{α β} \sum_{γ}^{γ \neq β} \int_{0}^{\infty} d τ C_{β γ γ β} (- τ) e^{- i ω_{β γ} τ} \\ + (1 - δ_{α β}) \int_{0}^{\infty} d τ C_{β α α β} (τ) e^{- i ω_{β α} τ} \\ + (1 - δ_{α β}) \int_{0}^{\infty} d τ C_{β α α β} (- τ) e^{- i ω_{α β} τ} . \end{align}

(14)

The rate constant of exciton relaxation between |α⟩ and |β⟩ is defined as

k_{α \to β} = \int_{- \infty}^{\infty} C_{α β β α} (τ) e^{- i ω_{β α} τ} d τ

⁠. In terms of the latter, K_αβ reads

K_{α β} = - δ_{α β} \sum_{γ}^{γ \neq β} k_{γ \to β} + (1 - δ_{α β}) k_{β \to α} .

(15)

By calculating the matrix exponential

e^{\hat{K} t^{'}}

⁠, the population transfer tensor is obtained. The final solution for the dynamics in the case of populations is

ρ_{α α}^{(n)} (t) = - i \int_{t_{0}}^{t} {[e^{\hat{K} (t - T)}]}_{α β} {[{\hat{H}}_{s - field} (T), {\hat{ρ}}^{(n - 1)} (T)]}_{β β} d T .

(16)

D. Third-order polarization and 2D spectrum

For a description of nonlinear spectroscopic techniques relying on the third-order polarization, three system–field interactions are required. Starting from an initial population of the electronic ground state, i.e.,

ρ_{g g}^{(0)} (t_{0}) = 1

⁠, the specific form of the transition dipole operator that facilitates transitions between ground state and singly excited states and between singly and doubly excited states leads to non-vanishing coherences resulting from Eq. , where the multiplication of the reduced density operator with the transition dipole operator from the left- or right-hand side in the commutator expression yields a coherence matrix element

ρ_{α g}^{(1)}

or

ρ_{g α}^{(1)}

⁠, respectively. As the commutator in the first-order expression does not result in any diagonal elements, all of the populations calculated according to Eq. are zero. However, in second order, a ground state population

ρ_{g g}^{(2)}

can be obtained for a multiplication of the transition dipole operator with

ρ_{α g}^{(1)}

or

ρ_{g α}^{(1)}

from the left- or right-hand side, respectively. A further possibility in this case is the generation of a coherence between a doubly excited state and the ground state, namely

ρ_{2 α g}^{(2)}

or

ρ_{g 2 α}^{(2)}

⁠. A population

ρ_{α α}^{(2)}

or a coherence

ρ_{α β}^{(2)}

or

ρ_{β α}^{(2)}

in the excited state subspace is obtained for a multiplication of the transition dipole operator with

ρ_{α g}^{(1)}

or

ρ_{g α}^{(1)}

from the right- or left-hand side, respectively. Likewise, the third interaction facilitates a transition with an increase or decrease in the excitation in the system, depending on the selection of the product expression between the transition dipole operator and the reduced density operator in the commutator expression in Eq. . Again, as in the case of the first interaction, only coherences can be generated. The third-order polarization is obtained by taking the trace with respect to the system over the product of the dipole operator with the third-order reduced density operator, as expressed by

\begin{align} {\vec{P}}_{mol}^{(3)} & = t r_{s} (\hat{\vec{μ}} {\hat{ρ}}^{(3)}) \\ = \sum_{α} ({\vec{μ}}_{g α} ρ_{α g}^{(3)} + {\vec{μ}}_{α g} ρ_{g α}^{(3)}) + \sum_{α} \sum_{2 β} ({\vec{μ}}_{α 2 β} ρ_{2 β α}^{(3)} + {\vec{μ}}_{2 β α} ρ_{α 2 β}^{(3)}) . \end{align}

(17)

We assume that the electric field, which interacts with the system, consists of non-overlapping pulses with central time T_i of the form

{\vec{E}}_{i} (x, t) = {\vec{e}}_{i} Z_{i} (x, t - T_{i}) .

(18)

Furthermore, we include an orientational average of the polarizations for different arrangements of the system and, hence, the transition dipole moments of the system relative to the field polarizations, resulting in

P (x, t) = {⟨{\vec{e}}_{4} \cdot {\vec{P}}_{mol}^{(3)} (x, t)⟩}_{orient}

⁠, which, according to our considerations above, can be written as

\begin{align} {⟨{\vec{e}}_{4} \cdot {\vec{P}}_{mol}⟩}_{orient} & = \int_{T_{3}}^{t} \int_{T_{2}}^{τ_{3}^{'}} \int_{T_{1}}^{τ_{2}^{'}} d τ_{1}^{'} d τ_{2}^{'} d τ_{3}^{'} S (t - τ_{3}^{'}, τ_{3}^{'} - τ_{2}^{'}, τ_{2}^{'} - τ_{1}^{'}) \\ \times Z_{3} (x, τ_{3}^{'} - T_{3}) Z_{2} (x, τ_{2}^{'} - T_{2}) Z_{1} (x, τ_{1}^{'} - T_{1}) . \end{align}

(19)

The involved third-order response

\begin{align} S (t - τ_{3}^{'}, τ_{3}^{'} - τ_{2}^{'}, τ_{2}^{'} - τ_{1}^{'}) & = i \sum_{n = 1}^{4} (R_{n} (t - τ_{3}^{'}, τ_{3}^{'} - τ_{2}^{'}, τ_{2}^{'} - τ_{1}^{'}) \\ - R_{n}^{*} (t - τ_{3}^{'}, τ_{3}^{'} - τ_{2}^{'}, τ_{2}^{'} - τ_{1}^{'})) \end{align}

(20)

consists of a sum of response functions, where each of them corresponds to a specific excitation pathway, depending on the selection of a combination of terms from the commutator expressions resulting from the system–field interaction. The respective response functions involve an orientational averaging, which will be specified later. With the substitution

\begin{align} τ (τ^{'}) & = \{\begin{cases} τ_{1} = τ_{2}^{'} - τ_{1}^{'}, \\ τ_{2} = τ_{3}^{'} - τ_{2}^{'}, \\ τ_{3} = t - τ_{3}^{'}, \end{cases} τ^{'} (τ) = \{\begin{cases} τ_{1}^{'} = t - τ_{3} - τ_{2} - τ_{1}, \\ τ_{2}^{'} = t - τ_{3} - τ_{2}, \\ τ_{3}^{'} = t - τ_{3}, \end{cases} \\ \det (\frac{\partial τ}{\partial τ^{'}}) = - 1, \end{align}

(21)

one arrives at

\begin{align} {⟨{\vec{e}}_{4} \cdot {\vec{P}}_{mol}⟩}_{orient} & = \int_{0}^{t - T_{3}} \int_{0}^{t - τ_{3} - T_{2}} \int_{0}^{t - τ_{3} - τ_{2} - T_{1}} d τ_{1} d τ_{2} d τ_{3} \\ \times S (τ_{3}, τ_{2}, τ_{1}) Z_{3} (x, t - τ_{3} - T_{3}) \\ \times Z_{2} (x, t - τ_{3} - τ_{2} - T_{2}) \\ \times Z_{1} (x, t - τ_{3} - τ_{2} - τ_{1} - T_{1}) . \end{align}

(22)

Because of the assumption of finite pulse widths, we may extend the upper limit of integration formally to ∞. Furthermore, we formally extend the lower limit of integration to −∞ by using the Heaviside theta function,

\begin{aligned} Θ (τ) = \{\begin{cases} 1 & for τ \geq 0, \\ 0, & else . \end{cases} \end{aligned}

(23)

In this way, we obtain

\begin{align} P (x, t) & = \int_{- \infty}^{\infty} \int_{- \infty}^{\infty} \int_{- \infty}^{\infty} d τ_{1} d τ_{2} d τ_{3} Θ (τ_{1}) Θ (τ_{2}) Θ (τ_{3}) \\ \times S (τ_{1}, τ_{2}, τ_{3}) Z_{3} (x, t - τ_{3} - T_{3}) \\ \times Z_{2} (x, t - τ_{3} - τ_{2} - T_{2}) Z_{1} (x, t - τ_{3} - τ_{2} - τ_{1} - T_{1}) . \end{align}

(24)

By assuming δ-shaped pulses with

Z_{j} (x, t^{'}) = δ (t^{'}) e^{i k_{j} x - i ω_{j} t^{'}} + c . c .

⁠, it turns out that the third-order polarization can be expressed by the response function with dependence on differences between the central times of the pulses T₁, T₂, and T₃ and the time variable t in terms of S(T₂ − T₁, T₃ − T₂, t − T₃). Furthermore, one finds that in the framework of the “rotating wave approximation” (RWA),⁴¹ it is required to select a contribution with appropriate sign in the complex exponential to compensate the oscillations of the pulses. Please note that the respective contributions depend on the considered response function. Furthermore, such a compensation is only possible for specific pathways. The selection of the complex exponentials from the pulses also determines the direction in which the signal is detected. This detection direction depends on the type of response function and results as −k₁ + k₂ + k₃ for rephasing and as k₁ − k₂ + k₃ for non-rephasing contributions.⁴³ Furthermore, there are double coherence contributions with detection direction k₁ + k₂ − k₃.⁴⁴ In the impulsive limit (δ-shaped pulses), pulse overlap effects are not relevant, but for finite pulse width, they can be taken into account as described in the supplementary material. For further evaluation, we use the definition of line shape functions,⁴⁵

g_{α β} (t) = \int_{0}^{t} \int_{0}^{\bar{t}} C_{α α β β} (τ) d τ d \bar{t},

(25)

where the complex conjugate can be determined by using the property of the correlation function

C_{α α β β}^{*} (τ) = C_{α α β β} (- τ)

⁠, and apply it to the tensor elements A_αβαβ, α ≠ β, accounting for coherences, thereby selecting correlation function contributions from Eq. where sequential pairs of equal indices appear in the subscript index pattern. These contributions are treated in a non-Markovian way, as indicated by Eq. , because of the involvement of diagonal elements of the system–bath coupling Hamiltonian in the respective correlation functions.³⁵ The remaining terms with a different index patterns lead to lifetime broadening and are specified in the supplementary material. For the treatment of populations, we apply Eqs. and . We furthermore replace the time arguments T₂ − T₁, T₃ − T₂, and t − T₃ in S(T₂ − T₁, T₃ − T₂, t − T₃) by τ, T, and t, respectively, and express the third-order polarization in terms of response function contributions, which remain after the application of the RWA, thereby indicating whether, in a selected pathway, besides the singly excited states, only the ground state (denoted as g) or also the doubly excited state (denoted as f) is involved. Furthermore, we distinguish between contributions involving either populations (indicated by subscript pop) or coherences (indicated by subscript coh) in the excited state. The contributions to the third-order polarization resulting from the different response functions lead to a reproduction of the respective response functions with modified time arguments, which will be formulated in the following, thereby denoting a factor from orientational averaging in the general form

{⟨({\vec{e}}_{4} \cdot {\vec{μ}}_{b}^{a}) ({\vec{e}}_{3} \cdot {\vec{μ}}_{d}^{c}) ({\vec{e}}_{2} \cdot {\vec{μ}}_{f}^{e}) ({\vec{e}}_{1} \cdot {\vec{μ}}_{h}^{g})⟩}_{orient}

⁠.

Depending on the excitation pathways, the response functions are associated with ground state bleaching (GSB), stimulated emission (SE), excited state absorption (ESA), and double coherence (DC) processes. The detailed form of the response functions is given in the supplementary material. We will reformulate the response functions after introducing some further definitions.

Note that the third-order polarization is a real quantity, as it also includes the complex conjugate counterparts of the specified response function expressions, which leads to a counter-propagating signal field. Furthermore, their contributions appear at negative-signed frequency positions when a two-dimensional Fourier transform is determined according to

\begin{align} σ_{2 D,R/NR} (ω_{τ}, T, ω_{t}) & = \frac{1}{{(2 π)}^{2}} \int_{- \infty}^{\infty} \int_{- \infty}^{\infty} d τ d t \\ \times i P_{R / N R}^{(3)} (τ, T, t) e^{\mp i ω_{τ} τ} e^{i ω_{t} t}, \end{align}

(26)

where the prefactor i accounts for generation of the signal field by the third-order polarization⁴⁶ and the sign in the Fourier transformation with respect to t₁ depends on whether a rephasing (R_2g,

R_{1 f}^{*}

⁠, R_3g) or non-rephasing (R_1g,

R_{2 f}^{*}

⁠, R_4g) contribution is considered, as indicated by the subscript R/NR. Double coherence contributions are not taken into account, as they only play a role for overlapping pulses not considered here. Between the response functions from the supplementary material and the corresponding expressions in Ref. 36, a connection can be drawn in view of the dependencies on frequency differences between electronic states, on decay constants, line shape functions, and population transfer tensor. However, the formulated response functions are only comparable with those from Ref. 36 if line shape functions with dependencies on different time arguments are disregarded in the latter ones. This finding reflects that in the present derivation of the response functions, a tracing over the bath was implicitly included in the quantum master equation for propagation over each time interval, equivalent to a thermalization before the next electronic transition takes place, whereas in Ref. 36, only a single tracing over the bath, comprising all time intervals for a given response function, was applied in the framework of an evaluation via second-order cumulant expansion. A comparison of the line shape function contributions with dependencies on single time arguments leads to the finding that the line shape functions depending on t and with both indices referring to a singly excited state in our formulation and in Ref. 36 are complex conjugate to each other. This finding indicates that the present formulation refers to an absorption process and the earlier refers to an emission process. As will be discussed later, the present expressions need to be corrected in order to take into account the nuclear relaxation in the excited states.

E. Treatment of dissipative dynamics relying on spectral density

In practice, the correlation function, from which the line shape functions and the transfer rates are obtained, is calculated from a spectral density J(ω). There exist different functional forms⁴⁷ and analytical treatments⁴⁸ of J(ω). For our purpose, the following type, originally extracted from fluorescence line narrowing spectra of B777 complexes,³⁵ was found to be suitable: J(ω) = SJ₀(ω) with the Huang–Rhys factor S and the normalized spectral density

J_{0} (ω) = \frac{1}{s_{1} + s_{2}} \sum_{i = 1, 2} π \frac{s_{i}}{7! 2 ω_{i}^{4}} ω^{5} \exp (- {(\frac{ω}{ω_{i}})}^{\frac{1}{2}}),

(27)

with the original parameters³⁵^, s₁ = 0.5025, s₂ = 0.4975, ω₁ = 0.557 cm⁻¹, and ω₂ = 1.94 cm⁻¹. Here, we use a Huang–Rhys factor S = 0.8, which has been adopted for WSCP from the temperature dependence of the linear absorption spectrum.¹⁹ Note also that, different from the earlier definition,³⁵ a factor of πω² of J₀(ω) is implicitly included in the present definition of the spectral density. Nevertheless, the description with the present formulation of the spectral density is equivalent to the one from Ref. 35, as the definitions of correlation functions and line shape functions depending on the spectral density are modified accordingly, resulting in formulations that are more common in the literature of 2D spectroscopy than those from Ref. 35. Further details, including a discussion of line narrowing spectra of WSCP^20,49 and alternative spectral densities, are given in the supplementary material.

The correlation function, which has been introduced in Sec. II B, is obtained from the spectral density as

C (t) = \frac{1}{π} \int_{0}^{\infty} d ω J (ω) (\cos (ω t) \coth (\frac{ω}{2 k_{B} T}) - i \sin (ω t)) .

(28)

Different from the definition in Eq. , in Ref. 35, the line shape function (yet without indices referring to the involved exciton states) has been introduced as

G (t) = \frac{1}{π} \int_{0}^{\infty} d ω \frac{J (ω)}{ω^{2}} (\cos (ω t) \coth (\frac{ω}{2 k_{B} T}) - i \sin (ω t)) .

(29)

Using a different formulation of the previously introduced line shape function⁵⁰ as

\begin{align} g (t) & = \frac{1}{π} \int_{0}^{\infty} d ω \frac{1 - \cos (ω t)}{ω^{2}} \coth (\frac{ω}{2 k_{B} T}) J (ω) \\ + \frac{i}{π} \int_{0}^{\infty} d ω \frac{\sin (ω t) - ω t}{ω^{2}} J (ω), \end{align}

(30)

one can easily convince oneself that the real part

R

of the expression from Eq. corresponds to the difference between

R (G (0))

and

R (G (t))

⁠. The imaginary part from Eq. corresponds to the difference between −iλt and

I (G (t))

⁠, where the reorganization energy λ can be calculated as

λ = \frac{1}{π} \int_{0}^{\infty} d ω \frac{1}{ω} J (ω) .

(31)

Altogether, one obtains

g (t) = G (0) - G (t) - i λ t .

(32)

The indices of the line shape functions can be attributed to products of transformation coefficients,³⁵^,

γ_{α β γ δ} = \sum_{m} \sum_{n} c_{m}^{(α)} c_{m}^{(β)} c_{n}^{(γ)} c_{n}^{(δ)} δ_{m n},

(33)

where δ_mn accounts for the assumption of uncorrelated fluctuations⁴⁰ and leads to non-vanishing line shape functions,

g_{α γ} (t) = γ_{α α γ γ} g (t) .

(34)

For the calculation of the rate constants,³⁵ which have been introduced in the context of Eq. , we identify the frequency-dependent correlation function as

\tilde{C} (ω) = \int_{0}^{\infty} d t \exp (i ω t) C (t) .

(35)

The real part of this function is obtained from the spectral density and from the Bose–Einstein distribution

n (ω) = \frac{1}{\exp (- \frac{ω}{k T}) - 1}

(36)

as

R (\tilde{C} (ω)) = (1 + n (ω)) J (ω) + n (- ω) J (- ω) .

(37)

The rate constant of exciton relaxation from state |α⟩ to state |β⟩ is given as

k_{α β} = 2 γ_{α β β α} R (\tilde{C} (ω_{α β})) .

(38)

The imaginary part of the complex-valued coherence dephasing rates specified in the supplementary material corresponds to a frequency shift, which, in combination with a shift due to the reorganization energy, leads to a modified transition frequency,^35,51

{\tilde{ω}}_{α_{g}} = ω_{α g} - γ_{α α α α} λ + \sum_{α \neq β} γ_{α β β α} I (\tilde{C} (ω_{α β})) .

(39)

The response functions given in the supplementary material can be reformulated using this modified definition of the frequencies. As the shifts related to the reorganization energy are found to be equivalent to the corresponding contributions entering the line shape functions according to Eqs. and , the replacement of line shape functions g(t) by G(0) − G(t) is required in such a reformulation. In the case of a dimer, the rate matrix containing the population transfer rates in the basis of the electronic states, including the ground state |g⟩, the exciton states |1⟩ and |2⟩, and the doubly excited state |d⟩, is given as

\hat{K} = (\begin{matrix} 0 & 0 & 0 & 0 \\ 0 & - k_{1 \to 2} & k_{2 \to 1} & 0 \\ 0 & k_{1 \to 2} & - k_{2 \to 1} & 0 \\ 0 & 0 & 0 & 0 \end{matrix}),

(40)

where |g⟩ and |d⟩ are not involved in exciton relaxation. The matrix exponential can be calculated analytically by transformation to a diagonal representation, where the matrix exponential is obtained by taking the exponentials of the diagonal elements, and subsequent back-transformation. In this way, one obtains

\begin{aligned} e^{\hat{K} t} = (\begin{matrix} 1 & 0 & 0 & 0 \\ 0 & \frac{k_{1 \to 2} f (t) + k_{2 \to 1}}{k_{1 \to 2} + k_{2 \to 1}} & - \frac{k_{2 \to 1} (f (t) - 1)}{k_{1 \to 2} + k_{2 \to 1}} & 0 \\ 0 & - \frac{k_{1 \to 2} (f (t) - 1)}{k_{1 \to 2} + k_{2 \to 1}} & \frac{k_{2 \to 1} f (t) + k_{1 \to 2}}{k_{1 \to 2} + k_{2 \to 1}} & 0 \\ 0 & 0 & 0 & 1 \end{matrix}), \end{aligned}

(41)

with

f (t) = e^{- (k_{1 \to 2} + k_{2 \to 1}) t}

⁠. Note that for a given initial population

\vec{P} (0)

⁠,

e^{\hat{K} t} \vec{P} (0)

will reveal the environment-induced population transfer leading to

\vec{P} (t)

⁠. The response functions are reformulated, thereby using the notation introduced in this section, in Appendix B. The appearance of only single time arguments of the line shape functions entering the derived response functions facilitates a separation of the respective expressions in terms of products containing absorption and emission spectra of the form⁵²^,

σ_{a b s / e m} = R \{\int_{0}^{\infty} d t \exp (i ω t) \sum_{α} | {\vec{μ}}_{α g} |^{2} C_{a b s / e m, α} (t)\} .

(42)

The correlation functions of absorption and emission of an exciton state |α⟩ are of the form

\begin{align} C_{abs, α} (t) & = \exp (- i {\tilde{ω}}_{α g} t) \exp (G_{α α} (t) - G_{α α} (0)) \\ \times \exp (- \frac{1}{2} \sum_{β}^{β \neq α} k_{α \to β} t) \end{align}

(43)

and

\begin{align} C_{em, α} (t) & = \exp (- i {\tilde{ω}}_{α g} t) \exp (G_{α α}^{*} (t) - G_{α α}^{*} (0)) \\ \times \exp (- \frac{1}{2} \sum_{β}^{β \neq α} k_{β \to α} t), \end{align}

(44)

respectively. Likewise, by introducing the magnetic transition dipole moment of state n in the site basis as

{\vec{m}}_{loc, n g} = - {\vec{R}}_{n} \times {\vec{μ}}_{loc, n g}

and by transforming it to the exciton basis, the circular dichroism (CD) spectrum can be expressed as⁵²

σ_{CD} = I \{\int_{0}^{\infty} d t \exp (i ω t) \sum_{α} {\vec{μ}}_{α g} \cdot {\vec{m}}_{α g} C_{abs, α} (t)\} .

(45)

Between the response functions in reformulated notation and the correlation functions of absorption and emission in Eqs. and , a connection can be drawn, which remains valid after calculating the corresponding spectra, provided that the absorption and emission spectra are considered complex-valued quantities instead of taking their real part as in Eq. . For an appropriate assignment of absorption- or emission-like character of the factors in the response functions with dependencies on τ and t, also consideration of the sign in the Fourier transformation according to Eq. is required, which depends on whether a response function is rephasing or non-rephasing. By drawing such connections, it turns out that the transitions that give rise to the dynamics during both τ and t can be either related to the complex-valued absorption spectrum or its complex conjugate counterpart, even though the final transition of each pathway is expected to correspond to an emission process. The absorption-like property of the transition associated with t stems from separate tracing over the unshifted bath oscillators in the QME for the evolution during each time interval, implying the assumption that at all stages of the excitation pathway, the system remains equilibrated with respect to the electronic ground state. This rather questionable assumption can be corrected by introducing thermal equilibration in excited-state populations in a phenomenological way, thereby applying a shift of the bath coordinate to the equilibrium position of an exciton state α once this exciton state is populated in a selected excitation pathway.³⁵ This situation is illustrated in Fig. 1. The influence of such a modification on the response functions is discussed in Appendix C, thereby addressing resulting modifications in the derivations described above. The influence of the assumption of thermal equilibration in the electronic ground state or in an exciton state is illustrated in the supplementary material by model calculations on Chl a WSCP (see Sec. III) for the specific case of response function R₂. These calculations reveal a slight (Stokes) shift of the peaks with respect to ω_t.

FIG. 1.

View large Download slide

Illustration of vibrational relaxation in exciton state prior to exciton relaxation, stimulated emission (SE), and excited state absorption (ESA). After optical excitation and exciton relaxation, a fast vibrational relaxation is assumed such that stimulated emission (SE) and excited state absorption (ESA) start from a vibrationally equilibrated excited state.

Altogether, in the response functions, absorption- and emission-type spectra associated with the evolution during τ and t, respectively, are connected to each other via T-dependent functions, which, in the case of excited-state populations, correspond to the transfer dynamics and, in the case of coherences, reflect the dephasing in terms of decaying oscillations. The response functions involving correction of thermal equilibration in the excited state are formulated in Eqs. (B1)–(B10). They also contain a factor to account for orientational averaging, which is further specified in Appendix D.

III. RESULTS

In the following, we will use the approach of calculating 2D spectra described in Secs. II D, II E, and Appendix B, thereby using the response functions given in Eqs. (B1)–(B10), to simulate measured 2D spectra²⁹ of WSCP of the plant species Brassica oleracea that were reconstituted with different variants of chlorophylls—either Chl a or Chl b, which differ in a single side chain of the porphyrin skeletal structure. A methyl group in the case of Chl a is replaced by a formyl group in Chl b. Four chlorophyll units are arranged in WSCP in such a way that two pairs of chlorophylls form dimers with a relatively small (center-to-center) distance of 10 Å between the monomer units, whereas the distance between the dimers is considerably larger (about 20 Å), as shown in Fig. 2. Therefore, the inter-dimer excitonic coupling is weak and substantial delocalization of exciton states only occurs in the dimers so that it is sufficient to rely on the model of a dimer instead of a tetramer, particularly in the calculation of linear spectra. A treatment of transfer processes due to inter-dimer coupling in the framework of generalized Förster theory⁵³ goes beyond the scope of the present article. For all calculations, we assume a temperature of 77 K, as in the experiment.²⁹ Inhomogeneous broadening is included by varying the site energies of the chlorophyll monomers according to a Gaussian distribution with a FWHM of 250 cm⁻¹. Please note that we had to increase the inhomogeneous width with respect to our earlier work (170 cm⁻¹)²⁰ in order to fit the optical spectra from Ref. 29, indicating differences in the sample preparation.

FIG. 2.

View large Download slide

Arrangement of chlorophyll pigments in WSCP.^12,15 Graphics prepared with VMD.

A. Chl a WSCP

In the case of Chl a, the monomer units of the two dimers in WSCP have the same site energies, as long as inhomogeneous broadening is not involved. The excitonic coupling between the Chl a pigments in each dimer is J = 83 cm⁻¹.²⁰ A comparison of the measured and calculated linear absorption and CD spectra is shown in Fig. 3. The low-energy region of the experimental spectra is reproduced quantitatively by the calculations.

FIG. 3.

View large Download slide

Comparison of experimental linear absorption spectrum of Chl a WSCP at 77 K and of experimental CD spectrum of Chl a WSCP at 300 K from Ref. 29 with results from the corresponding calculations (upper and lower panels, respectively).

The deviations between the calculated and measured absorption spectrum at higher energies are due to contributions from high-frequency intramolecular vibrations and the next higher electronic transition (Q_x) not taken into account in the present calculations. A quantitative description of this spectral region has recently been provided taking into account the non-adiabatic coupling between the Q_y and Q_x transitions.⁵² The deviations between the calculated and measured CD spectra in the high-frequency region in Fig. 3 are much less than in the absorption spectra due to the excitonic nature of the CD spectra and the small excitonic couplings involving these high-energy transitions.⁵² The measured and calculated rephasing contributions to the 2D spectra of Chl a WSCP complexes are shown in the upper and lower rows of Fig. 4, respectively. In the present case, the rephasing contributions contain all relevant features that are also found in the complete 2D spectra and in the non-rephasing contributions shown in the supplementary material (Figs. S2 and S3). Comparing the 2D spectra from experiment and simulation, one finds that substantial features are reproduced. In particular, the appearance of the main diagonal peak 1 at (ω_τ = 14 850 cm⁻¹, ω_t = 14 850 cm⁻¹) and its decay with increasing population time T is similar. At T = 20 fs, additionally, a negative-signed peak 2 at (ω_τ = 14 950 cm⁻¹, ω_t = 14 675 cm⁻¹) is observed, which also decays and has lost most of its intensity at larger values of T from 100 fs onward. The strong upper diagonal peak 1 reflects the large dipole strength of the upper exciton state, giving rise to GSB and SE contributions. As long as the upper exciton state is populated, there is ESA between this state and the two-exciton state appearing as a negative peak 2. Note that the GSB contribution of the upper exciton state survives exciton relaxation. Therefore, the intensity of peak 1 does not decay to zero for long T-times. By selecting a single point of each of the mentioned peaks, we can study the T-dependent decay, as displayed in Fig. 5, revealing good agreement between theory and experiment.²⁹ The deviations at small times could be due to overlap effects of the pulses that were not included in the calculations. A rough approximation of these overlap effects reveals an increase by a factor of 2 in the signal at T = 0 fs, as discussed in the supplementary material. Furthermore, the population dynamics reflected by the time dependence of the peak amplitudes could involve effects of higher order in the system–bath coupling, which are not captured by a second-order treatment. The influence of different orders in a perturbative treatment can, in principle, be extracted from the non-perturbative “Hierarchical Equations of Motion” (HEOM) method.^54,55 The deviations of the calculated results from the measured ones at long times in Fig. 5 indicate slight differences in relative contributions of ESA, GSB, and SE between theory and experiment at this particular point in the 2D spectrum. From the obtained curves, the decay time can be estimated to be about 60 fs, which is similar to the value inferred in Ref. 29 (∼100 fs). A calculated lifetime of 50–80 fs has been reported in our earlier work,¹⁹ and a shortest lifetime of 50 fs was found²⁴ in earlier 2D spectra on native WSCP containing a mixture of Chl a/Chl a and Chl b/Chl b homodimers and Chl a/Chl b heterodimers.

FIG. 4.

Experimental29 and calculated rephasing contribution to the 2D spectra of Chl a WSCP at 77 K at different population times, as indicated. The circle and square symbols indicate selected points of the upper diagonal peak (peak 1) and the cross-peak (peak 2), respectively, at which the dependence of the peak amplitude on the population time is determined. The coordinates of these points are (ωτ = 14 850 cm−1, ωt = 14 850 cm−1) (peak 1) and (ωτ = 14 950 cm−1, ωt = 14 675 cm−1) (peak 2). Experimental data of Ref. 29 were provided by Fresch et al.

View large Download slide

Experimental²⁹ and calculated rephasing contribution to the 2D spectra of Chl a WSCP at 77 K at different population times, as indicated. The circle and square symbols indicate selected points of the upper diagonal peak (peak 1) and the cross-peak (peak 2), respectively, at which the dependence of the peak amplitude on the population time is determined. The coordinates of these points are (ω_τ = 14 850 cm⁻¹, ω_t = 14 850 cm⁻¹) (peak 1) and (ω_τ = 14 950 cm⁻¹, ω_t = 14 675 cm⁻¹) (peak 2). Experimental data of Ref. 29 were provided by Fresch et al.

FIG. 5.

View large Download slide

Dependence of the amplitudes of the two main peaks of the measured and calculated rephasing contribution to the 2D spectra in Fig. 4 on the population time at the specified positions of peak 1 (black line: experiment;²⁹ green line: calculation) and peak 2 (red line: experiment; blue line: calculation). Experimental data of Ref. 29 were provided by Fresch et al.

FIG. 6.

View large Download slide

Squared linear absorption spectrum (red line) and diagonal cut through the 2D spectrum at T = 0 fs from Fig. 4 (black line) are compared for both calculated and measured²⁹ results. Experimental data of Ref. 29 were provided by Fresch et al.

In the sum of rephasing and non-rephasing contributions, Fourier transformation with respect to τ and t leads to an absorptive line shape⁵⁶ (in contrast to a dissipative one of the separate contributions) with an absorption-like spectral profile with respect to ω_τ and either an absorption- or emission-like spectral profile with respect to ω_t, depending on the selected pathway and on whether populations of the excited state are assumed as thermally equilibrated. In a diagonal cut through the 2D spectrum at T = 0 fs, the relative peak intensities are expected to be comparable to the square of the absorption spectrum. Further theoretical details of this comparison can be found in the supplementary material. As shown in Fig. 6, the respective curves, in comparison of both experimental and calculated results, are very similar, indeed.

In Fig. 7, the GSB-, SE- and ESA-type response functions with pathways resulting in rephasing contributions are displayed, which, in the present case, facilitate an interpretation of the features in the 2D spectrum in Fig. 4, which is a sum of these contributions. The corresponding energy diagrams including an illustration of the effect of exciton relaxation are shown in Fig. 8 in addition to the corresponding Feynman diagrams. The finding that the positive-signed diagonal peak in the sum spectrum in Fig. 4 still appears when the negative-signed cross-peak below has already vanished can be explained by the influence of the GSB contribution, which does not decay during T, as displayed in the first row of Fig. 7. This non-decaying appearance is related to the respective pathway, where the second transition is a de-excitation back to the electronic ground state, so that a constant ground-state population appears during T. Accordingly, the GSB process is independent of population transfer in the excited state, as illustrated in the first row of Fig. 8. While the GSB contribution can be considered as a constant background, SE and ESA contributions yield insight into the excited-state dynamics. In the second row of Fig. 7, the SE contribution R₂ is shown for different population times T. At early population time of 20 fs, its appearance is dominated by the diagonal peak at (ω_τ = 14 820 cm⁻¹, ω_t = 14 820 cm⁻¹), which results from excitation of the energetically higher exciton state. With increasing delay time, this peak decreases due to population decay, and simultaneously, a cross-peak at (ω_τ = 14 820 cm⁻¹, ω_t = 14 600 cm⁻¹) gains intensity because of population transfer to the energetically lower exciton state, as illustrated in the second row of Fig. 8. Since this low-energy exciton state has a small transition dipole strength, at large delay times, an additional diagonal peak at (ω_τ = 14 600 cm⁻¹, ω_t = 14 600 cm⁻¹) becomes visible that reflects initial excitation of the low-energy exciton state, from where also the actual SE process with de-excitation back to the electronic ground state takes place. Please note that this peak is also present at earlier times T but not visible because of the dominant diagonal peak at higher energies.

FIG. 7.

Selected rephasing response function contributions to the 2D spectrum of Chl a WSCP at 77 K, namely R3, R2, and R1* (first, second, and third rows, respectively) are shown for different population times T = 20 fs (left panel), T = 100 fs (middle panel), and T = 1 ps (right panel). The colors of the inserted points correspond to those of the arrows indicating the underlying transitions in Fig. 8.

View large Download slide

Selected rephasing response function contributions to the 2D spectrum of Chl a WSCP at 77 K, namely R₃, R₂, and $R_{1}^{*}$ (first, second, and third rows, respectively) are shown for different population times T = 20 fs (left panel), T = 100 fs (middle panel), and T = 1 ps (right panel). The colors of the inserted points correspond to those of the arrows indicating the underlying transitions in Fig. 8.

FIG. 8.

View large Download slide

Energy diagrams for different population times (left and middle columns) and Feynman diagrams (right column) for the rephasing response functions R₃ (first row), R₂ (second row), and $R_{1}^{*}$ (third row). The arrows in the energy diagrams correspond to optical transitions detected at different frequencies and times in the 2D spectra. The corresponding transitions are indicated in the response function plot (Fig. 7) as symbols with the respective color. The horizontal dashed lines in the Feynman diagrams denote exciton relaxation, whereas the horizontal solid lines indicate an interaction with the external fields.

In the case of ESA (third row of Fig. 7), the peaks have a negative sign, stemming from the subtraction of complex conjugate response functions in the third-order polarization according to Eq. (20) and indicating an absorption process. The redistribution of relative peak intensity with increasing T in the respective contribution to the 2D spectrum takes place in an opposite way compared to the case of SE. Initially, an excitation pathway involving the energetically higher exciton state with larger transition dipole strength leads to the intensive peak at (ω_τ = 14 820 cm⁻¹, ω_t = 14 600 cm⁻¹), which appears as a cross-peak because the energy shift due to the excitonic splitting has a complementary influence on the frequency of transitions involving ground- and doubly excited states. As the transition frequency between the ground state and an exciton state and between the other exciton state and the doubly excited state is expected to be similar, the peak position closer to the diagonal, namely at (ω_τ = 14 820 cm⁻¹, ω_t = 14 800 cm⁻¹), can be explained by population transfer from the upper to the lower exciton state, as illustrated in the bottom panel of Fig. 8. The peak resulting from population transfer gains intensity, while the initial peak decays until it is not recognizable anymore. The intensity of the former, however, is low, since the transition from the low-energy exciton state to the two-exciton state has the same low dipole strength as the one from the ground state to the low-energy exciton state.⁵⁷ Relying on these findings, one can conclude that the positive-signed diagonal peak in the 2D spectrum in Fig. 4, which is given as a sum of all response functions, results from the dominant GSB and SE contributions, where only the latter ones exhibit a decay during T. The negative-signed off-diagonal peak of the 2D spectrum in Fig. 4 stems from ESA and decays to zero with increasing population time. Even though only a single pronounced peak is recognizable in the 2D spectra of each response function initially, four peaks at different positions can appear. The missing appearance of some peaks expected from the excited-state electronic structure of the dimer can be explained by the involvement of a relatively small transition dipole moment between an exciton state and the ground- or doubly excited state in the underlying excitation pathway, which leads to a too low peak intensity to resolve the respective peak in the 2D spectrum. According to Eqs. (B1)–(B10) and (D9), the transition dipole moments enter to the fourth power in GSB and SE pathways involving only a single exciton state, whereas in the analogous ESA contributions, the transitions with involvement of ground- and doubly excited states are complementary so that the transition dipole moments involved in the initial and the final pair of excitations both enter squared. There is also the possibility that these peaks exhibit a vibrational substructure, in principle, but for the spectral density given in Eq. (27) and for the chosen parameters, it is not resolved. Note, however, that the waiting time dependence of the experimental peaks (Fig. 5) exhibits a small amplitude oscillation that could be related to a vibrational progression, which is neglected in our calculations.

B. Chl b WSCP

As mentioned further above, WSCP can be reconstituted with different types of chlorophylls. Now that we have investigated the 2D spectra of Chl a WSCP, we are going to study how the spectra change when Chl b is bound to WSCP instead. The electronic part of the excited-state Hamiltonian of Chl b slightly differs from the one of Chl a only by a 400 cm⁻¹ blueshift of the site energy, whereas the excitonic couplings are practically identical. The values of J = 83 cm⁻¹ and J = 82 cm⁻¹ have been inferred from a fit of the linear optical and hole-burning spectra for the excitonic coupling in Chl a WSCP and Chl b WSCP, respectively.²⁰ The difference of 1 cm⁻¹ cannot explain the extent of changes in the measured 2D spectra of Chl b WSCP. Different from Chl a, in the case of Chl b, the non-rephasing contribution is more informative and thus shown in Fig. 9, whereas the sum of rephasing and non-rephasing contributions and the separate rephasing contribution are shown in the supplementary material (Figs. S4 and S5). In particular, the peak at (ω_τ = 15 000 cm⁻¹, ω_t = 15 000 cm⁻¹) found in the experiment cannot be reproduced by calculations with our earlier²⁰ model parameters of Chl b WSCP. To get a better understanding of the appearing differences, we also compare the measured and calculated linear absorption and circular dichroism (CD) spectra in Fig. 10. Differences in the absorption spectra mainly appear in the frequency range above 15 500 cm⁻¹ because high-frequency intramolecular vibrations and the next higher electronic state Q_x (see, e.g., Ref. 52) have not been included in the present model, as noted already above. A closer inspection of the low-frequency region reveals that the experimental peak at E = 15 000 cm⁻¹ is slightly higher than in the calculations, and most importantly, at E = 14 600 cm⁻¹, there is an additional experimental peak that cannot be explained by the calculations. Such an additional peak in the Chl b WSCP absorption spectrum has been reported before²² and has been interpreted in terms of a chemically modified Chl b, termed Chl b derivate. The comparison of the experimental and calculated CD spectra in the lower panel of Fig. 10 shows that the peak from the energetically lower exciton state in the experiment has only about half of the intensity of the corresponding peak in the calculations, whereas the high-energy peaks agree nicely, apart from some deviations for energies larger than 15 600 cm⁻¹. Obviously, the experimental CD spectrum is non-conservative, a property that cannot be described by the present exciton theory. The difference between the conservative CD spectrum of Chl a WSCP and the non-conservative CD spectrum of Chl b WSCP has been related to the larger intrinsic CD signal of isolated Chl b as compared to Chl a.¹⁹ To overcome the remaining differences between the simulated and measured results, we extend our model for Chl b WSCP by assuming a mixture that contains some percentage of a derivate of Chl b.²² It turned out that the assumption of equal site energy shifts of a Chl b dimer subunit is appropriate in view of achieving a better agreement between the calculated spectra of the mixture and the measured ones. The mixture also contains a certain percentage of a Chl b derivate homodimer with a chosen energy shift such that the linear absorption spectrum shown in Fig. 11 has a side peak at E = 14 600 cm⁻¹. The energy shift between Chl b and the Chl b derivate was determined as ΔE = −360 cm⁻¹. Furthermore, the coupling was increased to J = 150 cm⁻¹, and the Huang–Rhys factor in the spectral density was multiplied with a factor of 2. Please note that the only critical parameter is the site energy shift of the derivate (supplementary material, Fig. S9).

FIG. 9.

Experimental29 and calculated non-rephasing contributions to the 2D spectra of Chl b WSCP at 77 K at different population times, as indicated. First row: measured 2D spectra; second row: calculated results for Chl b dimer;29 and third row: calculated results for Chl b dimer with involvement of derivate. The upper diagonal peak, the cross-peak, and the lower diagonal derivate peak are referred to as peak 1, peak 2, and peak 3, respectively. Experimental data of Ref. 29 were provided by Fresch et al.

View large Download slide

Experimental²⁹ and calculated non-rephasing contributions to the 2D spectra of Chl b WSCP at 77 K at different population times, as indicated. First row: measured 2D spectra; second row: calculated results for Chl b dimer;²⁹ and third row: calculated results for Chl b dimer with involvement of derivate. The upper diagonal peak, the cross-peak, and the lower diagonal derivate peak are referred to as peak 1, peak 2, and peak 3, respectively. Experimental data of Ref. 29 were provided by Fresch et al.

FIG. 10.

View large Download slide

Comparison of experimental linear absorption spectrum of Chl b WSCP at 77 K and of experimental CD spectrum of Chl b WSCP at 300 K from Ref. 29 with results from the corresponding calculations (upper and lower panels, respectively). Note that the simulation results in a conservative CD spectrum, while the experimental CD spectrum is non-conservative.

FIG. 11.

View large Download slide

Comparison of experimental linear absorption spectrum at 77 K (upper row) and of experimental CD spectrum at 300 K (lower row) from Ref. 29 with a superposition of a calculated spectrum of a Chl b homodimer and the spectrum of the derivate of which 16% is contained in the mixture.

The percentage of the mixture is determined from the fit of the experimental absorption spectrum, in particular from the relative intensities of the main peak at 15 250 cm⁻¹, the shoulder at 15 000 cm⁻¹, and the low-energy peak at 14 600 cm⁻¹. The best correspondence to the experimental data is obtained for a mixture containing 16% of the derivate. The combined absorption spectrum is displayed together with its components and with the measured one in the upper panel of Fig. 11. A further version of this figure, focusing on the low-energy region, is shown in the supplementary material (Fig. S8). The CD spectrum, calculated for this mixture, is compared in the lower panel of Fig. 11 with the experimental data. Unfortunately, the non-conservative nature of the experimental CD spectrum does not allow us to identify the low-energy peak of the derivate. As in the case of Chl a (Fig. 12), we again compare the diagonal cut through the 2D spectrum (sum of rephasing and non-rephasing contributions) at T = 0 fs with the squared absorption spectrum. The respective results from evaluation of both experiment and calculation are shown in Fig. 12. While the width of the most intensive peak is similar, as expected, the side bands at lower frequencies appear with a considerably lower intensity in the squared absorption spectra than in the diagonal cuts through the 2D spectrum both in the experimental and in the calculated data. This finding indicates that, indeed, separate components in the sense of distinguishable compounds are involved, as the sum of peaks of different components with a given percentage in the sample leads to a linear dependence of the relative intensities on this percentage in both absorption spectrum and 2D spectrum, so that in the squared absorption spectrum, the percentage enters quadratically. In contrast, if the peaks would result from states with different dipole strengths of the same component, the compared curves would be expected to become more similar. We used the parameters of the derivate, determined from the linear absorption spectrum, and calculated the 2D spectrum of the mixture (third row of Fig. 9). The additional peak at (ω_τ = 15 000 cm⁻¹, ω_t = 15 000 cm⁻¹) is successfully reproduced in the calculations, which becomes recognizable particularly in the non-rephasing contribution due to elongation of the peaks in antidiagonal direction and thus clearer separation of the peaks on the diagonal. In analogy to the time evolution of the peak amplitudes at selected positions in the case of Chl a WSCP, we also determined the peak amplitudes at selected positions of the measured 2D spectra of Chl b WSCP and the calculated ones of the model involving the derivate. The respective results are shown in Fig. S12 of the supplementary material.

FIG. 12.

View large Download slide

Squared linear absorption spectrum (red line) and diagonal cut through the 2D spectrum (sum of rephasing and non-rephasing contributions) at T = 0 fs from Fig. 9 (black line) are compared for both the calculated and mesured²⁹ results. Experimental data of Ref. 29 were provided by Fresch et al.

C. Treatment of Chl b WSCP with increased coupling as a tetramer

In the case of both Chl a and Chl b constituents, the structure of WSCP containing the chlorophyll dimer units in WSCP with inter-dimer coupling between them seems to require treatment as a tetramer at first glance. However, the inter-dimer coupling is considerably smaller than the reorganization energy of the exciton-vibrational coupling, so our description in the exciton basis would be of limited accuracy because of the Markov and secular treatment of the off-diagonal parts of the exciton-vibrational coupling. Because of these approximations, dynamic localization effects are disregarded in our approach and, hence, artificial delocalization of the exciton wave function in the tetramer could be a consequence. A non-perturbative treatment [e.g., with “Hierarchical Equations of Motion” (HEOM)], which is able to describe such dynamic localization effects, would go beyond the scope of the present work. In order to include these effects implicitly, we have approximated the tetramer by two non-interacting dimers so far. For Chl a WSCP and Chl b WSCP in non-derivate form, the electronic part of the tetramer Hamiltonian in the excited state with excitonic couplings obtained with the refined Poisson-TrEsp method²¹ is basically the same, namely (in units of cm⁻¹)

(46)

Using the given Hamiltonian for simulations in the tetramer configuration, the largest coupling strength between the two dimers is 24 cm⁻¹, which is less than the reorganization energy E_λ = ∫dωωJ(ω) of 63 cm⁻¹ resulting from the spectral density given in Eq. and the Huang–Rhys factor 0.8. Because of this weak coupling, the feasibility of a treatment of the dissipative dynamics in the framework of Redfield theory relying on a representation in the exciton basis is questionable and, therefore, the approximation of two separate dimers will lead to more reliable and meaningful results.

Different from the assumption that Chl b WSCP contains some percentage of a derivate,²² which has been investigated above, in Ref. 29, a different concept was proposed. The authors argue that the appearance of the energetically lower diagonal peaks in the 2D spectra of Chl b WSCP (Fig. 9) suggests an increased coupling between the dimer subunits so that the assumption of tetramer formation is supported. Their hypothesis is that for the given structural arrangement of the chlorophyll units in Chl b WSCP, the dipole strength of 3.83 D is insufficient and a dipole strength of at least 7 D is required to obtain a sufficiently large coupling via the dipole–dipole interaction to describe the measured data. It has been hypothesized that the formation of an H-bond between the formyl group of Chl b is responsible for this twofold enhancement of the dipole strength. While such a hypothesis sounds rather unexpected, the authors report support from an analysis of the overall oscillator strength of Chl a WSCP and Chl b WSCP. A remaining question is whether the assumption of the increased inter-dimer coupling actually improves the similarity of both the linear spectra and the 2D spectra and thus facilitates a mechanistic explanation of the different appearance of the spectra. To find out whether this is the case, we reformulate the model Hamiltonian with the increased coupling (with all matrix elements in units of cm⁻¹) as

(47)

This Hamiltonian is similar to the one from Eq. , just the excitonic couplings are upscaled by a factor

{(\frac{7 D}{3.83 D})}^{2} \approx 3.378

that takes into account the increased dipole strength. Please note that the dipole approximation is approximately valid for WSCP.²¹ The linear absorption and CD spectra for this model are shown in Fig. 13. The calculated linear absorption spectrum in the upper row of Fig. 13 lacks the experimental shoulder around 15 000 cm⁻¹ and exhibits a peak at 14 600 cm⁻¹, different from the measured absorption spectrum. The CD spectrum in the lower row of Fig. 13 differs qualitatively from the experimental data. In particular, the strong low-energy peak at 14 600 cm⁻¹ is not observed in the experiment. Substantial differences to the experimental findings also appear in the 2D spectra of the tetramer. Their non-rephasing contributions are displayed in the third row of Fig. 14.

FIG. 13.

View large Download slide

Upper row: comparison of linear absorption of Chl b WSCP tetramer with increased excitonic couplings described in Eq. (47) with experimental data at T = 77 K. Lower row: comparison of CD for Chl b WSCP tetramer described in Eq. (47) with experimental data at T = 300 K.

FIG. 14.

View large Download slide

Experimental and calculated non-rephasing contribution to the 2D spectra of Chl b WSCP at 77 K at different population times, as indicated. First row: measured 2D spectra; second row: calculated result for Chl b dimer with involvement of derivate; and third row: calculated result under the assumption of increased excitonic coupling and thus tetramer formation in Chl b WSCP.

A comparison with the rephasing contributions to the 2D spectra and the sum of both contributions shown in the supplementary material (Figs. S6 and S7) reveals that the negative-signed peak at (ω_τ = 15 200 cm⁻¹, ω_t = 15 350 cm⁻¹) appearing in the tetramer spectra in the third row of Fig. 14 contains a component from ESA apart from the negative-valued side wing due to the dissipative line shape. In contrast, the ESA peak observed at (ω_τ = 15 200 cm⁻¹, ω_t = 15 000 cm⁻¹) in the experimental 2D spectrum and in the calculated one of the derivate dimer does not yield a recognizable contribution in the calculated 2D spectrum of the tetramer. This ESA peak furthermore decays much faster than the one of the tetramer. The calculated tetramer spectra also fail to reproduce the additional peak at (ω_τ = 14 900 cm⁻¹, ω_t = 14 900 cm⁻¹) that is seen in the experimental spectra. Hence, we have to conclude that an H-bond-induced increase in the dipole strength of Chl b can neither explain the linear nor the 2D spectra of Chl b WSCP. Finally, we have investigated the Chl b WSCP Hamiltonian suggested in Ref. 29 that contains somewhat weaker excitonic couplings. In addition, for this Hamiltonian, we find qualitative deviations between the calculated and measured linear and nonlinear optical spectra (supplementary material, Fig. S12).

IV. CONCLUSIONS

We derived a theoretical description of 2D electronic spectroscopy in the framework of the “partial ordering prescription” (POP) cumulant expansion and applied it to the simulation of the measured spectra of WSCP reconstituted with different chlorophyll variants. It turned out that for Chl a WSCP, the features of 2D spectra from the experiment, which yield insight into energetic structure and transfer dynamics, could be reproduced well by the calculations. Even though Chl a WSCP and Chl b WSCP exhibit an almost negligible difference in the excitonic coupling according to our earlier estimates,²⁰ the measured 2D spectra differ substantially. To obtain a similarity of the measured and calculated 2D spectra in the case of Chl b WSCP, we have to assume that the sample consists of a mixture of Chl b and a derivate of it, where both of them form homodimers. Such a derivate fraction has been proposed previously in a hole-burning study.²²

An alternative explanation has been offered in terms of an increased dipole strength of Chl b, induced by a hydrogen bond, leading to increased excitonic couplings.^29,58 Already, the linear absorption and CD spectra from calculations with the increased coupling clearly differ from the measured ones. Furthermore, in the calculated 2D spectra of the tetramer, an additional ESA peak is observed, which does not appear in the experimental 2D spectra. These findings support our explanation that instead of an increased inter-dimer coupling, a shift of the site energies in the homodimer of a derivate of Chl b, which appears in a mixture together with Chl b WSCP, is responsible for the unexpected changes in the 2D spectra of Chl b WSCP compared to those of Chl a WSCP.

The chemical origin of this shift is an open question. A related question is whether this modification can also occur in light-harvesting complexes such as the LHC II of higher plants. In these complexes, however, the derivate will be much harder to identify since its absorption overlaps with that of Chl a, which is also bound in these complexes.

SUPPLEMENTARY MATERIAL

See the supplementary material for further information regarding analytical considerations about the influence of pulse overlap, formulation of response functions, influence of the assumption of excited-state equilibration to account for vibrational relaxation, connection between diagonal cut through 2D spectrum and squared absorption spectrum, composition of additional rephasing and non-rephasing components of the 2D spectra of Chl a and Chl b, low-energy region of absorption and CD spectra of Chl b WSCP with involvement of derivate, the role of the parameters entering the derivate model, time evolution of amplitudes of selected peaks in the non-rephasing part of the 2D spectra of Chl b WSCP, and the influence of modified couplings on spectra of tetramer model.

ACKNOWLEDGMENTS

The financial support by the Austrian Science Fund (FWF) (Grant Nos. P 33155-NBL and I 6484-N) is gratefully acknowledged. We thank Elisabetta Collini for stimulating discussions and for providing experimental data (linear absorption, circular dichroism, and 2D spectra).

AUTHOR DECLARATIONS

Conflict of Interest

The authors have no conflicts to disclose.

Author Contributions

Michael Riedl: Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Validation (equal); Visualization (equal); Writing – review & editing (equal). Thomas Renger: Conceptualization (equal); Investigation (equal); Funding (equal); Project administration (equal); Methodology (equal); Supervision (equal); Writing – review & editing (equal). Joachim Seibt: Conceptualization (equal); Investigation (equal); Methodology (equal); Software (equal); Supervision (equal); Writing – original draft (equal).

DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

APPENDIX A: DERIVATION OF QME BY USING A PROJECTION OPERATOR TECHNIQUE

A perturbative description of the dissipative dynamics is obtained by formulating the Nakajima–Zwanzig equation and by approximating it to second order in the system-bath coupling. For this purpose, the projection operator

\hat{P}

and its complement

\hat{Q}

with the properties

\begin{gathered} \hat{P} (\hat{O}) = {\hat{W}}_{eq} t r_{b} (\hat{O}), \hat{P} ({\hat{W}}_{sb} (t)) = {\hat{W}}_{eq} {\hat{ρ}}_{s} (t), \\ \hat{Q} (\hat{O}) = (1 - \hat{P}) (\hat{O}), \hat{Q} ({\hat{W}}_{sb} (t)) = δ {\hat{W}}_{sb} (t) \end{gathered}

(A1)

are introduced, where

δ {\hat{W}}_{sb} (t)

accounts for the deviation of the full density matrix from thermal equilibrium so that

δ {\hat{W}}_{sb} (t_{0}) = 0

taking into account initial equilibration. From the derivation, one obtains a description of the dissipative dynamics in terms of⁴¹

\begin{align} \partial_{t} {\hat{ρ}}_{s} (t) & = - i t r_{b} ([{\hat{H}}_{sb}^{s - b} (t), {\hat{W}}_{eq} {\hat{ρ}}_{s} (t)]) \\ - \int_{t_{0}}^{t} t r_{b} ([{\hat{H}}_{sb}^{s - b} (t), (1 - \hat{P}) [{\hat{H}}_{sb}^{s - b} (τ), {\hat{W}}_{eq} {\hat{ρ}}_{s} (τ)]]) d τ . \end{align}

(A2)

APPENDIX B: RESPONSE FUNCTIONS—REFORMULATED

In the following, we formulate the response functions, thereby using the notation introduced in Sec. II E and assuming a thermal equilibration of populations in the excited state, as described in Appendix C. The latter aspect leads to the appearance of modified rates of lifetime broadening

k_{μ ν}^{'}

and modified transition frequencies

{\tilde{ω}}_{μ ν}^{'}

⁠. The SE contributions (R₁ and R₂, involving only the ground state g, but not the doubly excited state f in their excitation pathway) and the ESA contributions (⁠

R_{1}^{*}

and

R_{2}^{*}

⁠, which also involve excitation to the doubly excited state f) with population of an exciton state from the singly excited state subspace during T are obtained as

\begin{align} R_{1 g, pop} (τ, T, t) & = \sum_{α, γ} e^{- i {\tilde{ω}}_{α g}^{'} t} e^{- i {\tilde{ω}}_{γ g} τ} e^{- \sum_{β}^{β \neq α} k_{β α}^{'} t} e^{- \sum_{β}^{β \neq α} k_{β γ} τ} \\ \times e^{G_{α α}^{*} (t) - G_{α α}^{*} (0)} {[e^{\hat{K} T}]}_{α γ} e^{G_{γ γ} (τ) - G_{γ γ} (0)} \\ \times {⟨({\vec{e}}_{4} \cdot {\vec{μ}}_{α}^{g}) ({\vec{e}}_{3} \cdot {\vec{μ}}_{g}^{α}) ({\vec{e}}_{2} \cdot {\vec{μ}}_{γ}^{g}) ({\vec{e}}_{1} \cdot {\vec{μ}}_{g}^{γ})⟩}_{orient}, \end{align}

(B1)

\begin{align} R_{2 g, pop} (τ, T, t) & = \sum_{α, γ} e^{- i {\tilde{ω}}_{α g}^{'} t} e^{i {\tilde{ω}}_{γ g} τ} e^{- \sum_{β}^{β \neq α} k_{β α}^{'} t} e^{- \sum_{β}^{β \neq α} k_{γ β} τ} \\ \times e^{G_{α α}^{*} (t) - G_{α α}^{*} (0)} {[e^{\hat{K} T}]}_{α γ} e^{G_{γ γ}^{*} (τ) - G_{γ γ}^{*} (0)} \\ \times {⟨({\vec{e}}_{4} \cdot {\vec{μ}}_{α}^{g}) ({\vec{e}}_{3} \cdot {\vec{μ}}_{g}^{α}) ({\vec{e}}_{2} \cdot {\vec{μ}}_{g}^{γ}) ({\vec{e}}_{1} \cdot {\vec{μ}}_{γ}^{g})⟩}_{orient}, \end{align}

(B2)

\begin{align} R_{1 f, pop}^{*} (τ, T, t) & = \sum_{α, 2 β, γ} e^{- i {\tilde{ω}}_{2 β α} t} e^{i {\tilde{ω}}_{γ g} τ} e^{- \sum_{2 δ}^{2 δ \neq 2 β} k_{2 δ 2 β}^{'} t - \sum_{β}^{β \neq α} k_{α β}^{'} t} e^{- \sum_{β}^{β \neq α} k_{γ β} τ} \\ \times e^{(G_{2 β 2 β} (t) - G_{2 β 2 β} (0)) + (G_{α α}^{*} (t) - G_{α α}^{*} (0))} \\ \times e^{- (G_{α 2 β} (t) - G_{α 2 β} (0)) - (G_{α 2 β}^{*} (t) - G_{α 2 β}^{*} (0))} \\ \times {[e^{\hat{K} T}]}_{α γ} e^{G_{γ γ}^{*} (τ) - G_{γ γ}^{*} (0)} \\ \times {⟨({\vec{e}}_{4} \cdot {\vec{μ}}_{2 β}^{α}) ({\vec{e}}_{3} \cdot {\vec{μ}}_{α}^{2 β}) ({\vec{e}}_{2} \cdot {\vec{μ}}_{g}^{γ}) ({\vec{e}}_{1} \cdot {\vec{μ}}_{γ}^{g})⟩}_{orient}, \end{align}

(B3)

and

\begin{align} R_{2 f, pop}^{*} (τ, T, t) & = \sum_{α, 2 β, γ} e^{- i {\tilde{ω}}_{2 β α} t} e^{- i {\tilde{ω}}_{γ g} τ} \\ \times e^{- \sum_{2 δ}^{2 δ \neq 2 β} k_{2 δ 2 β}^{'} t - \sum_{β}^{β \neq α} k_{α β}^{'} t} e^{- \sum_{β}^{β \neq α} k_{β γ} τ} \\ \times e^{(G_{2 β 2 β} (t) - G_{2 β 2 β} (0)) + (G_{α α}^{*} (t) - G_{α α}^{*} (0))} \\ \times e^{- (G_{α 2 β} (t) - G_{α 2 β} (0)) - (G_{α 2 β}^{*} (t) - G_{α 2 β}^{*} (0))} \\ \times {[e^{\hat{K} T}]}_{α γ} e^{(G_{γ γ} (τ) - G_{γ γ} (0))} \\ \times {⟨({\vec{e}}_{4} \cdot {\vec{μ}}_{2 β}^{α}) ({\vec{e}}_{3} \cdot {\vec{μ}}_{α}^{2 β}) ({\vec{e}}_{2} \cdot {\vec{μ}}_{γ}^{g}) ({\vec{e}}_{1} \cdot {\vec{μ}}_{g}^{γ})⟩}_{orient} . \end{align}

(B4)

The SE and ESA contributions with coherence in singly excited state read

\begin{align} R_{1 g, coh} (τ, T, t) \\ = \sum_{α, γ}^{α \neq γ} e^{- i {\tilde{ω}}_{α g} t} e^{- i {\tilde{ω}}_{α γ} T} e^{- i {\tilde{ω}}_{α g} τ} e^{- \sum_{β}^{β \neq α} k_{β α} t} e^{- \sum_{β}^{β \neq α} k_{β α} T - \sum_{c}^{c \neq γ} k_{γ c} T} e^{- \sum_{β}^{β \neq α} k_{β α} τ} \\ \times e^{G_{α α} (t) - G_{α α} (0)} e^{(G_{α α} (T) - G_{α α} (0)) + (G_{γ γ}^{*} (T) - G_{γ γ}^{*} (0))} \\ \times e^{- (G_{γ α} (T) - G_{γ α} (0)) - (G_{γ α}^{*} (T) - G_{γ α}^{*} (0))} e^{G_{α α} (τ) - G_{α α} (0)} \\ \times {⟨({\vec{e}}_{4} \cdot {\vec{μ}}_{α}^{g}) ({\vec{e}}_{3} \cdot {\vec{μ}}_{g}^{γ}) ({\vec{e}}_{2} \cdot {\vec{μ}}_{γ}^{g}) ({\vec{e}}_{1} \cdot {\vec{μ}}_{g}^{α})⟩}_{orient}, \end{align}

(B5)

\begin{align} R_{2 g, coh} (τ, T, t) \\ = \sum_{α, γ}^{α \neq γ} e^{- i {\tilde{ω}}_{α g} t} e^{- i {\tilde{ω}}_{α γ} T} e^{i {\tilde{ω}}_{γ g} τ} \\ \times e^{- \sum_{β}^{β \neq α} k_{β α} t} e^{- \sum_{β}^{β \neq α} k_{β α} T - \sum_{c}^{c \neq γ} k_{γ c} T} e^{- \sum_{β}^{β \neq α} k_{γ β} τ} \\ \times e^{(G_{α α} (t) - G_{α α} (0))} e^{(G_{α α} (T) - G_{α α} (0)) + (G_{γ γ}^{*} (T) - G_{γ γ}^{*} (0))} \\ \times e^{- (G_{γ α} (T) - G_{γ α} (0)) - (G_{γ α}^{*} (T) - G_{γ α}^{*} (0))} e^{G_{γ γ}^{*} (τ) - G_{γ γ}^{*} (0)} \\ \times {⟨({\vec{e}}_{4} \cdot {\vec{μ}}_{α}^{g}) ({\vec{e}}_{3} \cdot {\vec{μ}}_{g}^{γ}) ({\vec{e}}_{2} \cdot {\vec{μ}}_{g}^{α}) ({\vec{e}}_{1} \cdot {\vec{μ}}_{γ}^{g})⟩}_{orient} \end{align}

(B6)

and

\begin{align} R_{1 f, coh}^{*} (τ, T, t) & = \sum_{γ, 2 β, α}^{γ \neq α} e^{- i {\tilde{ω}}_{2 β γ} t} e^{- i {\tilde{ω}}_{α γ} T} e^{i {\tilde{ω}}_{γ g} τ} e^{- \sum_{2 δ}^{2 δ \neq 2 β} k_{2 δ 2 β} t - \sum_{β}^{β \neq α} k_{γ β} t} e^{- \sum_{β}^{β \neq α} k_{β α} T - \sum_{β}^{β \neq α} k_{γ β} T} e^{- \sum_{β}^{β \neq α} k_{γ β} τ} \\ \times e^{(G_{2 β 2 β} (t) - G_{2 β 2 β} (0)) + (G_{γ γ}^{*} (t) - G_{γ γ}^{*} (0))} e^{- (G_{γ 2 β} (t) - G_{γ 2 β} (0)) - (G_{γ 2 β}^{*} (t) - G_{γ 2 β}^{*} (0))} \\ \times e^{(G_{α α} (T) - G_{α α} (0)) + (G_{γ γ}^{*} (T) - G_{γ γ}^{*} (0))} e^{- (G_{γ α} (T) - G_{γ α} (0)) - (G_{γ α}^{*} (T) - G_{γ α}^{*} (0))} e^{G_{γ γ}^{*} (τ) - G_{γ γ}^{*} (0)} \\ \times {⟨({\vec{e}}_{4} \cdot {\vec{μ}}_{2 β}^{γ}) ({\vec{e}}_{3} \cdot {\vec{μ}}_{α}^{2 β}) ({\vec{e}}_{2} \cdot {\vec{μ}}_{g}^{α}) ({\vec{e}}_{1} \cdot {\vec{μ}}_{γ}^{g})⟩}_{orient}, \end{align}

(B7)

\begin{align} R_{2 f, coh}^{*} (τ, T, t) & = \sum_{γ, 2 β, α}^{γ \neq α} e^{- i {\tilde{ω}}_{2 β γ} t} e^{- i {\tilde{ω}}_{α γ} T} e^{- i {\tilde{ω}}_{α g} τ} e^{- \sum_{2 δ}^{2 δ \neq 2 β} k_{2 δ 2 β} t - \sum_{β}^{β \neq α} k_{γ β} t} e^{- \sum_{β}^{β \neq α} k_{β α} T - \sum_{β}^{β \neq α} k_{γ β} T} e^{- \sum_{β}^{β \neq α} k_{β α} τ} \\ \times e^{(G_{2 β 2 β} (t) - G_{2 β 2 β} (0)) + (G_{γ γ}^{*} (t) - G_{γ γ}^{*} (0))} e^{- (G_{γ 2 β} (t) - G_{γ 2 β} (0)) - (G_{γ 2 β}^{*} (t) - G_{γ 2 β}^{*} (0))} e^{(G_{α α} (T) - G_{α α} (0)) + (G_{γ γ}^{*} (T) - G_{γ γ}^{*} (0))} \\ \times e^{- (G_{γ α} (T) G_{γ α} (0)) - (G_{γ α}^{*} (T) - G_{γ α}^{*} (0))} e^{G_{α α} (τ) - G_{α α} (0)} {⟨({\vec{e}}_{4} \cdot {\vec{μ}}_{2 β}^{γ}) ({\vec{e}}_{3} \cdot {\vec{μ}}_{α}^{2 β}) ({\vec{e}}_{2} \cdot {\vec{μ}}_{γ}^{g}) ({\vec{e}}_{1} \cdot {\vec{μ}}_{g}^{α})⟩}_{orient}, \end{align}

(B8)

respectively. The GSB contributions are

\begin{align} R_{3 g} (τ, T, t) & = \sum_{α, γ} e^{- i {\tilde{ω}}_{α g} t} e^{i {\tilde{ω}}_{γ g} τ} e^{- \sum_{β}^{β \neq α} k_{β α} t} e^{- \sum_{β}^{β \neq α} k_{γ β} τ} \\ \times e^{G_{α α} (t) - G_{α α} (0)} e^{G_{γ γ}^{*} (τ) - G_{γ γ}^{*} (0)} \\ \times {⟨({\vec{e}}_{4} \cdot {\vec{μ}}_{α}^{g}) ({\vec{e}}_{3} \cdot {\vec{μ}}_{g}^{α}) ({\vec{e}}_{2} \cdot {\vec{μ}}_{g}^{γ}) ({\vec{e}}_{1} \cdot {\vec{μ}}_{γ}^{g})⟩}_{orient} \end{align}

(B9)

and

\begin{align} R_{4 g} (τ, T, t) & = \sum_{α, γ} e^{- i {\tilde{ω}}_{α g} t} e^{- i {\tilde{ω}}_{γ g} τ} e^{- \sum_{β}^{β \neq α} k_{β α} t} e^{- \sum_{β}^{β \neq α} k_{β γ} τ} \\ \times e^{G_{α α} (t) - G_{α α} (0)} e^{G_{γ γ} (τ) - G_{γ γ} (0)} \\ \times {⟨({\vec{e}}_{4} \cdot {\vec{μ}}_{α}^{g}) ({\vec{e}}_{3} \cdot {\vec{μ}}_{g}^{α}) ({\vec{e}}_{2} \cdot {\vec{μ}}_{γ}^{g}) ({\vec{e}}_{1} \cdot {\vec{μ}}_{g}^{γ})⟩}_{orient} . \end{align}

(B10)

Finally, the double coherence contributions are obtained as

\begin{align} R_{3 f}^{*} (τ, T, t) & = \sum_{α, 2 β, γ} e^{- i {\tilde{ω}}_{2 β α} t} e^{- i {\tilde{ω}}_{2 β g} T} e^{- i {\tilde{ω}}_{γ g} τ} \\ \times e^{- \sum_{2 γ}^{2 γ \neq 2 β} k_{2 γ 2 β} t - \sum_{β}^{β \neq α} k_{α β} t} e^{- \sum_{2 γ}^{2 γ \neq 2 β} k_{2 γ 2 β} T} \\ \times e^{- \sum_{β}^{β \neq α} k_{β γ} τ} e^{(G_{2 β 2 β} (t) - G_{2 β 2 β} (0)) + (G_{α α}^{*} (t) - G_{α α}^{*} (0))} \\ \times e^{- (G_{α 2 β} (t) - G_{α 2 β} (0)) - (G_{α 2 β}^{*} (t) - G_{α 2 β}^{*} (0))} \\ \times e^{G_{2 β 2 β} (T) - G_{2 β 2 β} (0)} e^{G_{γ γ} (τ) - G_{γ γ} (0)} \\ \times {⟨({\vec{e}}_{4} \cdot {\vec{μ}}_{2 β}^{α}) ({\vec{e}}_{3} \cdot {\vec{μ}}_{α}^{g}) ({\vec{e}}_{2} \cdot {\vec{μ}}_{γ}^{2 β}) ({\vec{e}}_{1} \cdot {\vec{μ}}_{g}^{γ})⟩}_{orient} \end{align}

(B11)

and

\begin{align} R_{4 f}^{*} (τ, T, t) & = \sum_{α, 2 β, γ} e^{- i {\tilde{ω}}_{α g} t} e^{- i {\tilde{ω}}_{2 β g} T} e^{- i {\tilde{ω}}_{γ g} τ} \\ \times e^{- \sum_{β}^{β \neq α} k_{β α} t} e^{- \sum_{2 γ}^{2 γ \neq 2 β} k_{2 γ 2 β} T} e^{- \sum_{β}^{β \neq α} k_{β γ} τ} \\ \times e^{G_{α α} (t) - G_{α α} (0)} e^{G_{2 β 2 β} (T) - G_{2 β 2 β} (0)} e^{G_{γ γ} (τ) - G_{γ γ} (0)} \\ \times {⟨({\vec{e}}_{4} \cdot {\vec{μ}}_{α}^{g}) ({\vec{e}}_{3} \cdot {\vec{μ}}_{2 β}^{α}) ({\vec{e}}_{2} \cdot {\vec{μ}}_{γ}^{2 β}) ({\vec{e}}_{1} \cdot {\vec{μ}}_{g}^{γ})⟩}_{orient} . \end{align}

(B12)

APPENDIX C: INFLUENCE OF ASSUMPTION OF EXCITED-STATE EQUILIBRATION

By expressing the position operator of a selected vibrational coordinate

{\hat{q}}_{i}

by the one of the shifted coordinate

{\hat{q}}_{i}^{'}

as

{\hat{q}}_{i} = {\hat{q}}_{i}^{'} - d_{i, α α}

in populations of excitonic states, thereby using

d_{i, α α} = | c_{i}^{(α)} |^{2} d_{i}

and

{\hat{u}}_{i}^{'} = ω_{vib, i}^{2} d_{i} {\hat{q}}_{i}^{'} = ω_{vib, i} \sqrt{2 λ_{i}} {\hat{q}}_{i}^{'}

with

d_{i} = \sqrt{\frac{2 S_{i}}{ω_{vib, i}}} = \frac{\sqrt{2 λ_{i}}}{ω_{vib, i}}

⁠, we obtain the modified system–bath coupling Hamiltonian,

\begin{align} {\hat{H}}_{sb}^{'} & = - \sum_{i} {\hat{u}}_{i}^{'} | c_{i}^{(α)} |^{2} | g 〉 〈 g | \\ + \sum_{β} \sum_{γ} \sum_{i} {\hat{u}}_{i}^{'} (c_{i}^{(β)} c_{i}^{(γ)} - δ_{β γ} | c_{i}^{(α)} |^{2}) | β 〉 〈 γ | \\ + \sum_{2 β} \sum_{2 γ} \sum_{i} {\hat{u}}_{i}^{'} (c_{i}^{(2 β)} c_{i}^{(2 γ)} - δ_{2 β 2 γ} | c_{i}^{(α)} |^{2}) | 2 β 〉 〈 2 γ |, \end{align}

(C1)

and the modified system Hamiltonian,

\begin{align} {\hat{H}}_{s}^{'} & = ε_{g}^{'} | g 〉 〈 g | + ε_{β}^{'} | β 〉 〈 β | + ε_{2 β}^{'} | 2 β 〉 〈 2 β | \\ = \sum_{i} λ_{i} | c_{i}^{(α)} |^{2} | g 〉 〈 g | + \sum_{β} (ε_{β} - \sum_{i} λ_{i} | c_{i}^{(β)} |^{4} \\ + \sum_{i} λ_{i} {(| c_{i}^{(β)} |^{2} - | c_{i}^{(α)} |^{2})}^{2}) | β 〉 〈 β | \\ + \sum_{2 β} (ε_{2 β} - \sum_{i} λ_{i} | c_{i}^{(2 β)} |^{4}) | 2 β 〉 〈 2 β | \\ + \sum_{2 β} \sum_{i} λ_{i} {(| c_{i}^{(2 β)} |^{2} - | c_{i}^{(α)} |^{2})}^{2} | 2 β 〉 〈 2 β |, \end{align}

(C2)

where |2β⟩ denotes a doubly excited state. The bath Hamiltonian remains unchanged. Please note that the corresponding contributions of the Hamiltonian in the exciton basis under the assumption of thermal equilibration in the electronic ground state,

{\hat{H}}_{sb}

⁠,

{\hat{H}}_{s}

⁠, and

{\hat{H}}_{b}

⁠, are obtained from Eqs. and by formally setting

c_{i}^{(α)}

to zero.

The assumption of thermal equilibration in an excited state population gains influence on the dissipative dynamics in the QME. By referring to the formulation in Eq. (11) with the auxiliary tensor given in Eq. (12), which we will refer to as $A_{α β γ δ}^{'}$ here, we will separately consider selected stages of an excitation pathway in the following:

For populations during T in SE and ESA contributions, one has α = β and γ = δ. The correlation function entering the auxiliary tensor can be generally expressed as $C_{κ λ μ ν}^{'} = t r_{b} {H_{sb, κ λ}^{'} (t), H_{sb, μ ν}^{'}}$ ⁠. One finds that by reformulating Eq. (12) with C′ in the case of α = γ, only terms with the index pattern ακκα for κ ≠ α remain, while those with index pattern αααα vanish because of vanishing diagonal elements of $H_{sb}^{'}$ in the singly excited state (please note that the different contributions with such an index pattern would compensate each other, anyway). For α ≠ γ, the contributions with index pattern γααγ remain. The remaining contributions are those with involvement of off-diagonal elements of the system–bath coupling in the singly excited state, which, according to Eq. (C2), are not influenced by the redefinition of the thermal equilibrium, as indicated by the appearance of δ_βγ. However, the frequency differences $ω_{α β}^{'} = ε_{α}^{'} - ε_{β}^{'}$ ⁠, which also appear in the auxiliary tensor, are modified due to the changes of the diagonal elements of the system Hamiltonian in Eq. (C2). Thus, the relaxation rates $k_{α β}^{'} = 2 γ_{α β β α} R (\tilde{C} (ω_{α β}^{'}))$ are obtained under the assumption of excited-state equilibration. Furthermore, different from Eq. (C3), one obtains the modified frequency shifts
${\tilde{ω}}_{α 0}^{'} = ω_{α 0} - γ_{α α α α} λ + \sum_{α \neq β} γ_{α β β α} I (\tilde{C} (ω_{α β}^{'})) .$
(C3)
We now consider a coherence during t, which is generated by a further transition to the ground state (SE) or to the doubly excited state (ESA), starting from a singly excited state population. In the case of SE, where a non-zero density matrix element ρ_α0 appears during t, the energy difference between singly excited state and ground state is $ω_{α 0}^{'} = ε_{α} - \sum_{i} 2 λ_{i} | c_{i}^{(α)} |^{4}$ ⁠. The real part of the dissipative terms $- \sum_{β}^{β \neq α} C_{α β β α}^{'} (τ) \exp (- i ω_{β α}^{'} τ)$ with $ω_{β α}^{'} = (ε_{β} - \sum_{i} λ_{i} | c_{i}^{(β)} |^{4} + \sum_{i} λ_{i} {(| c_{i}^{(β)} |^{2} - | c_{i}^{(α)} |^{2})}^{2}) - (ε_{α} - \sum_{i} 2 λ_{i} | c_{i}^{(α)} |^{4})$ leads to lifetime broadening, and the imaginary part results in a frequency shift. Furthermore, the term $- C_{gggg}^{'} (- τ)$ from the auxiliary tensor remains, which leads to a line shape function contribution according to Eq. (25). Please note that because of $- C_{gggg}^{'} (- τ) = - C_{gggg}^{' *} (τ)$ ⁠, the complex conjugate line shape function with positive time argument is obtained. The correlation function itself can be evaluated as $C_{gggg}^{'} = t r_{b} {H_{sb, g g}^{'} (t), H_{sb, g g}^{'}} = t r_{b} {- \sum_{i} {\hat{u}}_{i}^{'} (t) | c_{i}^{(α)} |^{2}, - \sum_{i} {\hat{u}}_{i}^{'} | c_{i}^{(α)} |^{2}} = | c_{i}^{(α)} |^{4} C^{'} (t)$ ⁠, which corresponds to the correlation function C_αααα in the description with thermal equilibrium in the electronic ground state. By considering the connection between the line shape function g(t) and G(t) in Eq. (34), one finds that the shifted frequency differences ${\tilde{ω}}_{α 0}^{'}$ for thermal equilibration in a population of exciton state α are analogous to those given in Eq. (C3) for thermal equilibration in the electronic ground state, as the contribution of the reorganization energy in g(t) leads to a cancellation of the respective differences when g(t) is replaced by G(t) in the SE response functions.
In the case of ESA, where a matrix element ρ_2βα appears during t, the energy difference between doubly and singly excited states is $ω_{2 β α}^{'} = (ε_{2 α} - \sum_{i} λ_{i} | c_{i}^{(2 β)} |^{4} + \sum_{i} λ_{i} {(| c_{i}^{(2 β)} |^{2} - | c_{i}^{(α)} |^{2})}^{2}) - (ε_{α} - \sum_{i} λ_{i} | c_{i}^{(α)} |^{4})$ ⁠. The lifetime broadening and frequency shift during t are the same as in the case of SE. Among the contributions of correlation functions with involvement of diagonal elements of $H_{sb}^{'}$ to the auxiliary tensor, only the one with all indices corresponding to 2β becomes non-zero. This correlation function can be identified as $C_{2 β 2 β 2 β 2 β}^{'} = t r_{b} {H_{sb, 2 β 2 β}^{'} (t), H_{sb, 2 β 2 β}^{'}} = t r_{b} {\sum_{i} {\hat{u}}_{i}^{'} (t) {(| c_{i}^{(2 β)} |^{2} - | c_{i}^{(α)} |^{2})}^{2}, \sum_{i} {\hat{u}}_{i}^{'} {(| c_{i}^{(2 β)} |^{2} - | c_{i}^{(α)} |^{2})}^{2}} = {(| c_{i}^{(2 β)} |^{2} - | c_{i}^{(α)} |^{2})}^{2} C (t)$ ⁠, where again a connection has been drawn to the correlation function of the system–bath coupling dynamics with thermal equilibration in the electronic ground state. Again, the reorganization energy term $\sum_{i} λ_{i} {(| c_{i}^{(2 β)} |^{2} - | c_{i}^{(α)} |^{2})}^{2}$ in the energy of the doubly excited state is compensated by reformulating the ESA response functions in terms of G(t) instead of g(t) so that ${\tilde{ω}}_{2 β α}^{'}$ can be evaluated in analogy to Eq. (C3) as
${\tilde{ω}}_{2 β α}^{'} = ω_{2 β α} - γ_{2 β 2 β 2 β 2 β}^{'} λ + \sum_{2 γ \neq 2 β} γ_{2 β 2 γ 2 γ 2 β}^{'} I (\tilde{C} (ω_{2 β 2 γ}^{'})),$
(C4)
where $γ_{2 α 2 β 2 γ 2 δ}^{'}$ contains the transformation coefficients that are part of the product expression $H_{sb, 2 α 2 β}^{'} (t) H_{sb, 2 γ 2 δ}^{'}$ ⁠. In the framework of different approaches, the assumption of thermal equilibration in an excited-state population is also possible by a polaron transformation, as discussed, e.g., in Ref. 59 in the framework of the non-perturbative “Hierarchical Equations of Motion” (HEOM), which can be reduced to a second-order description under certain approximations (for the connection to POP, see Ref. 42).

APPENDIX D: ORIENTATIONAL AVERAGE

By considering only one molecular complex, a connection can be drawn between the lab-system (X, Y, Z) where the experiment takes place and the internal coordinate system (x, y, z) of the one chosen complex. Vectors in the respective coordinate systems are related to each other via the transformation

{\vec{v}}_{lab} = Φ {\vec{v}}_{mol}

⁶⁰ with

Φ = (\begin{matrix} \cos (ϕ) \cos (θ) \cos (χ) - \sin (θ) \sin (χ) & \sin (ϕ) \cos (θ) \cos (χ) + \cos (ϕ) \sin (χ) & - \sin (θ) \cos (χ) \\ - \cos (ϕ) \cos (θ) \sin (χ) - \sin (θ) \cos (χ) & - \sin (ϕ) \cos (θ) \sin (χ) + \cos (ϕ) \cos (χ) & \sin (θ) \sin (χ) \\ \cos (ϕ) \sin (θ) & \sin (ϕ) \sin (θ) & \cos (θ) \end{matrix}) .

(D1)

The orientational average is defined as

〈 {\vec{f} 〉}_{orient} = \frac{1}{8 π^{2}} \int_{0}^{π} \sin (θ) d θ \int_{0}^{2 π} d ϕ \int_{0}^{2 π} d χ \vec{f} (θ, ϕ, χ) ρ (θ, ϕ, χ),

(D2)

where ρ(θ, ϕ, χ) is the density distribution with respect to the Euler angles. In the case of an isotropic sample, the distribution is uniform, i.e., ρ(θ, ϕ, χ) = 1. In the molecular polarization

{\vec{e}}_{4} \cdot {\vec{P}}_{mol}

⁠, an average is taken over terms of the form

({\vec{e}}_{4} \cdot {\vec{μ}}_{b}^{a}) ({\vec{e}}_{3} \cdot {\vec{μ}}_{d}^{c}) ({\vec{e}}_{2} \cdot {\vec{μ}}_{f}^{e}) ({\vec{e}}_{1} \cdot {\vec{μ}}_{h}^{g}) .

(D3)

We choose the lab-system such that

{\vec{e}}_{1} = {\vec{e}}_{Z} = (0, 0, 1)

⁠,

{\vec{e}}_{2} = (X_{2}, Y_{2}, Z_{2})

⁠,

{\vec{e}}_{3} = (X_{3}, Y_{3}, Z_{3})

⁠, and

{\vec{e}}_{4} = (X_{4}, Y_{4}, Z_{4})

⁠. Inserting these vectors in Eq. and transforming the molecular transition dipole vectors into the laboratory system, we get

\begin{align} ({\vec{e}}_{4} \cdot {\vec{μ}}_{b}^{a}) ({\vec{e}}_{3} \cdot {\vec{μ}}_{d}^{c}) ({\vec{e}}_{2} \cdot {\vec{μ}}_{f}^{e}) ({\vec{e}}_{1} \cdot {\vec{μ}}_{h}^{g}) \\ = (X_{4} \sum_{α}^{x, y, z} μ_{α b}^{a} Φ_{X}^{α} + Y_{4} \sum_{α}^{x, y, z} μ_{α b}^{a} Φ_{Y}^{α} + Z_{4} \sum_{α}^{x, y, z} μ_{α b}^{a} Φ_{Z}^{α}) \\ \times (X_{3} \sum_{α}^{x, y, z} μ_{α d}^{c} Φ_{X}^{α} + Y_{3} \sum_{α}^{x, y, z} μ_{α d}^{c} Φ_{Y}^{α} + Z_{3} \sum_{α}^{x, y, z} μ_{α d}^{c} Φ_{Z}^{α}) \\ \times (X_{2} \sum_{α}^{x, y, z} μ_{α f}^{e} Φ_{X}^{α} + Y_{2} \sum_{α}^{x, y, z} μ_{α f}^{e} Φ_{Y}^{α} + Z_{2} \sum_{α}^{x, y, z} μ_{α f}^{e} Φ_{Z}^{α}) \\ \times (\sum_{α}^{x, y, z} μ_{α h}^{g} Φ_{Z}^{α}) . \end{align}

(D4)

By multiplying this out, we get terms of the form

F_{4} G_{3} H_{2} μ_{α b}^{a} μ_{β d}^{c} μ_{γ f}^{e} μ_{ε h}^{g} Φ_{F}^{α} Φ_{G}^{β} Φ_{H}^{γ} Φ_{Z}^{ε},

(D5)

with F, G, H ∈ {X, Y, Z} in the lab-system. Now, we take the orientational average

μ_{α b}^{a} μ_{β d}^{c} μ_{γ f}^{e} μ_{ε h}^{g} {⟨Φ_{F}^{α} Φ_{G}^{β} Φ_{H}^{γ} Φ_{Z}^{ε}⟩}_{orient} .

(D6)

By calculating the orientational average of the product of transformation matrix elements for every possible combination of the different factors, one obtains⁶⁰

\begin{aligned} {⟨Φ_{Z}^{α} Φ_{Z}^{β} Φ_{Z}^{γ} Φ_{Z}^{ε}⟩}_{orient} & = \frac{δ_{α β} δ_{γ ε} + δ_{α γ} δ_{β ε} + δ_{α ε} δ_{β γ}}{15}, \\ {⟨Φ_{X}^{α} Φ_{X}^{β} Φ_{Z}^{γ} Φ_{Z}^{ε}⟩}_{orient} & = \frac{+ 4 δ_{α β} δ_{γ ε} - δ_{α γ} δ_{β ε} - δ_{α ε} δ_{β γ}}{30} = {⟨Φ_{Y}^{α} Φ_{Y}^{β} Φ_{Z}^{γ} Φ_{Z}^{ε}⟩}_{orient}, \\ {⟨Φ_{X}^{α} Φ_{Z}^{β} Φ_{X}^{γ} Φ_{Z}^{ε}⟩}_{orient} & = \frac{- δ_{α β} δ_{γ ε} + 4 δ_{α γ} δ_{β ε} - δ_{α ε} δ_{β γ}}{30} = {⟨Φ_{Y}^{α} Φ_{Z}^{β} Φ_{Y}^{γ} Φ_{Z}^{ε}⟩}_{orient}, \\ {⟨Φ_{Z}^{α} Φ_{X}^{β} Φ_{X}^{γ} Φ_{Z}^{ε}⟩}_{orient} & = \frac{- δ_{α β} δ_{γ ε} - δ_{α γ} δ_{β ε} + 4 δ_{α ε} δ_{β γ}}{30} = {⟨Φ_{Z}^{α} Φ_{Y}^{β} Φ_{Y}^{γ} Φ_{Z}^{ε}⟩}_{orient} . \end{aligned}

(D7)

Therefore, the final result for the orientational average is

\begin{align} {⟨ ({\vec{e}}_{4} \cdot {\vec{μ}}_{b}^{a}) ({\vec{e}}_{3} \cdot {\vec{μ}}_{d}^{c}) ({\vec{e}}_{2} \cdot {\vec{μ}}_{f}^{e}) ({\vec{e}}_{1} \cdot {\vec{μ}}_{h}^{g}) ⟩}_{orient} & = {⟨ μ_{4 b}^{a} μ_{3 d}^{c} μ_{2 f}^{e} μ_{1 h}^{g} ⟩}_{orient}, \\ = \frac{1}{30} (2 Z_{2} Z_{3} Z_{4} + 4 Z_{2} (X_{3} X_{4} + Y_{3} Y_{4}) - Z_{3} (X_{2} X_{4} + Y_{2} Y_{4}) - Z_{4} (X_{2} X_{3} + Y_{2} Y_{3})) \\ \times ({\vec{μ}}_{b}^{a} \cdot {\vec{μ}}_{d}^{c}) ({\vec{μ}}_{f}^{e} \cdot {\vec{μ}}_{h}^{g}) + \frac{1}{30} (2 Z_{2} Z_{3} Z_{4} + 4 Z_{3} (X_{2} X_{4} + Y_{2} Y_{4}) - Z_{4} (X_{2} X_{3} + Y_{2} Y_{3}) \\ - Z_{2} (X_{3} X_{4} + Y_{3} Y_{4})) ({\vec{μ}}_{b}^{a} \cdot {\vec{μ}}_{f}^{e}) ({\vec{μ}}_{d}^{c} \cdot {\vec{μ}}_{h}^{g}) + \frac{1}{30} (2 Z_{2} Z_{3} Z_{4} + 4 Z_{4} (X_{2} X_{3} + Y_{2} Y_{3}) \\ - Z_{2} (X_{3} X_{4} + Y_{3} Y_{4}) - Z_{3} (X_{2} X_{4} + Y_{2} Y_{4})) ({\vec{μ}}_{b}^{a} \cdot {\vec{μ}}_{h}^{g}) ({\vec{μ}}_{d}^{c} \cdot {\vec{μ}}_{f}^{e}) . \end{align}

(D8)

For parallel field-polarization vectors

{\vec{e}}_{1} = {\vec{e}}_{2} = {\vec{e}}_{3} = {\vec{e}}_{4} = {\vec{e}}_{Z}

⁠, this orientational average simplifies to

\begin{align} 〈 ({\vec{e}}_{4} \cdot {\vec{μ}}_{b}^{a}) ({\vec{e}}_{3} \cdot {\vec{μ}}_{d}^{c}) ({\vec{e}}_{2} \cdot {\vec{μ}}_{f}^{e}) {({\vec{e}}_{1} \cdot {\vec{μ}}_{h}^{g}) 〉}_{orient} \\ = \frac{1}{15} (({\vec{μ}}_{b}^{a} \cdot {\vec{μ}}_{d}^{c}) ({\vec{μ}}_{f}^{e} \cdot {\vec{μ}}_{h}^{g}) + ({\vec{μ}}_{b}^{a} \cdot {\vec{μ}}_{f}^{e}) ({\vec{μ}}_{d}^{c} \cdot {\vec{μ}}_{h}^{g}) \\ + ({\vec{μ}}_{b}^{a} \cdot {\vec{μ}}_{h}^{g}) ({\vec{μ}}_{d}^{c} \cdot {\vec{μ}}_{f}^{e})) . \end{align}

(D9)

REFERENCES

1.

G. D.

Scholes

,

G. R.

Fleming

,

A.

Olaya-Castro

, and

R.

van Grondelle

,

Nat. Chem.

3

,

763

(

2011

).

https://doi.org/10.1038/nchem.1145

Google Scholar

Crossref

PubMed

2.

G. J. S.

Fowler

,

R. W.

Visschers

,

G. G.

Grief

,

R.

van Grondelle

, and

C. N.

Hunter

,

Nature

355

,

848

(

1992

).

https://doi.org/10.1038/355848a0

Google Scholar

Crossref

PubMed

3.

H.

Witt

,

E.

Schlodder

,

C.

Teutloff

,

J.

Niklas

,

E.

Bordignon

,

D.

Carbonera

,

S.

Kohler

,

A.

Labahn

, and

W.

Lubitz

,

Biochemistry

41

,

8557

(

2002

).

https://doi.org/10.1021/bi025822i

Google Scholar

Crossref

PubMed

4.

R.

Saer

,

G. S.

Orf

,

X.

Lu

,

H.

Zhang

,

M. J.

Cuneo

,

D. A.

Myles

, and

R. E.

Blankenship

,

Biochim. Biophys. Acta, Bioenerg.

1857

,

1455

(

2016

).

https://doi.org/10.1016/j.bbabio.2016.04.007

Google Scholar

Crossref

5.

A.

Khmelnitskiy

,

R. G.

Saer

,

R. E.

Blankenship

, and

R.

Jankowiak

,

J. Phys. Chem. B

122

,

3734

(

2018

).

https://doi.org/10.1021/acs.jpcb.7b11763

Google Scholar

Crossref

6.

K.

Vrandecic

,

M.

Rätsep

,

L.

Wilk

,

L.

Rusevich

,

M.

Golub

,

M.

Reppert

,

K.-D.

Irrgang

,

W.

Kühlbrandt

, and

J.

Pieper

,

J. Phys. Chem. B

119

,

3920

(

2015

).

https://doi.org/10.1021/jp5112873

Google Scholar

Crossref

7.

M. J.

Llansola-Portoles

,

F.

Li

,

P.

Xu

,

S.

Streckaite

,

C.

Ilioaia

,

C.

Yang

,

A.

Gall

,

A. A.

Pascal

,

R.

Croce

, and

B.

Robert

,

Biochim. Biophys. Acta, Bioenerg.

1861

,

148078

(

2020

).

https://doi.org/10.1016/j.bbabio.2019.148078

Google Scholar

Crossref

PubMed

8.

J.

Adolphs

and

T.

Renger

,

Biophys. J.

91

,

2778

(

2006

).

https://doi.org/10.1529/biophysj.105.079483

Google Scholar

Crossref

PubMed

9.

F.

Müh

,

M. E.

Madjet

, and

T.

Renger

,

J. Phys. Chem. B

114

,

13517

(

2010

).

https://doi.org/10.1021/jp106323e

Google Scholar

Crossref

10.

L.

De Vico

,

A.

Anda

,

V. A.

Osipov

,

A. O.

Madsen

, and

T.

Hansen

,

Proc. Natl. Acad. Sci. U. S. A.

115

,

E9051

(

2018

).

https://doi.org/10.1073/pnas.1719355115

Google Scholar

Crossref

PubMed

11.

C.

Curutchet

,

J.

Kongsted

,

A.

Muñoz Losa

,

H.

Hossein-Nejad

,

G. D.

Scholes

, and

B.

Mennucci

,

J. Am. Chem. Soc.

133

,

3078

(

2011

).

https://doi.org/10.1021/ja110053y

Google Scholar

Crossref

12.

D.

Bednarczyk

,

O.

Dym

,

V.

Prabahar

,

Y.

Peleg

,

D. H.

Pike

, and

D.

Noy

,

Angew. Chem., Int. Ed.

55

,

6901

(

2016

).

https://doi.org/10.1002/anie.201512001

Google Scholar

Crossref

13.

J. N.

Sturgis

and

B.

Robert

,

J. Phys. Chem. B

101

,

7227

(

1997

).

https://doi.org/10.1021/jp963363n

Google Scholar

Crossref

14.

K.

McLuskey

,

S. M.

Prince

,

R. J.

Cogdell

, and

N. W.

Isaacs

,

Biochemistry

40

,

8783

(

2001

).

https://doi.org/10.1021/bi010309a

Google Scholar

Crossref

PubMed

15.

D.

Horigome

,

H.

Satoh

,

N.

Itoh

,

K.

Mitsunaga

,

I.

Oonishi

,

A.

Nakagawa

, and

A.

Uchida

,

J. Biol. Chem.

282

,

6525

(

2007

).

https://doi.org/10.1074/jbc.m609458200

Google Scholar

Crossref

16.

A.

Agostini

,

D. M.

Palm

,

F.-J.

Schmitt

,

M.

Albertini

,

M. D.

Valentin

,

H.

Paulsen

, and

D.

Carbonera

,

Sci. Rep.

7

,

7504

(

2017

).

https://doi.org/10.1038/s41598-017-07874-6

Google Scholar

Crossref

PubMed

17.

D. M.

Palm

,

A.

Agostini

,

A.-C.

Pohland

,

M.

Werwie

,

E.

Jaenicke

, and

H.

Paulsen

,

ACS Omega

4

,

7971

(

2019

).

https://doi.org/10.1021/acsomega.9b00054

Google Scholar

Crossref

PubMed

18.

J. L.

Hughes

,

R.

Razeghifard

,

M.

Logue

,

A.

Oakley

,

T.

Wydrzynski

, and

E.

Krausz

,

J. Am. Chem. Soc.

128

,

3649

(

2006

).

https://doi.org/10.1021/ja056576b

Google Scholar

Crossref

19.

T.

Renger

,

I.

Trostmann

,

C.

Theiss

,

M. E.

Madjet

,

M.

Richter

,

H.

Paulsen

,

H. J.

Eichler

,

A.

Knorr

, and

G.

Renger

,

J. Phys. Chem. B

111

,

10487

(

2007

).

https://doi.org/10.1021/jp0717241

Google Scholar

Crossref

20.

J.

Adolphs

,

M.

Berrer

, and

T.

Renger

,

J. Am. Chem. Soc.

138

,

2993

(

2016

).

https://doi.org/10.1021/jacs.5b08246

Google Scholar

Crossref

21.

C.

Friedl

,

D. G.

Fedorov

, and

T.

Renger

,

Phys. Chem. Chem. Phys.

24

,

5014

(

2022

).

https://doi.org/10.1039/d1cp03566e

Google Scholar

Crossref

PubMed

22.

J.

Pieper

,

M.

Rätsep

,

I.

Trostmann

,

H.

Paulsen

,

G.

Renger

, and

A.

Freiberg

,

J. Phys. Chem. B

115

,

4042

(

2011

).

https://doi.org/10.1021/jp111455g

Google Scholar

Crossref

23.

J.

Pieper

,

M.

Rätsep

,

I.

Trostmann

,

F.-J.

Schmitt

,

C.

Theiss

,

H.

Paulsen

,

H.

Eichler

,

A.

Freiberg

, and

G.

Renger

,

J. Phys. Chem. B

115

,

4053

(

2011

).

https://doi.org/10.1021/jp111457t

Google Scholar

Crossref

24.

J.

Alster

,

H.

Lokstein

,

J.

Dostal

,

A.

Uchida

, and

D.

Zigmantas

,

J. Phys. Chem. B

118

,

3524

(

2014

).

https://doi.org/10.1021/jp411174t

Google Scholar

Crossref

25.

A. M.

Rosnik

and

C.

Curutchet

,

J. Chem. Theory Comput.

11

,

5826

(

2015

).

https://doi.org/10.1021/acs.jctc.5b00891

Google Scholar

Crossref

26.

A.

Agostini

,

D. M.

Palm

,

H.

Paulsen

, and

D.

Carbonera

,

J. Phys. Chem. B

122

,

6156

(

2018

).

https://doi.org/10.1021/acs.jpcb.8b01906

Google Scholar

Crossref

27.

D. M.

Palm

,

A.

Agostini

,

V.

Averesch

,

P.

Girr

,

M.

Werwie

,

S.

Takahashi

,

H.

Satoh

,

E.

Jaenicke

, and

H.

Paulsen

,

Nat. Plants

4

,

920

(

2018

).

https://doi.org/10.1038/s41477-018-0273-z

Google Scholar

Crossref

PubMed

28.

V.

Prabahar

,

L.

Afriat-Jurnou

,

I.

Paluy

,

Y.

Peleg

, and

D.

Noy

,

FEBS J.

287

,

991

(

2020

).

https://doi.org/10.1111/febs.15068

Google Scholar

Crossref

PubMed

29.

E.

Fresch

,

E.

Meneghin

,

A.

Agostini

,

H.

Paulsen

,

D.

Carbonera

, and

E.

Collini

,

J. Phys. Chem. Lett.

11

,

1059

(

2020

).

https://doi.org/10.1021/acs.jpclett.9b03628

Google Scholar

Crossref

30.

Y.

Lahav

,

D.

Noy

, and

I.

Schapiro

,

Phys. Chem. Chem. Phys.

23

,

6544

(

2021

).

https://doi.org/10.1039/d0cp06582j

Google Scholar

Crossref

PubMed

31.

T.

Brixner

,

J.

Stenger

,

H. M.

Vaswani

,

M.

Cho

,

R. E.

Blankenship

, and

G. R.

Fleming

,

Nature

434

,

625

(

2005

).

https://doi.org/10.1038/nature03429

Google Scholar

Crossref

PubMed

32.

A.

Gelzinis

,

D.

Abramavicius

,

J. P.

Ogilvie

, and

L.

Valkunas

,

J. Chem. Phys.

147

,

115102

(

2017

).

https://doi.org/10.1063/1.4997527

Google Scholar

Crossref

33.

J.

Cao

,

R. J.

Cogdell

,

D. F.

Coker

,

H.-G.

Duan

,

J.

Hauer

,

U.

Kleinekathöfer

,

T. L. C.

Jansen

,

T.

Mančal

,

R. J. D.

Miller

,

J. P.

Ogilvie

,

V. I.

Prokhorenko

,

T.

Renger

,

H.-S.

Tan

,

R.

Tempelaar

,

M.

Thorwart

,

E.

Thyrhaug

,

S.

Westenhoff

, and

D.

Zigmantas

,

Sci. Adv.

6

,

eaaz4888

(

2020

).

https://doi.org/10.1126/sciadv.aaz4888

Google Scholar

Crossref

PubMed

34.

F.

Milota

,

V. I.

Prokhorenko

,

T.

Mancal

,

H.

von Berlepsch

,

O.

Bixner

,

H. F.

Kauffmann

, and

J.

Hauer

,

J. Phys. Chem. A

117

,

6007

(

2013

).

https://doi.org/10.1021/jp3119605

Google Scholar

Crossref

35.

T.

Renger

and

R. A.

Marcus

,

J. Chem. Phys.

116

,

9997

(

2002

).

https://doi.org/10.1063/1.1470200

Google Scholar

Crossref

36.

J.

Seibt

and

T.

Pullerits

,

J. Chem. Phys.

141

,

114106

(

2014

).

https://doi.org/10.1063/1.4897518

Google Scholar

Crossref

37.

J.

Seibt

and

T.

Mančal

,

Chem. Phys.

515

,

129

(

2018

).

https://doi.org/10.1016/j.chemphys.2018.08.026

Google Scholar

Crossref

38.

M.

Yang

and

G. R.

Fleming

,

Chem. Phys.

275

,

355

(

2002

).

https://doi.org/10.1016/s0301-0104(01)00540-7

Google Scholar

Crossref

39.

J.

Seibt

and

T.

Mančal

,

J. Chem. Phys.

146

,

174109

(

2017

).

https://doi.org/10.1063/1.4981523

Google Scholar

Crossref

40.

T.

Renger

,

A.

Klinger

,

F.

Steinecker

,

M.

Schmidt am Busch

,

J.

Numata

, and

F.

Müh

,

J. Phys. Chem. B

116

,

14565

(

2012

).

https://doi.org/10.1021/jp3094935

Google Scholar

Crossref

41.

V.

May

and

O.

Kühn

,

Charge and Energy Transfer Dynamics in Molecular Systems

(

Wiley-VCH

,

Weinheim

,

2011

).

Google Scholar

Crossref

42.

S.

Mukamel

,

I.

Oppenheim

, and

J.

Ross

,

Phys. Rev. A

17

,

1988

(

1978

).

https://doi.org/10.1103/physreva.17.1988

Google Scholar

Crossref

43.

D.

Abramavicius

,

B.

Palmieri

,

D. V.

Voronine

,

F.

Šanda

, and

S.

Mukamel

,

Chem. Rev.

109

,

2350

(

2009

).

https://doi.org/10.1021/cr800268n

Google Scholar

Crossref

PubMed

44.

J.

Kim

,

S.

Mukamel

, and

G. D.

Scholes

,

Acc. Chem. Res.

42

,

1375

(

2009

).

https://doi.org/10.1021/ar9000795

Google Scholar

Crossref

PubMed

45.

S.

Mukamel

,

Principles of Nonlinear Optical Spectroscopy

(

Oxford University Press

,

New York

,

1995

).

Google Scholar

46.

T.

Brixner

,

T.

Mančal

,

I. V.

Stiopkin

, and

G. R.

Fleming

,

J. Chem. Phys.

121

,

4221

(

2004

).

https://doi.org/10.1063/1.1776112

Google Scholar

Crossref

47.

A.

Kell

,

X.

Feng

,

M.

Reppert

, and

R.

Jankowiak

,

J. Phys. Chem. B

117

,

7317

(

2013

).

https://doi.org/10.1021/jp405094p

Google Scholar

Crossref

48.

M.

Toutounji

,

ChemPhysChem

25

,

e202300335

(

2024

).

https://doi.org/10.1002/cphc.202300335

Google Scholar

Crossref

PubMed

49.

J.

Adolphs

,

F.

Maier

, and

T.

Renger

,

J. Phys. Chem. B

122

,

8891

(

2018

).

https://doi.org/10.1021/acs.jpcb.8b08410

Google Scholar

Crossref

50.

J.

Seibt

,

T.

Hansen

, and

T.

Pullerits

,

J. Phys. Chem. B

117

,

11124

(

2013

).

https://doi.org/10.1021/jp4011444

Google Scholar

Crossref

51.

T.

Renger

and

V.

May

,

Phys. Rev. Lett.

84

,

5228

(

2000

).

https://doi.org/10.1103/physrevlett.84.5228

Google Scholar

Crossref

PubMed

52.

J.

Seibt

,

D.

Lindorfer

, and

T.

Renger

,

Photosynth. Res.

156

,

19

(

2023

).

https://doi.org/10.1007/s11120-022-00946-3

Google Scholar

Crossref

PubMed

53.

T.

Renger

,

M. E.

Madjet

,

A.

Knorr

, and

F.

Müh

,

J. Plant Physiol.

168

,

1497

(

2011

).

https://doi.org/10.1016/j.jplph.2011.01.004

Google Scholar

Crossref

54.

R.-X.

Xu

,

P.

Cui

,

X.-Q.

Li

,

Y.

Mo

, and

Y.

Yan

,

J. Chem. Phys.

122

,

041103

(

2005

).

https://doi.org/10.1063/1.1850899

Google Scholar

Crossref

55.

J.

Seibt

and

O.

Kühn

,

J. Chem. Phys.

153

,

194112

(

2020

).

https://doi.org/10.1063/5.0027373

Google Scholar

Crossref

56.

M.

Khalil

,

N.

Demirdöven

, and

A.

Tokmakoff

,

J. Phys. Chem. A

107

,

5258

(

2003

).

https://doi.org/10.1021/jp0219247

Google Scholar

Crossref

57.

T.

Renger

,

J.

Voigt

,

V.

May

, and

O.

Kühn

,

J. Phys. Chem.

100

,

15654

(

1996

).

https://doi.org/10.1021/jp960486z

Google Scholar

Crossref

58.

E.

Fresch

and

E.

Collini

,

Molecules

28

,

3553

(

2023

).

https://doi.org/10.3390/molecules28083553

Google Scholar

Crossref

PubMed

59.

J.

Seibt

and

O.

Kühn

,

J. Phys. Chem. A

125

,

7052

(

2021

).

https://doi.org/10.1021/acs.jpca.1c02684

Google Scholar

Crossref

60.

S. S.

Andrews

,

J. Chem. Educ.

81

,

877

(

2004

).

https://doi.org/10.1021/ed081p877

Google Scholar

Crossref

Theory of 2D electronic spectroscopy of water soluble chlorophyll-binding protein (WSCP): Signatures of Chl b derivate

Author Notes

I. INTRODUCTION

II. THEORETICAL BACKGROUND

A. Components of Hamiltonian

B. Time evolution and derivation of quantum master equation (QME)

C. Formulation of population and coherence dynamics

D. Third-order polarization and 2D spectrum

E. Treatment of dissipative dynamics relying on spectral density

III. RESULTS

A. Chl a WSCP

B. Chl b WSCP

C. Treatment of Chl b WSCP with increased coupling as a tetramer

IV. CONCLUSIONS

SUPPLEMENTARY MATERIAL

ACKNOWLEDGMENTS

AUTHOR DECLARATIONS

Conflict of Interest

Author Contributions

DATA AVAILABILITY

APPENDIX A: DERIVATION OF QME BY USING A PROJECTION OPERATOR TECHNIQUE

APPENDIX B: RESPONSE FUNCTIONS—REFORMULATED

APPENDIX C: INFLUENCE OF ASSUMPTION OF EXCITED-STATE EQUILIBRATION

APPENDIX D: ORIENTATIONAL AVERAGE

REFERENCES

Supplementary Material

Citing articles via

Submit your article

Sign up for alerts

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

Theory of 2D electronic spectroscopy of water soluble chlorophyll-binding protein (WSCP): Signatures of Chl b derivate

Author Notes

I. INTRODUCTION

II. THEORETICAL BACKGROUND

A. Components of Hamiltonian

B. Time evolution and derivation of quantum master equation (QME)

C. Formulation of population and coherence dynamics

D. Third-order polarization and 2D spectrum

E. Treatment of dissipative dynamics relying on spectral density

III. RESULTS

A. Chl a WSCP

B. Chl b WSCP

C. Treatment of Chl b WSCP with increased coupling as a tetramer

IV. CONCLUSIONS

SUPPLEMENTARY MATERIAL

ACKNOWLEDGMENTS

AUTHOR DECLARATIONS

Conflict of Interest

Author Contributions

DATA AVAILABILITY

APPENDIX A: DERIVATION OF QME BY USING A PROJECTION OPERATOR TECHNIQUE

APPENDIX B: RESPONSE FUNCTIONS—REFORMULATED

APPENDIX C: INFLUENCE OF ASSUMPTION OF EXCITED-STATE EQUILIBRATION

APPENDIX D: ORIENTATIONAL AVERAGE

REFERENCES

Supplementary Material

Citing articles via

Submit your article

Sign up for alerts

Related Content

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

This Feature Is Available To Subscribers Only