Comparison of methods to study excitation energy transfer in molecular multichromophoric systems

We compare theoretical methods for calculating excitation energy transfer rates in multichromophoric systems. The employed methods are the multichromophoric Förster resonance energy transfer (MC-FRET), the numerical integration of the Schrödinger equation (NISE), and the Haken-Strobl-Reineker (HSR) model. As a reference, we use the numerically exact Hierarchy of Equations of Motion (HEOM). We examine these methods in various system and bath parameter regimes including the regime relevant to biological light-harvesting systems. We apply them to a model system of a monomeric donor coupled to a multichromophoric acceptor ring of varying size in two limiting configurations, namely symmetric and asymmetric. We find that the symmetric case is more sensitive to the approximations of the methods studied. The NISE method gives the most reasonable results throughout the parameter regimes tested, while still being computationally tractable.


Introduction
The role of electronic excitation energy transfer (EET) is of utmost importance in photochemical and photophysical processes in biological and synthetic systems. In the last decades, there has been much interest in understanding the basic mechanisms of EET in systems consisting of organic molecules. Such interest stems from the aspiration not only to gain fundamental knowledge about nature, but also to comply with the technological demand for organic electronic materials. Basic mechanisms for controlled and efficient energy transfer in photosynthetic organisms can, for example, be adopted to design artificial antenna systems or molecular optoelectronic and photonic devices [1][2][3], photochromic switches [4][5][6], and molecular sensors [5,7]. A large number of theoretical studies have been performed to elucidate, explain, and characterize the energy transport [8][9][10][11]. In assemblies consisting of closely spaced molecules, the EET process is governed by the interplay of the intermolecular interactions in the system itself, giving rise to electronic excited states delocalized over multiple chromophores, and interactions with the noisy environment that inevitably affects the process. For such extended systems, the conventional Förster energy transfer theory does not apply [12] and computationally tractable methods providing reliable results are still not well established, despite the fact that a number of promising methods have been put forward [12][13][14][15][16]. The present work aims to perform a systematic study of the validity, efficiency, and performance of several popular methods to study EET in multichromophoric molecular systems.
The methods we want to evaluate need to contain essential features that have been found to be ubiquitous in natural light-harvesting systems and their synthetic analogues. The natural systems show a large variation in pigment composition, organization, and size, depending on their environments [9] or even light conditions [17,18]. They are typically comprised of tens of pigments, like in the bacterial light-harvesting complexes LH1 [19] and LH2 [20,21], hundreds of pigments, like in photosystems I [22] and II [23,24] of higher plants, or even thousands of pigments, as in the chlorosomes of green sulfur bacteria [25,26]. Synthetic analogues of light-harvesting systems obtained by self-assembly of organic dye molecules such as meso-tetra(4-sulfonatophenyl) porphyrin (TPPS4) [27], the cyanine dye 3,3'-bis(2-sulfopropyl)-5,5',6,6'-tetrachloro-1,1'-dioctylbenzimidacarbocyanine (C8S3) [28][29][30], and zinc chlorin (ZnChl) [31] are usually composed of a large number (many thousands) of pigments. Theoretical methods designed to describe EET within and between such systems, should be able to deal with many strongly coupled chromophores with good accuracy in a reasonable computational time.
Optically relevant excited states in molecular aggregates are typically described as Frenkel excitons, i.e., neutral collective electronic excitations which are superposition states of the excited states of individual molecules within the assembly. Due to strong intermolecular interactions, the electronic excitations are coherently shared between a number of molecules in the aggregate. In the absence of disorder, these excitons are delocalized throughout the whole aggregate. In the regime of weak interactions with vibrations, the EET is initially coherent in time, i.e., characterized by back and forth oscillations. The excitation thus propagates like a wave packet. In the regime of strong coupling to a bath of vibrations, phase relations decay rapidly in time, and the exciton motion takes place through an incoherent hopping process. A proper theoretical description should be able to account for both the coherent transport between strongly coupled chromophores and the incoherent transport between weakly coupled chromophores.
The theory of Förster resonance energy transfer (FRET) [32] has been successfully applied to EET between two weakly coupled chromophores -a donor and an acceptor -in the regime of incoherent energy transfer. It expresses the EET rate in terms of the overlap integral between the emission spectrum of the donor and the absorption spectrum of the acceptor. FRET, which is highly distance-dependent, became a popular tool to probe intermolecular distance in complex systems [33]. Despite its utility, the Förster theory fails when the donor and/or the acceptor system is an aggregate of strongly interacting chromophores and the donor-acceptor distance is comparable to the physical size of the aggregate. In this case, the donor-acceptor interaction can not be treated in the dipole approximation and EET may occur between optically forbidden exciton states. An extension of the Förster theory to multichromophoric donor and/or acceptor systems constitutes one set of methods that we will consider here. This is known as Multichromophoric Förster Resonance Energy Transfer (MC-FRET) theory and was first proposed by Sumi [12], and later further developed and evaluated by Jang et al. [34], Scholes and Fleming [35], Giorda et al. [36], and Cao et al. [37][38][39][40]. In this theory, distances between neighboring molecules within the donor and acceptor aggregates are assumed to be smaller than the closest separation between molecules of the two aggregates considered, but this latter separation is not necessarily larger than the size of the individual aggregates. Such arrangement gives rise to strong couplings within the donor and/or acceptor complexes, leading to delocalized excitonic states within these respective systems, while the weak donor-acceptor coupling still allows for a perturbative treatment of the EET process between the two complexes, albeit it not in the dipole approximation for the aggregates. An important advantage of the MC-FRET method is its ability to incorporate both memory effects and thermal relaxation in the bath description. As a consequence, the Stokes shift is included explicitly, which is important to properly describe the EET from the donor to the acceptor system. The main limitation of the MC-FRET is its restriction to the perturbative limit of the EET process, i.e., the donor-acceptor coupling must be weak compared to the interaction with the bath, and the EET process between the donor and acceptor aggregates must follow an incoherent hopping mechanism.
One of the most acclaimed models describing all regimes between coherent and incoherent motion in molecular aggregates is the Haken-Strobl-Reineker (HSR) model [13,41], which models the effect of the bath as classical fluctuations of the molecular transition energies and the intermolecular excitation transfer interactions. Its popularity originates in its simplicity. As MC-FRET, it is also capable of describing multichromophoric systems, but does not treat the donor-acceptor interaction perturbatively and is, therefore, not limited to the weak donor-acceptor coupling regime. The major limitation of the HSR model lies in its restriction to fast fluctuations of the bath as it assumes the white noise limit, and its intrinsic high-temperature limit. Importantly, the HSR model is able to describe the crossover between coherent and incoherent energy transfer.
In contrast to the HSR model, the Numerical Integration of the Schrödinger Equation (NISE) method [14,42], where the time-dependent Schrödinger equation is explicitly solved and averaged over explicit bath trajectories, properly incorporates the description of the bath from fast to slow modulation limits, yet preserving simplicity in its implementation. It is also capable of describing the multichromophoric nature of the aggregate system and has no restriction on the donoracceptor coupling. As the method neglects the effect of the system on the bath, it gives a thermal realization that corresponds to the high-temperature limit [43].
Finally, a formally exact numerical method describing coupled coherent-incoherent motion with a correct description of the effect of the system on the bath is the Hierarchy of Equations of Motion (HEOM) [44][45][46]. It is based on the stochastic Liouville approach where the system's density matrix is coupled to a hierarchy of auxiliary density matrices that accounts for a non-Markovian system-bath coupling. HEOM has become the "gold standard" against which other approximate methods can be compared. However, unfavorable computational scaling with the system size makes it unappealing for studying molecular systems larger than around 50 molecules [46].
In this paper, we compare the above theoretical methods for studying the EET in molecular aggregates and provide guidelines on the limits of their applicability and performance. We apply these methods to model systems of a monomeric donor coupled in two configurations to a multichromophoric molecular acceptor ring of varying size. We consider a wide range of parameters describing the system, the bath, and the system-bath coupling that includes both limiting cases and parameters comparable to those found in natural light-harvesting systems. We use HEOM as a reference for those parameter sets where the high-temperature approximation may fail.
The outline of the remainder of this paper is as follows. In Section 2, we present the model system, while in Section 3, we describe the pertinent details of the methods for calculating EET that we compare. In Section 4, we present and discuss the results of these methods applied to a monomeric donor coupled to a multichromophoric acceptor ring of varying size. Finally, in Section 5, we present our conclusions.

