Self-consistent approach to the dynamics of excitation energy transfer in multichromophoric systems
Veljko Janković Tomáš Mančal





 Institute of Physics Belgrade, University of Belgrade, Pregrevica 118, 11080 Belgrade, Serbia
 

 Faculty of Mathematics and Physics, Charles University, Ke Karlovu 5, 121 16 Prague 2, Czech Republic

§ ABSTRACT 

Computationally tractable and reliable, albeit approximate, methods for studying exciton transport in molecular aggregates immersed in structured bosonic environments have been actively developed.
Going beyond the lowest-order (Born) approximation for the memory kernel of the generalized quantum master equation typically results in complicated and possibly divergent expressions.
Starting from the memory kernel in the Born approximation, and recognizing the quantum master equation as the Dyson equation of the Green's functions theory, we formulate the self-consistent Born approximation to resum the memory-kernel perturbation series in powers of the exciton–environment interaction.
Our formulation is in the Liouville space and frequency domain, and handles arbitrary exciton–environment spectral densities.
In a molecular dimer coupled to an overdamped oscillator environment, we conclude that the self-consistent cycle significantly improves the Born-approximation energy-transfer dynamics.
The dynamics in the self-consistent Born approximation agree well with solutions of hierarchical equations of motion over a wide range of parameters, including the most challenging regimes of strong exciton–environment interactions, slow environments, and low temperatures.
This is rationalized by analytical considerations of coherence-dephasing dynamics in the pure-dephasing model.
We find that the self-consistent Born approximation is good (poor) at describing energy transfer modulated by an underdamped vibration resonant (off-resonant) with the exciton energy gap.
Nevertheless, it reasonably describes exciton dynamics in the seven-site model of the Fenna–Matthews–Olson complex in a realistic environment comprising both an overdamped continuum and underdamped vibrations.


Self-consistent approach to the dynamics of excitation energy transfer in multichromophoric systems
    Tomáš Mančal
    September 9, 2024
===================================================================================================





§ INTRODUCTION

The light absorption and thus initiated excitation energy transfer (EET) in molecular aggregates constitute the first steps of solar-energy conversion in both natural <cit.> and artificial <cit.> systems.
The EET takes place in a complex dynamic spatio–temporal landscape stemming from the competition of interactions promoting exciton delocalization (resonance coupling between molecules) and localization (static and dynamic disorder). <cit.>
As the energy scales of these counteracting interactions are typically comparable to one another, <cit.> theoretical descriptions of EET dynamics are quite challenging.

When approached from the perspective of the theory of open quantum systems, <cit.> the challenge transforms into describing non-Markovian quantum dynamics of excitons interacting with their environment.
As standard theories (such as Redfield <cit.> and Förster <cit.> theories) do not meet this challenge, <cit.> various numerically exact methods have been developed.
Two most common foundations of these are (i) the Feynman–Vernon influence functional theory, <cit.> and (ii) the Nakajima–Zwanzig [time-convolution (TC)] generalized quantum master equation (GQME). <cit.>
The approaches rooted in (i) include the hierarchical equations of motion (HEOM) <cit.> and a host of path-integral and process-tensor-based methods. <cit.>
Each of these approaches develops a different representation of the so-called exact reduced evolution superoperator (or the dynamical map) 𝒰(t), which becomes its central object.
On the other hand, the main aim of the approaches originating from (ii) is to evaluate the exact memory-kernel superoperator 𝒦(t). <cit.>
All the above-referenced approaches are in general computationally intensive.
Their applications to realistic EET models, which feature a larger number of chromophores and/or structured spectral densities (SDs) of the exciton–environment interaction extracted from experiments or atomistic simulations, <cit.> can thus be impractical.
A need for computationally less demanding and reliable, although approximate, approaches to study EET dynamics in realistic multichromophoric models cannot be overemphasized.

Recent years have witnessed the emergence of many such approaches, among others quantum–classical methods, <cit.> and related quantum-chemical approaches to nonadiabatic dynamics (e.g., the surface hopping). <cit.>
One can also devise approximations to 𝒰 and 𝒦 starting from the formally exact expressions of (i) and (ii), respectively.
For example, the former give rise to cumulant-expansion-based approaches, <cit.> while the latter result in various approximations to the exact memory kernel. <cit.>
The authors of Ref. PhysRevE.86.011915 noted that already the second-order TC (TC2) GQME, which retains only the lowest-order (second-order) term 𝒦^(2)(t) in the expansion of 𝒦 in powers of the exciton–environment interaction, can qualitatively reproduce the main features of EET dynamics in multichromophoric systems.
The most obvious improvement over the TC2 GQME, also known as the (second) Born approximation (BA), <cit.> is to retain some of the higher-order contributions to 𝒦. 
Such a route generally leads to involved expressions that might not always be convergent, and that can be practically evaluated only for the simplest models, where they show some improvements over the original theory. <cit.>
While already the lowest-order kernel of polaron-transformed QMEs <cit.> can significantly improve over TC2 GQME, decent results can also be obtained by partially resumming the perturbation expansion for 𝒦 using only low-order terms as the input.
Partial resummations are most commonly performed in the frequency domain, and one usually applies Padé or Landau–Zener resummation schemes to 𝒦^(2)(ω) and 𝒦^(4)(ω). <cit.>
Another possibility, the self-consistent resummation schemes strongly rooted in condensed-matter physics, <cit.> has received even less attention in this context.

Successful self-consistent improvements over the TC2 GQME have been reported in the context of quantum theory of electronic transport through molecular junctions. <cit.>
These approaches have mainly considered purely electronic quantum transport, when the role of the bath is played by electronic leads.
Their applications to the quantum transport in the presence of both electronic leads and molecular vibrations are much more recent and scarcer. <cit.>
A proposal for a self-consistent description of the dynamics of a general open quantum system has appeared only very recently. <cit.>
It leans on a diagrammatic representation of the perturbation series for 𝒰(t), which has indeed appeared in the context of open quantum dynamics, see, e.g., Refs. PhysRevA.41.6676,ChemPhys.347.185,JChemPhys.153.244122.
The novelty of the approach by Scarlatella and Schirò <cit.> lies in the subsequent formulation of the perturbation series for 𝒦(t), which is related to 𝒰(t) via the TC GQME.
In the context of the theory of Green's functions in quantum many-body systems, <cit.> the memory kernel 𝒦(t) is analogous to the self-energy, whereas the TC GQME plays the role of the Dyson equation.
Considering the zero-temperature dynamics of the spin–boson model, Scarlatella and Schirò <cit.> find that performing the self-consistent cycle starting from the BA memory kernel 𝒦^(2) produces promising results.
However, they consider only the zero temperature and SDs commonly used when studying the spin–boson model in its narrowest sense. <cit.>
Moreover, their time-domain formulation of the self-consistent Born approximation (SCBA) involves solving an integro-differential equation in each iteration of the cycle.

This study explores the applicability of the self-consistent memory-kernel resummation scheme introduced in Ref. scarlatella2023selfconsistent to EET dynamics in multichromophoric aggregates.
For the sake of completeness, we first reconsider the theoretical developments of Ref. scarlatella2023selfconsistent, and make them technically simpler by working in the Liouville space and replacing the above-mentioned integro-differential equation in the time domain by a matrix equation in the frequency domain.
In a molecular dimer coupled to an overdamped phonon environment, we find that the self-consistent cycle greatly improves the original BA (TC2 GQME) over a wide range of dimer parameters.
This success of the SCBA is rationalized by considering the analytically tractable example of coherence dephasing in the pure-dephasing model.
Remarkably, we find that the SCBA remains reliable even in the generally difficult regimes of strong exciton–environment interactions and/or slow environments, when multiply environmentally assisted processes dominate the dynamics, <cit.> as well as at low temperatures.
Considering exciton dynamics modulated by a single underdamped vibrational mode, we conclude that the SCBA delivers good results only when the vibrational energy is resonant with the exciton–energy gap.
Upon including overdamped phonon continuum on top of a number of underdamped modes, we find that the SCBA delivers decent results for EET dynamics in the seven-site model of the FMO complex interacting with the realistic environment extracted from atomistic simulations.

The paper is organized as follows.
Section <ref> introduces our theoretical framework, whose applicability to EET dynamics is assessed in Sec. <ref>.
Section <ref> summarizes our main findings and discusses prospects for future work.



§ THEORETICAL CONSIDERATIONS




 §.§ Model

We model EET dynamics in an aggregate composed of N chromophores using the Frenkel–Holstein Hamiltonian

    H=H_S+H_B+H_S-B.

We take into account only the ground and one excited state on each chromophore, so that the purely electronic part H_S reads as

    H_S=∑_nε_n|n⟩⟨ n|+∑_m≠ nJ_mn|m⟩⟨ n|.

