We simulate the long-range inter-complex electronic energy transfer in photosystem II—from the antenna complex, via a core complex, to the reaction center—using a non-Markovian (ZOFE) quantum master equation description that allows the electronic coherence involved in the energy transfer to be explicitly included at all length scales. This allows us to identify all locations where coherence is manifested and to further identify the pathways of the energy transfer in the full network of coupled chromophores using a description based on excitation probability currents. We investigate how the energy transfer depends on the initial excitation—localized, coherent initial excitation versus delocalized, incoherent initial excitation—and find that the overall energy transfer is remarkably robust with respect to such strong variations of the initial condition. To explore the importance of vibrationally enhanced transfer and to address the question of optimization in the system parameters, we systematically vary the strength of the coupling between the electronic and the vibrational degrees of freedom. We find that the natural parameters lie in a (broad) region that enables optimal transfer efficiency and that the overall long-range energy transfer on a ns time scale appears to be very robust with respect to variations in the vibronic coupling of up to an order of magnitude. Nevertheless, vibrationally enhanced transfer appears to be crucial to obtain a high transfer efficiency, with the latter falling sharply for couplings outside the optimal range. Comparison of our full quantum simulations to results obtained with a “classical” rate equation based on a modified-Redfield/generalized-Förster description previously used to simulate energy transfer dynamics in the entire photosystem II complex shows good agreement for the overall time scales of excitation energy transport.

## I. INTRODUCTION

The initial steps in photosynthesis of green plants occur when photosystem II (PSII) absorbs light and the nascent excitation is used to split water molecules into electrons, hydrogen ions, and molecular oxygen. The ensuing electrons and protons drive the subsequent photosynthetic reactions necessary for the production of adenosine tri-phosphate (ATP), the chemical currency of biological reactions.^{1} Greater than 80% of all excitations initiated in PSII will result in productive photochemistry, causing the creation of an irreversibly separated electron-hole pair contributing to water splitting.^{2} The mechanisms of efficient energy transport within PSII have been studied using spectroscopy^{3–14} and numerical simulations,^{15–18} with models based on spectroscopic data and structural information from X-ray crystal structures of the complexes.^{19–22} Recent work has focused on energy transfer within PSII supercomplexes, which represent the smallest photosynthetic unit that contains all of the relevant proteins for these first steps in photosynthesis.^{2,23,24} The largest PSII supercomplex isolated by Caffarri and co-workers^{24} (called C_{2}S_{2}M_{2}) contains 326 pigments bound into a protein assembly containing light harvesting antennae surrounding a pair of PSII reaction centers. The C_{2}S_{2}M_{2} supercomplex is a dimer, where each of the two monomers contains two antenna complexes (LHCII), three minor light harvesting complexes (CP24, CP26, and CP29), two core antenna complexes (CP43, CP47), and one reaction center (RC). The arrangement of the complexes is shown in Figure 1.

Absorption of sunlight by chlorophyll pigments in the antenna complexes of PSII creates electronic excitations that are transferred to other pigments located in the same protein complex and from there to pigments in neighboring complexes. In this way, the energy is transferred from the antenna complexes to the core complexes and from there to the reaction center. Because of these complicated dynamics, involving many degrees of freedom and a relatively large number of pigments, numerical simulations for PSII are very challenging. Consequently, simulations of energy transfer in PSII have been limited to smaller subcomplexes of PSII or to simpler rate-equation descriptions^{15,16,18,25} rather than the full quantum dynamics simulations that have been performed for smaller light-harvesting systems such as the Fenna-Matthews-Olson complex.^{26–31} Recently, in Ref. 16, simulations based on a combined modified-Redfield-generalized-Förster (MRGF) rate equation approach were carried out for the entire PSII supercomplex, providing new insight into the interplay of short-range transfer dynamics inside the individual subcomplexes (proteins) and long-range inter-complex dynamics within a rate-equation kinetic analysis. This MRGF formalism has subsequently been extended to also describe excitation transport across an intact thylakoid membrane containing approximately 40 000 chlorophyll molecules.^{32}

In recent years, spectroscopy experiments have provided evidence that there is coherence involved in excitonic energy transfer in light-harvesting complexes.^{33,34} This observation raises questions about the role of such coherence in the energy transfer process driving photosynthesis. In particular, is coherence important for the efficiency of the transfer? Is it relevant for design principles that aim to maximize this efficiency? These questions have been addressed in many recent papers.^{35–47} In a companion paper^{48} of the present paper, the role of electronic coherence in electronic excitation energy transfer is analyzed in the framework of excitation probability currents. This analysis shows that when the dynamics are described by a quantum master equation (QME), the probability currents between the pigments and thus the energy transfer are driven by the electronic coherence between the pigments. In the QME, the electronic coherence is explicitly taken into account as off-diagonal elements of the electronic density matrix. In contrast, in a “classical” master equation (rate equation), only the time dependence of the excitation probabilities of the pigments—not the coherence—is described explicitly. Nevertheless, coherence is still taken into account *implicitly* in the rate equation, where it is hidden in the “effective” rates. Without this, no transport would occur according to the more complete QME analysis, as we will show in Section II B 4 (see also the companion paper^{48}).

In the present paper, we investigate the energy transfer dynamics in PSII by means of a (non-Markovian) QME description that includes a full quantum description of the electronic coherence over all length scales. Since the QME simulation allows us to explicitly take into account the electronic coherence involved in the energy transfer both between and within individual subcomplexes, e.g., antenna and core complexes, it enables calculation of both long-range and short-range energy transfer pathways in the large network of coupled pigments in terms of excitation probability currents, using the full quantum framework described in Ref. 48. Using a non-Markovian QME allows us to account for the non-negligible memory times, i.e., the non-Markovianity, associated with both the intra-molecular vibrations of the pigments and the vibrations of the protein environment of the pigments to which the electronic degrees of freedom couple. The calculations presented here investigate the long-range inter-complex energy transfer from the antenna complex via a core complex to the reaction center, specifically, in a subsystem containing a LHCII monomer, the core complex CP43, and the reaction center (see Figure 1), which requires explicit incorporation of 33 pigments in the simulation. To be able to efficiently treat non-Markovian dynamics for this number of pigments, we use the non-Markovian ZOFE quantum master equation that has been developed in recent years,^{49} tested against exact results,^{30,31} and successfully applied to biological light-harvesting complexes. In addition to providing a computationally efficient approach that can handle tens of pigments coupled to a non-Markovian protein environment, the ZOFE approach allows calculations for wide parameter ranges of the model, as well as a flexible representation of the electron-vibration coupling. It also allows a detailed analysis of the role of vibrations—in particular, intra-molecular vibrations of the pigments and vibrations of the protein environment—that couple to the electronic degrees of freedom and that can enhance the electronic energy transfer through vibrationally enhanced transport, i.e., by creating resonances between pigments through vibronic energy levels of the pigments.

Since one of the goals of this work is to reveal the extent of explicit coherence between subcomplexes that is omitted from effective rate descriptions of long-range energy transfer given by MRGF and Förster theory, we use the excitonic Hamiltonian, phonon bath couplings, and spectral densities employed in the MRGF study of Ref. 16 and compare the simulation results of the ZOFE quantum master equation with results obtained using the MRGF rate equation. We find that while there are differences at intermediate time scales, the two approaches show good agreement of the total excitation population within the constituent subcomplexes for time scales longer than 200 ps, validating the use of the effective description provided by MRGF to study coarse-grained energy transport over long time and large length scales.^{16,32} The explicit calculation of coherences within the ZOFE calculations and the population-current analysis that this enables provide additional insight into the manifestations of coherence within and between subcomplexes that is not accessible with the MRGF framework.

In addition to the detailed investigation of how electronic coherence and non-Markovian exciton-phonon coupling are manifested in and affect the long-range energy transfer, we also address the role of the initial excitation conditions. There is an ongoing discussion about whether and how energy transfer in light-harvesting systems depends on the initial conditions—in particular, initial excitation in laser spectroscopy experiments can be coherent and localized, whereas in natural initial conditions through absorption of sunlight, delocalized, incoherent states are often assumed.^{38,41,50,51} Motivated by these discussions, in the present paper we simulate the energy transfer dynamics for very different initial states to ascertain how long-range energy transfer in PSII depends on the initial condition, finding that while at short times there is some dependence on initial conditions, this dependence disappears by ∼200 ps and the overall long time (ns) energy transfer rates are remarkably independent of the form of the initial excitation.

The remainder of the paper is structured as follows. In Section II, we describe the model that we use for the simulation of the energy transfer in PSII. The simulated energy transfer dynamics are then shown in Section III, for an initial state in which the initial excitation is localized in the antenna complex. In Section IV we compare to a simulation where the initial excitation is delocalized over all three complexes considered here. In Section IV, the dependency of the energy transfer on the coupling between the electronic and vibrational degrees of freedom is investigated. A summary of the results and concluding remarks are given in Section V.

## II. MODEL

Our goal is to investigate the long-range transport of excitation energy in the photosystem II supercomplex—from the surrounding antenna complexes, through the outer components of the supercomplex, to the reaction center—by means of a full quantum simulation. To this end, we model a subsection of the photosystem II supercomplex originally isolated by Croce *et al.*^{24} and subsequently used as a foundation for a structure based excitation energy transport model to explain the measured fluorescence lifetimes.^{2} Figure 1 shows the largest PSII supercomplex previously modeled.^{2,16} The colored pigments in LHCII, CP43, and the reaction center shown in Figure 1 represent the subsystem that we model in this work: it contains 33 chromophores distributed between a LHCII monomer (14 pigments), CP43 (13 pigments), and a truncated reaction center (6 pigments).

### A. System parameters

#### 1. Electronic degrees of freedom

Our model contains 33 chromophores in total, each of which we model as having two electronic states separated by the local excitation energy, referred to here as the “site energy,” which varies across different pigments depending upon their interactions with the local protein environment. We further describe the interaction between the chromophores by electronic matrix elements that couple electronic excitation of a chromophore to electronic excitation of other chromophores. In the subspace of single excitations that is relevant for energy transport, this is described by the electronic Hamiltonian,

containing *N* site energies (ε_{n}) and the *N* × *N* coupling matrix *V*. Here, |*n*〉 is the state in which only pigment *n* is excited and all others are in the ground state. All energies and coupling parameters are taken from Ref. 16, as will be described in detail below. The resulting site energies and coupling parameters are shown in Figure 2, together with the structure of the LHCII-CP43-RC truncated supercomplex.

#### 2. Coupling to vibrations

The electronic states of each chromophore are coupled to a large collection of vibrational modes that represent intra-molecular vibrational modes as well as vibrational modes of the surrounding protein scaffold. In the extreme limit of low-frequency vibrations, the long-time-scale protein conformation dynamics give rise to an inhomogeneous distribution of energies for each chromophore within any ensemble measurement. In contrast to Ref. 16, however, where ensemble averaging over static disorder in the site energies is performed, in the present work we do not describe ensemble averaging and instead use the electronic Hamiltonian with average site energies for our simulations. While ensemble averaging can be critical for comparing to experimental measurements, in this paper we investigate quantities that are primarily focused on the fundamental process of excitation energy transport within an individual PSII supercomplex. Here we may take the average energies instead of random energies of a single realization to create a “best-case scenario” for transport efficiency and emergence of coherence.

We describe here the intra-molecular vibrational modes of the pigments together with the vibrational modes of the protein environment by the single vibrational Hamiltonian