Model system
We will use the general framework of Frenkel excitons [47,48] to describe the energy transfer in molecular aggregates. For simplicity, we will limit our treatment in this paper to the energy transfer from one molecule (donor) to a complex of strongly interacting molecules (acceptor). All molecules are treated as two-level systems, with a ground state and an optically accessible excited state. The full system is then described by the following Hamiltonian: where the Hamiltonians H t ( ) D and H t ( ) A characterize the donor molecule and the acceptor complex, respectively, and J D A describes the interaction between them. In our study, H t ( ) D , is given by: where D is the electronic transition energy of the donor molecule and t ( ) denotes the Pauli creation (annihilation) operator for an excitation on acceptor molecule n. The summations over n and m run over all molecules ( … N 1, 2, , A ) in the acceptor complex. Furthermore, A is the average transition energy, which in this study is assumed equal for all acceptor molecules, and t ( ) n A indicates the fluctuation in the transition energy for molecule n caused by the interaction with the environment. These fluctuations are assumed uncorrelated for different molecules. Finally, J nm A is the electronic coupling between molecules n and m, which for the sake of simplicity is treated as a constant (no fluctuations); for several light-harvesting complexes this has been shown to be a reasonable approximation [49][50][51]. In the absence of fluctuations, the acceptor's one-quantum excited eigenstates are Frenkel excitons with energies E k A that result from diagonalizing H A with = 0 n A . The last term of Eq. (1) describes the electronic couplings between the donor and acceptor molecules: The electronic couplings in Eq. (3) and (4) can be calculated within the point-dipole approximation [52], multipole approximation [52], extended dipole approximation [53,54], or the transition charge electrostatic potential (TrEsp) method [55,56]. For our purpose of method comparison we will treat the couplings phenomenologically. In particular, within the acceptor ring to be considered, we restrict ourselves to nearest-neighbor interactions. We will consider two specific system geometries, one in which the donor is at the center of a circular acceptor aggregate and one in which the donor resides outside such a circular aggregate (see Section 4).