Here, |n⟩ is the collective singly excited state residing on chromophore n (all the other chromophores are unexcited), ε_n is the energy of the vertical transition from the ground state to the excited state on chromophore n, while J_mn is the resonance (typically dipole–dipole) coupling between chromophores m and n.
The aggregate is in contact with the environment modeled as collections of mutually independent harmonic oscillators associated with each chromophore:

    H_B=∑_nξω_nξb_nξ^† b_nξ.

Here, ω_nξ is the frequency of oscillator ξ on chromophore n, while bosonic operators b_nξ^† (b_nξ) create (annihilate) the corresponding oscillation quantum and obey [b_nξ,b_mξ']=0,[b_nξ,b_mξ'^†]=δ_nmδ_ξξ'.
The exciton–environment interaction is of the Holstein type, i.e., local and linear in oscillator displacements:

    H_S-B=∑_nξg_nξV_n(b_nξ+b_nξ^†),

where V_n=|n⟩⟨ n|.
Its strength, determined by the interaction constants g_nξ, is more conveniently described using the reorganization energy on chromophore n

    λ_n=∑_ξg_nξ^2/ω_nξ=∫_-∞^+∞dω/2π𝒥_n(ω)/ω,

where the SD of the exciton–environment interaction

    𝒥_n(ω)=π∑_ξ g_nξ^2[δ(ω-ω_nξ)-δ(ω+ω_nξ)]

is typically a continuous function of ω.
We focus on exciton dynamics starting from the factorized initial condition

    W(0)=ρ(0)ρ_B^eq,

where ρ(0) is the excitonic reduced density matrix (RDM) at the initial instant t=0, while

    ρ_B^eq=e^-β H_B/Tr_B e^-β H_B

represents the state of the environment (equilibrium at temperature T=β^-1) with no excitons present.
Within the Condon approximation, <cit.> this choice of the initial condition leads to the dynamics that can be probed in ultrafast nonlinear spectroscopies. <cit.>
As our formalism deals with the Green's function, we can still reconstruct the dynamics under arbitrary excitation condition as long as the light–matter interaction is weak and the Condon approximation is valid. <cit.>   
The only environmental quantity influencing the reduced excitonic dynamics starting from Eq. (<ref>) is the displacement autocorrelation function <cit.>

    C_n(t)   =∑_ξ g_nξ^2Tr_B{[b_nξ(t)+b_nξ^†(t)](b_nξ+b_nξ^†)ρ_B^eq}
       =∫_-∞^+∞dω/πe^-iω t𝒥_n(ω)/1-e^-βω.

The time dependence in Eq. (<ref>) is taken with respect to H_B, i.e., b_nξ(t)=e^-iω_nξtb_nξ.



 §.§ Real-time diagrammatic representation: Green's superoperator and self-energy superoperator

The assumption embodied in Eq. (<ref>) permits us to define the exact reduced evolution superoperator (dynamical map) 𝒰(t) by

    |ρ(t)=𝒰(t)|ρ(0).

Equation (<ref>) is formulated in the Liouville space, <cit.> in which the excitonic RDM at t>0, represented by the operator ρ(t) in the Hilbert space, becomes the N^2-component vector |ρ(t).
The reduced evolution superoperator is then a tetradic quantity comprising N^4 entries e'_2e'_1|𝒰(t)|e_2e_1, where {|e⟩} is an arbitrary basis in the single-exciton manifold. 
𝒰(t) is formally expressed as (see, e.g., Ref. PhysRevB.82.235307 and references therein)

    𝒰(t)=e^-iℒ_StTr_B{𝒯_te^-i∫_0^t ds ℒ_S-B^(I)(s)ρ_B^eq},

where the Liouvillian ℒ_a associated with the term H_a (a=S,B,S-B) in Eq. (<ref>) is defined by its action on an arbitrary Lioville-space vector |O corresponding to the Hilbert-space operator O, ℒ_a|O↔ [H_a,O], while the interaction-picture counterpart of ℒ_S-B is

    ℒ_S-B^(I)(t)=e^i(ℒ_S+ℒ_B)tℒ_S-Be^-i(ℒ_S+ℒ_B)t.

The time-ordering sign 𝒯_t imposes the chronological order (latest to the left) among the Liouvillians in the expansion of Eq. (<ref>) in powers of the exciton–environment interaction.
The Liouville-space approach adopted here is somewhat different from standard real-time approaches to the RDM of a particle in oscillator environment, <cit.> which consider the forward and backward Hilbert-space evolution operators separately.
Dealing with superoperators, we simplify the formalism, as we consider only the forward evolution in the Liouville space. <cit.> 
The average in Eq. (<ref>) can be performed term-by-term using the Wick's theorem, <cit.> and only even-order powers in ℒ_S-B remain.
The order 2k (k≥ 1) consists of (2k-1)!! terms.
Further manipulations usually proceed in two different manners.

The first possibility is to observe that, because of the 𝒯_t sign, all the terms appearing in any given order are mutually identical. <cit.>
Retaining the 𝒯_t sign, one obtains the well-known Feynman–Vernon expression <cit.>

    𝒰(t)=e^-iℒ_St𝒯_te^-Φ(t),

with

    Φ(t)=∑_n∫_0^t ds_2∫_0^s_2ds_1 V_n^(I)(s_2)^×
       [C_n^r(s_2-s_1)V_n^(I)(s_1)^×+iC_n^i(s_2-s_1)V_n^(I)(s_1)^∘].

In Eq. (<ref>), the superoperators V^× and V^∘ are defined by the correspondences V^×|O↔[V,O] and V^∘|O↔{V,O} (anticommutator), respectively.
Equations (<ref>) and (<ref>) serve as the starting point for HEOM and path integral-based approaches.



The possibility we opt for here is to make the time ordering explicit in each term, in which case the terms in any given order appear as different.
Further developments are facilitated by considering the (retarded) evolution superoperator [the (retarded) Green's function or Green's superoperator] <cit.>

    𝒢(t)=-iθ(t)𝒰(t)

instead of 𝒰(t).
Each term in order 2k of the expansion of 𝒢(t) can be represented by a diagram [Fig. <ref>(d)] comprising a total of 2k+2 (2 terminal and 2k internal) instants

    s_2k+1=t≥ s_2k≥…≥ s_1≥ 0=s_0,

and k chromophore indices n_k,…,n_1 associated with k pairs selected from the set {s_2k,…,s_1} by the application of the Wick's theorem.
The consecutive instants s_j+1 and s_j (0≤ j≤ 2k+1) of a diagram are connected by a straight line directed from s_j to s_j+1 [Fig. <ref>(d)] and representing the (retarded) free-exciton Green's superoperator [Fig. <ref>(b)]

    𝒢_S(t)=-iθ(t)e^-iℒ_St.

The directed double-line represents 𝒢(t) [Fig. <ref>(a)].
Dashed circumferences connect instants s_l and s_j (s_l≥ s_j) that are paired by the Wick's theorem, see Fig. <ref>(c).
These instants are accompanied by the same chromophore index n_jl [see Eq. (<ref>)] and the following superoperators [cf. Eq. (<ref>)]:

    s_l↔ V_n_jl^×,


    s_j↔ C_n_jl^r(s_l-s_j)V_n_jl^×+iC_n_jl^i(s_l-s_j)V_n_jl^∘.

Having placed all the superoperators in a chronologically ordered string (reading diagrams from left to right), one performs integrations ∫_0^t ds_2k…∫_0^s_2ds_1 over 2k internal instants, and summations over k independent chromophore indices.
Any diagram is either reducible or irreducible, depending on whether it can be cut in two by cutting a free Green's superoperator line or not. <cit.>
Examples of irreducible diagrams in Fig. <ref>(d) are diagrams (2), (4.2), and (4.3), while the diagram (4.1) is reducible.
Amputating the external lines corresponding to 𝒢_S(s_1) and 𝒢_S(t-s_2k) of irreducible diagrams, one obtains the diagrammatic representation of the so-called (retarded) self-energy superoperator, see Fig. <ref>(f).
The self-energy superoperator thus consists of all amputated diagrams that cannot be cut in two by cutting a free Green's superoperator line.
One can then derive the following Dyson equation: <cit.>

    𝒢(t)=𝒢_S(t)
       +∫_-∞^+∞ ds_2∫_-∞^+∞ds_1 𝒢_S(t-s_2)Σ(s_2-s_1)𝒢(s_1),

whose diagrammatic representation is shown in Fig. <ref>(e).
All the superoperators in Eq. (<ref>) are retarded, and the integrals in Eq. (<ref>) can run from -∞ to +∞, which facilitates the transition to the frequency domain.
Rewriting Eq. (<ref>) as an equation for 𝒰(t), inserting the result thus obtained into Eq. (<ref>), and differentiating with respect to t, one obtains the GQME

    ∂_t|ρ(t)=-iℒ_S|ρ(t)-∫_0^t ds 𝒦(t-s)|ρ(s),