where we approximate the vibrational modes to be harmonic. Here, *a*_{nλ} is the annihilation operator of mode *λ* with frequency *ω*_{nλ} that belongs to a pigment *n*. (Here and throughout the paper we set *ħ* = 1.) In this description, each pigment has its own separate bath of vibrational modes that does not directly couple to the modes of another pigment. Indirectly, however, such a coupling can come about through the electronic dipole-dipole interaction between two pigments.

We assume that electronic excitation of a pigment couples linearly to its own bath, such that the overall coupling between the electronic degrees of freedom and the vibrations is given by

where the system operators *L _{n}* are the projectors

*L*= − |

_{n}*n*〉〈

*n*|. The coupling constants

*κ*

_{nλ}that describe the strength of the coupling of electronic excitation of pigment

*n*to vibrational mode

*λ*of this pigment can be expressed through the vibrational spectral density of pigment

*n*,

Assuming that the spectral density is a continuous function, this leads to the vibrational correlation function

which depends on the temperature *T*. The Fourier transform

which we shall refer to as the “vibrational correlation spectrum,” will be useful later on for application of the ZOFE master equation. It can also be directly calculated from the spectral density *J _{n}*(

*ω*) (or vice versa) via the relation

^{52}

with the anti-symmetrized spectral density

Peaks in *C _{n}*(

*ω*) that are narrow compared to the relevant energies of the electronic degrees of freedom will lead to non-Markovian dynamics. When we apply the ZOFE master equation below, we will fit

*C*(

_{n}*ω*) using a sum of Lorentzians.

#### 3. Charge transfer in the reaction center

Up to this point, we have described a system where electronic excitations of individual pigments are coupled both to other pigments and to vibrational modes. Two other key features need to be incorporated into our model. These are (i) energy capture via charge separation at the reaction center, and (ii) loss of excitation via non-radiative or fluorescence processes. We include these dissipative processes by extending our description of the electronic degrees of freedom to include the phenomenological radical pair states (|RP1〉 and |RP2〉) as described in Ref. 16, and the ground state |0〉, where all pigments are in their ground electronic states.

In this truncated supercomplex model, energy is transferred from the excited states of three of the pigments in the reaction center to the first radical pair state RP1. On a slower time scale, some of this excitation is transported back to the excited states of the three pigments in the reaction center. From RP1, the energy is irreversibly and more slowly transferred to the second radical pair state RP2. The latter is the last step in our model and describes the irreversible trapping. We employ this simple phenomenological two-state rate model for the charge separation, introduced in Ref. 16, because our focus here is on the excitation transfer to the RC and the charge separation step is used primarily to measure the overall efficiency of this transfer. Using the same model as Ref. 16 also allows a consistent comparison of the ZOFE QME simulation to the MRGF approach used in Ref. 16. The parameterization of these transport processes between RC, RP1, and RP2 is given in Section II A 4 below.

The Hamiltonian term describing the radical pair states and the ground state is given by

where ε_{GS} is the ground state energy, and ε_{RP1} and ε_{RP2} are the energies of the RP1 and RP2 state. We shall describe the dynamical model employed to treat excitation transfer into the radical pairs and ground state of the pigment in more detail in Section II B 3.

#### 4. Choice of parameter values

To simulate excitation transfer and trapping in PSII supercomplexes we use five different kinds of parameters:

Pigment site energies,

spectral density of electron-phonon coupling at each site,

excitonic coupling between pigments,

phenomenological rate of excitation trapping at the reaction center,

rates of fluorescence and non-radiative decay.

For the site energies and spectral densities we use parameters extracted by matching generalized Förster/modified Redfield calculations with spectroscopic measurements on isolated LHCII proteins^{10,53,54} and core complexes.^{15} The work by Novoderezhkin on LHCII in Ref. 10 and Renger on CP43 in Ref. 15 represents the only Hamiltonians that have been demonstrated to reproduce transient absorption measurements, and thus excitation energy transport, on the isolated pigment-protein complexes. Given this, we have used a parameterization of the PSII supercomplex that contains the spectral density of Refs. 10, 53, and 54 for LHCII and that of Ref. 15 for CP43 and RC to describe the coupling of chlorophyll excited states in each of these complexes to their local vibrational baths. These spectral densities are plotted in Figs. 3(a) and 3(b), respectively (thin lines), together with fits described in Section II B 2 that are required to accommodate their use in the ZOFE calculations (thick lines).

We note that the spectral densities used for LHCII and for the core complex are quite different, with the LHCII spectral density showing a series of narrow peaks superposed on a broad distribution while the core complex spectral density shows only a broad distribution. This difference derives from the parameterizations used in their construction. The spectral density of chlorophyll in LHCII was generated in Ref. 53 by fitting to fluorescence line narrowing spectra with distinct vibronic peaks assigned to high frequency intramolecular modes and the effect of coupling to these represented by adding a number of underdamped oscillators to a single overdamped Brownian oscillator form of spectral density. In contrast, the spectral density of chlorophyll in the core complex was generated in Ref. 15 with incorporation only of low frequency protein modes and described instead using a smooth functional form containing no sharp features.^{55} The fit to the LHCII spectral density in Fig. 3(a) averages over the narrow peaks, resulting in a larger magnitude spectral density than the corresponding fitted core complex spectral density (Fig. 3(b)), which reflects the incorporation of coupling to intramolecular chlorophyll vibrations in addition to intermolecular and protein vibrations.

The spectral densities in Fig. 3 represent the bath spectral density for each complex in the site representation. However, since the most important effect of vibrations on excitonic energy transport is to induce relaxation between excitonic eigenstates, more insight into the relevant bath frequencies is gained by comparing the site spectral densities to exciton-vibration coupling density functions that reflect the energies of the transitions between exciton states and the relative strengths of the couplings of these transitions to the vibrations. This is explained in detail in Appendix D. The corresponding exciton-vibration coupling density functions for LHCII and the core complex are shown as bars at the top of each panel in Fig. 3, from which it is evident that the most relevant region of the spectral densities is frequencies below 200 cm^{−1}, with only a few excitonic transitions coupling to higher frequencies and none to frequencies above 500 cm^{−1}.

For item (iii), the excitonic couplings between pigments within each pigment-protein complex are taken from the same sources as the site energies and spectral densities. Excitonic coupling between pigments in different proteins is described by dipole-dipole coupling, as calculated previously in Ref. 16.

To describe excitation trapping at the RC, (iv), we use the characteristic transfer times that are found in Ref. 16, which are 0.64 ps from the excited states of the RC pigments to RP1, 160 ps from RP1 back to the excited states of the RC pigments, and 520 ps for the irreversible transfer from RP1 to RP2.

Electronic excitations present on any pigment are assumed to have equal probability of being lost through either fluorescent or non-radiative decay mechanisms, (v), with a characteristic decay time of 1.78 ns (combining the decay times of 2 ns for non-radiative decay and 16 ns for fluorescence used in Ref. 16). This site-basis representation of fluorescence does not correctly describe the change in fluorescence intensity as a function of time that results from super-radiance effects; however, this is a small influence on the overall population transfer since the probability of fluorescence is always much smaller than for non-radiative decay. Details of the model for excitation trapping at the RC and for radiative and non-radiative decay to the ground state are given in Section II B 3.

### B. Excitation energy transport

Based on Section II A, the total Hamiltonian of our model is now

Since this Hamiltonian comprises a large number of degrees of freedom, direct calculation of the excitation energy transport by explicitly solving the full quantum dynamics is not realistic. Therefore, we use instead the framework of the description of open quantum systems, and divide the total Hamiltonian into three components: a system (*H*_{sys}), an environment (*H*_{env}), and the interaction between system and environment (*H*_{sys–env}). The system part *H*_{sys} should contain the quantities relevant to the energy transport, that is, most importantly the probabilities of electronic excitation of the different pigments and the excitation of the radical pair states in the reaction center. On the other hand, we want to keep *H*_{sys} as small as possible to keep the calculation numerically manageable. Therefore, we choose this to consist only of the electronic degrees of freedom of the pigments and the radical pair states. The vibrations are then treated as the environment. That is, we have

It should be noted that this model, in which all vibrational degrees of freedom are in the environment part *H*_{env}, is equivalent to including vibrational modes explicitly in the system part *H*_{sys}, as is shown, e.g., in Ref. 56. At the level of the Hamiltonian, the division into system and environment is simply a matter of choice, based on what degrees of freedom one wants to observe and based on the capability of the method chosen to calculate the dynamics.

The excitation energy transport and the electronic coherence can be extracted from the reduced density matrix of the system part, i.e.,

which is obtained from the total density matrix by tracing out the degrees of freedom associated with the environment, i.e., the vibrations. In the basis of the states |*n*〉 of the local excitations of the pigments, the diagonal of the reduced density matrix gives the populations of the excited electronic states of the pigments, the radical pair states, and the ground state. The off-diagonal elements are the electronic coherences between the pigments, which provide excitation probability currents between the pigments and thereby constitute the energy transfer, as will be described below.

#### 1. ZOFE master equation

For quantum simulations of the energy transport we use the ZOFE master equation, which allows us to take into account the non-Markovian effects in the coupling between the electronic degrees of freedom and the vibrational environment. The ZOFE master equation is given by^{30,49}

where the auxiliary operator

captures the non-Markovian effects, and its time evolution is given by the auxiliary equation

with initial condition

Depending on the form of the environment correlation function *α _{n}*(

*t*−

*s*), a closed evolution equation for the auxiliary operator $ O \u0304 0 ( n ) ( t ) $ can be obtained.

^{30,49,57}For an environment correlation function that is a sum of damped oscillating terms,

such a closed evolution equation can be found.^{30} Writing the operator $ O \u0304 0 ( n ) ( t ) $ as the sum $ O \u0304 0 ( n ) = \u2211 j O \u0304 0 ( n j ) $, the closed evolution equation to obtain the operators $ O \u0304 0 ( n j ) ( t ) $ is given by^{30}

with initial condition $ O \u0304 0 ( n j ) ( t = 0 ) =0$ and $ O 0 ( n ) ( t , t ) = L n $. In the present work, we use these evolution equations with coupling operators *L _{n}* = − |

*n*〉〈

*n*| for each pigment

*n*. These coupling operators are the ones that appear in the electron-vibrational coupling Hamiltonian, Eq. (3). From the form of the coupling Hamiltonian, and from the discussion in Appendix D, one can see that the present model takes into account both vibronic relaxation and vibronic dephasing effects. This is also evident from the analysis of Ref. 56, where peaks in the environment spectral density are shown to be equivalent to damped vibrational modes that are included into the system part of the Hamiltonian and that couple to the electronic degrees of freedom and induce both vibronic relaxation and dephasing. The non-Markovian nature of the ZOFE-master-equation propagation allows to consistently incorporate these effects. In addition, the explicit non-Markovian formulation of the time-dependent quantum master equation includes the effect of system degrees of freedom “kicking” environment degrees of freedom and the environment “kicking back” the system at a later time.

#### 2. Environment spectral densities for ZOFE master equation

The form of *α _{n}*(

*t*−

*s*) in Equation (17) corresponds to a sum of Lorentzians centered at frequencies Ω

_{nj}and with weights Γ

_{nj}and widths

*γ*in the environment correlation spectrum

_{nj}*C*(

_{n}*ω*) (the Fourier transform of

*α*(

_{n}*t*−

*s*)),

We use this expression to fit the environment correlation spectra that are derived from experimentally determined spectral densities (see below), in order to obtain the parameters Ω_{nj}, Γ_{nj}, and *γ _{nj}* needed to solve the evolution equation for the auxiliary operator Eq. (18). These correlation function fits may then be transformed back to effective spectral densities

*J*(

_{n}*ω*) (Eq. (4)). The overall procedure is as follows:

Calculate the environment correlation spectra

*C*(_{n}*ω*) from the experimentally determined spectral densities*J*(_{n}*ω*) from Refs. 10 and 15, using Equation (7). The temperature entering in*C*(_{n}*ω*) is chosen as*T*= 300 K, for simulation of energy transport at room temperature.Fit these

*C*(_{n}*ω*) with Lorentzians according to Equation (19), to determine the parameters for the evolution equation (18) of the ZOFE auxiliary operators. In this work, we use two Lorentzians for the LHCII pigments and one Lorentzian for the CP43 and RC pigments, to keep the corresponding number of evolution equations small (see index*j*in Eq. (18)).

As a consistency check, we may convert these fitted sums of Lorentzians *C _{n}*(

*ω*) back to effective spectral densities

*J*(

_{n}*ω*) to compare with the original spectral densities. These are shown in Figure 3, together with the original, measurement-based spectral densities of Refs. 10 and 15. The effective spectral densities were constructed with the following parameter sets for Equation (19), where the same parameters were employed for all pigments in CP43 and RC, respectively, whereas in LHCII one set of parameters was used for all Chl A pigments and a second set for all Chl B pigments, respectively. For LHCII (Fig. 3(a)), both the Chl B spectral density (blue curve) and the Chl A spectral density (green curve), the locations of the two peaks are Ω

_{1}= 300 cm

^{−1}and Ω

_{2}= 1200 cm

^{−1}, and the widths of both peaks are given by

*γ*

_{1,2}= 150 cm

^{−1}. For Chl B, the weights of the peaks are Γ

_{1}= 36 000 (cm

^{−1})

^{2}for the low-energy peak and Γ

_{2}= 360 000 (cm

^{−1})

^{2}for the high-energy peak. For Chl A, the peaks are smaller, with Γ

_{1}= 27 000 (cm

^{−1})

^{2}for the low-energy peak and Γ

_{2}= 288 000 (cm

^{−1})

^{2}for the high-energy peak. For CP43 and the RC, Figure 3(b), we use a single peak in the fitted

*C*(

_{n}*ω*). The peak for CP43 (thick red curve) is located at Ω = 90 cm

^{−1}, with width

*γ*= 170 cm

^{−1}and weight Γ = 12 150 (cm

^{−1})

^{2}. The peak for RC (thick orange curve) has the parameters Ω = 100 cm

^{−1},

*γ*= 150 cm

^{−1}, and Γ = 14 000 (cm

^{−1})

^{2}.

As one can see in Figure 3(a), the effective spectral density used for the LHCII pigments in the ZOFE simulation smooths out the narrow peaks of the original spectral density of Ref. 10. This is necessary since the ZOFE approximation is quite limited in the degree of non-Markovianity that it can handle, i.e., the memory time of the vibrations cannot be too long compared to the time scale of the electronic dynamics. Narrow peaks would result in a long memory time, so it is necessary to approximate them by broad peaks that smooth out the sharp features in the coupling to the vibrations.

#### 3. ZOFE master equation: Excitation loss and radical pair states

To simulate the transfer of energy to the radical pair states in the reaction center, and to take radiative and non-radiative decay of the excited electronic states of the pigments into account, we add Lindblad-master-equation terms to the ZOFE master equation. An analogous treatment of radiative decay and trapping in the reaction center is employed in Ref. 29, where a non-Markovian master equation is also extended by the corresponding Lindblad terms, for the simulation of energy transport in the Fenna-Matthews-Olson complex. (See also Refs. 38 and 58 for related Lindblad models for trapping and decay.)

This Markovian Lindblad description of the trapping in the reaction center is a rough approximation of the full charge separation dynamics, which we use here in order to create a minimal model to capture the primary dynamics, and to be able to consistently compare the energy transport simulated with the ZOFE QME with the MRGF approach used in Ref. 16. We assume that the characteristic transfer times to the radical pair states that are found in Ref. 16 give a reasonable approximation and therefore describe the trapping part of our simulation with Lindblad terms based on these characteristic transfer times. This Markovian approximation for the interaction between charge transfer states and vibrational environment is sufficient for the present study, since our main focus is on the excitation transfer to the RC and the charge separation step is included here primarily to measure the efficiency of this excitation transfer.

In the Markov limit, where the environment correlation function *α*(*τ*) decays fast enough compared to all relevant time scales of the dynamics, the ZOFE master equation itself becomes a Lindblad master equation.^{49} It is therefore consistent to add Lindblad terms to the ZOFE master equation, since they are equivalent to additional ZOFE master equation terms with fast decaying correlation functions. In the Markov limit $ \alpha i ( \tau ) = \gamma \u02dc i \u2009\delta ( \tau ) $, the auxiliary operator of the ZOFE master equation becomes the constant operator $ O \u0304 0 ( i ) ( t ) = 1 2 \gamma \u02dc i L \u02dc i $ (see Eqs. (14) and (16)).^{49} (Here, we have replaced the index *n*, which is used above to refer to the pigments, by a more general index *i*, to refer to the respective process described by a system operator $ L \u02dc i $ and a corresponding coupling constant $ \gamma \u02dc i $; e.g., relaxation from a state |*f*〉 to a state |*t*〉 is described by $ L \u02dc f t =|t\u3009\u3008f|$.) Inserting this constant auxiliary operator into the ZOFE master equation (13) gives the well-known Lindblad master equation

with coupling (Lindblad) operators $ L \u02dc i $ and corresponding coupling constants $ \gamma \u02dc i $. We use such a Lindblad description for the following processes:

Radiative and non-radiative decay of the electronic excitation of all pigments (

*n*= 0…32) is described through Lindblad operators $ L \u02dc n decay =|0\u3009\u3008n|$ and the corresponding coupling constant $ \gamma \u02dc n decay $ is given by the decay rate. Here |0〉 denotes the electronic ground state. The decay rate is chosen such that it corresponds to a characteristic decay time of 1.78 ns (combining the two decay times of 2 ns for non-radiative decay and 16 ns for fluorescent decay assumed in Ref. 16).Transfer of excitation from the reaction center pigments 27 (P

_{D1}), 28 (P_{D2}), and 32 (Chl_{D1}) to radical pair state RP1 through $ L \u02dc n RP 1 =|RP1\u3009\u3008n|$ with rates $ \gamma \u02dc n RP 1 $ (where*n*= 27, 28, 32). From all three pigments the transfer to RP1 is assumed to occur with the same rate, which is set at the transfer time of 0.64 ps found in Ref. 16.Transfer from RP1 back to the reaction center pigments 27, 28, 32 through $ L \u02dc n RC =|n\u3009\u3008RP1|$ with rates $ \gamma \u02dc n RC $. It is assumed that equal amounts of population are transferred back to the three pigments, so that each of the three rates is one third of the overall rate corresponding to the characteristic transfer time of 160 ps found in Ref. 16.

Irreversible transfer from RP1 to radical pair state RP2 through $ L \u02dc RP 2 =|RP2\u3009\u3008RP1|$ with a rate $ \gamma \u02dc RP 2 $. This rate is chosen according to the characteristic transfer time of 520 ps found in Ref. 16.

For each of these processes, a corresponding Lindblad term is added to the ZOFE master equation, so that the complete master equation we solve in our simulations is

where the Lindblad operators $ L \u02dc i $ and the corresponding coupling constants (rates) $ \gamma \u02dc i $ in the Lindblad term in the second line of the equation belong to the electronic relaxation processes listed above, and the index *i* of the sum refers to the individual processes. In the ZOFE term in the first line, the ZOFE coupling operators *L _{n}* = − |

*n*〉〈

*n*| couple only electronic excitation of the pigments to the non-Markovian vibrational environment. The radical pair states |RP1〉 and |RP2〉 are not coupled to the non-Markovian vibrations.

The system density matrix *ρ*(*t*) now has components for the states |*n*〉, in which pigment *n* is excited, as well as components for the radical pair states |RP1〉 and |RP2〉 and for the ground state |0〉. *ρ*(*t*) is obtained by solving this master equation together with the evolution equation for the ZOFE auxiliary operators. For the truncated supercomplex considered here, *ρ*(*t*) is then a 36 dimensional square matrix. The system Hamiltonian does not contain coupling to or between RP1, RP2, and the ground state; all couplings to and between these states arise through the Lindblad electronic relaxation terms.

#### 4. Population currents and the contribution of coherence

The transfer of excitation energy between the pigments can be expressed in terms of currents of excitation probability. As shown in a companion paper,^{48} this population current description straightforwardly identifies the contribution of coherence to the energy transfer, and this description also clearly separates the respective contributions of unitary dynamics, dephasing, and relaxation to the population currents between the pigments.^{48} Another advantage of the population-current description is that it clearly shows the individual pathways that the energy transfer takes between the pigments. By integrating the currents over time, it is then possible to see the respective amounts of population that have been transported via each pathway. This insight is useful to identify and analyze the functionality of the light-harvesting supercomplex and its constituent complexes and pigments.

Because the total electronic excitation probability is conserved across the pigments—that is, the reduced electronic density matrix fulfills $ \u2211 n \rho n n ( t ) =1$ at all times—a continuity equation^{48}

holds, where *j _{mn}*(

*t*) is the (net) population current at time

*t*that transports population (i.e., excitation probability) from a pigment

*m*to another pigment

*n*. When

*j*(

_{mn}*t*) is positive, excitation is transported from pigment

*m*to pigment

*n*; when it is negative, the energy transfer goes in the other direction, i.e., from pigment

*n*to pigment

*m*.

Based on the continuity equation (22), for quantum master equations of the form used in the present work, Eq. (21), it is shown in Ref. 48 that the individual population currents *j _{mn}*(

*t*) between the pigments can be calculated through

where $ j m n unitary $, $ j m n dephas $, and $ j m n relax $ are the respective contributions of unitary dynamics, dephasing, and electronic relaxation, which are specified in the second line of the equation. Here, *ρ _{mn}*(

*t*) are the coherences between the sites,

*ρ*(

_{nn}*t*) are the populations of the sites, and

*V*are the inter-site couplings, and

_{mn}*γ*is the rate of electronic relaxation from a site (or state)

_{mn}*m*to a site

*n*.

When we have the density matrix *ρ*(*t*) at a time *t* from a numerical simulation as described above, we can use Equation (23) to calculate the currents between the sites. As described in detail in Ref. 48, the first term in Equation (23) stems from the unitary part of the dynamics—that is, in our model of PSII from the unitary transfer driven by the electronic coupling between the pigments. The second term stems from dephasing through the coupling to the vibrations and gives zero, that is, the dephasing does not contribute to the population currents. The third term stems from electronic relaxation processes, that is, it describes the transfer of population to the radical pair states enabling the charge separation in the reaction center, as well as the radiative and non-radiative decay to the ground state.

We see in Equation (23) that the unitary part of the currents comes entirely from the coherence between the sites—more precisely the imaginary part of the coherences. The relaxation part of the currents that describe the transfer to the RP states, on the other hand, comes entirely from the populations of the sites, as can be seen in the third term in Equation (23). This means that the energy transfer from LHCII to CP43 and to the RC derives entirely from the coherence between the sites, whereas the transfer to the RP states is entirely through the population. Thus, our calculations of the currents between the pigments that are presented below also show how much coherence is present between the pigments, or more precisely, how large the imaginary parts of the coherences between the pigments are. For instance, currents between distant pigments, i.e., currents that transport the excitation over a long distance through space, are based on the corresponding long-distance coherences between the pigments.