Methods for calculating the EET rate
In this section, we describe how in the methods considered, the bath is incorporated and the rate is obtained. In order to keep the paper reasonably self-contained and ensure proper translation of parameters between methods, for each method essential details are discussed. Their limitations will be summarized in Table 1.

General considerations
In the MC-FRET method, the EET rate constant between a donor molecule and an acceptor aggregate is directly calculated from the spectral overlap integrals, as will be described in detail later. By contrast, in the NISE, HSR, and HEOM methods, one calculates the probabilities for the excitation to reside on a particular (donor or acceptor) molecule as a function of time t, after placing it at the donor molecule at = t 0, and the EET rate constant is then extracted from the time dependence of these probabilities. In this subsection, we show how we do this.
In this study, we only consider the regime of incoherent energy transfer from a single donor molecule to the acceptor complex. In this case, initially coherent oscillations between the donor and acceptor states are damped fast on the scale of donor-acceptor EET, due to the interaction with the bath of vibrations and the EET may effectively be described through incoherent rate equations. In general, for the probability P t ( ) D for the donor molecule to be excited, this results in a shorttime plateau, possibly with small oscillations, followed by multi-exponential decay to the multichromophoric acceptor. In this paper, we are particularly interested in the situation where the resonance interactions within the acceptor complex are much larger than the donoracceptor interactions. In practice, this implies that equilibration of excitations within the acceptor complex occurs much faster than the EET from the donor to the acceptor complex. As a consequence, the acceptor may be described as one effective site, with excitation probability and the excitation probability is governed by one simple rate equation only: where k is the total EET from the donor molecule to the acceptor ring and k is the backtransfer rate. Eq. (5) gives rise to single-exponential decay, which indeed in all our simulations is observed (after the shorttime transient behavior). The rate k in Eq. (5) is the EET we are after in this study and the quantity that should be compared to the rate obtained from the MC-FRET method. In the following, we describe how k can be extracted from P t ( ) D as computed through the other methods. First, it is important to note that k and k are time-independent rates. Hence, they are related through the donor excitation probability in equilibrium, i.e., at long times. Thus, eq . Substituting this into Eq. (5) and solving for P t ( ) is the effective rate of the single-exponential decay observed for P t ( ) D . Thus, the procedure followed to obtain the rate k from the HSR, NISE, or HEOM methods, is to fit their respective results to Eq. (6), from which values for P D eq and k eff are obtained. From this, we calculate the total EET from the donor to the acceptor complex as = k P k (1 ) D eq eff . In the limit of high temperature, always applicable to the HSR and NISE methods, as these are intrinsically giving high-temperature equilibrium conditions, we have

Treatment of the thermal bath
The thermal noise of the bath is crucial for the EET process. In the following, we discuss how each of the methods for calculating the EET rate incorporates a noise correlation function. To limit the number of free parameters in the description of the bath, we will employ the overdamped Brownian oscillator model [57]. The bath vibrations are then modeled as a collection of independent harmonic oscillators and the system-bath coupling is assumed to be bilinear in system and bath frequency coordinates. The spectral density of this overdamped model is then given by [57]: where is the reorganization energy characterizing the coupling strength of the bath fluctuations to the electronic transition and is the inverse of the correlation time of the bath fluctuations: = 1/ c . The different EET methods employed here utilize different approximations for the bath. In essence, this is related to describing the effect of the bath on the optical line shape of the molecular systems considered at different levels of approximation. The interaction with the bath results in broadening of these line shapes, thus affecting the spectral overlap between the emission line shape of the donor and the absorption line shape of the acceptor complex. Below, we will specify these different approximations. The MC-FRET methods allow to use line shapes based on the full quantum time-correlation function. Within the fluctuation-dissipation theorem, which relates the response function to the correlation function, and using the second-order cumulant expansion, all the necessary information to calculate the optical response function is carried in the quantum time-correlation function C t ( ) q [57,58]. In general, it is complex and is expressed by [57,59,60]:

Table 1
Summary of the limitations of the methods considered in this work. "-" means no restrictions are imposed by the method, W is the exciton bandwidth, J D A is the largest donor-acceptor coupling.

Method
Temperature Bath fluctuations Coupling strength where D ( ) is the spectral density of the harmonic bath, which in our case is given by Eq. (7). The imaginary part of Eq. (8) accounts for the Stokes shift. The effect of the bath in the NISE method is accounted for by a stochastic model known as the overdamped Brownian oscillator. It assumes that the classical heat bath randomly pushes the system, thereby, affecting the molecular excitation energies of the system, expressed by the fluctuating terms t ( ) (2) and (3). Here, the influence of the system on the bath is neglected and as a result the Stokes shift is completely omitted. For the Brownian oscillators, the stochastic process is Gaussian with mean zero = t ( ( ) 0) and twotime correlation functions given by [57]: where the angular brackets … indicate averaging over the stochastic process and generically indicates the energy fluctuation on any of the molecules; correlations in these fluctuations between different molecules are neglected. Here, indicates the root-mean-squared amplitude (or magnitude) of the energy fluctuations, which may be related to the reorganization energy and temperature by taking the classical limit [57] This allows us to translate the reorganization energy and temperature of the MC-FRET method to the fluctuation amplitude in the NISE method.
The HSR model also uses the Gaussian stochastic process described above, but assumes the correlation time of the bath fluctuations, 1 , to be infinitely short. In other words, the bath fluctuations have a white noise spectrum. Physically, this implies that the frequency distribution of the vibrations is significantly broader than the exciton bandwidth. In this case, the correlation function of Eq. (9), C t ( ), reduces to where the ratio / 2 can be expressed as a single parameter, known as the homogeneous line width: This situation corresponds to the fast modulation limit in line shape theory [58], giving rise to a Lorentzian line shape of the optical spectrum of single molecules with a full-width at half-maximum (FWHM) of 2 . The value of corresponding to the bath used in MC-FRET and NISE can be determined from Eq. (12). A convenient dimensionless parameter, , is used for characterizing the speed of the dynamics of the bath [58]: In the fast modulation limit, the bath fluctuations are fast compared to the bath coupling strength, and thus 1. In the slow modulation limit, the opposite holds. The line shape of the optical spectrum evolves from a Lorentzian with no Stokes shift for 1 to a Gaussian with a Stokes shift of 2 for 1. In the intermediate regime, the line shape is described by a Voigt profile [57].
In the remainder of this section, we explain how in each method considered, the EET dynamics within the coupled donor-acceptor system described by the Hamiltonian in Eq. (1) and with the bath treatments given above, is evaluated.

Multichromophoric Förster resonance energy transfer (MC-FRET)
The equation for the MC-FRET rate constant derived in Ref. [12] has the form: where I ( ) Within the second-order cumulant approximation, the emission spectrum of the monomeric donor is given by [57]: where the line broadening function g t ( ) is derived from the quantum correlation function of Eq. (8), giving the equation [ with D ( ) the spectral density as defined in Eq. (7). For the multichromophoric acceptor system, in general, there is no simple expression for the absorption line shape. Below we consider three commonly used approximations that will define different ways to obtain the MC-FRET rates.
The first approximation considered uses the ideal absorption line shape (IAL) [59,60] and is labeled MC-FRET IAL . The expression for the IAL is based on the quantum master equation formalism and without quasistatic disorder has the form: where H A 0 is the unperturbed exciton Hamiltonian of the acceptor complex, i.e., Eq.
where E k A and c k n are the eigenvalues and the eigenvector coefficients of The second approximation considered is the so-called diagonal (secular) MC-FRET approximation [35,60,61], also referred to as the generalized Förster rate, labeled here MC-FRET diag . It is based on the assumption that the acceptor absorption density operators are diagonal in the eigenbasis. This means that the bath fluctuations are not large enough to considerably scramble the exciton basis. Therefore, the expression for the rate constant in Eq. (14) includes the acceptor-donor coupling in the eigenstate basis: 20) and the following expression for the diagonal elements of the absorption matrix: where g t ( ) is a single molecule line broadening function, which is assumed to be the same as for the donor molecule in Eq. (16).
The third approximation takes into account exchange narrowing in the acceptor complex, i.e., the fact that the electronic couplings reduce the line width of a coupled multichromophoric system compared to a single chromophore [62][63][64]. The exchange narrowing effect is well known for J-aggregates in the limit of static disorder [65]. Taking the exchange narrowing effect into account, we obtain the inverse participation ratio (IPR) MC-FRET approximation [66], labeled MC-FRET IPR , which uses the eigenstate absorption line shape of Eq. (21) scaled by the corresponding participation ratio of the eigenstates: where P = c | |