where the relation between the memory kernel 𝒦 and the retarded self-energy Σ is

    Σ(t)=-iθ(t)𝒦(t).

Equation (<ref>) differs from the standard GQME only by the representation of the memory kernel.
While the memory kernel is commonly expressed in terms of projection superoperators, here, we provide its representation as a perturbation expansion in the exciton–environment interaction.



 §.§ Frequency-domain description

Resummation techniques usually require transferring to the frequency space.
This is most conveniently done starting from Eq. (<ref>) and forming the dynamical equation for 𝒢(t)

    ∂_t𝒢(t)=-iδ(t)-iℒ_S𝒢(t)-i∫_-∞^+∞ds Σ(t-s)𝒢(s).

With the standard definition of the frequency-dependent quantities

    𝒢(t)=∫_-∞^+∞dω/2π e^-iω t𝒢(ω),

Eq. (<ref>) becomes the following algebraic equation:

    𝒢(ω)=[ω-ℒ_S-Σ(ω)]^-1.


The Hermitian property of the RDM implies that 𝒢(ω) satisfies

    e'_2 e'_1|𝒢(ω)|e_2e_1=- e'_1e'_2|𝒢(-ω)|e_1e_2^*.

This property shows that it is sufficient to compute 𝒢(ω) only for ω≥ 0, while the values for negative frequencies follow from this symmetry property.
We additionally note that the same symmetry property also characterizes the inverse Green's superoperator 𝒢^-1(ω), as well as the self-energy superoperator Σ(ω), which now carries the dimension of energy.

The frequency-domain diagrammatic representations of 𝒢(ω) and Σ(ω) appear the same as the time-domain representations in Figs. <ref>(d) and <ref>(f), respectively.
The general rules for translating diagrams into formulae can be inferred from the discussion in Sec. <ref>.



 §.§ Born and Redfield approximations


The lowest, second-order approximation to Σ is known as the (second) BA. <cit.>
The corresponding self-energy superoperator, shown in Fig. <ref>(a1), reads as (see also Refs. JChemPhys.132.194111,PhysRevA.96.032105)

    Σ_BA(t)=∑_nV_n^×𝒢_S(t)[C_n^r(t)V_n^×+iC_n^i(t)V_n^∘],

while its frequency-domain counterpart is

    Σ_BA(ω)=∑_n∫_-∞^+∞dν/2π𝒥_n(ω-ν)
       {(β(ω-ν)/2)V_n^×𝒢_S(ν)V_n^×+V_n^×𝒢_S(ν)V_n^∘}.


Upon instering 𝒦_BA(t)=iΣ_BA(t) into Eq. (<ref>), we obtain the well-known TC2 GQME <cit.>

    ∂_tρ(t)=-i[H_S,ρ(t)]
       -∑_n[V_n,∫_0^t ds C_n(t-s)e^-iH_S(t-s)V_nρ(s)e^iH_S(t-s)]
       +∑_n[V_n,∫_0^t ds C_n(t-s)^*e^-iH_S(t-s)ρ(s)V_ne^iH_S(t-s)].

For delta-correlated noise characterized by the dephasing rate Γ, C_n(t)=Γδ(t), Eq. (<ref>) assumes the Lindblad form

    ∂_tρ(t)=   -i[H_S,ρ(t)]
       -Γ∑_n(1/2{V_n^2,ρ(t)}-V_nρ(t)V_n).

We will use this result in Sec. <ref>.

The time-independent and nonsecular Redfield theory (see, e.g., Sec. 3.8.2 of Ref. May-Kuhn-book or Sec. 11.2 of Ref. Valkunas-Abramavicius-Mancal-book)

    ∂_t|ρ(t)=-iℒ_S|ρ(t)-ℛ|ρ(t),

is obtained by inserting Eq. (<ref>) into Eq. (<ref>), and then performing the Markov and adiabatic approximations.
These approximations result in the delta-like self-energy in the time domain, Σ_Red(t)=Σ_Redδ(t), the frequency-independent self-energy Σ_Red(ω)=Σ_Red, and the Redfield tensor ℛ=iΣ_Red [see Eq. (<ref>)], where 

    Σ_Red=∫_-∞^+∞ds Σ_BA(s) e^iℒ_Ss.

The Redfield theory is formulated in the exciton basis {|x⟩} defined through H_S|x⟩=ω_x|x⟩.
Equation (<ref>) then implies that the Redfield-tensor matrix elements read 

    x'_2x'_1|ℛ|x_2x_1=i x'_2x'_1|Σ_BA(ω_x_2-ω_x_1)|x_2x_1.

For delta-correlated noise, C_n(t)=Γδ(t), the Redfield equation [Eq. (<ref>)] also assumes the Lindblad form [Eq. (<ref>)]. <cit.>





 §.§ Self-consistent Born approximation

In this work, we improve upon the BA by replacing the free-exciton Green's superoperator 𝒢_S in Fig. <ref>(a1) or Eq. (<ref>) with the interacting Green's superoperator 𝒢, see Fig. <ref>(b1).
The resultant equation for the self-energy superoperator

    Σ(ω)=∑_n∫_-∞^+∞dν/2π𝒥_n(ω-ν)
       {(β(ω-ν)/2)V_n^×𝒢(ν)V_n^×+V_n^×𝒢(ν)V_n^∘}

is to be solved together with the Dyson equation [Eq. (<ref>)] in a self-consistent loop.
Namely, one starts from the free-exciton case, Σ^(0)(ω)=0, when Eq. (<ref>) gives the free-exciton Green's function (η→+0)

    𝒢^(0)(ω)=𝒢_S(ω)=[ω+iη-ℒ_S]^-1.

𝒢^(0)(ω) is then inserted into Eq. (<ref>) to yield Σ^(1)(ω)=Σ_BA(ω), which is then inserted into Eq. (<ref>) to yield 𝒢^(1)(ω), etc.
The above-described procedure is repeated until the difference between two consecutive iterations Σ^(k)(ω) and Σ^(k-1)(ω) for the self-energy (k≥ 1) becomes smaller than a prescribed numerical tolerance (see Sec. <ref> for more details).

Using either BA or SCBA for the self-energy superoperator, we do perform a partial resummation of the perturbation series for 𝒢 [Fig. <ref>(d)].
The resulting 𝒢_(SC)BA contains contributions from all orders in the exciton–environment interaction constant.
Nevertheless, the diagrammatic content of the SCBA is much richer than that of the BA, compare Fig. <ref>(b2) to <ref>(a2), suggesting that the SCBA is reliable in a much wider parameter range than the BA.
Still, the SCBA retains only a very limited subset of all possible diagrams appearing in Fig. <ref>(d), <cit.> and its reliability is to be carefully checked.



§ RESULTS




 §.§ Technical details


Equation (<ref>) features an infinitesimally small frequency η, which ensures the causality, i.e., 𝒢_S(t)=0 for t<0.
We always shift ω→ω+iη on the right-hand side of the Dyson equation [Eq. (<ref>)], which, apart from the causality, ensures that the matrix inversion in Eq. (<ref>) is numerically stable.
However, the numerical Fourier transformation of 𝒢(ω) thus obtained produces the exponentially damped Green's superoperator 𝒢(t)=e^-η t𝒢(t).
The results for the true Green's superoperator 𝒢(t)=e^η t𝒢(t) are thus the most reliable for t≪η^-1.
In all our computations, we set η=1 cm^-1, meaning that our results for exciton dynamics are bound to be reliable for t≪ 5 ps.
While we find that our results in the real-time domain are free of finite-η effects on time scales beyond η^-1, we always show only the initial 2–3 ps of exciton dynamics.

Generally, 𝒢(ω) slowly decays towards zero as |ω|→+∞.
The dominant component of the high-frequency tail of 𝒢(ω) can be inferred by performing one partial integration of 𝒢(ω)=∫_0^+∞dt e^i(ω+iη)t𝒢(t), which results in

    𝒢(ω)=1/ω+iη+𝒪(ω^-2),  |ω|→+∞.

Deriving Eq. (<ref>), we use 𝒢(t=0)=-i1, where 1 denotes the unit operator in the Liouville space.
The strongly pronounced high-frequency tail of 𝒢(ω) means that we have to consider many values of ω in order for the discrete Fourier transform to produce decent results in the time domain.
This is to be avoided as computing 𝒢(ω) and Σ(ω) involves inversion of an N^2× N^2 matrix [Eq. (<ref>)] and numerical integration [Eq. (<ref>)], respectively.
Defining 𝒯(ω)=1/ω+iη, and 𝒢^nt(ω)=𝒢(ω)-𝒯(ω), <cit.> the discrete (numerical) Fourier transformation of the non-tailed part 𝒢^nt(ω) produces 𝒢^nt(t), while the Fourier transformation of the high-frequency tail 𝒯(ω) can be performed analytically to yield 𝒯(t)=-iθ(t)e^-η t1.
Finally,

    𝒢(t)=-iθ(t)1+e^η t𝒢^nt(t).