When the coupling to the vibrations is strong compared to the electronic inter-pigment coupling, the resulting strong dephasing destroys a large part of the coherence between the pigments. We see from Equation (23) that this will decrease the energy transfer between the pigments, because the unitary contribution to the population current, which in our model amounts to all the energy transfer between the pigments, decreases since it is entirely described by the coherence. This suppression of the transport at strong coupling to the vibrations is well known and often explained in the picture of the quantum Zeno effect:^{59} the electronic excitation interacts with the vibrations—is “measured” by the vibrations—so strongly (i.e., at such a high rate) that the excitation is forced to stay in the original state and cannot move anymore.

To obtain the overall net currents between the different proteins, we first calculate the currents *j _{mn}*(

*t*) between individual sites

*m*and

*n*of the different proteins and then sum these currents up. That is, the current

*J*(

_{AB}*t*) between subcomplex

*A*and a subcomplex

*B*is given by

where when the current *J _{AB}*(

*t*) is positive, there is a net flow from

*A*to

*B*, and when

*J*(

_{AB}*t*) is negative, there is a net flow from

*B*to

*A*.

## III. ENERGY TRANSFER: SIMULATION RESULTS

In order to allow for a full quantum simulation of energy transfer, we construct a minimal model of light harvesting in PSII. In their natural environment (the thylakoid membrane), PSII supercomplexes are surrounded by light harvesting proteins (LHCII). Excitation energy transfer from the surrounding LHCII pigments into the supercomplex occurs via the light harvesting protein mechanically bound to the core complex (LHCII-S and LHCII-M). Our minimal model of the process of excitation energy transfer from the surrounding light harvesting proteins to the reaction center contains one LHCII monomer, the complete CP43 protein, and a slightly truncated RC. We use the ZOFE master equation to calculate the time-dependent density matrix in the space of the electronic states of the pigments, as described in Sec. II B. From the time-dependent density matrix, we obtain the time-dependent excitation probabilities of the pigments and, based on the coherence between the pigments, then calculate the excitation probability currents between the pigments and between the subcomplexes to analyze the energy transfer and its pathways.

It is important to differentiate between the characteristics of excitation energy transfer and the functional behavior of the light harvesting apparatus. In low light conditions, plant growth depends on efficient light harvesting—creating as many productive photochemical events in the RC as possible per photon absorbed.^{60} As such, in this paper, we shall equate the fraction of excitations that cause productive photochemistry (i.e., reach the RP2 state) with the function of the light harvesting antenna. A decrease in the fraction of excitations reaching RP2, then, is a decline in the functional behavior of the photosynthetic assembly. It follows that changes in the precise nature of the excitation dynamics do not necessarily result in changes to the functional behavior. Two different excitation dynamics can still give rise to equivalent quantities of excitation causing irreversible electron loss (reaching the RP2 state). As a result, to show that certain terms in an equation of motion are important to photosynthesis requires both that they change the excitation dynamics and that they increase (or decrease) the fraction of excitation that reaches RP2, the latter being a measure of how much they influence the function of photosynthesis. In particular, we use the amount of population that reaches the second radical pair state RP2 in the reaction center within 1 ns as a measure for overall efficiency of the transport. This measure of the amount of population at a fixed time (1 ns) also puts emphasis on the speed of transport. (We note that one could also use the equilibrated population on RP2 as a measure for transport efficiency, which would put more focus on the amount of population, rather than speed of transport.)

In our first study, we describe excitation energy transfer when initial excitation occurs in a pigment-protein complex outside the supercomplex and subsequent energy transfer occurs into the core complex via the attached LHCII trimer. We therefore initiate excitation on two Chl A molecules in the LHCII monomer (sites 7 and 10) to represent excitation that is entering the system from the neighboring LHCII monomers, since these sites are assumed to have a large rate of transfer with the neighboring LHCII monomers.^{16} Each of the sites 7 and 10 carry half of the initial population, and we choose the initial state to not contain any coherence between the sites, i.e., $\rho ( 0 ) = 1 2 |7\u3009\u30087|+ 1 2 |10\u3009\u300810|$. Thus, all off-diagonal elements of the initial density matrix are zero in the site basis. The dynamics following delocalized initial excitation will be discussed in Section IV A.

### A. Population

To investigate the (long-range) excitation energy transfer between the complexes and the transfer of energy to the charge separated states (RP1 and RP2), we consider the population of the electronic excitation in the individual complexes, $ P A = \u2211 m \u2208 A \rho m m ( t ) $, and the population of the RP states over time. Figure 4(a) shows these population dynamics. The solid lines in Figure 4(a) are the populations on each complex, *P _{A}*, for A = LHCIIC, CP43, RC, calculated with the ZOFE master equation. The population of site 7 in LHCII is shown separately (thin green line). Nearly all the excitation (population) is transferred from LHCII through CP43 and RC to the first radical pair state RP1 within about 200 ps. Irreversible transfer from RP1 to RP2 takes much longer because of the much weaker coupling between RP1 and RP2. After 1 ns, 80% of the population has been transferred to RP2. Within this time, only a small fraction of less than 5% is lost through radiative and non-radiative decay, in keeping with the relatively slow dynamics of excitation loss in a photosynthetic apparatus. As we shall see in the following discussion, this separation of time scales between excitation energy reaching the RC and the relatively slow loss mechanisms of fluorescence and non-radiative decay is a key feature of efficient light harvesting in PSII.

The dashed-dotted lines in Figure 4(a) show the population calculated with the modified-Redfield/generalized-Förster (MRGF) method of Ref. 16. In the MRGF calculations, pigments are grouped together into “domains.” The domain assignments used in this paper are identical to those previously described in Ref. 16. In brief, two pigments were included in the same domain if they contribute strongly to the same eigenstates of the electronic system Hamiltonian. Thus, *S*_{μ,γ} in Eq. (25) is a measure of how strongly pairs of pigments (*μ*, *γ*) contribute to the same exciton (*m*), normalized by each exciton’s inverse participation ratio $ ( \u2211 \xi U \u0303 \xi , m 4 ) $,

In keeping with previous work^{10} suggesting that exciton coupling below 15 cm^{−1} is insufficient to support delocalized excitons in PSII, all couplings below this value were set to zero prior to calculating $ U \u0303 \mu , m $, the matrix of electronic eigenvectors (*m*) in the site basis (*μ*). Two pigments were considered to be in the same domain when *S*_{μ,γ} > 0.1 for at least 50% of all realizations of disorder. Note that in the current paper disorder has not been considered, however we maintain the same domain definitions for the purpose of comparison. Within these domains, the dynamics are treated with modified Redfield theory. Between the domains, the transport is calculated based on transfer rates between the exciton states of one domain (electronic eigenstates of the domain) and the exciton states of another domain. These rates are determined using generalized Förster theory, i.e., employing the overlaps between emission and absorption spectra of the different domains.^{16} We note that in the MRGF calculation, the Chl 601b pigment in the LHCII-c monomer, present in the ZOFE calculation (pigment 0 in Figure 2), is excluded and instead its symmetric image in LHCII-b (Chl 601b’), which is more strongly coupled to LHCII-c,^{10} is included. This exchange of the outside pigments has little effect on the overall transport dynamics in Figure 4(a).

As we can see in Figure 4(a), at ns time scales the subcomplex and RP population dynamics from the MRGF calculation—especially the RP2 population—agree well with that of the ZOFE simulation, even though the former is a much simpler effective classical rate equation calculation. This is consistent with the findings of Ref. 38 that classical rate equations can be adequate for the calculation of *overall* transport efficiencies *provided that they are properly derived from the full quantum description*. At shorter time scales the MRGF subcomplex populations do nevertheless show some differences from the ZOFE populations, in particular MRGF overestimates transport rates between CP43 and the RC resulting in a faster growth of the RP1 population compared to ZOFE simulations. We can see the faster growth of RP1 by comparing the first time in the simulation that RP1 population reaches 40% of all excitation: 34 ps for MRGF compared to 45 ps for ZOFE. Of course, MRGF also provides no information on individual chromophore populations, thus giving no information on population currents (see Section III B below).

In Figure 4(a), we also show the result from a pure Förster calculation (dotted lines) for comparison, since this simple, perturbative method is often used to calculate energy transfer in biological systems. The pure Förster calculation is based on rates obtained from the overlaps between emission and absorption spectra of each pair of pigments (i.e., donor and acceptor pigments).^{61,62} This simple approximate method can give reasonable results in a regime where the coupling between donor and acceptor pigments is weak compared to the coupling to the vibrations. However, since this condition is not met for many pigment pairs in PSII, especially inside the more strongly coupled domains, we expect the pure Förster calculation to be less accurate than the ZOFE calculation as well as the MRGF calculation of Ref. 16. In Figure 4(a) we do indeed see a significant deviation of the pure Förster calculation (dotted lines) from the other simulations. Furthermore, the transport efficiency (population reaching RP2 within 1 ns) calculated with pure Förster is lower.

The observation that the MRGF and ZOFE calculations are in close agreement with each other, while the pure Förster calculation gives a lower transport efficiency, supports the key role of exciton delocalization in excitation transport through PSII. This dynamic effect, well known from the early theoretical literature,^{63–67} has recently been rediscovered and referred to as “supertransfer,” since it occurs when strongly coupled pigments allow an excitation to coherently delocalize within the donor domain, thereby enhancing the rate of transport by up to a factor equal to the number of pigments in the donor domain (see, e.g., Ref. 41 and references therein). This dynamical enhancing effect of excitonic delocalization within a subcomplex or domain is naturally included in the MRGF and ZOFE calculations, but it is not taken into account in the pure Förster description, which thus leads to a lower transport efficiency.

### B. Population currents

Figure 4(b) shows the population currents *j*(*t*) between proteins and from the reaction center to the radical pair states. These currents are calculated from the time-dependent density matrix using Equations (23) and (24).

From Equation (23), we know that the energy transfer from LHCII to CP43 and on to the RC is entirely through the coherence between the pigments, and the transfer to the RP states is entirely through the population. This is combined in Figure 4(b), where the solid lines are the total currents, and the dashed and dotted lines show the coherence and population contribution, respectively.

We see that all currents between complexes are positive over the entire time range. This means that the net flow of energy is directed towards the RP2 state at all times and no temporary intercomplex backflow occurs. Because of the high rate of transport between the RC and the RP1 state, in Figure 4(b) the current from RC to RP1 almost completely follows the current from CP43 to the RC; there is only a very small delay of the current to RP1. This means the population is almost passed right through the RC to RP1. The current to the RP2 state has a smaller magnitude but stretches over a longer time range. Therefore, the overall population that is transferred to RP2 during this time span, which is simply the current integrated over time, is not much smaller than the population transferred to the RC, namely 80% versus ≳95%. As expected, the *direct* current between LHCII and the RC (cyan curve in Figure 4(b)) is zero, due to the large distance between the LHCII and RC pigments.

To see the individual pathways of the energy transfer, we integrate the population currents between the individual pigments over time. The result is shown in Figure 5(a). We see in Figure 5(a) that there is a complex network of pathways rather than just one or two dominant pathways. Previous simulations using the MRGF model have shown that transport through PSII does not occur by a single, obligate pathway.^{16} Our current simulations demonstrate that there are multiple transport pathways from LHCII to CP43, and from CP43 to the RC. However, particularly within complexes, there are a few pathways that transport a particularly large amount of excitation energy, identified in Figure 5(a).