Numerical integration of the Schrödinger equation (NISE)
In the NISE method [14,42], the dynamics is modeled explicitly averaging over trajectories using wave function theory. The dynamics of the system is determined by the time-dependent Schrödinger equation: where t ( ) is the exciton wave function at time t. This equation is solved numerically by dividing time into small intervals, t (with t 1 ), and solving the Schrödinger equation for each interval assuming the Hamiltonian to be constant over it: where , ) is the time-evolution matrix: Here, H t ( ) is the Hamiltonian which parametrically depends on the time-evolution of the classical bath coordinates. The time-evolution for longer times is obtained through multiplying consecutive time-evolution matrices, leading to the expression: The effect of the bath on the system is accounted for by the stochastic fluctuations t ( ) D and t ( ) n A in Eqs. (2) and (3), respectively, characterized by Eq. (9). To this end, we construct the Hamiltonian trajectories in the following way [68]: where G ( ) is a random number taken from a Gaussian distribution with width and mean zero. At the first time step the fluctuating molecular frequencies are initialized as random numbers from a Gaussian distribution with width and mean zero.
In the NISE method, we monitor the population transfer between the donor molecule and the acceptor complex through the decay of the population of the donor: where D is the state in which the donor is excited. To obtain the results reported in Section 4, we averaged over 50 000 random realizations of the Hamiltonian trajectories. The rate constant can be obtained by analyzing the population transfer, as discussed in Section 3.1 (Eq. (6)).
The main advantages of the NISE method are its simplicity in implementation and the ability to properly incorporate colored noise, i.e., to use arbitrary values of . The most important disadvantages of the NISE method are its intrinsic high-temperature limit (k T with the system size, because the most timeconsuming part of the calculation is determining the matrix exponentials in Eq. (25). A comparison of alternative propagation schemes can be found in Ref. [69]. Another drawback is that the calculations have to be repeated for many disorder realizations in order to obtain an average over the bath fluctuations.

Haken-Strobl-Reineker (HSR) model
The HSR model [13,41] is based on the stochastic Liouville equation (SLE) description of coupled coherent and incoherent exciton dynamics [41]. The coherent part of the dynamics is characterized by the unperturbed Hamiltonian H 0 , i.e., Eq. (1) without fluctuations in the molecular transition energies. The incoherent part results from these fluctuations which, as stated above, are treated as Gauss-Markov stochastic processes described by Eqs. (11) and (12). The equation of motion for the density matrix element nm is then given by [13]: where n is a generic label for a molecule (either donor or any acceptor molecule). In this approach, the excited-state phase relations between two molecules decay with rate 2 , which destroys the back and forth oscillations in the occupation probabilities nn of individual molecules, either on a time step fast compared to EET between both molecules involved (incoherent limit, J nm ) or slow (J nm ). The main advantage of the HSR model is its simplicity in implementation. However, like the NISE method, a major limitation of the HSR model is its high-temperature limit: k T B and k T W B . Furthermore, the HSR model is naturally limited to the fast fluctuation limit, 1. The computation effort for the HSR model scales approximately as O N ( ) A 6 with the system size, because it requires propagating the density matrix. However, in contrast to the NISE method, all the dynamic disorder is automatically included and no additional averaging over realizations of the fluctuations is needed as long as only dynamic disorder is considered.

Hierarchy of Equations of Motion (HEOM)
The HEOM method is based on the exact equations of motion for a system with bilinear coupling to a quantum-mechanical harmonic oscillator bath. A description of the implementation of this method for a bath of overdamped Brownian oscillators can be found in Ref. [46]. The method is numerically exact as long as the bath can be described by a spectral density of independent quantum harmonic oscillators.
We use the PHI software package [46,70] to calculate HEOM exciton dynamics and obtain P t ( ) D , from which the EET rate constant was extracted (see Section 3.1). For the system parameters studied here, it turned out not necessary to include Matsubara frequencies [15]. We use a hierarchy depth of 11, which was tested for convergence by comparing to results with higher hierarchy depth.

Results and discussion
In this study, the acceptor complex is a ring with the number of molecules, N A , ranging from 1 to 20 and the molecular transition dipoles oriented perpendicularly to the plane of the ring (Fig. 1). Such configuration of transition dipoles corresponds to an H-type aggregate. Within this ring, only nearest-neighbor electronic interactions, > J ( 0) A , are considered. We consider two locations for the molecular donor: (a) "symmetric system", where the donor is placed in the center of the ring (Fig. 1a), resulting in equal couplings for all donor-acceptor pairs, J n D A = J; (b) "asymmetric system", where the donor is placed outside the ring (Fig. 1b) and is coupled only to the nearest acceptor molecule, again with coupling strength denoted as J. These configurations resemble typical chromophore arrangements as found for transfer within and between natural light-harvesting aggregates, while preserving a high degree of simplicity to ease the interpretation and facilitate comparison of the methods considered. The parameters are chosen to ensure a strongly delocalized acceptor system and a weak coupling between the donor and the acceptor system. All energetic quantities are expressed in units of the electronic donor-acceptor coupling J. Furthermore, we either use = In the following subsections, we present and discuss the results for different regimes of the system and bath parameters, all of which are summarized in Table 2. ) and fast modulation dynamics ( 1). In this range, all methods are expected to work well and thus agree with each other. We examine this for the example of the symmetric unbiased donor-acceptor system with = N 4 A (Fig. 1a). In Fig. 2, the rate constant, calculated with the different methods, are presented as a function of . As can be seen from the plot, all methods apart from the MC-FRET diag, IPR agree. The discrepancy for the MC-FRET diag, IPR is attrib-

