Signatures of intramolecular vibrational and vibronic Q\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$_{\mathrm{x}}$$\end{document}x–Q\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$_{\mathrm{y}}$$\end{document}y coupling effects in absorption and CD spectra of chlorophyll dimers

An electron-vibrational coupling model that includes the vibronic (non-adiabatic) coupling between the Q\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$_{\mathrm{y}}$$\end{document}y and Q\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$_{\mathrm{x}}$$\end{document}x transitions of chlorophyll (Chl), created by Reimers and coworkers (Scientific Rep. 3, 2761, 2013) is extended here to chlorophyll dimers with interchlorophyll excitonic coupling. The model is applied to a Chl a dimer of the water-soluble chlorophyll binding protein (WSCP). As for isolated chlorophyll, the vibronic coupling is found to have a strong influence on the high-frequency vibrational sideband in the absorption spectrum, giving rise to a band splitting. In contrast, in the CD spectrum the interplay of vibronic coupling and static disorder leads to a strong suppression of the vibrational sideband in excellent agreement with the experimental data. The conservative nature of the CD spectrum in the low-energy region is found to be caused by a delicate balance of the intermonomer excitonic coupling between the purely electronic Q\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$_{\mathrm{y}}$$\end{document}y transition and the Q\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$_{\mathrm{y}}$$\end{document}y transition involving intramolecular vibrational excitations on one hand and the coupling to higher-energy electronic transitions on the other hand. Supplementary Information The online version contains supplementary material available at 10.1007/s11120-022-00946-3.