Remarkably, in Figure 5(a) we see that there are a number of relatively strong intra-complex currents that run in a circular fashion, where energy is transported in a loop between pigments. For instance, there is a loop in LHCII from pigment 7 to 13 to 2 to 3, and back to 7. In CP43, there are several such loops. For instance, from 25 to 26 to 23, and back to 25. And there is a longer one from 25 to 26 to 18 to 22 to 23, and back to 25. There is also one from 16 to 21 to 24, and back to 16.

In the reaction center, there is a strong circular current from 32 (Chl_{D1}) to 27 (P_{D1}) to 28 (P_{D2}), and back to 32. These ring currents show significant correlation with the domains of delocalized excitation used in the MRGF calculation. In particular, we find that the strong ring currents in CP43 and the RC are each confined to within one MRGF domain. Only in LHCII, where the electronic energy gaps between the pigments are larger (see bars in Fig. 2(b)) and the domains are smaller, do the ring currents run across different domains. In this situation the ring currents are weaker, consistent with the weaker excitonic coupling between chromophores in different domains. Thus the ring currents reflect excitonic delocalization.

We also see in Figure 5(a) that in the reaction center overall the dominant net flow to the RP1 state is via pigment 32, which in turn obtains most of the population from pigment 30 (Pheo_{D1}) which is the most important connector for inter-complex transport between CP43 and RC. Even though pigments 27 (P_{D1}) and 28 (P_{D2}) are part of the strong circular current, they do not appear to contribute significantly to the net transport to RP1. There is even a net backflow from the RP1 state to the pigments 27 and 28, which acts to decrease the efficiency of the transport. Further, pigments 27 and 28 only receive a small amount of population from CP43 pigments. In the light of these pathways, it appears that the efficient excitation energy transfer to and in the RC relies mainly on the two pigments 30 (Pheo_{D1}) and 32 (Chl_{D1}), and that therefore other reaction center pigments may rather be needed for their role in the electron transport (described by the RP states in our model) than for the excitation transfer. To test this hypothesis, we ran additional simulations where—except for the pigments 30 and 32—we decoupled all RC pigments from all other pigments, inhibiting electronic excitation transfer to all RC pigments but 30 and 32, while leaving the coupling to the RP states unmodified. We found the resulting efficiency of the energy transport to the RP2 state in this restricted model to be indeed the same as in the original model where the excitation transport to the other RC pigments was allowed. This shows that only the two pigments 30 (Pheo_{D1}) and 32 (Chl_{D1}) may be needed for the excitation energy transport to the RP states, and the role of other RC pigments appears to be mainly in the electron transport (described by the RP states).

## IV. ROBUSTNESS OF ENERGY TRANSPORT AND FUNCTION IN PSII

In Sec. III we have described excitation energy transfer when initial excitation occurs in a pigment-protein complex on the periphery of the supercomplex and subsequent energy transfer occurs into the core complex via the attached LHCII trimer. In this section, we shall explore how robust the excitation energy dynamics and the functional behavior of PSII are to changes in the initial state and the electron-phonon coupling.

### A. Delocalized initial excitation

There is an ongoing discussion about the initial conditions of light absorption in light-harvesting systems, and how they affect the subsequent transport dynamics.^{38,41,50} In particular, initial conditions in nature—excitation by sunlight—versus initial conditions in laser spectroscopy experiments—excitation through coherent laser light—have been discussed by a number of authors.^{38,41,50}