As we assume local exciton–environment interaction [Eq. (<ref>)], 𝒢(ω) and Σ(ω) are most conveniently represented in the site basis {|n⟩}.
In all the examples to be discussed, we assume for simplicity that the environments of individual chromophores are identical.
The integration in Eq. (<ref>) expressing  Σ in terms of 𝒢 is performed using the ordinary trapezoidal rule.
Because of the symmetry property in Eq. (<ref>), we perform numerical integration of Eq. (<ref>) from 0 to ω_max with frequency step Δω.
We take care to correctly treat the contribution around ν=ω, which is finite for Ohmic SDs.
We stop the self-consistent cycle once we achieve δΣ^(k)≤ε_tol, where (k≥ 1)

    δΣ^(k)=max_n'_2n'_1
    n_2n_1|2 n'_2n'_1|Σ^(k)(ω)-Σ^(k-1)(ω)|n_2n_1/ n'_2n'_1|Σ^(k)(ω)+Σ^(k-1)(ω)|n_2n_1|,

while ε_tol is the desired (relative) accuracy.
We summarize the values of adjustable numerical parameters (η,ω_max,Δω,ε_tol) involved in our computations in Table <ref>.
Figures <ref>(a) and <ref>(b) respectively show δΣ^(k) as a function k when the SCBA is used on the dimer in the overdamped phonon continuum [the examples analyzed in Figs. <ref>(a)–<ref>(d)] and on the seven-site FMO model in the realistic environment [the examples analyzed in Figs. <ref>(b1) and <ref>(b2)].
Notably, the convergence of the self-consistent algorithm is achieved in a couple of tens of steps, even in a multichromophoric system immersed in a structured environment.
After its initial increase with k, starting from the value of 2 [see the text above Eq. (<ref>)], δΣ^(k) decreases in a power-law fashion for sufficiently large k. 





In both BA and SCBA, the trace of the RDM is preserved because of the outermost commutator in self-energies in Eqs. (<ref>) and (<ref>).
The RDM positivity is a much subtler issue, but we observe that, whenever SCBA improves over the BA, the results of both approximations conform to the positivity requirement on the time scales analyzed.

As we are primarily interested in examining the reliability of the BA and SCBA, we use

    |ρ(0)=1/N∑_n_2n_1|n_2n_1

as the initial condition in all numerical computations.
Apart from its simplicity, this initial condition provides a fairer assessment of the approximation performance than the widely used initial condition |n_0n_0, in which the exciton is placed at chromophore n_0.
In more detail, exciton dynamics in an arbitrary basis {|e⟩} is

    e_2e_1|ρ(t)=∑_n_2n_1
    n'_2n'_1⟨ e_2|n_2⟩⟨ n_1|e_1⟩×
        n_2n_1|𝒢(t)|n'_2n'_1 n'_2n'_1|ρ(0).

If the exciton is initially placed at chromophore n_0, its subsequent dynamics is determined by only N^2 matrix elements n_2n_1|𝒢(t)|n_0n_0 of 𝒢 out of the total of N^4 elements.
Quite generally, <cit.> the quality of approximate dynamics is different for different matrix elements of 𝒢, that is, for different starting chromophores n_0. 
Inserting the initial condition of Eq. (<ref>) into Eq. (<ref>), we can assess the overall approximation performance, which is effectively "averaged" over different matrix elements of 𝒢.



 §.§ Reliability of the SCBA: Asymmetric dimer


The advantages and shortcomings of the SCBA are most transparently identified on the simplest model relevant for EET, the molecular dimer.
In the site basis {|1_n⟩,|2_n⟩}, the exciton Hamiltonian H_S is represented by the matrix 

    H_S=[ Δε  J;  J  0 ],

where Δε is the site-energy gap, while J is the resonance coupling. 
The exciton state of lower (higher) energy is denoted as |1_x⟩ (|2_x⟩).

In Secs. <ref>–<ref>, we consider exciton dynamics in the featureless phonon environment described by the overdamped Brownian oscillator (OBO) SD <cit.> 

    𝒥_ph(ω)=2λ_phωγ_ph/ω^2+γ_ph^2,

where λ_ph is the reorganization energy, while γ_ph^-1 determines the environment-reorganization time scale.
Table <ref> summarizes the default values of model parameters, which are broadly representative of photosynthetic aggregates.
The performance of our approximations upon varying these parameters is studied in Secs. <ref> and <ref>. 


In Sec. <ref>, we assume the exciton interacts with a single underdamped vibrational mode, so that the interaction can be modeled by the underdamped Brownian oscillator SD <cit.>

    𝒥_vib(ω)=S_0ω_0[ωγ_0/(ω-ω_0)^2+γ_0^2+ωγ_0/(ω+ω_0)^2+γ_0^2],

where ω_0 is the vibrational frequency, S_0 the Huang–Rhys factor, while γ_0 is the relaxation rate.
Quite generally, γ_0≪ω_0 and S_0≪ 1.
The values used in benchmarks are taken from Ref. JChemPhys.157.095103 and listed in Sec. <ref>.



  §.§.§ Overdamped phonon continuum. Variations in excitonic parameters





In this section, we fix the exciton–environment interaction parameters λ_ph and γ_ph, and study the quality of SCBA and BA for different excitonic parameters J and Δε.
Figure <ref> provides an overall assessment of the performance of BA [Fig. <ref>(a)] and SCBA [Fig. <ref>(b)], and of the improvement of the SCBA over the BA introduced by the self-consistent cycle [Fig. <ref>(c)].
A convenient performance measure is the trace distance between the approximate [ρ_(SC)BA(t)] and numerically exact [ρ_HEOM(t)] RDM.
As the trace distance is time-dependent, while we also limit ourselves to the short-time dynamics (see Sec. <ref>), we quantify the performance of our approximations using the maximum trace distance over the time window [0,t_max]:

    𝒟((SC)BA|HEOM)=
       max_0≤ t≤ t_max1/2∑_k=1^2 |r_k^(SC)BA-HEOM(t)|,

where r_k^(SC)BA-HEOM(t) is the k-th eigenvalue of the operator ρ_(SC)BA(t)-ρ_HEOM(t).
We set t_max=2 ps in Eq. (<ref>).
The colorbar ranges in Figs. <ref>(a) and <ref>(b) suggest that the SCBA is generally a better approximation to the exact dynamics than the BA.
To quantify the improvement of the SCBA over the BA, in Fig. <ref>(c) we plot the ratio 𝒟(BA|HEOM)/𝒟(SCBA|HEOM), which we find to be greater than or equal to unity for all the pairs (J,Δε) examined.
The larger the ratio, the more pronounced the improvement of the SCBA over the BA.

Figures <ref>(a) and <ref>(b) show that the reliability of both BA and SCBA generally improves with increasing J and/or decreasing Δε.
This suggests that the approximations are best suited for relatively delocalized excitons, when the mixing angle θ∈[0,π/4], defined as tan(2θ)=2J/Δε, significantly deviates from zero.
Still, for Δε=0, when excitons are perfectly delocalized, the quality of both BA and SCBA increases with increasing J, i.e., decreasing λ_ph/J.
The impact of the exciton–environment interaction on the approximation reliability will be analyzed in detail in Sec. <ref>. 
Figure <ref>(c) reveals that the improvement of SCBA over BA is the most pronounced in the region of moderate resonance coupling 50 cm^-1≲ J≲ 150 cm^-1 and large site-energy gap Δε≳ 200 cm^-1.
The improvement is also appreciable for small resonance coupling, irrespective of the value of Δε, when both BA and SCBA perform relatively poorly, see Figs. <ref>(a) and <ref>(b).



  §.§.§ Overdamped phonon continuum. Analytical insights into the pure-dephasing model



For J=0, the model reduces to the pure-dephasing model, in which there is no population dynamics, but only dephasing of the initially present interexciton coherences.
The superoperators entering Eq. (<ref>) are then time independent, and combining Eqs. (<ref>)–(<ref>) we readily obtain the following exact expression for the reduced evolution superoperator: 

    n'_2n'_1|𝒢(t)|n_2n_1=-iθ(t)δ_n'_2n_2δ_n'_1n_1{δ_n_2n_1+ . 
    .
        (1-δ_n_2n_1)e^-iε_n_2n_1te^-2g^r(t)},