Introduction
Non-adiabatic couplings between different excited electronic states of chlorophyll (Chl) and related pigments in photosynthetic antennae give rise to a fast internal conversion of excitation energy to the first excited state S 1 , from where excitation energy transfer to the photosynthetic reaction center (RC) starts. In this way, the absorption spectrum of the RC is increased spectrally and spatially. Internal conversion timescales in the 100 fs range have been reported both experimentally, using time-resolve spectroscopy (Shi et al. 2005;Bricker et al. 2015;Meneghin et al. 2017;Song et al. 2019;Do et al. 2022), and theoretically from non-adiabatic excited state molecular dynamics simulations (Bricker et al. 2015;Zheng et al. 2017;Fortino et al. 2021).
For the two lowest excited states, S 1 and S 2 , the non-adiabatic coupling is so strong that a quantum mechanic mixing between the S 0 → S 1 (termed Q y ) and the S 0 → S 2 (termed Q x ) transitions can occur. This mixing is particularly strong if the 0-1 transition of Q y and the 0-0 transition of Q x are in near resonance. Here, "0-1" refers to an electronic transition that is accompanied by the excitation of the first excited state of a high-frequency intramolecular vibrational mode and "0-0" is a purely electronic transition. Such a resonance condition is met, e.g., in Chl a and to a much lesser extent in bacteriochlorophyll a (Reimers et al. 2013).
In a landmark study, Reimers and coworkers (Reimers et al. 2013) provided compelling evidence for this mixing and thereby solved a long-standing puzzle concerning the assignment of the Q y vibrational sideband and the Q x transition (see Fig. 1). Using a simple model that includes a non-adiabatic coupling between the S 1 and S 2 states, which is assumed to be linear in the coordinate of a vibronic coupling mode, and a remaining set of intramolecular modes that couple separately to the Q y and Q x transitions, they were able to explain the absorption, fluorescence, and magnetic circular dichroism spectra of 32 different chlorophyllides in various environments. We follow Reimers et al. and will call this non-adiabatic mixing between Q y and Q x "vibronic coupling". Independent evidence for this vibronic coupling has been reported from polarized two-dimensional electronic 1 3 spectroscopy on Chl a (Song et al. 2019) and Chl c (Bukartė et al. 2020) in solution.
In the present study, the vibronic coupling model of Reimers et al. (2013) is extended to describe optical spectra of chlorophyll dimers with strong excitonic coupling. A suitable model system to study the interplay of intrachromophore vibronic and interchromophore excitonic couplings is the water-soluble chlorophyll binding protein (WSCP). It is a tetramer with quasi D2 symmetry (Horigome et al. 2007;Bednarczyk et al. 2016) containing 4 chlorophyll a pigments that are arranged in two dimers with strong intra-and weak inter-dimer excitonic couplings. WSCP has been used as a simple model system for the study of pigment-pigment and pigment-protein interactions (Hughes et al. 2006;Renger et al. 2007;Adolphs et al. 2016;Friedl et al. 2022;Pieper et al. 2011a, b;Alster et al. 2014;Rosnik and Curutchet 2015;Bednarczyk et al. 2016;Agostini et al. 2018;Palm et al. 2018;Prabahar et al. 2020;Fresch et al. 2020;Lahav et al. 2021). Due to their quasi-symmetric protein environment, all four Chls have the same mean local transition energy (site energy). Another simplifying aspect of WSCP is the fact that wavefunction overlap between the Chls is sufficiently low, such that short-range contribution to the site energies and excitonic couplings can safely be neglected. On the other hand, the Chls in the dimers are close enough for strong excitonic coupling, giving rise to a delocalization of excited states. The latter leads to a strong redistribution of oscillator strength between the two lowest exciton states of the dimer (Hughes et al. 2006;Renger et al. 2007). Due to the open sandwich geometry of local transition dipole moments, the high-energy exciton state gets most of the oscillator strength of the Q y 0-0 transitions of the two Chls.
Whereas this low-energy region of the optical spectra of WSCP has been thoroughly investigated (Hughes et al. 2006;Renger et al. 2007;Renger 2015, 2016;Adolphs et al. 2016;Friedl et al. 2022;Pieper et al. 2011b), the spectrum at higher energies is less understood theoretically. In this region, we expect a mixing between the Q x and Q y transitions of the Chls in the dimer that are coupled vibronically in the monomers and excitonically between the monomers. Interestingly, there are relatively strong experimental signals in the absorption spectrum at higher energies, but barely any signals in the circular dichroism (CD) spectrum (see Fig. 2). Since the intrinsic CD of Chls is small, the CD spectrum is determined by exciton contributions that vanish for localized excited states. In contrast, in the absorption spectrum the excitonic coupling can only redistribute the oscillator strength without changing the overall intensity of the spectrum. Hence, our working hypothesis for the present study is that a localization of excited states involving excited intramolecular vibrational modes occurs that leads to a suppression of the CD signal at high energies.
The remaining parts of this work are organized as follows: We start by summarizing the monomer Hamiltonian of Reimers et al. (2013) containing the vibronic coupling between the Q y and Q x transitions. Next, this Hamiltonian is extended by introducing the intermonomer excitonic coupling. Afterward, the low-frequency part of the spectral density of the local exciton-vibrational coupling of the monomers is used to introduce a system-bath interaction  (Reimers et al. 2013) (black solid line) at T = 300K and calculation (blue dashed line) using a vibronic coupling model adopted from Reimers et al. (2013), as described in the theory part.
In the calculation a 0-0 transition energy of 14910 cm −1 ( 670. Hamiltonian that is applied afterward to derive line-shape functions of optical transitions. Next, we summarize the parameterization of the monomer and dimer Hamiltonians that are applied to describe optical spectra of isolated Chl a in ether and of a Chl a dimer in WSCP. Our theoretical analysis comprises a step-by-step investigation of the influence of different parts of the Hamiltonian on the optical spectra presented in Figs. 1 and 2. Finally we summarize our findings concerning the structure of the Q y -Q x vibrational side band in the absorption spectrum, the suppression of this band in the CD spectrum and the conservative nature of the latter.

Monomer Hamiltonian
Based on the model proposed in Reimers et al. (2013), we start from the Hamiltonian of the monomer with electronic g round st ate | S 0 ⟩ and singly excited st ates | m⟩ ∈ {| S 1 ⟩, | S 2 ⟩, | S 3 ⟩, | S 4 ⟩} . While the states | S 1 ⟩ and | S 2 ⟩ are vibronically coupled, this is not the case for the states | S 3 ⟩ and | S 4 ⟩ , which are taken into account additionally compared to Reimers et al. (2013). Electronic transitions from | S 0 ⟩ to | S 3 ⟩ and to | S 4 ⟩ , denoted as B y and B x , respectively, are distinct components of the Soret band and have been characterized by quantum chemical calculations in an earlier work (Lindorfer et al. 2017). For simplicity they are assumed to be purely electronic states without vibrational substructure in our model. The corresponding eigenenergies are 0 = 0 and m . Electronic excitation from | 0⟩ to | m⟩ goes along with displacement of the potentials of N m vibrational oscillator modes, where each of them is addressed by the index i. We denote their respective vibrational frequencies, position coordinate operators, displacements, and momentum operators as i , q i , d m,i and p i , respectively. While we assume that the displacement of each vibrational mode depends on the electronic state, the vibrational frequency is assumed to be equal in all electronic states. Furthermore, we separately specify the corresponding parameters and operators of the vibronic coupling mode as VC , q VC , d VC and p VC , which are all independent of the electronic state. The vibronic coupling mode can be associated with the coupling mode of a conical intersection, while among the further vibrational modes at least those with different Huang-Rhys factors in Q y and Q x excitation share the properties of tuning modes (Robb 2011). Different from Reimers et al. (2013), we do not apply a scaling of the position coordinate by the square root of the inverse vibrational frequency, q = Reimers et al. (2013) can be identified with ℏ 2 2 q 2 . We furthermore set ℏ = 1 . The scaling of the vibrational coordinate also influences the vibronic coupling constant in terms of ̃V C = √ 1 VC VC . Altogether the monomer Hamiltonian can be written as Please note that d m,i = 0 for m ∈ {S 3 , S 4 } , that is we do not take into account vibrational excitations of the B y (S 0 → S 3 ) and B x (S 0 → S 4 ) transitions. We select the vibronic coupling mode and the high-frequency intramolecular vibrational oscillators to be treated as part of the relevant system, while the low-frequency intermolecular ones enter in terms of contributions to an environment. Note that the harmonic oscillator potentials of the vibronic coupling mode in the diagonal elements are not displaced, but the linear dependence of the off-diagonal elements of the above Hamiltonian on the position operator of the vibronic coupling mode leads to an effective displacement d VC,ef f and accordingly to an effective Huang-Rhys factor S VC,ef f = 1 2 VC d 2 VC,ef f if a diagonalization is applied. The basis representation of the monomer Hamiltonian is described in Appendix A. While this basis representation is formulated without restriction to a certain number of vibrational excitation quanta, we only take a small number of the possible vibrational excitations into account, as described in more detail in "Model assumptions and parameters" Section. However, we want to mention explicitly that simultaneous excitation of an intramolecular vibrational mode and a vibronic coupling mode are accounted for in our model.

Dimer Hamiltonian
The Hamiltonian of an excitonic dimer is composed of Hamiltonians of the monomer units given in Eq. (1), which in the following will be referred to by introducing an additional superscript index, and excitonic coupling contributions between electronic transitions of the different monomer units. Thus, it is of the form where ā and b count the pigments and m and n the electronic states of pigment ā and b , respectively. For the basis (1) (2)

3
representation we use a product of the monomer bases with one of the monomer units being excited and the other one being de-excited (i.e., in the electronic ground state) in the formulation of singly excited states of the dimer. To reduce the size of the basis, we apply the so-called one-particle approximation (OPA) (Spano 2002(Spano , 2010Hestand and Spano 2018) and only take a single product state composed of the lowest vibrational eigenfunctions of the explicitly treated intramolecular vibrational modes into account if the respective monomer unit is de-excited, while the complementary monomer unit in the excited electronic states also exhibits vibrational excitations. The basis representation of the dimer Hamiltonian is described in Appendix B.

System-bath coupling Hamiltonian
From N m = N m,explicit + N m,bath vibrational modes in the monomer Hamiltonian in Eq. (1), N m,explicit high-frequency modes are treated explicitly by including them in the system Hamiltonian, and the remaining N m,bath low-frequency modes enter as contributions to a thermal bath in the framework of the concept of an open quantum system (May and Kühn 2011). The contribution of a bath mode with index i to the system-bath coupling Hamiltonian Ĥ SB is obtained by selecting the sum of the terms with linear dependence on q i from expansion of the squared expression (q i − d m,i ) 2 , whereas the terms q 2 i and d 2 m,i are attributed to the bath Hamiltonian Ĥ B and the system Hamiltonian Ĥ S , respectively. In the bath Hamiltonian also the squared momentum operators p 2 i enter. The system Hamiltonian corresponds to the Hamiltonian of the dimer involving the explicitly treated modes, as specified above, with an energy correction of the basis states by the reorganization energies of the bath modes introduced by the term containing the squared displacement. The combination of the system-bath coupling contributions of the monomer units in a dimer system leads to The basis representation of this operator is formulated in Appendix C, where also the treatment of the dissipative dynamics is described.

Calculation of linear optical spectra
in direction . We furthermore introduce a time evolution operator Û (t) , which comprises the evolution of coherent and dissipative dynamics and is evaluated by time propagation of the QME, in practice. Initially, all elements of the system density matrix are zero, apart from the ground state population. We denote the initial system density matrix as ̂0 ,g and assume an initial equilibration of the thermal bath with density matrix ̂B ,eq . With these definitions we can write the dipole-dipole correlation function of absorption as where ⟨⟨•⟩⟩ = Tr S � Tr B � •̂B ,eq �� denotes the trace over both system and bath. As we take the transition dipole operators as vectorial quantities, the contributions of the directional components are separated from each other. Therefore, an orientational average is obtained by multiplication with a global factor of 1 3 . Furthermore, to take inhomogeneous broadening into account, we vary the electronic excitation energies of the monomer units (site energies) independent of each other by adding random numbers from a Gaussian distribution. This approach is equivalent to an independent incremental variation of the site energies and subsequent weighting of the resulting spectra by the respective Gaussian distribution function.
In the calculation of circular dichroism (CD) spectra not only the electronic transition dipole moments, which under the assumption of localized excitation of electronic state n of a single monomer unit a from the electronic ground state are given as ⃗ (a) n0 =

∑ (a)
,n0 ⃗ e , but also the magnetic transition dipole moments enter. The latter are imaginary and proportional to the cross product of the position vector ⃗ R a of monomer unit a, where the position is identified with the center of mass, and the corresponding electronic transition dipole moment: (Lindorfer and Renger 2018;Weigang 1965). Note that we neglect the small intrinsic CD of Chl a (Lindorfer et al. 2017). In analogy to the case of the electronic transition dipole moment an operator of the magnetic transition dipole moment in the basis of the electronic states for each directional component can be defined as m (ā) = ∑̄n m (ā) ,n0 |n (ā) ⟩⟨0 (ā) | +c.c . Accordingly, after transformation to the exciton basis one obtains the lower and upper triangular part of the vectorial magnetic transition dipole operators ⃗ m exc,+ = ∑ ∑ ⃗ e m exc, , 0 | ⟩⟨0 | and ⃗ m exc,− = ∑ ∑ ⃗ e m exc, ,0 | 0⟩⟨ | , respectively. In analogy to Eq. (4) we formulate the dipole-dipole correlation function of CD as The absorption and CD spectrum is obtained via Fourier transformation as (4) C abs (t) = ⟨⟨⃗ −Û (t)⃗ +̂0,g ⟩⟩, (5) C CD (t) = ⟨⟨⃗ m −Û (t)⃗ +̂0,g ⟩⟩. and respectively.
The correlation functions C abs (t) and C CD (t) are obtained as (see Appendix 9) and with which contains the renormalized excitation energy of exciton state | ⟩, and with Here, is given as ∑ m ∑ nÂ mÂmÂ nÂn and thus corresponds to products of coefficients entering in the transformation from the localized basis to the exciton basis, C ( ) is the Fourier transform of the system-bath correlation function introduced in Appendix C and is the reorganization energy of a local optical excitation. .

Model assumptions and parameters
The model system for the monomer which we use for our calculations is adopted from Reimers et al. (2013). However, we apply some modifications, such as a reduction of the number of intramolecular vibrational modes attributed to the system (see Table 1) and additional involvement of modes attributed to a thermal bath. Furthermore, we extend the model in such a way that it is also capable of describing an excitonic dimer in the framework of the OPA. As we are aiming at a validation of the model by comparison with measured spectra of Chl a from the literature, we adopt the model parameters of the monomer specified in Reimers et al. (2013) for Chl a dissolved in ether. Accordingly, we assume an energy gap between electronic excitation of Q x and Q y as ΔE Q x Q y = Q x − Q y = 1640 cm −1 , a ratio of the squared optical transition dipole moments for electronic excitation of Q x and Q y as f Q x Q y = 0.1 and the full width at half maximum (FWHM) for the inhomogeneous broadening of Q x and Q y as FWHM(Q x ) = 720 cm −1 and FWHM(Q y ) = 360 cm −1 , respectively. To compensate the influence of the bath modes, which are not taken into account in Reimers et al. (2013), we take the liberty of adjusting the inhomogeneous broadening width of Q y to FWHM(Q y ) = 240 cm −1 . For the description of a dimer of Chl a bound to WSCP the inhomogeneous broadening is modified anyway, and values of FWHM(Q x ) = 340 cm −1 and FWHM(Q y ) = 170 cm −1 are assumed. The latter value has been determined previously (Renger et al. 2007), whereas for the former we assume that the relative magnitude with respect to the width of the Q y transition remains the same as in the solvent described above. The B y (S 0 → S 3 ) and B x (S 0 → S 4 ) transitions, which are decoupled from Q y and Q x transitions of the same monomer, are energetically shifted with respect to Q y by ΔE B y Q y = 7570 cm −1 and ΔE B x Q y = 8740 cm −1 , and the ratios of the transition dipole strengths of B y and B x and of Q y are f B y Q y = 2.52 and f B x Q y = 2.43 , respectively, as determined previously (Lindorfer et al. 2017). Note, however, that the exact energies of the B x and B y states are not so critical, since we want to study the influence of these transitions on Table 1 Explicitly treated intramolecular vibrational modes i with their respective frequencies i and Huang-Rhys factors S i . These modes were obtained by separation of sets of five (or six) modes with similar frequencies from the 51 modes given in Reimers et al. (2013), thereby calculating a summed Huang-Rhys factor and an averaged frequency, where in the calculation of the latter the frequency contribution of each mode in a selected set is weighted by the relative contribution to the Huang-Rhys factor of the combined mode. the low-energy (Q y and Q x ) region of the spectra. According to Reimers et al. (2013), we assume a vibronic coupling constant of ̃V C = 750 cm −1 and a frequency of the vibronic coupling mode of VC = 1500 cm −1 . However, different from Reimers et al. (2013), we further reduce the number of intramolecular vibrational modes included in the model. Starting from the 51 effective vibrational modes identified there, we compose nine groups of five modes and one group with six modes by appropriately partitioning the modes sorted by their frequencies and replace each group by a single mode with summed Huang-Rhys factor and averaged frequency, where in the calculation of the frequency average a weighting by the Huang-Rhys factors of the respective modes is applied. The summed Huang-Rhys factors of the ten resulting modes are then attributed to the involvement of the respective modes in the Q y transition. The frequencies and Huang-Rhys factors of the resulting modes are composed in Table 1.
In the case of the Q x transition the Huang-Rhys factors of the lowest-frequency modes ( i = 1, 2, 3 ) are rescaled in such a way that the displacement of these modes in the Q x transition is three times larger than in the Q y transition, as suggested in Reimers et al. (2013), leading to a factor of 9 for the Huang-Rhys factors. This rescaling in the low-frequency range leads to an increase of the total reorganization energy from 262 cm −1 for Q y to 380 cm −1 for Q x . The systembath coupling is determined by the spectral density given in Eq. (C18), where in comparison with Renger and Marcus (2002) the Huang-Rhys factors are rescaled in such a way that the total Huang-Rhys factor becomes equal to 0.8, as determined previously from the temperature dependence of the absorption spectrum of WSCP (Renger et al. 2007). This rescaling leads to the values s 1 = 0.402 and s 2 = 0.398 , while the corresponding frequencies 1 = 0.557 cm −1 and 2 = 1.94 cm −1 remain unchanged. To determine the directions of the transition dipole moments we rely on structural data for the WSCP (Horigome et al. 2007) and determine the vectors connecting the nitrogen atoms at opposite sites of the Chl a molecules. While the transition dipole moment ⃗ Q y is supposed to be aligned with the connecting vector of the nitrogen atoms identified as N B and N D , the transition dipole moment ⃗ Q x is associated with the connecting vector of the nitrogen atoms identified as N A and N C . In fact, this rule of thumb merely gives some orientation. In previous works it turned out that a more accurate description is obtained by applying a rotation of the transition dipole moment in the plane spanned by the connecting vectors of the nitrogen atoms. More precisely, in the case of ⃗ Q y a rotation by −7 • has turned out to be appropriate to reproduce the measured rotational strength of the dimer in WSCP by calculations (Renger et al. 2009). In the present work also rotations of ⃗ Q x , ⃗ B y and ⃗ B x by 20 • , 20 • and −20 • have been applied, respectively, as these values have turned out to be the most appropriate ones for simulation of the measured spectra (see Supporting Information (SI), Section III, Figs. S9-S16). Note, however, that these rotations are not critical for the optical spectra. Qualitatively similar spectra are obtained if ⃗ Q x , ⃗ B y and ⃗ B x are assumed to be oriented along the x-, y-, and x-axes, respectively (SI, Fig. S17). For the excitonic couplings between the electronic transitions of the different monomer units we used the following estimates: For the coupling between the Q y 0-0 transitions we use J Q y Q y = 83 cm −1 . This value has been obtained from a fit of optical spectra (Adolphs et al. 2016) and, recently from quantum chemical/electrostatic calculations (Friedl et al. 2022). In dipole-dipole approximation the excitonic coupling is given as | is the center-to-center distance between monomers a and b and the additional factor (a,b) , respectively. For the present system, the dipoledipole approximation gives a deviation of about 10% for the coupling between the Q y transitions J Q y Q y (Friedl et al. 2022). This deviation is small enough to estimate the excitonic couplings between the remaining transitions as is the ratio of transition dipole magnitudes between the k-th and the Q y transition of monomer c. The resulting values of the excitonic coupling in the subspace of Q x and Q y excitation are .5 cm −1 and J Q y B y = 136.2 cm −1 are obtained. If we are aiming at a comparison of calculated spectra with measured ones, we adjust the energetic position of Q y accordingly, otherwise we set it to zero.
As mentioned in Section "Monomer Hamiltonian" we only take a very limited number of vibrational excitations into account. While only a single intramolecular vibration can be excited to its first vibrational eigenstate, excitation of up to four vibrational quanta of the vibronic coupling mode is possible in the framework of our treatment. The applied restriction with respect to the number of intramolecular vibrational excitations can be justified by the small Huang-Rhys factors of the respective modes. Even though there is no displacement of the vibronic coupling mode upon electronic excitation, the vibronic coupling leads to an effective displacement in the diagonalized Hamiltonian. As the corresponding Huang-Rhys factor is estimated to be larger than the Huang-Rhys factors of the individual intramolecular vibrational modes, we take four excited vibrational eigenstates into account to ensure convergence. Please note that simultaneous excitation of an intramolecular vibrational mode and of the vibronic coupling mode is feasible in our treatment. Such combined excitations gain relevance in the side band region of absorption and CD spectra. However, including them drastically increases the numerical effort, making convergence tests with respect to the number of vibrational modes, included in the model, and with respect to restrictions regarding the number of vibrationally excited basis states difficult. We therefore applied such tests only at an earlier stage of the development of our model, where simultaneous excitation of an intramolecular vibrational mode and of the vibronic coupling mode had not been included yet. Some of the respective results are shown in Section IB of the SI.
It is worth mentioning that the vibronic coupling relies on the concept of a conical intersection. For the description of conical intersections separable tuning and a coupling modes can be identified (Robb 2011), where only the latter actually couples the involved electronic states and thus corresponds to the vibronic coupling mode in our description. For a tuning mode different Huang-Rhys factors in the involved electronic states are required. This criterion is fulfilled for the first three modes with the lowest frequencies from the ten recombined modes in our description. As we treat the vibronic coupling mode independent of the specific vibrational modes with the character of tuning modes (and also independent of the remaining vibrational modes), thereby accounting for simultaneous excitations of these different types of modes in the product basis of the vibrational eigenstates, our description is compatible with the concept of a conical intersection which implies the independence of coupling and tuning mode.

Results and discussion
To illustrate the influence of different aspects of the model on the optical spectra of Chl a monomers in ether ( Fig. 1) and of Chl a dimers in WSCP (Fig. 2), we start with a minimal model and successively extend it until all relevant aspects are included.

Monomer absorption spectra
The homogeneous and inhomogeneous absorption spectra of the Chl a monomer resulting from various approximations are displayed in the upper and lower panel of Fig. 3, respectively. First we consider a very simplified situation, where both the vibronic coupling mode and the intramolecular vibrational modes are disregarded by setting the vibronic coupling constant and the assigned Huang-Rhys factors to zero, respectively, so that the influence of vibrational dynamics only enters in terms of components associated with the thermal bath. Under the assumption that only Q y is excited without involvement of its vibrational structure, only a single peak appears in the absorption spectrum displayed as a black line. This peak is centered at zero on the frequency axis, which corresponds to the electronic excitation energy of the Q y transition. If the coupling of the intramolecular vibrations to the Q y excitation is taken into account by assuming the Huang-Rhys factors in Table 1, the peak centered at zero, which stems from the 0-0 transition, decreases, and due to the involvement of vibrational excitation of the different intramolecular modes a side band appears, as indicated by the red line. Additional involvement of the Q x transition with its vibrational structure leads to an enhancement of the absorption profile in the energetic region of Q x (green and red line). Furthermore, if vibronic coupling between the Q x and Q y transitions is taken into account,  Fig. 3 Absorption spectra of Chl a monomer at T = 300 K without and with inhomogeneous broadening are displayed in upper and lower panel, respectively; black line: without intramolecular vibrations, excitation of Q y only; red line: with intramolecular vibrations, excitation of Q y only; green line: with intramolecular vibrations, excitation of Q y and Q x ; blue line: with additional vibronic coupling between the Q x and Q y transitions the resulting absorption spectrum, displayed as a blue line, exhibits a slightly redshifted 0-0 line and substantial changes in the vibrational structure. Due to the vibronic coupling a band splitting occurs with new bands around 1000 cm −1 and 2000 cm −1 . The splitting appears because of the coupling between the first excited vibrational eigenstate of the vibronic coupling mode in S 1 and the lowest vibrational eigenstate of the vibronic coupling mode in S 2 , which are energetically close to each other and give rise to a quantum mechanical mixing between the 0-1 Q y transition and the 0-0 Q x transition. The origin of the resulting peaks in the monomer absorption spectrum is discussed Section IA in the SI. As we take simultaneous excitation of the vibronic coupling mode and of a single intramolecular mode to its first vibrational eigenstate into account, also the intramolecular modes get involved the vibronic coupling, even though the influence of simultaneous excitation turns out to be negligible from a comparison of Figs. S1 and S2 from the SI. A minor contribution to the peak in the region above 2000 cm −1 and the small peak above 3500 cm −1 can be associated with the Q y transition under the influence of vibronic coupling. If inhomogeneous broadening is involved, the absorption spectra become smoother, and the peaks are spread over a broader frequency range, as displayed in the lower panel of Fig. 3. From the comparison of the measured absorption spectrum for Chl a in ether Reimers et al. (2013) and the calculated result in Fig. 1 (involving adjustment of the energetic position of the zero phonon line to the actual transition frequency of Q y excitation) it becomes recognizable that the involvement of vibronic coupling is required to reproduce the progression of the side bands. Note that the involvement of all 51 modes from Reimers et al. (2013) instead of the ten recombined modes did not lead to recognizable changes in the absorption spectra (see Section IB and Fig. S3 in the SI). In Section IB of the SI we also discussed whether taking into account only the first excited state of a single intramolecular vibrational mode is sufficient. The excitation of more than one vibrational quantum does not lead to qualitative changes in the shape of the spectra (SI, Section IB).

Dimer absorption spectra
Next we study the properties of the dimer model with excitonic coupling between the monomer units, aiming at an understanding of the measured (Palm et al. 2017) and calculated spectra in Fig. 2. We disregard inhomogeneous broadening at first and successively take into account the involvement of different aspects of the model in analogy to the corresponding discussion in the monomer case (Fig. 3). The resulting homogeneous absorption and circular dichroism spectra are shown in the upper and lower left half of Fig. 4, respectively. The corresponding inhomogeneous spectra are shown in the right half of Fig. 4. It turns out that the changes of the dimer absorption spectrum in the region of intramolecular vibrational sidebands are similar to those of the monomer for the considered cases. This finding can be explained by the relatively small values of the excitonic couplings compared to the frequencies of the vibronic coupling mode and of the explicitly treated intramolecular vibrations. Recognizable differences between the absorption spectra of dimer and monomer consist in the less smooth vibrational side bands in the dimer case due to underlying splitting of the excitonically coupled states with a single intramolecular vibrational excitation. The relative amplitude of the peak in the 0-0 Q y region with respect to that of the high-frequency vibrational sidebands (and the Q x transition) is smaller in the case of the dimer (Figs. 2 and 4) than in the case of the monomer (Figs. 1 and 3). This effect is caused by the excitonic splitting between the exciton states formed by the 0-0 Q y transitions and the strong lifetime broadening of the upper exciton state (Renger et al. 2007). Obviously, the lifetime broadening of the (mixed Q y -Q x ) excited states is smaller.
In all of the considered cases in the homogeneous CD spectrum (lower left panel of Fig. 4) a pair of oppositesigned peaks appears close to the zero position associated with the electronic excitation energy of the 0-0 Q y transitions. In the green curve, which is obtained under the assumption of involvement of Q x excitation with its vibrational structure, a peak combination with a sign change appears at the electronic excitation energy of Q x (around 1500 cm −1 ). These peak structures in the energetic region of the Q y and Q x transitions can be attributed to the influence of the excitonic couplings J Q y Q y and J Q x Q x , respectively. From the point of view that J Q x Q x is more than a magnitude smaller than J Q y Q y it seems surprising that the integral over the absolute value of the peak structures in the energetic region of Q x still corresponds to about 1 3 of the integral over the peak structures in the energetic region of the 0-0 Q y transitions. Obviously, the influence of the arrangement of the transition dipole moments involved in the respective transitions overcompensates the tendencies expected for the relative peak intensities by relying on the absolute values of the excitonic coupling matrix elements. Switching on the vibronic coupling between the Q y and Q x transitions leads to a slight redshift of the CD peaks in the Q y (0, 0) region and to a splitting of the Q x (0, 0) CD peaks with new peaks at around 1000 cm −1 and 2000 cm −1 . The rotational strength in the vicinity of the two new peaks is reduced with respect to that of the Q x (0, 0) transition observed in the absence of vibronic coupling. Because of the mixing of Q x and Q y the aspect of the different geometric arrangement of the assigned transition dipole moments gains importance, resulting in a decrease of the integrated CD spectrum due to the involvement of Q y by its vibronic coupling to Q x .
By including inhomogeneous broadening the spectra become considerably smoother (right hand side of Fig. 4). It becomes recognizable that in the energetic regions where side bands with considerable intensity appear in the absorption spectra, the relative intensity with respect to the spectral region of the 0-0 Q y transitions is much lower in the CD spectra. Obviously the sign changes in the peaks, which were observed in the CD spectrum without inhomogeneous broadening, result in cancelation when inhomogeneous broadening gets involved. In addition, the site energy disorder leads to a localization of the exciton states involving excited intramolecular vibrations, because of their small Franck-Condon factors that give rise to a small excitonic coupling. These localized excited states contribute only to the absorption, but not to the circular dichroism spectrum.
Interestingly, the coupling of the Q y transition to intramolecular modes leads to a non-conservative shape of the CD spectrum in the 0-0 spectral region of the Q y transition, which means that integration over the selected spectral region does not lead to compensation of contributions from positive-and negative-signed peaks (Lindorfer et al. 2017;Georgakopoulou et al. 2002Georgakopoulou et al. , 2006Georgakopoulou et al. , 2007Lindorfer et al. 2021). This property becomes recognizable from the colored lines in the lower right panel of Fig. 4. There is redistribution of negative rotational strength of the 0-0 transition to the high-frequency vibrational sideband region that is hardly visible in Fig. 4 because of the large spectral width of the sideband region. Since the experimental CD spectrum of Chl a WSCP is conservative in the 0-0 Q y spectral region (Fig. 2), an additional mechanism seems to be active in the experiment. Such a mechanism is provided by the coupling of Q x and Q y to B x and B y . In order to quantify this effect, we have rescaled all inhomogeneous CD spectra from the lower right panel of Fig. 4, such that their positive peaks get equal amplitude. The resulting spectra are compared in Fig. 5 with a calculation that includes the excitonic couplings to high-energy B x and B y transitions (the orange line in Fig. 5). While the involvement of B x and B y plays a minor role in the absorption spectra (SI, Fig. S4 Fig. 4 are rescaled in such a way that the positive-signed peak has the same amplitude. Different from Fig. 4, also a case with additional involvement of B x and B y is displayed as an orange dashed line a substantial influence on the intensities of the positive-and negative-signed peak component in the energetic region of Q y in the CD spectrum (SI, Fig. S4, lower half). As a measure of the conservativity of the CD spectrum in the lowenergy region we can take the spectrum resulting by including only the 0-0 transitions of the monomers (the black solid line in Fig. 5). Including the intramolecular vibrations of the Q y transition (red line) leads to a redistribution of negative rotational strength from the main peak around 0 to the vibrational sideband around 1000 cm −1 . This redistribution persists upon including the Q x transition and the vibronic coupling between Q x and Q y . Including the B x and B y transitions, the main CD peaks finally get conservative again, in agreement with the experimental data in Fig. 2. In the SI the corresponding absorption and CD spectra are displayed over a wider frequency range which also captures the features resulting from excitation of B x and B y (SI, Fig. S4). Here we concentrate on the low-energy region.
To get further insight into this mechanism which determines the appearance of the dimer spectra, we investigate the role of excitonic coupling involving vibrational excitation.
To study the role of vibrational excitation in excitonic coupling, we include the intramolecular vibrational sidebands of the monomers, but include only the excitonic coupling between the 0-0 transitions of the monomers. The resulting CD spectra are composed in Fig. 6 together with the results of calculations that include all excitonic couplings which have already been discussed in the context of Figs. 4 and 5.
The Q x transition and its vibronic coupling to Q y is included in all calculations and we investigate the involvement of the B x and B y transitions. In the absorption spectrum (upper part of Fig. 6) we find a slight decrease of the vibrational sideband around 1000 cm −1 when excitonic coupling involving the 0-1 Q y transitions is disregarded. Otherwise the vibrational sideband stays practically the same, independent of excitonic coupling to the high-energy B x and B y transitions. The small changes arising from involvement of the 0-1 Q y transitions in the excitonic coupling indicate that the vibrational sideband is that of a monomer (including the vibronic mixing with the Q x transition). This result is in agreement with a recent numerical study by Reppert (2020) on the character of excited states of molecular dimers and with a recent comparison of perturbative and numerically exact lineshape theories (Caycedo-Soler et al. 2022). It supports our earlier treatments Müh and Renger 2012), where intramolecular excitations that include excitonic couplings as well as high-frequency vibrational excitations of pigment-protein complexes were treated as localized transitions, in contrast to alternative treatments that include highfrequency intramolecular modes into the spectral density and use perturbative line-shape theory that cannot distinguish between the different degrees of delocalization of the 0-0 and the 0-1 transitions of the chromophores (Gelzinis et al. 2017;Novoderezhkin et al. 2005).
The lower part in Fig. 6 contains the corresponding CD spectra (As in Fig. 5, for better comparison these spectra were rescaled to yield the same height of their positive peak). The results from calculations with and without contributions of vibrationally excited states in the excitonic coupling, both of them involving B x and B y , are displayed with orange and red line color. Furthermore, the corresponding results for the case without involvement of B x and B y are displayed as a blue and a violet curve, respectively. As described above, there is redistribution of negative rotational strength between the 0-0 Q y transition and the intramolecular vibrational side bands by the excitonic coupling between the 0-0 Q y transition of one monomer and the 0-1 Q y transition of the other monomer, where the increase of relative intensity of the side band is larger than in the absorption spectrum. Similar magnitudes of the negative peak in the CD spectrum are obtained from a full calculation (the orange line) and from a calculation where the coupling between the 0-0 Q y transitions of one monomer and the B x and B y transitions as well as the

Fig. 6
Absorption and CD spectra of dimer at T = 300 K with inhomogeneous broadening ( FWHM(Q y ) = 240 cm −1 , FWHM(Q x ) = 720 cm −1 ), excitation of Q y and Q x and vibronic coupling are displayed in upper and lower panel, respectively; if all matrix elements of the excitonic coupling are taken into account, the blue and the orange line are obtained from calculations without and with involvement of B y and B x ( FWHM(B y ) = FWHM(B x ) = 1800 cm −1 ), respectively; if only those coupling matrix elements with involvement of lowest vibrational eigenfunctions are taken into account, the violet and the red line are obtained from calculations without and with involvement of B y and B x , respectively. Please note that the absorption and CD spectra were rescaled to obtain a common maximum of the 0-0 line (the lowenergy peak). 0-1 Q y transitions of the other monomer are neglected (violet line). The coupling to B x and B y leads to a more negative rotational strength in the Q y 0-0 region and the coupling to Q y 0-1 shows opposite behavior. Hence, an error compensation occurs if both types of couplings are neglected. Now that we have studied the influence of the different parts of our model Hamiltonian on the optical spectra of the Chl a dimer in WSCP, we are ready to discuss the comparison of our calculated spectra with the experimental data in Fig. 2 in detail. The experimental absorption spectrum shows a main peak at around 15000 cm −1 and two broad peaks at 16000 cm −1 and 17000 cm −1 . Our calculations show that the main peak is formed by two exciton transitions with dominant contributions from the 0-0 Q y transitions of the two monomers. Due to the large homogeneous and inhomogeneous broadening these transitions appear as one peak at 300 K . At cryogenic temperatures the lowest exciton transition results in a low-energy shoulder of the main peak (see SI, Fig. S5). As in the monomer spectrum, the two broad high-energy peaks are mixed Q y (0, 1) -Q x (0, 0) transitions. The excitonic coupling between the Q y (0, 1) transitions as well as between the Q x transitions (and those between Q y (0, 1) and Q x ) of the two monomers is found to have a negligible effect on the intensity and shape of the absorption spectrum. Overall, the agreement between theory and experiment is somewhat better for the absorption spectrum of the monomer (Fig. 1) than for that of the dimer (Fig. 2). This result is not surprising since the parameters of the vibronic coupling Hamiltonian of the monomer were fitted to the experiment (here Chl a in ether) by Reimers et al. (2013). We used these parameters also for Chl a in WSCP. The overall Huang-Rhys factor S = ∑ i S i = 0.278 (Table 1) is indeed close to the S = 0.23 estimated recently (Friedl et al. 2022) from fluorescence line narrowing data of WSCP (Pieper et al. 2011a), taking into account the effect of the excitonic coupling on the redistribution of oscillator strength. Concerning the value of the vibronic coupling strength we do not have any independent estimate for Chl a in WSCP, but note that the values estimated by Reimers et al. for Chl a can vary by roughly 100% depending on the solvent (see Table S9 from Reimers et al. (2013)). In the SI we analyze the dependence of the spectra on the variations of Huang-Rhys factors and vibronic coupling constants (SI Section II, Figs. S6-S8). The analysis suggests that reducing the vibronic coupling and increasing the Huang-Rhys factors somewhat improves the fit of the experiment. Our dimer spectra were calculated in the framework of the OPA, as discussed above. To explain the remaining differences between measured and calculated dimer spectra an extension of the model by also including vibrations in the electronic ground state of the de-excited monomer unit in the framework of the so-called two-particle approximation (TPA) (Spano 2002(Spano , 2010Hestand and Spano 2018) could be revealing. Additional vibrational levels in the electronic ground state facilitate additional resonances between electronic de-excitation and excitation involving the vibrational substructure of excited states and ground state, which contribute to the excitonic coupling and would be expected to particularly influence the side band region. However, such an extension of the model would substantially increase the computational effort and would go beyond the scope of the present study. Preliminary calculations with application of the TPA only with respect to the vibronic coupling mode and without the possibility of simultaneous excitation of an intramolecular vibrational mode did not result in substantial differences compared to calculations with OPA, as discussed in Section IB of the SI.
Excellent agreement between theory and experiment is obtained for the CD spectrum (Fig. 2, lower part). In particular, the calculations explain the conservative nature of the spectrum and the absent signal in the high-energy spectral region. Our calculations show that the conservative nature of the CD spectrum rests on a delicate balance between different excitonic couplings between the Q y (0, 0) transition and higher electronic (B x , B y ) and vibronic (Q y (0, 1) -Q x (0, 0) ) transitions. In earlier studies the focus has been on the coupling of the Q y (0, 0) transition of Chl and bacteriochlorophyll (BChl) with higher electronic states of Chls, BChls and carotenoids (Lindorfer et al. 2017;Georgakopoulou et al. 2002Georgakopoulou et al. , 2006Georgakopoulou et al. , 2007. Since these transitions have a much larger dipole strength than the Q y (0, 1) or the Q x (0, 0) transition, it could be expected that their influence is larger. However, the present study shows that at least for Chl a these low-intensity transitions can compensate their small dipole strength by a much smaller energy difference with respect to the Q y (0, 0) transition. Often in the earlier studies only part of the non-conservativity could be explained and the missing part was assumed to originate from the intrinsic circular dichroism of the pigments. It could well be that for different dipole geometries the coupling to high-energy electronic transitions and that to Q y (0, 1) and Q x (0, 0) give nonconservative contributions in the Q y (0, 0) region that add up instead of compensating each other as in the present study. The present calculations also explain why there is practically no CD signal outside the low-energy region that is due to the 0-0 Q y transition. The Q x transitions would have a favorable dipole geometry to create some vibrational strength in the high-energy region. However, due to the vibronic coupling to the Q y transition with unfavorable geometry, the rotational strength is reduced. The remaining high-energy features in the CD spectrum are removed by the inhomogeneous distribution of site energies giving rise to a localization of excited states and overlapping positive and negative bands between the different homogeneous spectra that cancel each other in the inhomogeneous spectrum.
Please note that the aspect of assuming ten modes with summed Huang-Rhys factors and averaged frequency of the 51 modes included in Reimers et al. (2013) is more critical in the case of the dimer than in the case of the monomer. Vibronic effects might be exaggerated, since coherence is maintained much longer in a basis with few stronger vibrations compared to more weaker vibrations with slightly different frequencies. While we found some changes of distinct features in the dimer spectra without inhomogeneous broadening when we included 51 modes, the differences were much less pronounced in the corresponding inhomogeneously broadened spectra (SI, Section IB).

Conclusion
We have developed a description of absorption and CD spectra of a Chl a dimer in WSCP, relying on the model proposed in Reimers et al. (2013) for the Chl a monomer units. While we reduced the number of explicitly treated intramolecular vibrations by combining vibrational modes with similar vibrational frequencies, we extended the model by including additional modes attributed to a thermal bath. It turned out that particularly the involvement of a vibronic coupling mode is of importance for the appropriateness of the model for description of both monomer and dimer spectra. As in the monomer case investigated earlier (Reimers et al. 2013), in the absorption spectrum of the dimer we find strong signatures of the intra-monomer vibronic coupling between the 0-1 Q y transition and the 0-0 Q x transition. The inter-monomer excitonic couplings involving these transitions are rather weak and hidden under the inhomogeneous broadening. The latter leads to an almost complete cancelation of the CD signal outside the spectral region of the 0-0 Q y transition. The conservative nature of the CD signal in the latter region is perturbed by the inter-monomer excitonic coupling between the 0-0 Q y transition and the mixed 0-1 Q y / 0-0 Q x transitions, an effect reported here for the first time. For the present Chl a dimer in WSCP this perturbation is compensated by the excitonic coupling between the 0-0 Q y and the high-energy B x and B y transitions. In related systems with a non-conservative CD spectrum the above compensation might be incomplete.

Appendix A Basis representation of Monomer Hamiltonian
For a representation of vibronic coupling in the basis of the respective vibrational eigenfunctions, we express the position and momentum operators of the vibronic coupling mode in terms of bosonic creation and annihilation operators b † VC and b VC as q VC = √ 1 2 VC (b † VC +b VC ) and . We denote the K-th vibrational eigenfunction of the unshifted oscillator in the electronic states | S 1 ⟩ and | S 2 ⟩ as | S 1 ,VC,K ⟩ and | S 2 ,VC,K ⟩ , respectively, with K ≥ 0 . Using the orthogonality of the vibrational eigenfunctions, ⟨ m,VC,K | n,VC,L ⟩ = K,L , and the properties of t h e c r e a t i o n a n d a For the treatment of the other vibrational modes which are attributed to the system we also introduce a basis of vibrational eigenfunctions, however with respect to the shifted harmonic oscillators in the excited electronic states. We denote these vibrational eigenfunctions as | m,i,k (d m,i )⟩ , where the electronic state to which they are assigned is specified by m, the index i refers to the selected vibrational mode and k is the quantum number of the selected vibrational eigenfunction. The property of the assigned harmonic oscillator to be shifted with respect to the electronic ground state is indicated by the dependence on the respective displacement d m,i , which can formally be expressed by applying a shift operator to the corresponding vibrational eigenfunction of the unshifted oscillator in terms of | m,i,l (d m,i )⟩ = exp(ip m,i d m,i ) | m,i,l ⟩ . Note that, for convenience, we keep referring to the equilibrium position of the electronic ground state as the reference position, different from Reimers et al. (2013), where the equilibrium position in S 1 was taken as a reference. Because of the adjustment of the basis of the vibrational eigenstates to the displacement of the harmonic oscillator, the vibrational basis representation in the subspace of each electronic state m becomes diagonal. The respective matrix element can be specified as In the vibrational basis representation of subspaces of an operator involving different electronic states non-zero matrix elements appear if the displacements of the harmonic oscillator potentials in these electronic states are different. This is the case for non-diagonal subspaces of the Hamiltonian in the electronic basis, such as | S 2 ⟩⟨S 1 | and | S 1 ⟩⟨S 2 | from Eq. (A1), for which the matrix elements in the basis of vibrational eigenstates of the vibronic coupling mode are given in Eq. (A2). For each pair of vibrational basis states of an intramolecular vibrational mode i with vibrational quantum numbers k and l in the system Hamiltonian these matrix elements are multiplied by overlap integrals ⟨ S 2 ,i,k (d S 2 ,i ) | S 1 ,i,l (d S 1 ,i )⟩ and ⟨ S 1 ,i,k (d S 1 ,i ) | S 2 ,i,l (d S 2 ,i )⟩ , respectively. Also in the matrix elements of an electronic transition dipole operator ̂= ∑ m m0 | m⟩⟨0 | +c.c overlap integrals appear, so that the transition dipole strength m0 is multiplied with which for a nonzero displacement d m,i leads to non-vanishing contributions also for different quantum numbers k and l, whereas the corresponding overlap factors of the vibronic coupling mode ⟨ m,VC,K | 0,VC,L ⟩ = K,L require equal quantum numbers K and L for a non-vanishing contribution. The overlap factors of the intramolecular vibrational modes can be calculated by using the recursion schemes from Manneback (1951) and where the absolute value of the expression √ i 2 (d n,i − d m,i ) can be identified with the square root of a Huang-Rhys factor S mn,i , which in the case that m is the electronic ground state results as S n,i = S n0,i = S nm,i = 1 2 i d 2 n,i , because of d m,i = d 0,i = 0 . If a vibrational quantum number k or l in the recursion scheme becomes equal to zero, no further recursion steps are taken with respect to this vibrational quantum number. In the case of l = 0 the right hand side of Eq. (A4) and (A5) is zero, anyway. By using the recursion, any overlap integral can be drawn back to the overlap between the lowest vibrational eigenfunctions, . Therefore, the factor exp − S nm,i 2 also appears in overlap factors involving higher vibrational eigenfunctions. It accounts for the normalization of the vector of basis states and can be adjusted in view of reducing the dimension of the basis of vibrational eigenstates to an extent which still facilitates sufficiently accurate description of relevant aspects of the vibrational dynamics (in our case the vibrational features in a selected region of the absorption spectrum are of interest, as we will discuss later). In principle, the basis states with respect to the N n,explicit intramolecular vibrational modes which are explicitly taken into account in a selected state n of the system Hamiltonian consist of a product of the vibrational eigenfunctions with index i involving all possible combinations of vibrational quantum numbers k i . We specify a vector l of dimension N l,explicit which contains l i as its i-th component to formulate a single product state in terms of | n,l ⟩ = ∏ j | n,j,l j (d n,j )⟩ . The vibrational basis in the subspace of state n can be expressed as a vector of such product states | n,l ⟩ , where each vector component is a product state with a specific combination of vibrational quantum numbers. For the normalization we determine the overlap of the product state of the lowest vibrational eigenfunctions in electronic state m, namely | m,0 ⟩ = ∏ i | m,i,0 ⟩ , with each vector component of | n,l ⟩ . If we, for example, restrict our basis to product states of the lowest vibrational eigenfunctions, except for a single vibrational eigenfunction at most, the norm of the resulting vector can be determined by introducing normalization factors mn,i instead of exp − To arrive at an analogous notation of the basis states as in Reimers et al. (2013), we introduce a product basis of electronic state | m⟩ , of the product state of the vibrational eigenfunction of the explicitly treated intramolecular vibrational modes | m,k ⟩ and of the vibrational eigenfunctions of the vibronic coupling mode | m,VC,K ⟩ in terms of | m,k,K ⟩ =| m⟩ | m,k ⟩ | m,VC,K ⟩ . In this product basis in the 1 3 subspace of the singly excited states with {m, n} ∈ {S 1 , S 2 } the matrix elements and are obtained, where max(K,L) gives the larger of the two vibrational quantum numbers.

Appendix B Basis representation of Dimer Hamiltonian
In the general formulation we obtain the basis representation of the dimer Hamiltonian The first term on the right hand side of Eq. (B8) results as For the evaluation of the third term from Eq. (B8) the coupling contribution Ĵ from the dimer Hamiltonian given in Eq.
(2) is expressed by the sum of its components J (ā)(b) mn |m (ā) ⟩⟨n (b) | , where each component yields a contribution The matrix elements for a coupling contribution with opposite assignment of excited and de-excited monomer unit, i.e., for can be obtained by taking the complex conjugate of Eq. (B8), thereby switching the bra-and ket side of the basis vectors. As a result, one finds a symmetry under reassignment of excited and de-excited monomer unit.
If a given coupling constant is associated with the excitonic coupling between 0-0 transitions, this aspect can be accounted for by division by the squared rescaling factor of the vibrational product states 2 m0 to compensate the appearance of this factor in Eq. (B11).
The transition dipole operator of monomer unit ā with a selected directional component ∈ {x, y, z} is given as ,m0 |m (ā) ⟩⟨0 (ā) | +c.c , where a selected sum component in the basis specified above leads to matrix elements An exciton basis representation of operators given in the basis introduced above with localization of the excitation at a single monomer unit is obtained by applying a transformation with the matrix (B10) where the transformation coefficients c ;

Appendix C System-bath Hamiltonian and dissipative dynamics
The system part of the system-bath coupling contribution of mode q (a) i can be identified as The representation of such an operator in the basis with excitation of a single monomer unit involving the eigenstates of the assigned vibrational modes is obtained in analogy to the corresponding basis representation of the dimer Hamiltonian in Eq. (B8) as and These results lead to the finding that in the given basis the subspace of the excited monomer unit is represented as a unity matrix, whereas the subspace of the de-excited monomer unit is represented by a matrix containing only elements equal to zero. In the following we assume that all operators attributed to the system comprising the electronic states, the explicitly treated intramolecular vibrational modes and the vibronic coupling mode are already represented in the basis introduced above. For each bath mode i the bath contribution of the system-bath coupling, 2 i d m,iqi , enters in a component of the correlation function. We assume that the fluctuations in different electronic states m and n are uncorrelated, leading to the correlation function where the time dependence in q (a) i (t) stems from a representation in the interaction picture with respect to the bath Hamiltonian. For the description of the frequency-dependent oscillator strength for a continuum of bath oscillators covering a selected frequency range we introduce the spectral density with parameters s and , where in comparison with the definition from Renger and Marcus (2002) a factor of 2 is included. From the spectral density the correlation function is obtained as Aiming at a description in the framework of Redfield theory (May and Kühn 2011), we choose a representation in the exciton basis. Accordingly, the transformed system Hamiltonian Ĥ S,exc =Â †Ĥ SÂ corresponds to a diagonal matrix. The coherent part in the quantum master equation (QME) of the transformed system density matrix ̂e xc =Â †̂Â is Kühn (2011), the upper integration border in Eq. (C21) is taken as variable. Otherwise, if the upper integration border is set to infinity, one obtains matrix elements corresponding to the Fourier-transformed correlation function at the difference frequencies associated with the time evolution of the respective matrix elements. The assumption of a variable upper integration border in Λ (a) i,exc (t) (and in its adjoint operator) leads to a description of coherence dynamics at the same level as with the cumulant expansion technique (Mukamel 1995;Renger and Marcus 2002). In the case with variable upper integration border the dissipative contribution to the Liouville-von-Neumann equation in Markov approximation (here meant in the sense that the time integration is drawn into the definition of Λ (a) i,exc (t) instead of involving a convolution with the density matrix) can then be formulated as where Λ †(a) i,exc (t) contains a convolution with the complex conjugate correlation function. In the framework of the secular approximation the population dynamics is separated from the coherence dynamics (Pisliakov et al. 2006). By using subscript indices to denote matrix elements with respect to states and in the exciton basis (involving the vibrational structure of the monomer units of the dimer), the rate expressions are obtained, which enter in the rate equation of population transfer, and in those for coherence dephasing,

Appendix D Lineshape functions
The description of the dissipative dynamics in the framework of Redfield theory with secular approximation turned out to be equivalent to a treatment relying on the partial ordering description (POP) cumulant expansion used in Renger and Marcus (2002). The formulation of this approach involves the frequency-dependent correlation function The real part of this frequency-dependent correlation function can also be obtained directly from the spectral density and from the Bose-Einstein distribution as thereby keeping in mind that different from Renger and Marcus (2002) a factor of 2 has been taken into account in the spectral density given in Eq. (C18). This aspect also plays a role if the line-shape function is compared with the corresponding expression from Renger and Marcus (2002). The connection between this definition of the line-shape function and the definition where it is calculated from the correlation function (Mukamel 1995) as is the following: The analytical expression for the relation between the line-shape function and the spectral density, resulting from the combination of Eq. (D32) with Eq. (C19), can be expressed as (D32) One can easily convince oneself that the real part of the expression from Eq. (D33) corresponds to the difference between ℜ(G(0)) and ℜ(G(t)) . The imaginary part from Eq. (D33) corresponds to the difference between −i t and ℑ(G(t)) , where the reorganization energy is identified as Altogether one obtains Due to the transformation of the system component of the system-bath coupling contributions products of transformation coefficients are multiplied with the line-shape function of the bath component, leading to This definition with appearance of an index pattern implies that only diagonal elements of the system-bath coupling in the exciton basis are selected. Thus, they can be related to selected terms from the lifetime broadening contribution given in the context of Eq. (C25), namely to those involving only diagonal elements of operators attributed to system-bath coupling. Note that integration of the latter, resulting from time evolution of the QME, is required to obtain a line-shape function expression. According to Renger and Marcus (2002), in the treatment of the dissipative dynamics using the partial ordering prescription (POP) cumulant expansion it is justified to treat terms involving off-diagonal matrix elements of the system part of systembath coupling contributions in a different way than terms involving the respective diagonal elements if either the condition >>| | is fulfilled or if the dissipative dynamics takes place on a much longer timescale than the one given by 1∕ . Then the upper integration border in Eq. (C21) is replaced by infinity, leading to an evaluation of the Fourier-transformed correlation function at the frequency difference between the eigenenergies of the involved states of the system.