Table 2
System and bath parameters used in the calculations: the nearest-neighbor coupling in the acceptor (J A ), the donor-acceptor coupling (J D A ), the number of molecules in the acceptor ring (N A ), the temperature (k T B ), the inverse of the correlation time ( ), and the reorganization energy ( ).

Intermediate regime (Section 4.3)
Real system parameters ( Fig. 9a and 10a) 300 20 1-20 200 333.6 100 High temperature (Fig. 9b and 10b) 300 20 1-20 5000 333.6 100 High temperature, small ( Fig. 9c and 10c) 300 20 1-20 5000 333. 6 4 uted to the fact that the effect of the exchange narrowing, taken into account by a simple scaling of the line shape using the participation ratio, is considerably overestimated by using the IPR for the ordered acceptor system. While it was shown that the effect of exchange narrowing, often seen in the static-disorder limit, is also found in the fast fluctuation limit [71], the effect is not as prominent as it is in the staticdisorder limit. Closer inspection of the graph reveals that the rate constant calculated with the HSR model (green curve) deviates slightly as decreases. This behavior is expected, because smaller corresponds to slower dynamics, i.e., moving out of the validity regime of the HSR model. We note that the current high-temperature, fast-fluctuation limit case serves as a good validation of the implementation of the different methods.
The fact that the maximum rate constant is reached at finite ( 7) can be explained in terms of the optimal overlap between the emission line shape of the donor molecule and the line shape corresponding to the absorption matrix of the acceptor complex. When is too large, the line widths become too narrow, resulting in partial overlap. As decreases, the width gets larger, as does the overlap, thereby resulting in a large rate constant. This overlap will reach an optimum however, because when the line shapes become too broad, the absolute value of the overlap gets smaller again, because the peak heights decrease.
We now consider the effect of the size of the acceptor system on the rate constant for one fixed value of = 11.2, again for the symmetric unbiased donor-acceptor system, and compare all methods to the HEOM standard (Fig. 3). As can be seen, all methods apart from the MC-FRET diag, IPR , are in excellent agreement. The deviation for the MC-FRET diag, IPR , as discussed previously, is attributed to the overestimation of the effect of exchange narrowing. Due to the symmetry, the accepting state of the acceptor complex considered here is predominantly the totally symmetric ( = k 0) exciton state, which for an H-aggregate is the state with the highest energy. For the unbiased donor-acceptor system, the overlap between the donor emission line shape and the absorption line shape for = N 1 A is maximal, because the monomer energy of the donor and acceptor is the same and there is no Stokes shift. The overlap decreases for = N 2, 3 A due to the lower overlap between the monomer donor line shape and the excitonic acceptor states, because the latter split and the accepting state lies higher in energy than the donor excited state. When N A increases beyond 4, the energy of the excitonic acceptor state that overlaps with the donor state hardly changes anymore and the rate constant only depends on the donor-acceptor coupling which increases linearly with N A . This results in the linear behavior for the rate constant beyond = N 4 A observed in Fig. 3.