where g^r(t)=∫_0^t ds_2∫_0^s_2ds_1 C^r(s_1) is the real part of the lineshape function, whereas ε_n_2n_1=ε_n_2-ε_n_1 (for the dimer, ε_n_2n_1=±Δε).
The derivation of Eq. (<ref>), in which only g^r appears, crucially relies on our assumption that individual-chromophore environments are identical. 
While the exact coherence dynamics, which follows from Eqs. (<ref>), (<ref>), and (<ref>), can be recovered from the time-convolutionless second-order QME, <cit.> the corresponding BA and SCBA results remain only approximations to the exact solution, as both involve an explicit convolution in the time domain.
Still, relevant analytical insights concerning the (SC)BA can be obtained for the pure-dephasing model in the high-temperature limit 2π T≫γ_ph.
In Appendix <ref>, we derive that the matrix elements of the exact self-energy superoperator n_2n_1|Σ(ω)|n_2n_1 for n_2≠ n_1 have the following continued-fraction expansion (CFE):


    n_2n_1|Σ(ω)|n_2n_1=
        4λ_phTω-ε_n_2n_1+iγ_ph-2· 4λ_phTω-ε_n_2n_1+2iγ_ph-3· 4λ_phTω-ε_n_2n_1+3iγ_ph-⋯.


Truncating the CFE in the first layer, we obtain the BA result (n_2≠ n_1)

    n_2n_1|Σ_BA(ω)|n_2n_1=4λ_phT/ω-ε_n_2n_1+iγ_ph.

The same result follows from Eq. (<ref>) upon approximating coth(β(ω-ν)/2)≈2T/ω-ν, as appropriate in the high-temperature limit, and performing a contour integration by closing the contour in the upper half-plane.
In the pure-dephasing model, we can use Eqs. (<ref>) and (<ref>) to obtain the matrix elements of the self-energy superoperator in the Redfield theory

    n_2n_1|Σ_Red|n_2n_1=-i4λ_phT/γ_ph.

Appendix <ref> also demonstrates that the CFE of the SCBA self-energy reads


    n_2n_1|Σ_SCBA(ω)|n_2n_1=4λ_phTω-ε_n_2n_1+iγ_ph-4λ_phTω-ε_n_2n_1+2iγ_ph-4λ_phTω-ε_n_2n_1+3iγ_ph-⋯.



Comparing Eq. (<ref>) to Eq. (<ref>), we find that the CFE of the SCBA result can be obtained from the exact result by changing all the coefficients multiplying 4λ_phT in the CFE numerators to unity.
While this could suggest that the SCBA is, in general, a poor approximation to the exact solution, <cit.> we note that the broadening factors are the same at each CFE denominator in both the SCBA and exact self-energy.
One can then expect that the SCBA is superior to both the BA and the Redfield theory at reproducing the timescale of coherence dephasing.
This expectation is confirmed in Fig. <ref>(a) comparing different approximations to the coherence dynamics in a symmetric (Δε=0) pure-dephasing dimer.
The Redfield theory predicts an excessively fast coherence dephasing whose analytical form reads as [see Appendix <ref> and the circles in Fig. <ref>(a)]

    1_n2_n|ρ_Red(t)=1/2exp(-4λ_phT/γ_pht).

On the contrary, the BA predicts an excessively slow coherence dephasing that can be reasonably described by [see Appendix <ref> and the down-triangles in Fig. <ref>(a)]

    1_n2_n|ρ_BA(t)≈1/2cos(2√(λ_phT) t)e^-γ_pht/2.

The exact coherence-dephasing timescale is in between the results of the Redfield theory and BA, while the corresponding exact dynamics can be reasonably approximated by [see Appendix <ref> and the up-triangles in Fig. <ref>(a)]

    1_n2_n|ρ(t)≈1/2exp(-2λ_phTt^2).

The SCBA indeed reproduces the correct order of magnitude of the dephasing timescale, although it displays oscillations similar to that predicted by Eq. (<ref>).
The performance of different approximations can also be inferred from (the imaginary part of) the corresponding self-energy profile presented in Fig. <ref>(a).
The self-energy within the Redfield theory [Eq. (<ref>)] does not bear even a qualitative resemblance to the exact result.
The BA, SCBA, and the exact result all display a peak centered around ω=0.
The BA peak is the narrowest and highest and has a Lorentzian shape whose full width at half-maximum is determined by γ_ph only, see Eq. (<ref>).
The exact peak [Eq. (<ref>)] is much broader than the BA peak, while the width of the SCBA peak [Eq. (<ref>)] is somewhat smaller than, yet comparable to, the width of the exact peak.





  §.§.§ Overdamped phonon continuum. Variations in exciton–environment interaction parameters and temperature





Here, we fix exciton parameters Δε and J, and vary λ_ph and γ_ph to examine the reliability of BA [Fig. <ref>(a)] and SCBA [Fig. <ref>(b)], as well as the improvement over the BA brought about by the self-consistent cycle [Fig. <ref>(c)].
Overall, we find that the quality of both BA and SCBA improves with decreasing the reorganization energy and/or shortening the environment-reorganization timescale, see Figs. <ref>(a) and <ref>(b).
Figure <ref>(c) shows that the SCBA is always better at approximating the exact dynamics than the BA.
Within the range of values of λ_ph and γ_ph^-1 typically used in models of photosynthetic EET (25 cm^-1≲λ_ph≲ 250 cm^-1 and 50 fs≲γ_ph^-1≲ 200 fs), <cit.> the self-consistent cycle significantly improves BA results, see the dashed-line rectangle in Fig. <ref>(c).
Notably, for very fast environments (γ_ph^-1≲ 10 fs), we find that the quality of BA and SCBA is virtually the same, i.e., 𝒟(SCBA|HEOM)≈𝒟(BA|HEOM), for all reorganization energies examined.
To understand this, we observe that when γ_ph→∞ and T→∞ so that γ_ph/T→ 0, our model reduces to the well-known Haken–Strobl–Reineker white-noise model, <cit.> within which C_n(t)=Γδ(t), with Γ=2λ_phT/γ_ph.
Inserting this C_n(t) into Eq. (<ref>) for the exact evolution superoperator 𝒰(t), we conclude that the exact equation of motion for ρ(t) coincides with the Lindblad-like Eq. (<ref>). 
To reach this conclusion, we use that the white-noise assumption renders the time-ordering sign upon differentiating Eq. (<ref>) effective only on the superoperator e^-Φ(t) and ineffective on the superoperator ∂_tΦ(t).
Our discussion in Sec. <ref> then implies that the BA, and even the Redfield theory, become exact in the white-noise limit.
The exactness of the BA implies that the self-consistent cycle cannot improve on BA any further, meaning that the SCBA result is also exact in the white-noise limit.



Figure <ref>(c) shows that the improvement of the SCBA over the BA is the most pronounced for sluggish environments or large reorganization energies. 
Both observations are ultimately rooted in the richer diagrammatic content of the SCBA compared to that of the BA, see Fig. <ref>, and are remarkable as computing exciton dynamics modulated by strong exciton–environment interactions and/or slow environments has to take into account higher-order environmentally assisted processes. <cit.>
Similar challenges are encountered when studying exciton dynamics at low temperatures.
Additional insights into the quality of different approximations can be gained from Fig. <ref>, which compares approximate and numerically exact exciton dynamics for the default values of model parameters [Figs. <ref>(a1) and <ref>(a2)], slow environment [γ_ph=5 cm^-1, Figs. <ref>(b1) and <ref>(b2)], strong exciton–environment interaction [λ_ph=200 cm^-1, Figs. <ref>(c1) and <ref>(c2)], and low temperature [T=77 K, Figs. <ref>(d1) and <ref>(d2)]. 
The overall performance of approximate methods is essentially as discussed in Sec. <ref>.
The Redfield theory shows pronounced deviations from the exact result already on shortest time scales, whereas the BA and SCBA reproduce the very initial stages (t≲ 50 fs) of the exact dynamics quite well.
While the subsequent dynamics within the BA generally exhibits oscillatory features that are damped relatively slowly, cf. Fig. <ref>(a), the SCBA dynamics of exciton populations and the interexciton coherence follows the corresponding HEOM results very reasonably.
The SCBA approximates the true exciton-population dynamics in both realistically slow [Fig. <ref>(a1)] and excessively slow [Fig. <ref>(b1)] environments quite well, both on subpicosecond and somewhat longer time scales.
On the other hand, Figs. <ref>(a2) and <ref>(b2), as well as the insets of Figs. <ref>(a1) and <ref>(b1), suggest that the SCBA is not that good at reproducing the true interexciton-coherence dynamics.
Still, it does capture the correct long-time behavior (the imaginary part of the interexciton coherence tends to zero, and the real part tends to a non-zero value).
For strong interactions and at low temperatures, the subpicosecond dynamics of exciton populations within the SCBA is quite close to the true dynamics, see Figs. <ref>(c1) and <ref>(d1), while some quantitative differences between them appear on longer time scales (these are more pronounced for the higher-energy exciton state).
On the other hand, the SCBA dynamics of the interexciton coherence agrees well with the corresponding numerically exact dynamics, especially at lower temperatures, see Figs. <ref>(c2) and <ref>(d2) and the insets of Figs. <ref>(c1) and <ref>(d1).
Finally, the convergence of the self-consistent algorithm [Fig. <ref>(a)] generally slows down with increasing the interaction energy or the memory time scale, while its pace appears to be relatively weakly affected by temperature variations.



  §.§.§ Overdamped phonon continuum. Self-energy superoperator