Motivated by this, we now investigate the influence of the initial state on the long-range energy transfer in PSII. So far, we have considered initially localized excitation of only two pigments in the LHCII, with initial coherence in the exciton basis but not in the site basis. Now, we contrast this by considering an initial state that is a completely delocalized superposition of all exciton states of all three proteins. We choose this initial state to be an (equally weighted) *incoherent* superposition of all exciton states |Φ_{k}〉 of the electronic system Hamiltonian, i.e., no initial coherence in the exciton basis. Thus, we take an initial density matrix, where in the exciton basis all diagonal elements (populations) are 1/*N*_{sites} (*N*_{sites} is the number of sites), except for the population of the radical pair states and the ground state, which are zero, and with all off-diagonal elements (coherences) set to zero, i.e., $\rho ( 0 ) = \u2211 k = 0 32 1 33 | \Phi k \u3009\u3008 \Phi k |$.

Since the exciton states are the eigenstates of the system Hamiltonian, this state would persist in the course of purely unitary system time propagation. But the coupling to the vibrations, the coupling to the radical pair states, and the radiative and non-radiative decay to the ground state cause the initial state to evolve in a non-unitary manner, driving the excitation transport through PSII. These dynamics lead to the creation of excitonic coherence in the course of time evolution, even though it is not present initially. The transport dynamics that result from this incoherent, delocalized initial state are summarized in Figure 6.

The population dynamics are shown in Figure 6(a) and are compared to the dynamics resulting from the localized initial state (Figure 6(a)). We see that for the delocalized initial state, the population dynamics of the pigments now look quite different. In particular, 40% of the initial excitation is in the CP43 protein and almost 20% is in the reaction center pigments. As a result, 60% excitation is spatially closer to the reaction center than in the previous simulation. This results in a much earlier rise in the RP1 population (black line) than seen in our previous simulation (black, dashed line). Some of this excitation, however, is initiated in higher energy states allowing for back-transfer from CP43 to LHCII as can be seen in Figure 6(b). The green curve, representing the flow of excitation from LHCII to CP43, is negative now for an initial period of time, showing that there is a flow of excitation from CP43 back to LHCII. That is to say, after initial excitation there is a net excitation transport in the opposite direction from what was seen before for the localized initial state, namely away from the reaction center.

Despite the changes in excitation transport dynamics, the functional behavior of PSII is invariant to the change in initial excitation. In particular, the overall transport of population to the radical pair state RP2 is almost identical for both initial excitation conditions, as can be seen in Figure 6(a). (We note that because of the logarithmic time axis, the earlier deviation between the RP1 populations for the different initial states appears to be much larger than the deviation between the RP2 populations at later times, but in fact these deviations are of comparable magnitude.) This shows that changing the initial state from the previously considered localized state to a delocalized state does not change the efficiency of the energy transport. The stability of the functional behavior with respect to changes in excitation energy dynamics can be explained by the separation of time scales between the excitation energy transport and the loss mechanisms. Excitation is lost in PSII on a time scale of ∼2 ns through a combination of fluorescence and non-radiative decay. Excitation transport, however, occurs on a ∼100 ps time scale. If we approximate the excitation transport to the radical pair states by a single rate constant, then the overall behavior is described by a two-rate (three-state) model—one rate for excitation to drive charge separation and a second rate of excitation undergoing non-radiative or radiative decay. Within this simplified model, the fraction of excitation that performs charge separation (as opposed to non-radiative or radiative decay) as time goes to infinity is given by

—the population reaching the RP2 state in the long-time limit—where *k*_{exc→RP} is the rate of population transfer to the RP states and *k*_{decay} is the rate of loss through non-radiative and radiative decay. We will use the phrase long-time efficiency to refer to this quantity *P*_{RP2}(*t* → ∞) and differentiate it from our 1 ns measure of efficiency, i.e., *P*_{RP2}(*t* = 1 ns), used in the remainder of the paper. In Figure 6(a) we see that at 100 ps the total excitation left on the pigments amounts to a population of *P*_{pig}(*t* = 0.1 ns) ≈ 0.2. Inserting this value in $ k exc \u2192 RP =\u2212 1 t ln ( P pig ( t ) \u2009 e k decay t ) $, we estimate the rate of excitation transport from the LHCII initial state (dashed lines) to the RP states as *k*_{exc→RP} ≈ 15 ns^{−1}, which, using Eq. (26), gives a long-time efficiency of 97%. To substantially change this value, the rate of excitation transport to the RP states would need to slow enough to become comparable to the 0.5 ns^{−1} rate of excitation loss. For example, slowing the rate of transport to the RP states by a factor of 2 only decreases the long-time efficiency to 93%. As a result, small modulations in the rate of transport to the RP states have only a minimal influence on the long-time efficiency within our model, and thus will not substantially shift the functional behavior of PSII.

Figure 5(b) shows the time-integrated population currents between the individual pigments for the delocalized initial state, and should be compared to Figure 5(a) for the localized initial state. The time-integrated currents inside LHCII and from LHCII to CP43 are smaller now than before. This is expected, since now only part of the initial population is in LHCII. Aside from these changes, however, the overall pattern of the pathways is very similar to what we have seen for the localized initial state before. The currents have the same directions as before, and even their relative magnitudes seem very similar. This is a remarkable result, since we know in this initial condition, every pigment has an equal initial excitation. Nevertheless, both the short-range (intra-complex) and long-range (inter-complex) energy transfer are very similar for both initial conditions. Interestingly, this suggests that while there is not a single obligate pathway of excitation energy transport within PSII—there is a select subset of PSII pigments that are used to create an excitation transport network. This observation may allow future work to differentiate between pigments that are involved in long-range transport from those that are predominately only associated with absorbing light and then transferring excitation locally (and hence showing little to no net flux in either of our initial conditions).

We have also performed simulations for other initial states, including localized, delocalized superpositions with coherence in the site or exciton basis, and for localized initial excitations with specific phase relations that were previously found to lead to directed transport for models of chains of coupled sites.^{68} In all these cases, we do not find any significant change in the overall efficiency of the transport to the RP2 state, showing a remarkable robustness of the overall long-range energy transport with respect to the initial conditions.

### B. Influence of the coupling to vibrations

We now investigate the influence of the coupling between electronic and vibrational degrees of freedom on the energy transfer. The coupling to the vibrations is treated approximately in our model so that it is important to consider how sensitive the energy transfer is to changes in these parameters. Furthermore, we expect the coupling to the vibrations to be an important component of the energy transfer dynamics, since it can enhance the transport, as is well known and has been studied for several smaller scale (single complex) systems (see, e.g., Refs. 59, 69, and 70).

In the following, we investigate these relationships by varying the strength of the coupling to the vibrations, that is, the magnitude of the environment spectral densities of the pigments Eq. (4). For all pigments, we simply multiply their spectral densities *J _{n}*(

*ω*) by a constant global factor, i.e.,

where the factor *c* is the same for all pigments *n* in all proteins. As the initial state, we take again the state where the excitation is localized on only the two pigments 7 and 10 in LHCII that we have used in Section III.

To investigate the effect of decreasing coupling to the vibrations on the energy transport, we decrease the coupling in stages. Figure 7 shows how the energy transfer changes when we decrease the coupling to the vibrations by factors of ten, i.e., *c* = 0.1, 0.01, 0.001, 0.0001. For each of these cases, Figure 7(a) shows the population of the RP2 state (blue, thick bar) and RP1 state and the total population in each of the complexes LHCII, CP43, and the RC (colored, thin bars) at a time of 1 ns. We see that when we decrease the vibrational coupling, the transport efficiency—the population of the RP2 state at 1 ns—for a decrease of the coupling by factor 0.1 first increases slightly. Then, as we further decrease the coupling to the vibrations, the efficiency of the transport decreases monotonically. We can also see that when the vibrational coupling is smaller and smaller, there is more and more population trapped in LHCII. It is remarkable that even when we decrease the coupling by a factor of hundred (i.e., *c* = 0.01), the amount of population transported to RP2 within 1 ns is still higher than 60% (compared to about 80% in the case of original coupling strength), showing that the energy transport in PSII appears to be very robust with respect to changes of the coupling to the vibrations. Only when we decrease the coupling further do we find that the energy transfer efficiency decreases significantly. In the extreme case where we decrease the coupling to the vibrations by a factor of thousand (i.e., *c* = 0.001), it goes down to less than 20% population on RP2 after 1 ns. Decreasing the coupling further (*c* = 0.0001) causes most of the population to be trapped in LHCII and only a very small amount of population (∼2%) reaches RP2. We can see in Figure 7(a) that the total length of the bars, i.e., the total population of the RP states and of the complexes together at 1 ns, is less than 1 (100%). The difference between the total length of the bars and 100% is the amount of population that is lost due to radiative and non-radiative decay to the ground state during the time of 1 ns. We see that this loss of population increases as the coupling to the vibrations is decreased, because the less population reaches the RP states, the more population remains in the excited electronic states of the pigments where it is subject to the decay.

For the two cases of *c* = 0.1 and *c* = 0.001 (the second and fourth bars from the top in Figure 7(a)), we take a closer look at the energy transport dynamics in Figures 7(b) and 7(c) (for *c* = 0.1) and Figures 7(d) and 7(e) (for *c* = 0.001). Figure 7(b) shows how for *c* = 0.1 the excitation is transferred slightly faster from LHCII to CP43 and the RC than for the original coupling (shown as dashed lines for comparison), leading to the slightly larger amount of population reaching RP2 within 1 ns that is observed in the bar diagram Figure 7(a). The population of pigment 7, which carries half of the initial population, shows oscillations now, indicating that part of the population is transported back and forth between pigments. But this back-and-forth oscillation of a relatively small amount of population does not make the overall transport less efficient. Also in the current from LHCII to CP43, shown in Figure 7(c), there are now oscillations in the magnitude of the current. But the current still is always directed from the LHCII towards CP43, i.e., there is no backflow to LHCII, and the current is positive at all times. On the other hand, the maxima of this current are more than twice as large as in the case of the original coupling to the vibrations. (See Figure 4(b) for comparison, and note the different scales on the y axes). The current from CP43 to the RC and from the RC to RP1 looks qualitatively similar to the original case. There are no strong fluctuations in the magnitude of these currents. But they are larger now compared to the original situation where the coupling to the vibrations is stronger.

The transport dynamics for the case of *c* = 0.001 are shown in Figures 7(d) and 7(e). In panel (d), aside from the much smaller amount of population transported to RP2, we can also see that it takes much longer for the population to be transported away from the LHCII. The reason for that becomes apparent when we look at the population of pigment 7 in LHCII: the population oscillates back and forth between the LHCII pigments many times, and is not able to leave the LHCII. It is trapped in the LHCII for quite a long time due to (unitary) localization inside the complex because of non-uniform pigment energies. In panel (e), we see that the current between LHCII and CP43 oscillates fast between positive and negative values. That means that excitation is flowing back and forth between these two proteins, quickly changing direction. The other currents are almost zero, because only a small amount of population can pass through CP43 and reach the reaction center.

Why the strong decrease of the vibrational coupling leads to such a trapping can be understood with the notion of vibrationally enhanced transport: the coupling to vibrations leads to transport through vibronic resonance and its absence leads to trapping. This concept is depicted in Figure 8: between the excited electronic states of the pigments there are energy gaps that, when they are large compared to the dipole-dipole interaction strength between the respective pigments, cause the excitation to be trapped on one or a few pigments, because of the off-resonance of the electronic transitions. But if the electronic excitation of a pigment couples to vibrational modes with vibrational energies large enough and with a coupling strength large enough compared to the energy gaps, then the vibrational energy levels can fill the gaps and make the transport of excitation from one pigment to another resonant, thereby enhancing the transport. Furthermore, the coupling to the continuum of vibrational modes leads to the relaxation down the energy manifold that is needed for directing the transfer to the RP states. Figure 8 depicts a situation in the limit where the electronic excitation is mostly localized on single pigments. We should note that one can of course also consider the vibrational enhancement of the transport in a picture in which the formation of delocalized, mixed vibronic states leads to resonances and stronger effective couplings between the delocalized vibronic states. The latter picture arises naturally in our model when one considers a delocalized basis instead of the localized representation shown in Figure 8.

For a delocalized initial state, in the case of very weak coupling to the vibrations, the trapping in the LHCII and CP43 will have smaller consequences for the transfer efficiency, because there is also initial excitation in the reaction center that can be transferred to the RP states. Further simulations have shown that for a delocalized initial condition the population that is transferred to RP2—in the case of the very weak coupling to the vibrations—can be three times higher than for the localized initial state in Figure 7(d). However, the population transferred to RP2 is still lower by more than a factor of two than when the coupling to the vibrations has the original strength. On the other hand, we have found in further simulations that for the original vibrational coupling *c* = 1, and when *c* = 0.1, 0.01, changing the initial state has little influence on the efficiency of the energy transfer. Only when the vibrational coupling is decreased beyond this, *c* = 10^{−3}–10^{−4}, does the initial state have a significant impact on the overall efficiency. This shows again a high under natural conditions robustness of the energy transfer to variations of the initial condition.

### C. Full Markovian Lindblad simulation: Larger range of coupling strength

So far, we have only *decreased* the coupling to the vibrations in relation to our original model. It would also be instructive to investigate the energy transfer at stronger coupling as well. Stronger coupling, however, leads to problems in the ZOFE calculation, because of the limited range of non-Markovianity for which the ZOFE approximation is valid.^{49,57,71} To still be able to gain an insight into the behavior at stronger coupling and—at least roughly—assess the energy transfer efficiency, we therefore use a full Markovian Lindblad calculation, where the Markov approximation is applied to the electron-vibration coupling. Based on the reasoning of Ref. 38, an assessment of the transport efficiency can be reasonable with an appropriate Markovian description. As a check, we also calculate the cases of weaker vibrational coupling again with the Markovian Lindblad simulation, so that we can compare with the previous ZOFE simulation results. A description of the Lindblad calculation is given in Appendix A.

Figure 9 summarizes the effect of stronger couplings with such a Lindblad calculation, where the strength of the coupling to the vibrational environment is varied by a factor between 0.0001 and 50. Here, factor 1 corresponds to the coupling strength in the original model. We see that as before for the ZOFE calculation, the Lindblad results also show the highest transport efficiency in the case where the vibrational coupling is decreased by a factor of ten with respect to the original model. In this case, also the amount of integrated population of the RP2 state is very similar to the ZOFE result. However, the RP2 populations for the original coupling and for the cases of weakest coupling, calculated with Lindblad, are in less good agreement with the ZOFE results.

Let us now look at the cases where the vibrational coupling is increased compared to the original model. In Figure 9, we see that as the coupling is increased, the transport efficiency becomes lower and lower, until for a fifty-fold increase, it is almost zero. At the same time the amount of population that stays in the LHCII becomes larger and larger. This suppression of the transport when the coupling to the vibrations is stronger may be understood in terms of the quantum Zeno effect,^{59} where the electronic excitation interacts with the vibrations—i.e., it is “measured” by the vibrations in the site basis—so strongly, that is, at such a high rate, that the excitation is forced to stay in the original state and cannot move anymore. In terms of the expression for the population currents, Equation (23), this effect can be understood from the fact that the stronger coupling to the vibrations causes stronger dephasing, destroying the coherence between the sites that drives the energy transfer, and thereby diminishes the energy transfer. Because of this effect at stronger coupling and the trapping of the excitation at weak coupling, the highest transport efficiency is obtained in the intermediate regime. We should emphasize, however, that the stronger coupling to the vibrations also increases the non-Markovian effects, causing the ZOFE approximation to eventually fail in this regime. Thus, we should expect that the Lindblad results, which neglect these non-Markovian influences altogether, are also increasingly inaccurate in this strong-coupling regime and can serve merely to gain a rough picture or “tendency” of what happens at stronger coupling. Potential resonance effects caused by individual strongly coupled vibrational modes at specific frequencies are neglected by the averaging character of the Lindblad description. In the intermediate regime, the magnitude of the coupling to the vibrations is such that it gives the optimal balance between vibrationally enhanced transport (as described in Figure 8) and the quantum Zeno effect. These results of the Lindblad calculation therefore provide more evidence that the strength of the coupling to the vibrations that we used in our original model lies roughly in the region of optimal coupling strength. However, remarkably, these results again show that the energy transfer efficiency could be even (slightly) higher at a factor ∼10 weaker coupling to the vibrations.

## V. SUMMARY AND CONCLUSIONS

In this work, we simulated the long-range inter-complex energy transfer in the PSII supercomplex—from LHCII via the core antenna CP43 to the reaction center, involving 33 pigments—by means of a full quantum calculation that uses a non-Markovian ZOFE quantum master equation description. We found that nearly all excitation reaches the reaction center within about 200 ps, and that after 1 ns, about 80% of the energy is transformed into charge separation in the reaction center. Less than 5% of excitation is lost through radiative and non-radiative decay within this time.

These findings for the overall (ns time scale) energy transfer dynamics are found to be in good agreement with the results of Ref. 16, where a modified-Redfield-generalized-Förster (MRGF) rate-equation method was applied to calculate the long-range energy transfer in the entire PSII supercomplex. In the MRGF method, strongly coupled pigments are assigned to domains treated with a modified-Redfield description, whereas excitation transfer between the domains, where the couplings are weaker, is treated with a generalized-Förster description.

Comparison of our ZOFE quantum-master-equation results for the excitation probabilities within each subcomplex of the truncated PSII supercomplex LHCII-CP43-RC with the corresponding results obtained with the MRGF rate-equation description of Ref. 16 showed differences in short time dynamics, but very good agreement for times beyond ∼200 ps. At shorter times, MRGF still reproduces the relative time scales of excitation energy transfer found with ZOFE but the subcomplex population dynamics show differences that reflect the more complete description of electronic coherence, both intra-complex and inter-complex, in the ZOFE description. On long time scales (beyond 200 ps), both ZOFE and MRGF dynamics are dominated by relatively slow irreversible transfer between radical pair states, the kinetic effect of which is to reduce the effect of differences at earlier times. This agreement of MRGF and ZOFE results for overall long time energy transfer indicates that the underlying PSII model possesses a separation of effective interaction strengths between the pigments that is sufficient for a reasonable subdivision into weakly coupled domains composed of strongly coupled pigments with enhanced inter-domain transfer rates resulting from coherent intra-domain exciton delocalization.^{41,63–67} This separation of interaction strengths is the essential ingredient needed for MRGF to work well as an effective approach to describe long range energy transport, and is equivalent to the separation of time scales that is explicitly observed in the full quantum calculations of the excitation transport dynamics over the three subcomplexes of PSII included in the current simulations.

We have further shown that excitation energy transfer is significantly slower within a pure Förster calculation, producing time scales of excitation transport that can differ by almost an order of magnitude compared to ZOFE, which further supports the importance of exciton delocalization and the coherence deriving from this for rapid excitation transport within PSII. The improved agreement between MRGF and ZOFE, compared to the agreement achieved with pure Förster calculations, is consistent with the findings of Ref. 38, where a full quantum description of electronic excitation transfer is compared to that of a classical rate equation, and it is shown that, provided this is properly derived from the full quantum description, the rate equation description is capable of capturing the overall features of the energy transfer.

In addition to calculating the usual time-dependent populations of electronically excited pigments, we also analyzed the energy transfer dynamics in terms of excitation population *currents* between the pigments. These currents are proportional to the coherences between the pigments. By employing the ZOFE quantum master equation we were then able to directly calculate the currents from the elements of the time-dependent electronic density matrix that represent the coherence. This analysis revealed the pathways and directions of the energy transfer between the individual pigments and the net flow of energy between the different complexes. We found that when initially all excitation is in the LHCII antenna complex, the net flow of energy is directed towards the reaction center at all times and no temporary inter-complex backflow occurs. Integrating the population currents over time showed that there is a complex network of pathways rather than just one or two dominant pathways, in accordance with the findings of Ref. 16. Interestingly, instead of being uniformly directed towards the reaction center, often excitation energy is transported in loops between pigments, i.e., excitation is transferred in a closed circuit through a series of neighboring pigments. Such loops were found in all three subcomplexes that we considered: LHCII, CP43, and the reaction center, and were shown to correlate with the domains of delocalized excitation used in the MRGF calculations. This raises the interesting question as to whether such ring currents might be detectable experimentally in the future and provide a new route to characterization of excitonic delocalization.

In the reaction center, the efficient excitation energy transport to the charge separated states appeared to rely mainly on only two pigments (30 (Pheo_{D1}) and 32 (Chl_{D1})), providing evidence that the rest of the reaction center pigments are needed primarily for their role in the electron transport rather than for the excitation transfer; some of these reaction center pigments appear to constitute a strong loop excitation current, but do not contribute to the net transport of excitation energy to the charge separated states. To test this conclusion, which was drawn from the analysis of the pathways of integrated population currents, we ran additional simulations with all reaction center pigments—except for 30 and 32—decoupled from the rest of the pigments so that transfer of excitation to and from these decoupled pigments was inhibited, but leaving the coupling to the charge separated states unmodified. The resulting efficiency of the energy transport to the charge separated states in this restricted model was found to be the same as in the original model where excitation transport to the other reaction center pigments was allowed, showing that the excitation energy transport to and inside the reaction center appears to rely indeed only on the two pigments 30 (Pheo_{D1}) and 32 (Chl_{D1}). This confirmed the conclusions drawn from the analysis of the integrated population currents and showed that this analysis is suited to identify transport pathways and assess their importance, providing a useful tool for the investigation of design-function relationships in light-harvesting systems.

To find out how the energy transfer in PSII depends on the initial condition—which is determined by the specific features of the light used for the excitation—we ran simulations of the energy transfer for different initial states; first we considered initial excitation that is localized on only two LHCII pigments; then we started from an initial state in which the excitation is completely delocalized over the LHCII, CP43, and reaction center pigments. We found that even though for the delocalized initial state the energy transfer dynamics of course look different from that for the localized initial state, the overall efficiency of the energy transfer to the charge separation in the reaction center is nevertheless the same in both cases. Also the overall pattern of the energy transfer pathways between the individual pigments is very similar for the different initial conditions—for both short-range (intra-complex) and long-range (inter-complex) energy transfer. Similar results were observed for a variety of other initial conditions. This shows that the energy transfer in PSII is remarkably robust with respect to the initial conditions of excitation.

To learn how sensitive the energy transfer dynamics are to variations in the coupling between the electronic and vibrational degrees of freedom, we varied the strength of this coupling to the vibrations in our simulations: decreasing the coupling to the vibrations for all pigments by a factor of ten compared to the original situation did not change the energy transfer efficiency much—it even slightly increased the efficiency, showing that the energy transfer is very robust to variations of the coupling to the vibrations. Even when the coupling to the vibrations was decreased by a factor of 100, the amount of energy transferred to the charge separation within 1 ns decreased only by about 10%. We assign this to a decrease of the well known effect of vibrationally enhanced transport, which helps the transfer of energy through more resonant energy levels when coupling to vibrations (vibrational energy levels) is present. Only when the electron-vibration coupling was decreased by a factor of 1000 does the amount of energy transferred decrease significantly (by up to a factor 4 compared to the original situation), due to trapping of the excitation in the LHCII complex. This trapping is characteristic of the predominantly unitary dynamics occurring in the absence of vibrationally enhanced transport. These observations show that for PSII the strength of the coupling to the vibrations that is needed for the effect of vibrational enhanced transport to direct the energy transfer efficiently to the reaction center is rather small.

In contrast to the Fenna-Matthews-Olson (FMO) complex, where calculations have revealed that the transport efficiency does not fall under 70% when the coupling to the vibrational environment is decreased, even to extremely small coupling strength,^{42} for PSII we found that the efficiency can become very low—almost zero—for extremely small coupling to the vibrations. This shows that unlike the situation in FMO, for PSII the coupling to the vibrations is crucial for achieving efficient transport. (We note that in Ref. 42 the long-time limit of the trapped population is used to measure the transport efficiency, whereas in the present work we measured the efficiency by the amount of population that is trapped in the reaction center after 1 ns, i.e., when the trapped population is still increasing.)

To assess how the energy transfer changes when the coupling to the vibrations is increased compared to the original situation, we performed Markovian Lindblad simulations, instead of the ZOFE calculations because the latter were not applicable in this parameter range. Here, we found that if the coupling to the vibrations is stronger than in the original situation and is further increased step by step, the efficiency of the energy transfer decreases monotonically, until for a fifty-fold increase of the coupling, almost no energy is transferred to the reaction center anymore (within the reference time of 1 ns) due to the strong dephasing that destroys the coherence needed for the energy transfer.^{48} These approximate Lindblad simulations, together with the ZOFE simulations, thus provide evidence that the original strength of the coupling to the vibrations lies in the parameter region for optimal energy transfer efficiency. We showed, however, that for a coupling that is about ten times weaker than in the original situation, the efficiency could be even slightly higher. This could be an important feature to reflect upon when considering the possible (evolutionary) optimization of natural light-harvesting systems like PSII.

The ability of the ZOFE approach to characterize electronic coherence in a photosynthetic system with very inhomogeneous excitonic couplings and finely tuned electron-phonon interaction characterized by non-Markovian effects demonstrated in this work suggests further application to more detailed modeling of the charge separation and charge transfer steps that follow electronic excitation energy transfer to the reaction center. The present work employed the simplified reaction center model from Ref. 16, which contains only two phenomenological, effective radical-pair states connected by incoherent rates, and does not explicitly include mixing of the radical-pair states with excitons.^{34} In future work, it will be important to improve this modeling of the reaction center with more detailed models and accurate parameters for the charge separation and transport steps and for their coupling to the prior exciton dynamics, in order to provide more detailed insight into the coherent dynamics of charge separation and charge transport subsequent to the excitonic energy transport studied here.

## Acknowledgments

This material was supported by DARPA under Award No. N66001-09-1-2026 and in part by the U.S. Department of Energy (DOE) through the SciDAC program at Lawrence Berkeley National Laboratory. J. J. R. thanks the Whaley group for helpful discussions. This work was posted to the Cornell preprint archive on January 27, 2015, with identifier arXiv:1501.06674. A related work focusing on the role of reorganization energy in energy transfer efficiency was subsequently posted at arXiv:1502.02657.

### APPENDIX A: PURE LINDBLAD SIMULATION OF EXCITATION ENERGY TRANSFER

To be able to roughly assess the efficiency of excitation energy transfer in our truncated PSII model for a range of the electron-vibration coupling that is extended significantly to higher coupling strengths, in Section IV C we used a pure Lindblad simulation for the excitation energy transfer, applying the Markov approximation to the electron-vibration coupling. Here, we briefly describe the calculation and how the corresponding parameters are obtained.

For the simulation, we employ a full Lindblad master equation, where we now use Markovian dephasing terms for the electron-vibrational coupling of the pigments, instead of the non-Markovian ZOFE terms used previously. The new dephasing terms describing the electron-vibrational coupling contain Lindblad operators *L _{n}* = |

*n*〉〈

*n*| and coupling parameters $ \gamma n Lind $ for each pigment

*n*, so that—instead of Eq. (21) for the ZOFE simulation—we now have the pure Lindblad master equation

where the terms in the second line again describe electronic relaxation processes (Section II B 3). In the dephasing term in the first line, the electronic excitation of each pigment *n* is coupled to vibrations with respective coupling strength $ \gamma n Lind $. To calculate the excitation energy transfer from the full Lindblad description (A1), we need to assume values for the coupling parameters $ \gamma n Lind $. In the following, we show how we obtain these parameters through a simple approximation from the parameters that we used to describe the non-Markovian electron-vibration coupling within the ZOFE description. To this end, let us again consider the special form of the environment correlation function (17) used in the ZOFE description

The pure Lindblad description (A1) is obtained in the Markov limit, where

i.e., *α _{n}*(

*t*−

*s*) decays fast compared to the time scales of the system dynamics. To express the parameters $ \gamma n Lind $ through the parameters of Eq. (A2) in this limit, we consider the Dirac series $ \delta \u03f5 ( x ) = 1 2 \u03f5 exp ( \u2212 | x | / \u03f5 ) $ for

*ϵ*→ 0, which lets us write

since $ O 0 ( n ) ( t , t ) = L n $ (see Eq. (16)). Since we know that in the Markov limit $ \alpha n ( t \u2212 s ) = \gamma n Lind \delta ( t \u2212 s ) $, the auxiliary operator of the ZOFE master equation becomes the constant operator $ O \u0304 0 ( n ) ( t ) = 1 2 \gamma n Lind L n $ (see Eqs. (14) and (16)),^{49} from comparison with Eq. (A6) it follows that we can approximate the needed coupling parameters for the full Lindblad descriptions as

i.e., in the case where the *γ _{nj}* are large compared to the system energies. We note that in order for the $ \gamma n Lind $ to stay finite in the

*γ*→ ∞ limit, not only the

_{nj}*γ*, but also at least one of the parameters Γ

_{nj}_{nj}has to go to infinity. Inserting the parameters for the electron-vibration coupling of the ZOFE simulation (according to Section II B 2) into Eq. (A7), we find the following values for the parameters $ \gamma n Lind $ that we used for our full Lindblad simulation in Section IV C:

The parameters *γ _{nj}* for the ZOFE electron-vibration coupling (Section II B 2) that we used here to calculate the Lindblad parameters $ \gamma n Lind $ are about

*γ*≈ 150 cm

_{nj}^{−1}. We can see from the fact that these

*γ*are not at all large compared to the system energies (the differences of the eigenenergies of the system Hamiltonian are in the range of up to 1300 cm

_{nj}^{−1}), that the Markov approximation is not well justified here, and that therefore our full Lindblad simulations should be considered to be a rather crude approximation of the non-Markovian dynamics.

To demonstrate this explicitly, in Figure 10 we show the transport dynamics resulting from the Lindblad calculation described above, together with the ZOFE and pure Förster results from Figure 4(a). It is evident that there is considerable deviation of the Lindblad results (dashed lines) from the ZOFE results (solid lines), illustrating the inaccuracy of the Lindblad method in this parameter regime. For the time-dependent populations of LHCII and CP43, the deviation of Lindblad from ZOFE results is even larger than the deviation of the pure Förster (dotted lines) from ZOFE results. Remarkably however, the populations of the two radical pair states RP1 and RP2 that result from the Lindblad calculation agree very closely with the populations from the pure Förster results, despite the deviations between Lindblad and pure Förster in the LHCII and CP43 populations. This shows once again that because of a separation of time scales, here between the fast transport to the reaction center and the slower subsequent transport to the RP1, RP2 states, the populations of RP1, RP2 can be quite robust with respect to variations in the preceding exciton dynamics.

### APPENDIX B: AMOUNT OF COHERENCE INVOLVED IN THE ENERGY TRANSFER

According to the population current description of Equation (23), the *imaginary component* of coherence between pigments is responsible for the population currents through PSII. In the ongoing discussion about the role of coherence in the energy transfer in biological systems, coherence and its contribution to the energy transfer have so far been quantified according to different metrics that generally do not look separately at the imaginary and real parts of the coherence (see, e.g., Refs. 41 and 72 for quantification and discussion of coherence). However, as described in detail in Ref. 48, the imaginary and the real part of the coherence have very different effects—the imaginary part drives the transfer, whereas the real part is related to the transfer-diminishing localization effects in the presence of energy gaps between the pigments. Therefore, we shall not only quantify the total coherence, but also look separately at its imaginary and real parts.

Solving the ZOFE master equation for the full quantum dynamics gives the time-dependent electronic density matrix and therefore allows direct quantification of the coherences, the off-diagonal elements of the density matrix in the respective basis of interest. We first quantify the total amount of coherence by the sum

over the absolute values of all off-diagonal elements of the density matrix *ρ*, which provides an intuitive measure for the overall amount of coherence. This “*l*_{1}-norm of coherence” has been widely used in previous studies and is reviewed in Ref. 72.

Aside from the coherence in the site and exciton bases, we are also interested in the coherence in the domain-exciton basis of Ref. 16. This is of interest since in the MRGF calculation of Ref. 16, a classical rate equation is used to describe long-range intercomplex transfer that explicitly contains only the populations in the domain-exciton basis, but not the coherences. Thus, coherence in this basis is not explicitly taken into account in the calculation giving the population dynamics in Figure 4(a) (dashed-dotted lines), yet good agreement with the ZOFE population dynamics is nevertheless achieved. This is consistent with Ref. 38, where it is discussed that the overall population dynamics, or at least the energy transfer efficiency, that is, the population transferred to the radical pair states, should be retained under a proper transition from the full quantum description to the classical rate description where coherence is not explicitly taken into account.

In Figure 11(a), we show the total amount of coherence measured by Equation (B1) for the site basis, the exciton basis, and the domain-exciton basis. The domain-exciton coherence shows the coherence not (explicitly) taken into account in the MRGF calculations of Ref. 16.

Since the initial state is an incoherent state in the site basis, at time zero there is no coherence in the site basis. Similarly, in the domain-exciton basis, the coherence is initially zero. (Sites 7 and 10, which carry the initial excitation, belong to the same domain that consists of only these two sites.)

In the exciton basis, on the other hand, there is initial coherence at time zero, because the excitation localized on sites 7 and 10 is expressed through a coherent superposition of exciton states that each have part of their excitation on sites 7 and 10, but are also delocalized over other sites, mainly sites of the LHCII antenna protein.

In both site and domain-exciton basis, the amount of coherence rises initially, because the dipole-dipole interaction between the pigments drives the initially localized state towards a more delocalized state. During the transport dynamics, the total electronic coherence decreases due to interaction with the vibrations, and later due to the accumulation of population in the radical pair states.

Next, we look at the time-integrated coherences between the individual pigments and also separately at the imaginary and real parts of the coherence, because of their different role in the energy transfer. To this end, we time-integrate the absolute value of the imaginary and real parts of the elements of the density matrix in the site basis during the first 1 ns of its propagation (at which point, as we have shown in Figure 4(a), the majority of the excitation has been irreversibly trapped at RP2). In Figure 11(b), the resulting site-basis elements of the time-integrated density matrix are shown as a color matrix. The time-integrated pigment populations are the diagonal elements. One can see that there is an accumulation—a “damming-up”—of population on the energetically lowest sites of each protein. Such a damming-up effect occurs because the transport of excitation that is fast inside the complexes, where the couplings between the pigments are relatively strong, is slowed down between the complexes. (It also depends on the coupling to the vibrations and also on the gap between the “donor” protein’s pigment energies and the “acceptor” protein’s pigment energies.) The upper triangle in Figure 11(b) shows the (absolute values of) the imaginary parts of the coherences between the sites integrated over time. That is, it shows ∫|Im(*ρ*_{n,m≠n})| *dt*. As seen in Equation (23), the imaginary parts of the coherences give the unitary contribution to the population currents between the sites. That is, the imaginary parts of the coherence constitute the whole of the energy transfer, except for the transfer from the reaction center to the radical pair states, which is described by relaxation.

The lower triangle in Figure 11(b) shows the real parts of the coherence, that is, ∫|Re(*ρ*_{n,m≠n})| *dt*. We see that between many sites, the real part of the coherence is larger than the imaginary part. In the analysis in Ref. 48, it is shown that the real part of the coherence is related to localization effects in the energy transfer that stem from the electronic energy gaps between the pigments, adding to localization stemming from other influences, e.g., coupling to vibrations. The smaller the real part of the coherence, the smaller are these localization effects that can hinder the energy transfer. (Elements that are very small are ignored in the matrix shown in Figure 11(b).) In contrast to the real parts, the imaginary parts of the coherences directly induce transport between the pigments, as described by the first term of Equation (23) and discussed in more detail in Ref. 48. There are relatively strong imaginary parts of coherences visible in the blocks that show the coherence between LHCII sites and CP43 sites, and also between CP43 sites and RC sites. These are the coherences that are responsible for the currents between these complexes, shown by the green and red curves in Figure 4(b). In Figure 11(b), we see that there are even non-zero coherences between LHCII sites and reaction center sites. However, these are very weak and the resulting net population current that they cause between LHCII and the RC appears to be close to zero on the scale of Figure 4(b).

### APPENDIX C: EFFECT OF COUPLING TO THE RADICAL PAIR STATES

To see the influence of the radical pair states on the transport dynamics, we now look at the dynamics with the coupling to the RP states switched off. That is, there is no transfer of population to the radical pair states.

Figure 12 shows the resulting dynamics for the same quantities as in Figure 4, but with the coupling to the RP states switched off. We see in Figure 12(a) that up to a certain time, the populations in the proteins are the same, with or without this coupling. The further a protein is away from the RP states, the later the population without the coupling deviates from the population where the coupling is present. Instead of falling off fast, now—without the coupling to the RP states—after around 200 ps the populations decline much slower; this residual, slower decline is due to the radiative and non-radiative decay of the pigment excited states. After 200 ps, the residual, average population per pigment is highest in the RC, and it is lowest in LHCII. This is because the average energy in the RC is lowest and that in LHCII is highest, providing the “energy funnel” that directs the transport to the RC.

Let us next look at the coherence. When we compare the total coherence in Figure 12(b), where the coupling to the RP states is switched off, with that in Figure 11(a), where the coupling to the RP states is present, we see that up to about 10 ps the coherence with and without coupling is equal. After that, however, the coherence without coupling to the RP states goes to finite values that decrease only slowly, whereas in the case when the coupling is present, the coherence falls off to zero. This behavior can be explained with the Cauchy-Schwarz inequality

that holds for the coherences *ρ _{ij}*(

*t*) and the populations

*ρ*(

_{ii}*t*),

*ρ*(

_{jj}*t*) of any density matrix

*ρ*(

*t*). It says that coherences can only ever arise between states that are populated. The larger the populations of two given states are, the larger can the coherence between these states potentially be. In the case where there is transfer to the RP states, population is taken away from the excited states of the pigments by the RP states, and consequently, the coherence, being limited by the population through the Cauchy-Schwarz inequality, has to fall off to zero. On the other hand, when the transfer to the RP states is switched off, the population remains in the excited states of the pigments and the coherence stays finite.

### APPENDIX D: DETERMINATION OF RELEVANT ENERGY RANGE OF BATH SPECTRAL DENSITY FOR EXCITONIC TRANSPORT

The model for energy transport described in Section II employs an open quantum system approach in which the coupling between the electronic degrees of freedom (system) and the vibrational degrees of freedom (environment/bath) at each pigment is described by a bath spectral density. To identify the energy region of the bath spectral densities that is relevant for the transport dynamics, we analyze the coupling of the transitions between eigenstates of the system Hamiltonian (exciton states) to the environment/bath (vibrational degrees of freedom).

In the total Hamiltonian *H*_{tot} = *H*_{sys} + *H*_{sys–env} + *H*_{env} of Section II, the system-bath interaction is written in the standard form

describing linear coupling between system operators *L _{n}* = − |

*n*〉〈

*n*| and bath operators $ B n \u2261 \u2211 \lambda \kappa n \lambda ( a n \lambda \u2020 + a n \lambda ) $ for each pigment

*n*. The system operator

*L*represents electronic excitation of pigment

_{n}*n*and the bath operator

*B*represents a weighted sum of the linear coupling of the vibrational modes

_{n}*λ*to the electronic excitation at this site, where the individual weights are given by the parameters

*κ*

_{nλ}. These parameters determine the bath spectral density of pigment

*n*, $ J n ( \omega ) = \u2211 \lambda | \kappa n \lambda | 2 \delta ( \omega \u2212 \omega n \lambda ) $ (see Eq. (4)).

We now express the system operators *L _{n}* in terms of exciton states |Φ

_{k}〉, i.e.,

*H*

_{sys}|Φ

_{k}〉 =

*E*|Φ

_{k}_{k}〉, which allows the system-bath interaction Hamiltonian to be written as

with new system operators $ L \u0303 l k \u2261| \Phi k \u3009\u3008 \Phi l |$ and factors $ f l k n \u2261\u3008 \Phi k |n\u3009\u3008n| \Phi l \u3009$. The latter can be assumed to be real, i.e., $ f l k n = f k l n $, since *H*_{sys} is real and symmetric. This representation of the system-bath interaction Hamiltonian provides important insight into the influence of the electron-vibrational coupling on the energy transport. It shows that each transition between a pair of exciton states |Φ_{l}〉 and |Φ_{k}〉, described by the operator $ L \u0303 l k $, couples to the vibrational bath of each pigment *n* with a relative coupling strength (weighting factor) $ f l k n $ that is given by the overlap of these two exciton states with the localized excitation state |*n*〉. It is this coupling that induces the relaxation of excitation in the excitonic manifold that is essential for energy transfer towards the reaction center. Together with the coupling parameters *κ*_{nλ}, contained in *B _{n}*, i.e., together with the bath spectral density, the relative weights $ f l k n $ determine the overall exciton-vibrational coupling for each excitonic transition involved in relaxation. Each excitonic transition from |Φ

_{l}〉 to |Φ

_{k}〉 has a different energy Δ

_{lk}=

*E*−

_{l}*E*. Taking these energies together with the corresponding relative coupling strengths $ f l k n $ for the excitonic transitions reveals the energy regions in which the bath spectral densities couple significantly to the excitonic transitions, and thus shows the energy range of the bath spectral density that is relevant for the energy transfer dynamics.

_{k}To quantify this, we construct the following measure of the exciton-vibration coupling strength density as a function of excitonic transition energy *ϵ _{i}* = Δ

_{lk}:

where the spacing between the consecutive transition energies, Δ*ϵ _{i}* ≡ (

*ϵ*

_{i+1}−

*ϵ*

_{i−1})/2 (and Δ

*ϵ*

_{1}=

*ϵ*

_{2}−

*ϵ*

_{1}, Δ

*ϵ*=

_{N}*ϵ*−

_{N}*ϵ*

_{N−1}at the boundaries), takes into account the variable density of excitonic transitions in energy. The resulting exciton-vibration coupling density functions are shown as bars with height

*F*(

_{n}*ϵ*) at the top of each panel in Figure 3. These bars are the collection of the

_{i}*F*(

_{n}*ϵ*) for all pigments

_{i}*n*, where pigments with the same bath spectral density have bars of the same color (blue: Chl B pigments in LHCII, green: Chl A pigments in LHCII, red: pigments in CP43, and orange: pigments in RC). One sees that the excitonic transitions carrying significant coupling strength are only present in the energy region up to 500 cm

^{−1}, which indicates that above this energy the bath spectral density is not important for the energy transfer.

## REFERENCES

We employ a slightly modified set of parameters supplied by V. I. Novodorezhkin to D. I. G. Bennett that is reproduced in the supplementary material of Ref. 16.