High-temperature and slow-modulation limit
Next, we consider the regime of high temperature and move out of the fast-fluctuation regime ( of the order of unity or smaller). The rate constant, calculated with the different methods applied to the symmetric system and compared to the HEOM method, are presented as a function of in Fig. 4. As can be seen from the plot, the methods deviate significantly from each other in this regime. First, when approaching the static limit, i.e., , 0, the rate constant for the HSR method vanishes. This is in contrast to all the other methods and is a clear sign of the limitation of the HSR model to the white noise limit. The MC-FRET diag, IPR consistently overestimates the rate predicted by the other methods. This can once more be attributed to overestimating the effect of exchange narrowing, resulting in too narrow lines, and thus higher rate constant. The NISE and MC-FRET diag methods agree with the HEOM method over the whole range of , indicating that these methods both work well in this regime. While the MC-FRET IAL gives results close to NISE and MC-FRET diag , it gives an unreasonable absorption line shape in the slow dynamics limit: instead of giving a Gaussian with a FWHM , which would be expected in the slow-modulation limit, Eq. 17 results in two Lorentzians of FWHM when decreases. This discrepancy is illustrated in Fig. 5, where the MC-FRET IAL absorption line shape of the tetrameric acceptor is plotted together with the MC-FRET diag absorption line shape and the emission line shape of the  As the dynamics gets slow, and decrease. In the static limit, the line width is proportional to and thus independent of (see Eq. 10), and therefore it has a finite value. That is why the rate constant for all methods except the HSR model approaches a constant when approaches zero. For the HSR model, the rate approaches zero in the static limit, because the line shape in this method is given by a Lorentzian of width 2 . When approaches zero, gets very large (see Eq. 12) resulting in vanishing overlap, which in turn leads to a vanishing rate constant.
We now consider the effect of the size of the acceptor system on the rate constant calculated for biased symmetric systems and for one fixed value of , namely = 0.4. The results are presented in Fig. 6. Since = 0.4 corresponds to the regime of slow dynamics, HSR is not expected to be valid, and indeed a large deviation from HEOM is observed for all the graphs. From Fig. 6a ( = E 0), it can be seen that the absolute deviation from HEOM for all methods is increasing with growing N A . However, the relative error is independent of N A (see Fig. S1 of SI). Having an energy bias between the donor and the acceptor molecules results in the shift of the excitonic energy levels relative to the donor excited state and hence influences the possible resonances. Therefore, by tuning the energy bias, the energy transfer can be either enhanced or suppressed. The energy bias also affects the agreement between the methods as can be seen especially for the case of the largest bias (Fig. 6c). The results can be understood from the fact that the absorption line shape differs for the considered methods: with no energy bias the overlap between the absorption and emission line shapes is largest due to the resonance (effect), hence discrepancies between the methods are largest; with a large energy bias, the overlap between the two line shapes is only partial, hence discrepancies are smaller, since the differences in the absorption line shapes are not as crucial as in the resonant case.
Next, we increase the nearest-neighbor interactions in the acceptor ring to = J J 150 A compared to the donor-acceptor coupling, = J J D A , resulting in a much larger shift of the highest excitonic energy level for the acceptor. The other parameters are the same as in Fig. 6. From  Fig. 7, it can be seen that the larger value of J A generally results in a smaller rate constant than in Fig. 6, which directly results from the bad match between the energy of the donor with the acceptor's exciton band. The bump observed for = N 2 A comes from the resonance effect discussed previously. This effect has more impact in systems of smaller size for larger coupling in the acceptor ring.
Finally, we test the methods for the asymmetric system, where the donor molecule is coupled to only one of the acceptor molecules. We use the same parameters as for the symmetric system. Fig. 8 shows that all methods are in good agreement, except the HSR model, which fails to describe the slow dynamics of the bath. Thus, either of the MC-FRET or the NISE methods can be used in this case. It should be noted that unlike the symmetric case where the EET is dominated by transfer to the totally symmetric ( = k 0) exciton state, in the asymmetric system all states within the exciton band participate. For small N A (<10), resonance effects appearing and disappearing due to the discrete positions of the exciton levels in the band, give rise to the observed  oscillations. By contrast, for the symmetric system such oscillations with N A do not occur, as the predominant accepting state is the symmetric state, at fixed energy . For large N A , the rate constant almost does not change, because the addition of new exciton states does not change much the density of states and, therefore, hardly influences the overlap between the acceptor's absorption line shape and the donor's emission line shape. Furthermore, in the symmetric case, there are N A couplings of magnitude J, resulting in the linear increase of the rate constant with the size of the acceptor's system. This is different from the asymmetric case, where there is only one coupling of magnitude J, and as a result the rate constant does not change with the increase of N A .