We formulate our approximate approaches in the frequency domain, with the self-energy (memory-kernel) superoperator as their central quantity, and Figs. <ref>(a)–<ref>(d) discuss the reflections of the above-summarized time-domain observations on the frequency domain.
We choose the slow-environment regime analyzed in Fig. <ref>(b) and concentrate on the matrix elements 2_x2_x|Σ(ω)|2_x2_x and 1_x2_x|Σ(ω)|2_x2_x describing respectively the population flux out of the higher-energy exciton state and the population-to-coherence transfer from that state.
The results in Figs. <ref>(a)–<ref>(d) are obtained setting the artificial-broadening parameter to η=1 cm^-1, and we have checked that varying η over the range [0.5,5] cm^-1 does not qualitatively (and to a large extent quantitatively) affect the results presented here. 
Within the BA, the imaginary parts of both self-energy matrix elements display very narrow peaks centered around the exciton energy gap (at ω=±Δε_X), see Figs. <ref>(a) and <ref>(c), which is compatible with the oscillatory features of the BA in Fig. <ref>(b).
On the other hand, the peaks of the numerically exact profiles are much wider, and their centers are somewhat shifted from Δε_X.
While the SCBA profiles overall reasonably reproduce the numerically exact ones both in terms of peak positions and shapes, we note the tendency of the SCBA towards somewhat excessive peak shifting and narrowing.
The relative advantage of the SCBA over the BA is the most obvious for the population-to-coherence transfer, when the BA completely misses the peak appearing somewhat below Δε_X in both HEOM and SCBA profiles, see Fig. <ref>(c).

In Figs. <ref>(e) and <ref>(f), we compare the BA and SCBA results for the self-energy superoperator in the time domain to the corresponding numerically exact result.
As discussed in Sec. <ref>, these time-domain results do not depend on η.
The general tendencies observed in the exciton dynamics in Fig. <ref>(b) are also seen on the level of the time-dependent memory kernel.
Namely, the BA result for the time-dependent memory kernel displays pronounced weakly damped oscillatory features, while the SCBA reproduces the exact result reasonably well.
The quantitative agreement between the SCBA and HEOM results for the matrix element connected to population-to-population transfer is somewhat better than in the case of population-to-coherence transfer, compare Fig. <ref>(e) to Fig. <ref>(f).  
We finally note that obtaining the memory-kernel superoperator in the time domain is generally difficult. <cit.>
On the contrary, our time-domain results are readily obtained by the numerical Fourier transformation of the corresponding frequency-domain self-energies.





  §.§.§ A single underdamped vibrational mode




Here, using the values of J,Δε, and T summarized in Table <ref>, we benchmark the (SC)BA when the SD of the exciton–environment interaction is modeled using Eq. (<ref>).
We set γ_0=3 cm^-1, corresponding to the relaxation timescale of γ_0^-1=1.77 ps.
Keeping in mind that the interaction with an individual intrachromophore mode is, in general, relatively weak (S_0∼ 0.01, irrespective of the vibrational frequency ω_0), <cit.> one could expect that already the BA recovers the numerically exact exciton dynamics.
Figures <ref>(a) and <ref>(b) reveal that this is indeed the case when the vibrational mode is not resonant with the excitonic energy gap (ω_0≠Δε_X).
The true dynamics then exhibits weakly damped oscillations in both exciton populations and interexciton coherence, which reflect the relatively slow and inefficient interchromophore population transfer discussed in the literature. <cit.>
While the BA reproduces the oscillatory behavior fairly well on the time scales we focus on, the SCBA produces an excessively fast oscillation damping and an overall incorrect dynamics of exciton populations.
One may then expect that the general tendency of the self-consistent cycle towards fast equilibration of the excitonic subsystem could render the SCBA a reliable approximation when the vibrational frequency is nearly resonant with the exciton-energy gap (ω_0≈Δε_X).
At resonance, the existing results show that the damping of the interexciton coherence and the concomitant interchromophore population transfer are particularly fast. <cit.> 
Figures <ref>(a) and <ref>(b) show that the SCBA is indeed sufficiently good at reproducing this rapid equilibration of both exciton populations and coherences.

In the off-resonant case, the exciton populations and interexciton coherence predominantly oscillate at frequencies |ω_0-Δε_X| and Δε_X, respectively, see the most pronounced features of the spectra in Figs. <ref>(c) and <ref>(d).
The spectra in Figs. <ref>(c), <ref>(d), <ref>(c), and <ref>(d) originate from our frequency-domain computations, so that they are broadened with parameter η. 
The incorrect SCBA population dynamics in Fig. <ref>(a) is reflected on the frequency domain as the spurious low-frequency feature in Fig. <ref>(c).
In the resonant case, the dynamics of interexciton coherence displays beats most probably stemming from oscillations at two similar frequencies, see the most intensive features in Fig. <ref>(d).
In contrast to the off-resonant case, the relation between the beating frequency or the frequency of population oscillations in Fig. <ref>(c) to the inherent energy scales of the problem (Δε,J,ω_0) is not obvious.   
To establish such a relation, we note that even the weak exciton–vibration interaction can induce appreciable mixing of the vibrational and excitonic levels, and thus render the description in terms of vibronic-exciton states more appropriate.
To discuss our observations in terms of these hybrid states, we assume, for simplicity, that the vibrational mode is undamped, i.e., γ_0=0.
Then, the Hamiltonian is analogous to the Jaynes–Cummings Hamiltonian of quantum optics. <cit.>
When S_0≪ 1, it is sufficient to consider at most a single vibrational quantum, <cit.> and in Appendix <ref> we conclude that the first two states above the lowest-lying vibronic-exciton state |1_x,v_0=0⟩ are linear combinations of the states |1_x,v_0=1⟩ and |2_x,v_0=0⟩.
In the nearly resonant case (Δε_X≈ω_0), the energies corresponding to the oscillatory features in exciton populations

    Ω_pop^res≈ω_0√(2S_0)|sin(2θ)|,

and the interexciton coherence

    Ω_coh,±^res≈Δε_X±ω_0√(S_0/2)|sin(2θ)|

are proportional to √(S_0) and agree reasonably well with the corresponding numerically exact and SCBA results, see Figs. <ref>(c) and <ref>(d).
Interestingly, the frequency of the beats in the interexciton-coherence dynamics virtually coincides with the frequency of population oscillations, i.e., Ω_pop^res≈Ω_coh,+^res-Ω_coh,-^res.  
In the off-resonant case Δε_X≠ω_0, we find that the frequency of population oscillations is shifted from |Δε_X-ω_0| by an amount proportional to S_0, i.e.,

    Ω_pop^off-|Δε_X-ω_0|≈ S_0sin^2(2θ)ω_0^2/|Δε_X-ω_0|,

which agrees very well with the HEOM result, and not so well with BA and SCBA results, see Fig. <ref>(c).
The frequency shift of the most pronounced component of exciton-coherence oscillations from Δε_X is also linear in S_0,

    Ω_coh,-^off-Δε_X≈-1/2S_0sin^2(2θ)ω_0^2/|Δε_X-ω_0|,

in good agreement with HEOM and BA results in Fig. <ref>(d).
We mention that the much less intensive feature of the coherence-oscillation spectrum appearing around the vibration energy is shifted from ω_0 by 

    Ω_coh,+^off-ω_0≈1/2S_0sin^2(2θ)ω_0^2/|Δε_X-ω_0|.




 §.§ Reliability of the SCBA: Seven-site model of the FMO complex






In Figs. <ref> and <ref>, we benchmark the SCBA and BA on the widely studied seven-site model of the FMO complex found in green sulphur bacteria.
Detailed benchmarks of our approximations for chromophore populations (Fig. <ref>) and interchromophore coherences (Fig. <ref>) against numerically exact results are possible for the OBO SD [Eq. (<ref>)].
In Figs. <ref>(a1), <ref>(b1), <ref>(a1), and <ref>(b1), we use λ_ph=35 cm^-1 and γ_ph^-1=50 fs (γ_ph=106.2 cm^-1).
To explore the viability of our methodology, we also apply it to the model using the structured SD emerging from atomistic simulations performed in Ref. JPhysChemLett.7.3171.
As we are not aware of any numerically exact results for the dynamics modulated by this structured bath, Figs. <ref>(a2), <ref>(b2), <ref>(a2), and <ref>(b2) compare SCBA and BA results.
Details on the excitonic Hamiltonian H_S and the structured SD used are summarized in Appendix <ref>.

Using the OBO SD, already the BA provides a good approximation to the subpicosecond dynamics of chromophore populations [Figs. <ref>(a1) and <ref>(b1)], in agreement with findings of Ref. PhysRevE.86.011915.
Nevertheless, on a picosecond timescale, the BA overestimates (underestimates) populations of low-energy (high-energy) states, and this effect becomes more pronounced with decreasing the temperature, compare Fig. <ref>(a1) to Fig. <ref>(b1).
While the SCBA suffers from the same deficiency, its predictions for chromophore populations are systematically closer to the numerically exact results throughout the time window examined.
At higher temperatures, Figs. <ref>(a1) and <ref>(a1) suggest that the SCBA is better at approximating population dynamics than coherence dynamics.
On the other hand, Figs. <ref>(b1) and <ref>(b1) suggest that the reverse is true at lower temperatures.
Figures <ref>(a1) and <ref>(b1) do not show BA results, which, at all temperatures, display oscillatory features lasting much longer than the numerically exact method predicts.  

Comparing panels (a2) and (a1) [(b2) and (b1)] in Figs. <ref> and <ref>, we conclude that our approximations deliver reasonable results for exciton dynamics in the structured environment, as it is overall similar to the dynamics in the featureless environment.
The problems with longer-time population dynamics of extremal-energy states are exacerbated in the structured environment at lower temperatures, when SCBA (and also BA) predicts nonphysical (greater than 1 or negative) populations of such states, see Fig. <ref>(b2).
A possible origin of these problems can be understood from our analysis of the dynamics modulated by an underdamped vibrational mode, see Sec. <ref>.
There, we find that the resonance condition between the exciton-energy gap and the vibrational energy quantum is crucial to the success of the SCBA.
Here, exciton-energy gaps fluctuate due to low-frequency phonon modes, their fluctuations becoming more pronounced with increasing temperature.
Therefore, at higher temperatures, satisfying resonance conditions between exciton-energy gaps and weakly damped vibrational modes is more probable than at lower temperatures.
In other words, the SCBA is expected to work better at higher temperatures.
The SCBA coherence dynamics in Figs. <ref>(a2) and <ref>(b2) is physically sensible, while the lifetime of interchromophore coherences in the structured environment is somewhat shorter than in the featureless environment at both temperatures examined.



§ SUMMARY AND OUTLOOK


We have developed and benchmarked the self-consistent Born approximation for studying the dynamics of EET through a multichromophoric aggregate linearly interacting with a bosonic environment.
We start from the lowest-order approximation for the memory kernel of the GQME, and improve it in the self-consistent cycle based on the GQME represented in the Liouville space.
We find that the SCBA reproduces the exact exciton dynamics modulated by an overdamped phonon continuum very well, even in the generally difficult regimes of strong exciton–environment interaction, slow environmental reorganization, and low temperature.
This success of the SCBA can be understood from the analytically tractable example—coherence-dephasing dynamics in the pure-dephasing model.
We conclude that the reliability of the SCBA when the dynamics is modulated by an underdamped vibrational mode leans on the resonance between the mode frequency and the exciton-energy gap.
Nevertheless, in a structured environment comprising both an overdamped phonon continuum and a number of underdamped vibrational modes, the SCBA reasonably describes exciton dynamics through the seven-site model of the FMO complex.

Importantly, our method does not introduce any assumptions on the form of the exciton–environment SD, making it a strong candidate for studying exciton dynamics modulated by structured environments whose properties are extracted from experiments or atomistic simulations.
Moreover, the method can, in principle, be systematically improved by performing the self-consistent cycle starting from a higher-order approximation for the memory kernel.
The next member of the family of self-consistent approximations thus obtained is the so-called one-crossing approximation, in which the starting memory kernel is the sum of the first two diagrams in Fig. <ref>(f).
The final memory kernel then contains all diagrams in which the lines representing the environmental assistance cross at most once.
However, the one-crossing approximation is computationally much more demanding than the SCBA, as each iteration involves a double integral over frequency.
The practical applicability of the one-crossing approximation is thus determined by the balance between its computational requirements and the improvements it offers over the SCBA.
Finally, our Liouville-space frequency-domain formulation of the SCBA suggests that it might be used as a computationally efficient and reasonably accurate approach to compute experimentally accessible nonlinear response functions.
To test such a possibility, one should generalize the present formulation so that memory kernels in different excited-state sectors can be computed. 



§ AUTHOR CONTRIBUTIONS

V.J. developed the method and performed all analytical derivations.
V.J. and T.M. numerically implemented the method and conceived the examples on which the method was tested.
V.J. performed all numerical computations and prepared the initial version of the manuscript.
Both authors contributed to the final version of the manuscript.

Work in Prague is funded by the Czech Science Foundation (GAČR), Grant No. 22-26376S.
Work in Belgrade is supported by the Institute of Physics Belgrade, through the grant by the Ministry of Science, Technological Development, and Innovation of the Republic of Serbia.
Numerical computations were performed on the PARADOX-IV supercomputing facility in the Scientific Computing Laboratory, National Center of Excellence for the Study of Complex Systems, Institute of Physics Belgrade.



§ COHERENCE DEPHASING IN THE PURE-DEPHASING MODEL (OBO SD AND HIGH-TEMPERATURE LIMIT)

We first show that the CFE of the exact Σ(ω) is given by Eq. (<ref>).
In the high-temperature limit, the real part of the lineshape function can be approximated as <cit.>

    g^r(t)=2λ_phT/γ_ph^2(e^-γ_pht+γ_pht-1).

Inserting Eq. (<ref>) into Eq. (<ref>) and performing the Fourier transformation of the latter gives (n_2≠ n_1)

    n_2n_1|𝒢(ω)|n_2n_1=e^Λ_ph∑_k=0^+∞(-Λ_ph)^k/k!1/ω-ε_n_2n_1+iΛ_phγ_ph+ikγ_ph,

where we introduce the dimensionless parameter Λ_ph=4λ_phT/γ_ph^2.
It is known that the zero-temperature absorption lineshape (the excitation-addition spectral function) of a two-level system whose energy gap ε_eg between the ground and excited state is modulated by an undamped vibrational mode (frequency ω_0 and HR factor S_0) is proportional to the negative imaginary part of the retarded Green's function (see, e.g., Ch. 8 of Ref. Mukamel-book or Ch. 4 of Ref. Mahanbook):

    G^R(ω)=e^-S_0∑_k=0^+∞S_0^k/k!1/ω-ε_eg+S_0ω_0-kω_0+iη.

The CFE of Eq. (<ref>) reads as <cit.>

    G^R(ω)=1ω-ε_eg-(S_0ω_0)^2ω-ε_eg-ω_0-2(S_0ω_0)^2ω-ε_eg-2ω_0-⋯.

The CFE of n_2n_1|𝒢(ω)|n_2n_1 is then obtained by substituting ω-ε_eg→ω-ε_n_2n_1, S_0→-Λ_ph, and ω_0→-iγ_ph in Eq. (<ref>), which follows from comparing Eq. (<ref>) with Eq. (<ref>).
Keeping in mind that n_2n_1|Σ(ω)|n_2n_1=ω-ε_n_2n_1- n_2n_1|𝒢(ω)|n_2n_1^-1, see Eq. (<ref>), one immediately obtains Eq. (<ref>).

We now argue that the CFE of the SCBA self-energy is given by Eq. (<ref>).
Our starting point is the BA propagator in the frequency domain (n_2≠ n_1)

    n_2n_1|𝒢_BA(ω)|n_2n_1=1ω-ε_n_2n_1-4λ_phTω-ε_n_2n_1+iγ_ph,

which is obtained by inserting the BA self-energy [Eq. (<ref>)] into the Dyson equation [Eq. (<ref>)].
The BA propagator is then inserted into Eq. (<ref>) to obtain the following self-energy Σ^(2) in the next iteration:

    n_2n_1|Σ^(2)(ω)|n_2n_1=∫_-∞^+∞dν/2π4λ_phγ_phT/(ω-ν)^2+γ_ph^2 n_2n_1|𝒢_BA(ω)|n_2n_1.