Intermediate regime
The intermediate regime corresponds to the bath dynamics of many known light-harvesting systems, like, for example, LH2 [20,21], or selfassembled C8S3 nanotubes [28,30]. We chose the parameters close to those used previously for the LH2 system [51,72,73] (see Table 2) and apply them first to the symmetric system. Note that in this subsection, we use absolute units (cm −1 and ps −1 ). The results for the rate constant in this parameter regime as a function of N A are presented in Fig. 9a   methods predict too low rates for monomer to monomer transfer. MC-FRET diag and MC-FRET IAL give similar results and their predictions for the rate lie in between those for the HSR and the MC-FRET IAL methods. For better understanding the predicted results, we calculate the rate constant for the same system, except that we use the high-temperature limit ( = k T 5000 B cm −1 ), corresponding to = 0.33 (Fig. 9b), or alternatively use the high-temperature limit ( = k T 5000 B cm −1 ) and a low reorganization energy ( = 4 cm −1 , corresponding to = 1.67) (Fig. 9c). The results show that NISE is in good agreement with HEOM in both cases and has the lowest relative error compared to other methods (see SI). This confirms that its deviation from HEOM in the regime of parameters relevant to LH2 mainly comes from the hightemperature limitation. Here, MC-FRET diag and MC-FRET IAL give a rate constant close to that of HEOM, while MC-FRET diag, IPR overestimates the rate constant in Fig. 9b and underestimates it in Fig. 9c. This unpredictability is explained by the narrowness of the DOS line shape scaled by the IPR for the acceptor system compared to the line shape predicted by the other MC-FRET methods.
Again a resonance effect is observed around = N 3 A and the same linear growth for large N A , as discussed previously. Finally, we consider the asymmetric system using the same parameters as for the symmetric system. The results are shown in Fig. 10. The observed resonance effect for the smaller acceptor size result from participation of all exciton states in the process of EET. Moreover, the overall rate constant for this configuration is lower than in the symmetric one, since in the latter there are more couplings, which increase the rate, as discussed previously.
The agreement between the methods is the same as for the symmetric system. However, the relative error is much lower. In this sense, varies from 1 to 20; = 333.6 cm −1 . (a) Regime of room temperature; (b) high-temperature limit; (c) high-temperature limit and small reorganization energy. Fig. 10. Donor-acceptor rate constant for the energy transfer in the asymmetric system. The parameters are the same as in Fig. 9.
the HSR model becomes a good approximation for approximately > N 10 A in the parameter regime of Fig. 10a and 10b. The good performance of the NISE and HSR models found for room temperature for the symmetric and asymmetric systems should be expected to deteriorate if either or is increased significantly compared to the temperature, as shown in Table 1. Therefore as temperature is lowered, these methods are expected to get worse both in the limits of fast and slow fluctuations. Similar effects will emerge in the MC-FRET based methods unless absorption and emission matrices correctly including these effects are used. Furthermore, if one were to generalize the methods to multichromophoric donor complexes, at low temperature the emission from the donor system will be collected in the lowest levels, breaking the high-temperature approximation of these methods.

Conclusions
In this paper, we compared different theoretical methods to study EET in molecular aggregates. The methods considered were the multichromophoric Förster resonance energy transfer (MC-FRET) approach (with three different approximations used for calculating the absorption line shape of the acceptor complex), the Numerical Integration of the Schrödinger Equation (NISE), and the Haken-Strobl-Reineker (HSR) model. We applied them to a system consisting of a monomeric donor coupled to a multichromophoric acceptor ring in a symmetric or asymmetric way. We tested the methods in a number of regimes and for certain parameter ranges we compared them against the "gold standard" Hierarchy of Equations of Motion (HEOM) method.
First, we verified that in the high-temperature and fast-modulation limit, all methods performed accurately, apart from MC-FRET diag, IPR that failed to give correct results due to its overestimation of the exchange narrowing effect.
In the high-temperature and slow dynamics limit, the methods deviated considerably from each other. HSR then was clearly seen to fail as a result of the white-noise limit. MC-FRET diag, IPR significantly overestimated the rate predicted by the other methods, again due to an overly pronounced exchange narrowing effect. MC-FRET IAL , which uses the expression for the ideal absorption line shape based on the quantum master equation formalism [59,60], gave an unreasonable absorption line shape despite giving reasonable results for the EET rate. We found that MC-FRET diag , which uses the diagonal approximation for the absorption line shape, agrees with NISE over the whole range of the inverse of the correlation time of the bath fluctuations, . This suggests that MC-FRET diag works well in this limit, since the NISE method has no limitations in this limit (its only limitation is its inherent high-temperature assumption).
In the intermediate regime, with bath dynamics comparable to that of many known natural light-harvesting systems, NISE gives the most reasonable results, though it systematically overstimates the rates predicted by HEOM. In this parameter regime, when looking at the high temperature limit, either MC-FRET diag , MC-FRET IAL , or NISE gives correct results. In the intermediate regime of the bath fluctuations, an accurate prediction of the line shape is at least as important as temperature effects. This explains the variability of the performance of different MC-FRET methods.
In the symmetric case, the donor predominantly couples to a single state of the acceptor complex -a totally symmetric exciton state ( = k 0). This makes the description of the line shape of this state very critical: small changes in this line shape result in large variations of the rate constant. This explains the large sensitivity of the rate constant to the method used. In contrast, in the asymmetric case the donor couples to the whole exciton band. As a consequence, the line shape (of a single state) is not as critical anymore. This results in lower relative errors and hence less crucial dependency of the rate constant on the chosen method.
The MC-FRET methods showed to be by far the most computationally efficient. However, one should take into account that the description of the line shape is crucial here. Studying exciton dynamics and energy transfer from a large multichromophoric donor to a large multichromophoric acceptor remains a big challenge. This can only be solved by implementing an accurate and efficient way of calculating the emission line shape function of the involved multichromophoric donor complexes. On the other hand, NISE being computationally more demanding, shows the most reasonable results despite its high-temperature limitation. In contrast to HEOM, it is straightforward to efficiently parallelize the NISE method.