Here, we have approximated coth(β(ω-ν)/2)≈2T/ω-ν.
The last integral is solved by a integrating along a contour that is closed in the upper half-plane to ensure causality.
As the poles of the BA propagator are in the lower half-plane by construction, the only pole of the integrand in the upper half-plane is at ν=ω+iγ_ph.
Evaluating the corresponding residue, we obtain

    n_2n_1|Σ^(2)(ω)|n_2n_1=4λ_phTω-ε_n_2n_1+iγ_ph-4λ_phTω-ε_n_2n_1+2iγ_ph.

Repeating the above-described procedure ad infinitum, we obtain Eq. (<ref>).

We end this section by presenting analytical results for the matrix elements of the Green's superoperator that determine the dynamics of coherence dephasing within the BA and the Redfield theory, see Eqs. (<ref>) and (<ref>).
Equation (<ref>) implies that the BA propagator in the time domain can be expressed as

    n_2n_1|𝒢_BA(t)|n_2n_1=e^-iε_n_2n_1t∫_-∞^+∞dΩ/2πe^-iΩ tΩ+iγ_ph/(Ω-Ω_+)(Ω-Ω_-),

where

    Ω_±=γ_ph/2[±√(16λ_phT/γ_ph^2-1)-i]≈± 2√(λ_phT)-iγ_ph/2

are the roots of the quadratic equation Ω^2+iγ_phΩ-4λ_phT=0, both of which lie in the lower half-plane.
We use 2π T/γ_ph≫ 1 to obtain Ω_± in the high-temperature limit.
The integral in Eq. (<ref>) is solved using the contour integration with the final result

    n_2n_1|𝒢_BA(t)|n_2n_1=-iθ(t)[Ω_++iγ_ph/Ω_+-Ω_-e^-iΩ_+t-Ω_-+iγ_ph/Ω_+-Ω_-e^-iΩ_-t].

Retaining the contributions that are the most dominant in the high-temperature limit 2π T/γ_ph≫ 1, Eq. (<ref>) can be further simplified to

    n_2n_1|𝒢_BA(t)|n_2n_1=-iθ(t)e^-iε_n_2n_1tcos(2√(λ_phT) t)e^-γ_pht/2.

Equation (<ref>) implies that the Redfield propagator in the time domain can be expressed as

    n_2n_1|𝒢_Red(t)|n_2n_1   =e^-iε_n_2n_1t∫_-∞^+∞dΩ/2πe^-iΩ t/Ω+i4λ_phT/γ_ph
       =-iθ(t)e^-iε_n_2n_1texp(-4λ_phT/γ_pht).

The approximation to the exact coherence-dephasing dynamics embodied in Eq. (<ref>) is obtained by inserting the short-time approximation g^r(t)≈λ_phTt^2 to the lineshape function [Eq. (<ref>)] into Eq. (<ref>).



§ VIBRONIC-EXCITON MODEL


We consider the model dimer, Eq. (<ref>), in which the excitons interact with an undamped vibrational mode of frequency ω_0 and HR factor S_0 [we set γ_0=0 in Eq. (<ref>)].
We assume that S_0≪ 1, so that the reorganization energy S_0ω_0 is much smaller than all the other energy scales in the problem (Δε,J).
It is known that the center-of-mass motion of the intrachromophore vibrations, described by B_+=(b_1+b_2)/√(2), does not affect the single-exciton dynamics. <cit.>
The energies of the exciton states, as well as transitions between them, are then modulated by the relative motion of intrachromophore vibrations, which is described by B_-=(b_1-b_2)/√(2).
The single-exciton dynamics is governed by the Hamiltonian H=H_X+H_B_-+H_X-B_-, where

    H_X=∑_k=1^2|k_x⟩[ε_k_x+(-1)^kcos(2θ)ω_0√(S_0/2)(B_-^†+B_-)]⟨ k_x|,


    H_B_-=ω_0 B_-^† B_-,


    H_X-B_-=-sin(2θ)ω_0√(S_0/2)(|2_x⟩⟨ 1_x|+|1_x⟩⟨ 2_x|)(B_-^†+B_-).

The energy of the bare exciton state k=1,2 is ε_k_x=Δε/2[1+(-1)^k√(1+tan^2(2θ))].
The dynamic modulation of the exciton energy can be taken into account exactly by transferring to the polaron frame, H=UHU^†, using the polaron transformation

    U=∑_k=1^2|k_x⟩⟨ k_x|e^(-1)^kS_θ,   S_θ=cos(2θ)√(S_0/2)(B_-^†-B_-).

The bare exciton energies ε_k_x are then shifted to ε_k_x=ε_k_x-S_0ω_0/2cos^2(2θ), while the eigenstates |k_x,v_0⟩ of the transformed Hamiltonian H_X+H_B_- are enumerated by the exciton number k and the number v_0=0,1,… of excited vibrational quanta.
In more detail,

    |k_x,v_0⟩=|k_x⟩(B_-^†)^v_0/√(v_0!)|∅⟩,  ε_k,v_0=ε_k_x-S_0ω_0/2cos^2(2θ)+v_0ω_0.

As we assume that S_0≪ 1, the leading contribution (proportional to √(S_0)) to the transformed interaction term H_X-B_- in powers of S_0 is identical to H_X-B_-.
To proceed further, we additionally perform the rotating-wave approximation, which amounts to

    H_X-B_-≈ -ω_0√(S_0/2)(|2_x⟩⟨ 1_x| B_-+|1_x⟩⟨ 2_x|B_-^†).

We now limit ourselves to the subspace containing at most a single vibrational excitation, where we find that H_X-B_- mixes the states |1_x,v_0=1⟩ and |2_x,v_0=0⟩.
In the resonant case ω_0≈Δε_X, these two states are nearly degenerate, and H_X-B_- lifts this degeneracy.
While there is no degeneracy in the off-resonant case ω_0≠Δε_X, the vibronic mixing still affects the energy difference between these two states.
The energies of vibronically mixed states measured with respect to the energy of the lowest-lying state |1_x,v_0=0⟩ then read as

    E_v_0=1,±-E_v_0=0=1/2[Δε_X+ω_0±√((Δε_X-ω_0)^2+2S_0ω_0^2sin^2(2θ))].

The exciton populations then exhibit oscillatory features of frequency

    Ω_pop=E_v_0=1,+-E_v_0=1,-=√((Δε_X-ω_0)^2+2S_0ω_0^2sin^2(2θ)),

while the oscillations in interexciton coherence have the frequencies

    Ω_coh,±=E_v_0=1,±-E_v_0=0.


In the nearly resonant case (Δε_X-ω_0≈ 0), the frequencies of oscillatory features in exciton-population and interexciton-coherence dynamics are given by Eqs. (<ref>) and <ref>, respectively. 
To be consistent with the simplifications made up to now, our considerations in the off-resonant case have to keep only the lowest-order term in small S_0, so that

    Ω_pop^off≈|Δε_X-ω_0|[1+S_0sin^2(2θ)ω_0^2/(Δε_X-ω_0)^2],


    Ω_coh,±^off=1/2[Δε_X+ω_0±|Δε_X-ω_0|± S_0sin^2(2θ)ω_0^2/|Δε_X-ω_0|].

Equations (<ref>), (<ref>), and (<ref>) are then readily obtained from Eqs. (<ref>) and (<ref>).


§ COMPUTATIONS ON THE SEVEN-SITE MODEL OF THE FMO COMPLEX


The excitonic Hamiltonian H_S is taken from Ref. BiophysJ.91.2778, and the values of interchromophore couplings and average chromophore energies are summarized in Table <ref>.


We take the structured SD of the exciton–environment interaction from the spreadsheet  appearing in the Supporting Information to Ref. JPhysChemLett.7.3171.
More specifically, we use the data from the sheet , which reports the total SD (comprising both the interchromophore and intrachromophore contributions) for individual BChls in one of the FMO subunits.
In our computations, we assume that the SDs for all BChls are identical, and thus use the arithmetic average of the SD data for BChl1,…,BChl7.
The SD thus obtained is plotted in Fig. <ref>.
The frequency grid is equidistant, with spacing Δω=0.53 cm^-1, which we find sufficiently fine for the numerical integration of Eq. (<ref>) using the ordinary trapezoidal rule.
The SD is Ohmic, and we obtain its linear behavior around ω=0 by fitting the points (0,0),(Δω,𝒥(Δω)),…,(5Δω,𝒥(5Δω)) to a straight line.
As the SD is available up to ω_avail=2650.8 cm^-1, our numerical computations use ω_max=2ω_avail, and assume 𝒥(ω)=0 for ω_avail<ω<ω_max.



While we take the excitonic and the exciton–environment interaction parameters from different studies and different (seven-site <cit.> and eight-site <cit.>) FMO models, we emphasize that our main goal is to examine the applicability of our approximate methods to a relatively realistic model of a multichromophoric complex.
The results in Figs. <ref> and <ref> indeed suggest that our approximate methods can provide decent results on exciton dynamics in realistic models.   


