Longrange coherent energy transport in Photosystem II
Abstract
We simulate the longrange intercomplex electronic energy transfer in Photosystem II – from the antenna complex, via a core complex, to the reaction center – using a nonMarkovian (ZOFE) quantum master equation description that allows us to quantify the electronic coherence involved in the energy transfer. We identify the pathways of the energy transfer in the 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 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 vary the strength of the coupling between the electronic and the vibrational degrees of freedom. We find that the original parameters lie in a (broad) region that enables optimal transfer efficiency, and that the energy transfer appears to be very robust with respect to variations in the vibronic coupling. Nevertheless, vibrationally enhanced transfer appears to be crucial to obtain a high transfer efficiency. We compare our quantum simulation to a “classical” rate equation based on a modifiedRedfield/generalizedFörster description that was previously used to simulate energy transfer dynamics in the entire Photosystem II complex, and find very good agreement between quantum and rateequation simulation of the overall energy transfer dynamics.
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 triphosphate (ATP), the “chemical currency” of biological reactions Knowles [1980]. Greater than 80% of all excitations initiated in PSII will result in productive photochemistry, causing the creation of an irreversibly separated electronhole pair contributing to water splitting Caffarri et al. [2011]. The mechanisms of efficient energy transport within PSII have been studied using spectroscopy Ruban and Murchie [2012], Belgio et al. [2012], Marin et al. [2012], Krüger et al. [2012, 2010], SchlauCohen et al. [2014], Novoderezhkin et al. [2011a, b], Fuller et al. [2014], Wells et al. [2014], SchlauCohen et al. [2012], Broess et al. [2006] and numerical simulations Raszewski and Renger [2008], Bennett et al. [2013], Kreisbeck et al. [2014], Shibata et al. [2013], with models based on spectroscopic data and structural information from Xray crystal structures of the complexes Umena et al. [2011], Liu et al. [2004], Pan et al. [2011], Müh et al. [2014]. 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 Caffarri et al. [2011], van Oort et al. [2010], Caffarri et al. [2009]. The largest PSII supercomplex isolated by Caffarri and coworkers Caffarri et al. [2009] (called CSM) contains 326 pigments bound into an assembly of light harvesting proteins surrounding a PSII reaction center. The CSM 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 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 rateequation descriptions Raszewski and Renger [2008], Bennett et al. [2013], Shibata et al. [2013] rather than the full quantum dynamics simulations that have been performed for smaller lightharvesting systems such as the FennaMatthewsOlson complex Rebentrost et al. [2009a], Ishizaki and Fleming [2009], Chin et al. [2010], Kreisbeck et al. [2011], Ritschel et al. [2011a, b]. Recently, in Ref. Bennett et al. [2013], simulations based on a combined modifiedRedfieldgeneralizedFörster (MRGF) rate equation approach were carried out for the entire PSII supercomplex, providing new insight into the interplay of shortrange transfer dynamics inside the individual subcomplexes (proteins) and longrange intercomplex dynamics within a rateequation kinetic analysis.
In recent years, spectroscopy experiments have provided evidence that there is coherence involved in excitonic energy transfer in lightharvesting complexes Engel et al. [2007], Romero et al. [2014]. This observation raises questions about the role of such coherence in the energy transfer. 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 OlayaCastro et al. [2008], Scholes et al. [2011], Hoyer et al. [2012], Jesenko and Žnidarič [2013], Chenu et al. [2014], Scholes et al. [2012], Kassal et al. [2013], Rebentrost et al. [2009b], Hoyer et al. [2010], Shim et al. [2012], Abramavicius and Mukamel [2010], Sarovar and Whaley [2013], Chenu and Scholes [2015]. In a companion paper Roden and Whaley [2015] of the present paper, the role of electronic coherence in electronic excitation energy transfer is analyzed in the framework of excitation probability currents.
In the present paper, we investigate the energy transfer dynamics in PSII by means of a full quantum simulation, using a nonMarkovian quantum master equation description. This not only allows us to quantify the electronic coherence involved in the energy transfer over length scales including several subcomplexes, but also to calculate its contribution to the energy transfer in terms of excitation probability currents, using the framework of Ref. Roden and Whaley [2015]. The analysis of the probability currents also provides insight into the pathways of the energy transfer in the large network of coupled pigments. Using a nonMarkovian quantum master equation allows us to account for the nonnegligible memory times, i.e., the nonMarkovianity, associated with both the intramolecular vibrations of the pigments and the vibrations of the protein environment of the pigments to which the electronic degrees of freedom couple. Specifically, we investigate the longrange intercomplex energy transfer from the antenna complex via a core complex to the reaction center. To this end, we simulate the energy transfer in a subsystem containing a LHCII monomer, the core complex CP43, and the reaction center (see Figure 1). This requires incorporation of pigments in the simulation. To be able to efficiently treat nonMarkovian dynamics for this number of pigments, we use the nonMarkovian ZOFE quantum master equation that has been developed in recent years Strunz and Yu [2004], successfully applied to biological lightharvesting complexes, and also tested against exact results Ritschel et al. [2011a, b]. This approach allows calculations for wide parameter ranges of the model. To compare the simulation results of the ZOFE quantum master equation and the MRGF rate equation of Ref. Bennett et al. [2013], we analyze here the results of both methods for excitonic energy transfer in the subsystem of PSII described above.
There is an ongoing discussion about whether and how energy transfer in lightharvesting 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 Kassal et al. [2013], Jesenko and Žnidarič [2013], Han et al. [2013], Dorfman et al. [2013]. Motivated by these discussions, in the present paper we simulate the energy transfer dynamics for very different initial states to ascertain how longrange energy transfer in PSII depends on the initial condition. We consider first initial excitation that is localized in the antenna complex, then look at how the energy transfer changes when all pigments in the antenna complex, core complex, and the reaction center carry the initial electronic excitation, i.e., when the initial state is completely delocalized over all three complexes. We also analyze the role of vibrations – in particular, intramolecular 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. In the present work, we vary the coupling between the electronic and the vibrational degrees of freedom to investigate how the energy transfer depends on this coupling to the vibrations.
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
We aim to investigate the longrange 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. Caffarri et al. [2009] and subsequently used as a foundation for a structure based excitation energy transport model to explain the measured fluorescence lifetimes Caffarri et al. [2011]. Figure 1 shows the largest PSII supercomplex previously modeled Caffarri et al. [2011], Bennett et al. [2013]. 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 reation center (6 pigments).
ii.1 System parameters
ii.1.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,
(1) 
containing site energies () and the coupling matrix . Here, is the state in which only pigment is excited and all others are in the ground state. The construction of an electronic Hamiltonian for the PSII supercomplex has been described previously Bennett et al. [2013]. The couplings between pigments contained within the same protein have been extracted from the literature for each complex Bennett et al. [2013], where they were constructed to reproduce the available linear and nonlinear spectroscopy data for ensembles of the isolated pigmentprotein complex. The coupling between pigments in different proteins has been calculated assuming a dipoledipole coupling Bennett et al. [2013]. All energies and coupling parameters are taken from Ref. Bennett et al. [2013], and are shown in Figure 2, together with the structure of the LHCIICP43RC truncated supercomplex.
ii.1.2 Coupling to vibrations
The electronic states of each chromophore are coupled to a large collection of vibrational modes that represent intramolecular vibrational modes as well as vibrational modes of the surrounding protein scaffold. In the extreme limit of lowfrequency vibrations, the longtimescale protein conformation dynamics give rise to an inhomogeneous distribution of energies for each chromophore within any ensemble measurement. In contrast to Ref. Bennett et al. [2013], 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 “bestcase scenario” for transport efficiency and emergence of coherence.
We describe here the intramolecular vibrational modes of the pigments together with the vibrational modes of the protein environment by the single vibrational Hamiltonian
(2) 
where we approximate the vibrational modes to be harmonic. Here, is the annihilation operator of mode with frequency that belongs to a pigment . (Here and throughout the paper we set ). 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 dipoledipole interaction between two pigments.
We assume that electronic excitation of a pigment couples linearly to its own vibrational modes, such that the overall coupling between the electronic degrees of freedom and the vibrations is given by
(3) 
where the system operators are the projectors . The coupling constants that describe the strength of the coupling of electronic excitation of pigment to vibrational mode of this pigment, can be expressed through the vibrational spectral density of pigment
(4) 
Assuming that the spectral density is a continuous function, this leads to the vibrational correlation function
(5) 
which depends on the temperature . The Fourier transform
(6) 
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 (or vice versa) via the relation Mukamel [1982]
(7) 
with the antisymmetrized spectral density
(8) 
Peaks in that are narrow compared to the relevant energies of the electronic degrees of freedom will lead to nonMarkovian dynamics. When we apply the ZOFE master equation below, we will fit using a sum of Lorentzians.
ii.1.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 to make it a reasonable model of the initial stages of photosynthesis in PSII. These are, i) energy capture via charge separation at the reaction center, and ii) loss of excitation via nonradiative 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 ( and ), as described in Ref. Bennett et al. [2013], and the ground state , 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. The characteristic transfer times that are found in Ref. Bennett et al. [2013] 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 nonradiative decay mechanisms, with a characteristic decay time of 1.78 ns (combining the decay times of 2 ns for nonradiative decay and 16 ns for fluorescence used in Ref. Bennett et al. [2013]).
The Hamiltonian term describing the radical pair states and the ground state is given by
(9) 
where is the ground state energy, and and are the energies of the RP1 and RP2 state.
We shall describe in detail the dynamical model employed to treat excitation transfer into either the radical pair or ground state of the pigment in Section II.2.3.
ii.2 Excitation energy transport
Based on the previous section (Section II.1), the total Hamiltonian of our model is then
(10) 
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 (), an environment (), and the interaction between system and environment (). The system part 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 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
(11) 
The excitation energy transport and the electronic coherence can be extracted from the reduced density matrix of the system part, i.e.,
(12) 
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 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 offdiagonal elements are the electronic coherences between the pigments, which constitute a major focus of the present work.
ii.2.1 ZOFE master equation
For the simulations of the dynamics of energy transport and coherence we use the ZOFE master equation, which allows us to take into account the nonMarkovian effects in the coupling between the electronic degrees of freedom and the vibrational environment. The ZOFE master equation is given by Strunz and Yu [2004], Ritschel et al. [2011a]
(13) 
where the auxiliary operator
(14) 
captures the nonMarkovian effects, and with initial condition
(15) 
Depending on the form of the environment correlation function , a closed evolution equation for the auxiliary operator can be obtained Yu et al. [1999], Strunz and Yu [2004], Ritschel et al. [2011a]. For an environment correlation function that is a sum of damped oscillating terms,
(16) 
such a closed evolution equation can be found Ritschel et al. [2011a]. Writing the operator as the sum , the closed evolution equation to obtain the operators is given by Ritschel et al. [2011a]
(17) 
with initial condition and . In the present work, we use these evolution equations with coupling operators for each pigment .
ii.2.2 Approximate environment spectral density for ZOFE master equation
The form of in Equation (16) corresponds to a sum of Lorentzians centered at frequencies and with weights and widths in the environment correlation spectrum (the Fourier transform of ):
(18) 
We use this expression to fit the environment correlation spectra of the pigments to experimentbased spectra, in order to obtain the parameters , , and that we need when we solve the evolution equation for the auxiliary operator Eq. (17). The resulting spectral densities (Eq. (4)) that we use in the ZOFE simulation are shown in Figure 3, together with the original, measurementbased spectral densities of Ref. Bennett et al. [2013].
As one can see in Figure 3A, the fitted spectral density that we use for the LHCII pigments in the ZOFE simulation is very broad compared to the narrow peaks of the original spectral density. We approximate the original spectral density with these broad peaks, because the ZOFE approximation that we use in the simulation is quite limited in the degree of nonMarkovianity that it can handle. Thus, the memory time of the vibrations can not be too long compared to the time scale of the electronic dynamics. Narrow peaks would result in a long memory time, so we are limited to broad peaks that capture only the general features of the coupling to the vibrations.
The approximate spectral densites in Figure 3A also neglect the highenergy contributions of the original spectral density around and above 1,500 cm. Since these peaks are offresonant with the system energies, which are marked by the bars at the upper edge, we assume that they do not have a large influence on the system dynamics and it is reasonable to neglect them. In Ref. Ritschel et al. [2011a], it was shown that the influence of offresonant parts of the spectral density on excitation energy transport is indeed negligible.
In Figure 3B, the spectral densities for CP43 and the RC that we use for the ZOFE simulation are shown. They are only in rough agreement with the original, measurementbased spectral densities from Ref. Bennett et al. [2013]. This difference arises because these spectral densities are obtained not by directly fitting the spectral densities , but by fitting the environment correlation spectra with Lorentzians, according to Equation (18). To summarize, we apply the following procedure:

Fit these with Lorentzians according to Equation (18), to determine the parameters for the evolution equations (17) for the ZOFE auxiliary operators. In this work, we use only two Lorentzians for the LHCII pigments and one Lorentzian for the CP43 and RC pigments, to keep the corresponding number of evolution equatons small (see index in Eq. (17)).
As a consistency check, we convert these fitted sums of Lorentzians back to spectral densities to compare with the original spectral densities, as is done in Figure 3.
ii.2.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 nonradiative decay of the excited electronic states of the pigments into account, we add Lindbladmasterequation terms to the ZOFE master equation. An analogous treatment of radiative decay and trapping in the reaction center is employed in Ref. Kreisbeck et al. [2011], where a nonMarkovian master equation is also extended by corresponding Lindblad terms, for the simulation of energy transport in the FennaMatthewsOlson complex. (See also Refs. Caruso et al. [2009], Jesenko and Žnidarič [2013] 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. We assume that the characteristic transfer times to the radical pair states that are found in Ref. Bennett et al. [2013] give a reasonable approximation and therefore describe the trapping part of our simulation with Lindblad terms based on these characteristic transfer times.
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 Strunz and Yu [2004]. 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 , the auxiliary operator of the ZOFE master equation becomes the constant operator (see Eqs. (14) and (15)) Strunz and Yu [2004]. (Here, we have replaced the index , which is used above to refer to the pigments, by a more general index , to refer to the respective process described by a system operator and a corresponding coupling constant ; e.g. relaxation from a state to a state is described by .) Inserting this constant auxiliary operator into the ZOFE master equation (13) gives the wellknown Lindblad master equation
(19) 
with coupling (Lindblad) operators and corresponding coupling constants . We use such a Lindblad description for the following processes:

Radiative and nonradiative decay of the electronic excitation of all pigments () is described through Lindblad operators and the corresponding coupling constant is given by the decay rate. Here, 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 nonradiative decay and 16 ns for fluorescent decay assumed in Ref. Bennett et al. [2013]).

Transfer of excitation from the reaction center pigments 27, 28, and 32 to radical pair state RP1 through with rates (where ). 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. Bennett et al. [2013].

Transfer from RP1 back to the reaction center pigments 27, 28, 32 through with rates . 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. Bennett et al. [2013].

Irreversible transfer from RP1 to radical pair state RP2 through with a rate . This rate is chosen according to the characteristic transfer time of 520 ps found in Ref. Bennett et al. [2013].
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
(20) 
where the Lindblad operators and corresponding coupling constants (rates) in the Lindblad term in the second line of the equation belong to the respective processes that are listed above, and the index of the sum refers to the individual processes. In the ZOFE term in the first line, the ZOFE coupling operators couple only electronic excitation of the pigments to the nonMarkovian vibrational environment. The radical pair states and are not coupled to the nonMarkovian vibrations.
The system density matrix now has components for the states , in which pigment is excited, as well as components for the radical pair states and and for the ground state . is obtained by solving this master equation together with the evolution equation for the ZOFE auxiliary operators. For the truncated supercomplex considered here, is then a 36 dimensional square matrix. The system Hamiltonian does not contain coupling to or between RP1, RP2, and the ground state; all coupling to and between these states arises through the Lindblad terms.
ii.2.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 between the pigments. As shown in a companion paper Roden and Whaley [2015], 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 Roden and Whaley [2015]. Another advantage of the populationcurrent 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 lightharvesting 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 at all times – a continuity equation Roden and Whaley [2015]
(21) 
holds, where is the (net) population current at time that transports population (i.e., excitation probability) from a pigment to another pigment . When is positive, excitation is transported from pigment to pigment ; when it is negative, the energy transfer goes in the other direction, i.e., from pigment to pigment .
Based on the continuity equation (21), for quantum master equations of the form used in the present work Eq. (20), it is shown in Ref. Roden and Whaley [2015] that the individual population currents between the pigments can be calculated through
(22) 
where , , and are the respective contributions of unitary dynamics, dephasing, and electronic relaxation, which are specified in the second line of the equation. Here, are the coherences between the sites, are the populations of the sites, and are the intersite couplings and is the rate of electronic relaxation from a site (or state) to a site .
When we have the density matrix at a time from a numerical simulation as described above, we can use Equation (II.2.4) to calculate the currents between the sites. As described in detail in Ref. Roden and Whaley [2015], the first term in Equation (II.2.4) 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 nonradiative decay to the ground state.
We see in Equation (II.2.4) 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 (II.2.4). 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.
When the coupling to the vibrations is strong compared to the electronic interpigment coupling, the resulting strong dephasing destroys a large part of the coherence between the pigments. We see from Equation (II.2.4) 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 Rebentrost et al. [2009c]: 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 between individual sites and of the different proteins and then sum these currents up. That is, the current between subcomplex and a subcomplex is given by
(23) 
where when the current is positive, there is a net flow from to , and when is negative, there is a net flow from to .
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 (LHCIIS and LHCIIM). Our minimal model of the process of excitation energy transfer from the surrounding light harvesting proteins to the reactions center contains one LHCII monomer, the complete CP43 protein and a slightly truncated RC. We use the ZOFE master equation to calculate the timedependent density matrix in the space of the electronic states of the pigments, as described in the previous section. This enables us to analyze the transfer of the electronic excitation energy and the electronic coherence between the pigments.
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 Blankenship [2014]. As such, in this paper, we will 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 efficiency of the transport, since this state undergoes irreversible charge separation in our model. (We note, however, that fluorescencelifetime measurements on PSII are only sensitive to the total population remaining in chlorophyll excitation. As a result, the experiment cannot directly differentiate the kinetics of electron transport, and the assignments to RP1 and RP2 must be considered phenomenological Bennett et al. [2013].)
In our first study, we describe excitation energy transfer when initial excitation occurs in a pigmentprotein complex outside the supercomplex and subsequent energy transfer occurs into the core complex via the attached LHCII trimer. We therefore initiate excitation on two chlorophyll 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 Bennett et al. [2013]. 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. Thus, all offdiagonal elements of the initial density matrix are zero in the site basis. The dynamics following delocalized initial excitation will be discussed in Section IV.1.
iii.1 Population
To investigate the (longrange) 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 and the population of the RP states over time. Figure 4A shows these population dynamics.
The solid lines in Figure 4A are the populations on each complex calculated with the ZOFE master equation. The population of site 7 in LHCII is shown separately (thin green line). Nearly all 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 nonradiative 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 timescales between excitation energy reaching the RC and the relatively slow loss mechanisms of fluorescence and nonradiative decay is a key feature of efficient light harvesting in PSII.
The dashed lines in Figure 4A show the population calculated with the modifiedRedfield/generalizedFörster (MRGF) method of Ref. Bennett et al. [2013]. In the MRGF calculations, pigments are grouped together into “domains”. The assignment of pigments to different domains is determined by a threshold electronic coupling and the distribution of site energies. 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 Bennett et al. [2013].
As we can see in Figure 4A, the population dynamics from the MRGF calculation – especially the RP2 population – agree fairly well with that of the ZOFE simulation, even though the former is a much simpler classical rate equation calculation. This is consistent with the findings of Ref. Jesenko and Žnidarič [2013] that classical rate equations are adequate for the calculation of transport efficiencies when they are properly derived from the full quantum description.
In Figure 4A, 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) Förster [1948], May and Kühn [2008]. 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. Bennett et al. [2013].
In Figure 4A 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. The dynamic effect is sometimes 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. Kassal et al. [2013] and references therein). This enhancing effect of supertransfer due to 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.
iii.2 Population currents
Figure 4B shows the population currents between the proteins and from the reaction center to the radical pair states. These currents are calculated from the timedependent density matrix using Equations (II.2.4) and (23).
From Equation (II.2.4), 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 by Figure 4B, 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. That 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 4B 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. That 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 . As expected, the direct current between LHCII and the RC (cyan curve in Figure 4B) 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.
We see in Figure 5 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 Bennett et al. [2013]. 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.
Remarkably, in Figure 5 we see that there are a number of relatively strong 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 to 27 to 28, and back to 32. We also see in Figure 5 that in the reaction center overall the dominant net flow to the RP1 state is via pigment 30 and 32. Even though pigments 27 and 28 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. Pigments 27 and 28 receive a small amount of population from CP43 pigments, but the backflow from RP1 to these two pigments is larger and therefore seems to outweigh their contribution to the transport to RP1. 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 and 32, 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 the rest of the 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 and 32 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).
iii.3 Amount of coherence
According to the population current description of Equation (II.2.4), 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 has 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. Baumgratz et al. [2013], Kassal et al. [2013] for quantification and discussion of coherence). However, as described in detail in Ref. Roden and Whaley [2015], 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 transferdiminishing 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 timedependent electronic density matrix and therefore allows direct quantification of the coherences, the offdiagonal elements of the density matrix in the respective basis of interest. We first quantify the total amount of coherence by the sum
(24) 
over the absolute values of all offdiagonal elements of the density matrix , which provides an intuitive measure for the overall amount of coherence. This “norm of coherence” has been widely used in previous studies and is reviewed in Ref. Baumgratz et al. [2013].
Aside from the coherence in the site and exciton bases, we are also interested in the coherence in the domainexcition basis of Ref. Bennett et al. [2013]. This is of interest since in the MRGF calculation of Ref. Bennett et al. [2013], a classical rate equation is used to describe longrange intercomplex transfer that explicitly contains only the populations in the domainexciton 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 4A (dashed lines), nevertheless good agreement with the ZOFE population dynamics is achieved. This is consistent with Ref. Jesenko and Žnidarič [2013], 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 4C, we show the total amount of coherence measured by Equaiton (24) for the site basis, the exciton basis, and the domainexcition basis. The domainexciton coherence shows the coherence not (explicitly) taken into account in the MRGF calculations of Ref. Bennett et al. [2013].
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 domainexciton 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 domainexciton basis, the amount of coherence rises initially, because the dipoledipole 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 timeintegrated 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 timeintegrate 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 above, the majority of the excitation has been irreversibly trapped at RP2). In Figure 4D, the resulting sitebasis elements of the timeintegrated density matrix are shown as a color matrix. The timeintegrated pigment populations are the diagonal elements. One can see that there is an accumulation – a “dammingup” – of population on the energetically lowest sites of each protein. Such a dammingup 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 4D shows the (absolute values of) the imaginary parts of the coherences between the sites integrated over time. That is, it shows . As seen in Equation (II.2.4), 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 decribed by relaxation.
The lower triangle in Figure 4D shows the real parts of the coherence, that is, . We see that between many sites, the real part of the coherence is larger than the imaginary part. As discussed in Ref. Roden and Whaley [2015], the real part of the coherence is related to localization effects in the energy transfer due to the energy gaps between the pigments. 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 in Figure 4D).
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 4B. In Figure 4D, we see that there are even nonzero 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 in Figure 4B.
Iv Robustness of energy transport and function in PSII
In the previous section we have described excitation energy transfer when initial excitation occurs in a pigmentprotein 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 electronphonon coupling.
iv.1 Delocalized initial excitation
There is an ongoing discussion about the initial conditions of light absorption in lightharvesting systems, and how they affect the subsequent transport dynamics Kassal et al. [2013], Jesenko and Žnidarič [2013], Han et al. [2013]. 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 Kassal et al. [2013], Jesenko and Žnidarič [2013], Han et al. [2013].
Motivated by this, we now investigate the influence of the initial state on the longrange 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 – 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 ( is the number of sites), except for the population of the radical pair states and the ground state, which are zero, and with all offdiagonal elements (coherences) set to zero.
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 nonradiative decay to the ground state cause the initial state to evolve in a nonunitary 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 6A and are compared to the dynamics resulting from the localized initial state (Figure 4A). 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 backtransfer from CP43 to LHCII as can be seen in Figure 6B. 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 as seen before for the localized initial state, namely away from the reaction center.
Figure 7 shows the timeintegrated population currents between the individual pigments for the delocalized initial state, and should be compared to Figure 5 for the localized initial state.
The 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 shortrange (intracomplex) and longrange (intercomplex) 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 longrange 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).
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 6A. (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 timescales between the excitation energy transport and the loss mechanisms. Excitation is lost in PSII on a timescale of ns through a combination of fluorescence and nonradiative decay. Excitation transport, however, occurs on a ps timescale. If we approximate the excitation transport to the radical pair states by a single rate constant, then the overall behavior is described by a tworate (threestate) model – one rate for excitation to drive charge separation and a second rate of excitation undergoing nonradiative or radiative decay. Within this simplified model, the fraction of excitation that performs charge separation (as opposed to nonradiative or radiative decay) as time goes to infinity is given by
(25) 
– the population reaching the RP2 state in the longtime limit – where is the rate of population transfer to the RP states and is the rate of loss through nonradiative and radiative decay. We will use the phrase longtime efficiency to refer to this quantity and differentiate it from our 1 ns measure of efficiency, i.e. , used in the remainder of the paper. In Figure 6A we see that at 100 ps the total excitation left on the pigments amounts to a population of . Inserting this value in , we estimate the rate of excitation transport from the LHCII initial state (dashed lines) to the RP states as , which, using Eq. (25), gives a longtime 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 rate of excitation loss. For example, slowing the rate of transport to the RP states by a factor of 2 only decreases the longtime efficiency to 93%. As a result, small modulations in the rate of transport to the RP states have only a minimal influence on the longtime efficiency within our model, and thus will not substantially shift the functional behavior of PSII.
We now analyze the dependence of the amount of coherence on the initial state. In Figure 6C, we see that for the delocalized initial state there is a significantly smaller amount of excitonic coherence involved in the transfer dynamics compared to the amount for the localized state in Figure 4C (note the different range of the y axis of Figure 6C compared to Figure 4C). This is mainly due to the fact that there is no initial excitonic coherence present in the case of the delocalized initial state. In Figure 6D, the pattern of the timeintegrated coherence between the individual sites appears overall similar to that for the localized initial state in Figure 4D. However, looking closer, there are differences, which are reflected in the large difference between the currents resulting from the delocalized and localized initial states. We note also, that since in Figure 6C and D absolute values are summed up, changes of the signs of the individual terms that would correspond to opposite directions of the corresponding population currents are not visible. For example currents that lead to a flow of excitation back to LHCII in Figure 6B, which have negative signs, cannot be identified in Figure 6D, because the absolute value is taken before integrating the coherences. But we see in Figure 6D that just as the current between LHCII and CP43 is smaller for the delocalized initial state, so the imaginary parts of the coherences between LHCII and CP43 pigments are also weaker than before in Figure 4D.
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 Eisfeld [2011]. 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 longrange energy transport with respect to the initial conditions.
iv.2 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 with respect to changes in this part of the model. 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. Rebentrost et al. [2009c], Kolli et al. [2012], Chin et al. [2012]).
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 by a global factor that is the same for all of them, i.e.,
(26) 
where the constant factor is the same for all pigments 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 8 shows how the energy transfer changes when we decrease the coupling to the vibrations by factors of ten, i.e., .
For each of these cases, Figure 8A 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. ), 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. ), it goes down to less than 20% population on RP2 after 1 ns. Decreasing the coupling further (), causes most of the population to be trapped in LHCII and only a very small amount of population () reaches RP2. We can see in Figure 8A 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 nonradiative 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 and (the second and fourth bars from the top in Figure 8A), we take a closer look at the energy transport dynamics in Figures 8B, C (for ) and Figures 8D, E (for ). Figure 8B shows how for 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, observed in the bar diagram Figure 8A. 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 backandforth 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 8C, 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 4B 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 are shown in Figures 8D, 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 a quite long time due to (unitary) localization inside the complex because of nonuniform 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 9:
between the excited electronic states of the pigments there are energy gaps that, when they are large compared to the dipoledipole interaction strength between the respective pigments, cause the excitation to be trapped on one or a few pigments, because of the offresonance 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 relaxation down the energy manifold that is needed for the energy funnel, directing the transfer to the RP states, to work.
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 tranfer 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 8D. However, the population transferred to RP2 is still by more than a factor of two lower 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 , and when , , changing the initial state has little influence on the efficiency of the energy transfer. Only when the vibrational coupling is decreased beyond this, –, does the initial state have a significant impact on the overall efficiency. This shows again a high robustness of the energy transfer to variations of the initial condition.
Next, let us consider how decreasing the coupling to the vibrations affects the coherence involved in the transport dynamics. Figure 10 shows the magnitude of coherence for the case of and should be compared to Figures 4C, D, for the original vibrational coupling. For the weaker coupling to the vibrations (Figure 10), there is more coherence between the sites now, because the dephasing through the vibrations is weaker. The excitonic coherence, by contrast, decreases (see panel A compared to Figure 4C) which is expected because part of it is created by the coupling to the vibrations. Remarkably, these changes in the coherence do not affect the overall transport efficiency much, as we have seen in Figure 8.
iv.3 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 tranfer at stronger coupling as well. Stronger coupling, however, leads to problems in the ZOFE calculation, because of the limited range of nonMarkovianity for which the ZOFE approximation is valid Strunz and Yu [2004], Yu et al. [1999], Roden et al. [2011]. 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 electronvibration coupling. Based on the reasoning of Ref. Jesenko and Žnidarič [2013], 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 11 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 11, we see that as the coupling is increased, the transport efficiency becomes lower and lower, until for a fiftyfold 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 Rebentrost et al. [2009c], 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 (II.2.4), 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. There, the magnitude of the coupling to the vibrations is such that it gives the optimal balance between vibrationally enhanced transport (as described in Figure 9) 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 lie 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 weaker coupling to the vibrations.
V Summary and conclusions
In this work, we simulated the longrange intercomplex energy transfer in the PSII supercomplex – from LHCII via the core complex CP43 to the reaction center – by means of a full quantum calculation, using a nonMarkovian ZOFE quantum master equation description. We found that nearly all energy from the initial electronic excitation is transported to 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% is lost through radiative and nonradiative decay within this time.
These findings for the overall energy transfer dynamics are in good agreement with the results of Ref. Bennett et al. [2013], where a modifiedRedfieldgeneralizedFörster (MRGF) rateequation method was applied to calculate the longrange energy transfer in the entire PSII supercomplex. In the MRGF method, strongly coupled pigments are assigned to domains treated with a modifiedRedfield description, whereas excitation transfer between the domains, where the couplings are weaker, is treated with a generalizedFörster description. In the present work, to compare our ZOFE quantummasterequation results to the MRGF rateequation description of Ref. Bennett et al. [2013], we used the MRGF method to calculate the energy transfer in the truncated supercomplex LHCIICP43RC of PSII that we considered. The resulting energy transfer dynamics are in very good quantitative agreement with that calculated with the ZOFE master equation, showing that the simpler MRGF rate equation provides an adequate description of the quantum energy transfer dynamics and energy transfer efficiency of PSII. The observation that the MRGF results agree well with the ZOFE results shows that the subdivision into domains, which is related to the notion of supertransfer, i.e., coherent delocalization inside domains leads to enhanced transfer rates between the domains Kassal et al. [2013], is a reasonable approach to describe energy transfer in large lightharvesting complexes, where the coupling between the pigments can be very inhomogeneous. We have further shown that excitation energy transfer is significantly slower within a pure Förster calculation, which further supports the importance of supertransfer for rapid excitation transport within PSII. The close agreement between ZOFE and MRGF is also consistent with the findings of Ref. Jesenko and Žnidarič [2013], where a full quantum description of electronic excitation transfer is compared to that of a classical rate equation, and it is shown that the rate equation description is capable of capturing the overall features of the energy transfer when properly derived from the full quantum description.
In addition to the commonly used description in terms of timedependent populations of the excited electronic states of the pigments, we also analyzed the energy transfer dynamics in terms of excitation population currents between the pigments. This analysis revealed the pathways and directions of the energy transfer between the individual pigments and the net fow 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 intercomplex 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. Bennett et al. [2013]. The energy transfer between different proteins was found to rely on several different pathways. Interestingly, instead of being uniformly directed towards the reaction center, often excitation energy is transported in loops between pigments, i.e., excitation is transferred from one pigment to a neighboring pigment and from there to the next, etc., in such a way that it goes around in a closed circuit of pigments and in the end comes back to the initial pigment. Such loops were found in all three subcomplexes that we considered, LHCII, CP43, and the reaction center. In the reaction center, the efficient excitation energy transport to the charge separated states appeared to rely mainly on only two pigments (30 and 32), 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 and 32. 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 designfunction relationships in lightharvesting systems.
By simulating the energy transfer with the ZOFE master equation, we were able to quantify the amount of electronic coherence that is involved in the energy transfer dynamics. Motivated by the description of the energy transfer in terms of population currents and the different roles of imaginary and real part of the coherence presented in a companion paper Roden and Whaley [2015], we quantified the real and imaginary components of the coherence between the pigments separately. We found that the real components, which are related to localization of excitation due to energy gaps between the pigments Roden and Whaley [2015], are overall larger than the imaginary components, which are proportional to the population currents between the pigments and drive the excitation energy transfer in the PSII supercomplex.
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 shortrange (intracomplex) and longrange (intercomplex) energy transfer. This shows that the energy transfer in PSII is remarkably robust with respect to the initial conditions.
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 electronvibration coupling was decreased by a factor of 1,000, 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 occuring 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.
However, in contrast to the FennaMatthewsOlson (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 Rebentrost et al. [2009b], for PSII we found that the efficiency can become very low – almost zero – for extremely small coupling to the vibrations. This shows that contrary to FMO, for PSII the coupling to the vibrations is crucial for achieving efficient transport.
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 fiftyfold 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 Roden and Whaley [2015]. 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, in the light of the discussion about the possible (evolutionary) optimization of natural lightharvesting systems like PSII.
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 an extended range of the electronvibration coupling, in Section IV.3 we used a pure Lindblad simulation for the excitation energy transfer, applying the Markov approximation to the electronvibration coupling. Here, we briefly describe the calculation and how the corresponding parameters are obtained.
For the simulation, we employ a Lindblad master equation (19), where we now include Lindblad terms describing the coupling of the electronic excitation to the vibrations of the pigments and protein environment – in addition to the Lindblad terms for transfer to the RP1 and RP2 states and radiative and nonradiative decay to the ground state (relaxation terms) that were also included in the ZOFE simulations. The new Lindblad terms describing the electronvibrational coupling contain Lindblad operators and coupling parameters for each pigment , so that – instead of Eq. (20) for the ZOFE simulation – we now have the pure Lindblad master equation
(27) 
where the terms in the second line again describe transfer to the RP1 and RP2 states and radiative and nonradiative decay. In the Lindblad term in the first line, the electronic excitation of each pigment is coupled to vibrations with respective coupling strength . To calculate the excitation energy transfer from the full Lindblad description (27), we need to assume values for the coupling parameters . In the following, we show how we obtain these parameters through a simple approximation from the parameters that we used to describe the nonMarkovian electronvibration coupling within the ZOFE description. To this end, let us again consider the special form of the environment correlation function (16) used in the ZOFE description
(28) 
The pure Lindblad description (27) is obtained in the Markov limit, where
(29) 
i.e. decays fast compared to the timescales of the system dynamics. To express the parameters through the parameters of Eq. (28) in this limit, we consider the Dirac series for , which lets us write
(30) 
Inserting Eq. (30) in Eq. (28) we have
(31) 
Inserting (31) in the Eq. (14) for the ZOFE auxiliary operator, we get
(32) 
since (see Eq. (15)). Since we know that in the Markov limit , the auxiliary operator of the ZOFE master equation becomes the constant operator (see Eqs. (14) and (15)) Strunz and Yu [2004], from comparison with Eq. (32) it follows that we can approximate the needed coupling parameters for the full Lindblad descriptions as
(33) 
i.e. in the case where the are large compared to the system energies.
We note that in order for the to stay finite in the limit, not only the , but also at least one of the parameters has to go to infinity.
Inserting the parameters for the electronvibration coupling of the ZOFE simulation (according to Section II.2.2) into Eq. (33), we find the following values for the parameters that we used for our full Lindblad simulation in Section IV.3:
for the Chl A pigments in LHCII,
for the Chl B pigments in LHCII,
for the pigments in CP43,
and for the pigments in the reaction center.
We can see from the fact that these values 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 1,300 cm) 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 nonMarkovian dynamics.
Appendix B 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 12A that up to a certain time, the populations in the proteins with and without coupling are equal. 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 nonradiative 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 12C, where the coupling to the RP states is switched off, with that in Figure 4C, 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 CauchySchwarz inequality
(34) 
that holds for the coherences and the populations , of any density matrix . 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 CauchySchwarz 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.
Acknowledgements.
This material was supported by DARPA under Award No. N660010912026. J. R. thanks Brendan Abolins, Daniel Freeman, Siva Darbha, Jon Aytac, Ty Volkoff, Felix Motzoi, Loren Greenman, Donghyun Lee, and Aleksey Kocherzhenko for helpful discussions.References
 Knowles [1980] Jeremy R Knowles. Enzymecatalyzed phosphoryl transfer reactions. Annual review of biochemistry, 49(1):877–919, 1980.
 Caffarri et al. [2011] Stefano Caffarri, Koen Broess, Roberta Croce, and Herbert van Amerongen. Excitation energy transfer and trapping in higher plant photosystem ii complexes with different antenna sizes. Biophysical journal, 100(9):2094–2103, 2011.
 Ruban and Murchie [2012] Alexander V Ruban and Erik H Murchie. Assessing the photoprotective effectiveness of nonphotochemical chlorophyll fluorescence quenching: a new approach. Biochimica et Biophysica Acta (BBA)Bioenergetics, 1817(7):977–982, 2012.
 Belgio et al. [2012] Erica Belgio, Matthew P Johnson, Snježana Jurić, and Alexander V Ruban. Higher plant photosystem ii lightharvesting antenna, not the reaction center, determines the excitedstate lifetime—both the maximum and the nonphotochemically quenched. Biophysical journal, 102(12):2761–2771, 2012.
 Marin et al. [2012] Alessandro Marin, Ivo HM van Stokkum, Vladimir I Novoderezhkin, and Rienk van Grondelle. Excitationinduced polarization decay in the plant lightharvesting complex lhcii. Journal of Photochemistry and Photobiology A: Chemistry, 234:91–99, 2012.
 Krüger et al. [2012] Tjaart PJ Krüger, Cristian Ilioaia, Matthew P Johnson, Alexander V Ruban, Emmanouil Papagiannakis, Peter Horton, and Rienk van Grondelle. Controlled disorder in plant lightharvesting complex ii explains its photoprotective role. Biophysical journal, 102(11):2669–2676, 2012.
 Krüger et al. [2010] Tjaart PJ Krüger, Vladimir I Novoderezhkin, Cristian Ilioaia, and Rienk Van Grondelle. Fluorescence spectral dynamics of single lhcii trimers. Biophysical journal, 98(12):3093–3101, 2010.
 SchlauCohen et al. [2014] Gabriela S SchlauCohen, Samuel Bockenhauer, Quan Wang, and WE Moerner. Singlemolecule spectroscopy of photosynthetic proteins in solution: exploration of structure–function relationships. Chemical Science, 2014.
 Novoderezhkin et al. [2011a] Vladimir I Novoderezhkin, Elisabet Romero, Jan P Dekker, and Rienk van Grondelle. Multiple chargeseparation pathways in photosystem ii: Modeling of transient absorption kinetics. ChemPhysChem, 12(3):681–688, 2011a.
 Novoderezhkin et al. [2011b] Vladimir Novoderezhkin, Alessandro Marin, and Rienk van Grondelle. Intraand intermonomeric transfers in the light harvesting lhcii complex: the redfield–förster picture. Physical Chemistry Chemical Physics, 13(38):17093–17103, 2011b.
 Fuller et al. [2014] Franklin D Fuller, Jie Pan, Andrius Gelzinis, Vytautas Butkus, S Seckin Senlik, Daniel E Wilcox, Charles F Yocum, Leonas Valkunas, Darius Abramavicius, and Jennifer P Ogilvie. Vibronic coherence in oxygenic photosynthesis. Nature chemistry, 6(8):706–711, 2014.
 Wells et al. [2014] Kym L Wells, Petar H Lambrev, Zhengyang Zhang, Gyözö Garab, and HoweSiang Tan. Pathways of energy transfer in lhcii revealed by roomtemperature 2d electronic spectroscopy. Physical Chemistry Chemical Physics, 16(23):11640–11646, 2014.
 SchlauCohen et al. [2012] Gabriela S SchlauCohen, Akihito Ishizaki, Tessa R Calhoun, Naomi S Ginsberg, Matteo Ballottari, Roberto Bassi, and Graham R Fleming. Elucidation of the timescales and origins of quantum electronic coherence in lhcii. Nature chemistry, 4(5):389–395, 2012.
 Broess et al. [2006] Koen Broess, Gediminas Trinkunas, Chantal D van der Weijde Wit, Jan P Dekker, Arie van Hoek, and Herbert van Amerongen. Excitation energy transfer and charge separation in photosystem ii membranes revisited. Biophysical journal, 91(10):3776–3786, 2006.
 Raszewski and Renger [2008] Grzegorz Raszewski and Thomas Renger. Light harvesting in photosystem ii core complexes is limited by the transfer to the trap: can the core complex turn into a photoprotective mode? Journal of the American Chemical Society, 130(13):4431–4446, 2008.
 Bennett et al. [2013] Doran IG Bennett, Kapil Amarnath, and Graham R Fleming. A structurebased model of energy transfer reveals the principles of light harvesting in photosystem ii supercomplexes. Journal of the American Chemical Society, 135(24):9164–9173, 2013.
 Kreisbeck et al. [2014] Christoph Kreisbeck, Tobias Kramer, and Alán AspuruGuzik. Scalable highperformance algorithm for the simulation of exciton dynamics. application to the lightharvesting complex ii in the presence of resonant vibrational modes. Journal of Chemical Theory and Computation, 10(9):4045–4054, 2014.
 Shibata et al. [2013] Yutaka Shibata, Shunsuke Nishi, Keisuke Kawakami, JianRen Shen, and Thomas Renger. Photosystem ii does not possess a simple excitation energy funnel: timeresolved fluorescence spectroscopy meets theory. Journal of the American Chemical Society, 135(18):6903–6914, 2013.
 Umena et al. [2011] Yasufumi Umena, Keisuke Kawakami, JianRen Shen, and Nobuo Kamiya. Crystal structure of oxygenevolving photosystem ii at a resolution of 1.9 å. Nature, 473(7345):55–60, 2011.
 Liu et al. [2004] Zhenfeng Liu, Hanchi Yan, Kebin Wang, Tingyun Kuang, Jiping Zhang, Lulu Gui, Xiaomin An, and Wenrui Chang. Crystal structure of spinach major lightharvesting complex at 2.72 å resolution. Nature, 428(6980):287–292, 2004.
 Pan et al. [2011] Xiaowei Pan, Mei Li, Tao Wan, Longfei Wang, Chenjun Jia, Zhiqiang Hou, Xuelin Zhao, Jiping Zhang, and Wenrui Chang. Structural insights into energy regulation of lightharvesting complex cp29 from spinach. Nature structural & molecular biology, 18(3):309–315, 2011.
 Müh et al. [2014] Frank Müh, Dominik Lindorfer, Marcel Schmidt am Busch, and Thomas Renger. Towards a structurebased exciton hamiltonian for the cp29 antenna of photosystem ii. Physical Chemistry Chemical Physics, 16(24):11848–11863, 2014.
 van Oort et al. [2010] Bart van Oort, Marieke Alberts, Silvia de Bianchi, Luca Dall’Osto, Roberto Bassi, Gediminas Trinkunas, Roberta Croce, and Herbert van Amerongen. Effect of antennadepletion in photosystem ii on excitation energy transfer in arabidopsis thaliana. Biophysical journal, 98(5):922–931, 2010.
 Caffarri et al. [2009] Stefano Caffarri, Roman Kouřil, Sami Kereïche, Egbert J Boekema, and Roberta Croce. Functional architecture of higher plant photosystem ii supercomplexes. The EMBO journal, 28(19):3052–3063, 2009.
 Rebentrost et al. [2009a] Patrick Rebentrost, Rupak Chakraborty, and Alán AspuruGuzik. Nonmarkovian quantum jumps in excitonic energy transfer. The Journal of chemical physics, 131(18):184102, 2009a.
 Ishizaki and Fleming [2009] Akihito Ishizaki and Graham R Fleming. Theoretical examination of quantum coherence in a photosynthetic system at physiological temperature. Proceedings of the National Academy of Sciences, 106(41):17255–17260, 2009.
 Chin et al. [2010] Alex W Chin, Animesh Datta, Filippo Caruso, Susana F Huelga, and Martin B Plenio. Noiseassisted energy transfer in quantum networks and lightharvesting complexes. New Journal of Physics, 12(6):065002, 2010.
 Kreisbeck et al. [2011] Christoph Kreisbeck, Tobias Kramer, Mirta Rodríguez, and Birgit Hein. Highperformance solution of hierarchical equations of motion for studying energy transfer in lightharvesting complexes. Journal of Chemical Theory and Computation, 7(7):2166–2174, 2011.
 Ritschel et al. [2011a] Gerhard Ritschel, Jan Roden, Walter T Strunz, and Alexander Eisfeld. An efficient method to calculate excitation energy transfer in lightharvesting systems: application to the fenna–matthews–olson complex. New Journal of Physics, 13(11):113034, 2011a.
 Ritschel et al. [2011b] Gerhard Ritschel, Jan Roden, Walter T Strunz, Alán AspuruGuzik, and Alexander Eisfeld. Absence of quantum oscillations and dependence on site energies in electronic excitation transfer in the fenna–matthews–olson trimer. The Journal of Physical Chemistry Letters, 2(22):2912–2917, 2011b.
 Engel et al. [2007] Gregory S Engel, Tessa R Calhoun, Elizabeth L Read, TaeKyu Ahn, Tomáš Mančal, YuanChung Cheng, Robert E Blankenship, and Graham R Fleming. Evidence for wavelike energy transfer through quantum coherence in photosynthetic systems. Nature, 446(7137):782–786, 2007.
 Romero et al. [2014] Elisabet Romero, Ramunas Augulis, Vladimir I Novoderezhkin, Marco Ferretti, Jos Thieme, Donatas Zigmantas, and Rienk Van Grondelle. Quantum coherence in photosynthesis for efficient solarenergy conversion. Nature Physics, 10(9):676–682, 2014.
 OlayaCastro et al. [2008] Alexandra OlayaCastro, Chiu Fan Lee, Francesca Fassioli Olsen, and Neil F Johnson. Efficiency of energy transfer in a lightharvesting system under quantum coherence. Physical Review B, 78(8):085115, 2008.
 Scholes et al. [2011] Gregory D Scholes, Graham R Fleming, Alexandra OlayaCastro, and Rienk van Grondelle. Lessons from nature about solar light harvesting. Nature chemistry, 3(10):763–774, 2011.
 Hoyer et al. [2012] Stephan Hoyer, Akihito Ishizaki, and K Birgitta Whaley. Spatial propagation of excitonic coherence enables ratcheted energy transfer. Physical Review E, 86(4):041911, 2012.
 Jesenko and Žnidarič [2013] Simon Jesenko and Marko Žnidarič. Excitation energy transfer efficiency: Equivalence of transient and stationary setting and the absence of nonmarkovian effects. The Journal of chemical physics, 138(17):174103, 2013.
 Chenu et al. [2014] Aurélia Chenu, Pavel Malỳ, and Tomáš Mančal. Dynamic coherence in excitonic molecular complexes under various excitation conditions. Chemical Physics, 2014.
 Scholes et al. [2012] Gregory D Scholes, Tihana Mirkovic, Daniel B Turner, Francesca Fassioli, and Andreas Buchleitner. Solar light harvesting by energy transfer: from ecology to coherence. Energy & Environmental Science, 5(11):9374–9393, 2012.
 Kassal et al. [2013] Ivan Kassal, Joel YuenZhou, and Saleh RahimiKeshari. Does coherence enhance transport in photosynthesis? The Journal of Physical Chemistry Letters, 4(3):362–367, 2013.
 Rebentrost et al. [2009b] Patrick Rebentrost, Masoud Mohseni, and Alán AspuruGuzik. Role of quantum coherence and environmental fluctuations in chromophoric energy transport. The Journal of Physical Chemistry B, 113(29):9942–9947, 2009b.
 Hoyer et al. [2010] Stephan Hoyer, Mohan Sarovar, and K Birgitta Whaley. Limits of quantum speedup in photosynthetic light harvesting. New Journal of Physics, 12(6):065041, 2010.
 Shim et al. [2012] Sangwoo Shim, Patrick Rebentrost, Stéphanie Valleau, and Alán AspuruGuzik. Atomistic study of the longlived quantum coherences in the fennamatthewsolson complex. Biophysical journal, 102(3):649–660, 2012.
 Abramavicius and Mukamel [2010] Darius Abramavicius and Shaul Mukamel. Quantum oscillatory exciton migration in photosynthetic reaction centers. The Journal of chemical physics, 133(6):064510, 2010.
 Sarovar and Whaley [2013] Mohan Sarovar and K Birgitta Whaley. Design principles and fundamental tradeoffs in biomimetic light harvesting. New Journal of Physics, 15(1):013030, 2013.
 Chenu and Scholes [2015] Aurélia Chenu and Gregory D Scholes. Coherence in energy transfer and photosynthesis. Annual review of physical chemistry, pages 69–96, 2015.
 Roden and Whaley [2015] Jan JJ Roden and K Birgitta Whaley. A probability current analysis of energy transport in open quantum systems. arXiv preprint arXiv:1501.06090, 2015.
 Strunz and Yu [2004] Walter T Strunz and Ting Yu. Convolutionless nonmarkovian master equations and quantum trajectories: Brownian motion. Physical Review A, 69(5):052115, 2004.
 Han et al. [2013] Alex C Han, Moshe Shapiro, and Paul Brumer. Nature of quantum states created by one photon absorption: Pulsed coherent vs pulsed incoherent light. The Journal of Physical Chemistry A, 117(34):8199–8204, 2013.
 Dorfman et al. [2013] Konstantin E Dorfman, Dmitri V Voronine, Shaul Mukamel, and Marlan O Scully. Photosynthetic reaction center as a quantum heat engine. Proceedings of the National Academy of Sciences, 110(8):2746–2751, 2013.
 Mukamel [1982] S Mukamel. Principles of nonlinear optical spectroscopy, 1995. Oxford University Press, New York, Mukamel S., J. Chem. Phys, 77:173, 1982.
 Yu et al. [1999] Ting Yu, Lajos Diósi, Nicolas Gisin, and Walter T Strunz. Nonmarkovian quantumstate diffusion: Perturbation approach. Physical Review A, 60(1):91, 1999.
 Caruso et al. [2009] Filippo Caruso, Alex W Chin, Animesh Datta, Susana F Huelga, and Martin B Plenio. Highly efficient energy excitation transfer in lightharvesting complexes: The fundamental role of noiseassisted transport. The Journal of Chemical Physics, 131(10):105106, 2009.
 Rebentrost et al. [2009c] Patrick Rebentrost, Masoud Mohseni, Ivan Kassal, Seth Lloyd, and Alán AspuruGuzik. Environmentassisted quantum transport. New Journal of Physics, 11(3):033003, 2009c.
 Blankenship [2014] Robert E Blankenship. Molecular mechanisms of photosynthesis. John Wiley & Sons, 2014.
 Förster [1948] T Förster. Energy transfer and fluorescence between molecules. Ann Phys, 437:55–75, 1948.
 May and Kühn [2008] Volkhard May and Oliver Kühn. Charge and energy transfer dynamics in molecular systems. John Wiley & Sons, 2008.
 Baumgratz et al. [2013] T Baumgratz, M Cramer, and MB Plenio. Quantifying coherence. arXiv preprint arXiv:1311.0275, 2013.
 Eisfeld [2011] Alexander Eisfeld. Phase directed excitonic transport and its limitations due to environmental influence. Chemical Physics, 379(1):33–38, 2011.
 Kolli et al. [2012] Avinash Kolli, Edward J O’Reilly, Gregory D Scholes, and Alexandra OlayaCastro. The fundamental role of quantized vibrations in coherent light harvesting by cryptophyte algae. The Journal of chemical physics, 137(17):174109, 2012.
 Chin et al. [2012] AW Chin, SF Huelga, and MB Plenio. Coherence and decoherence in biological systems: principles of noiseassisted transport and the origin of longlived coherences. Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 370(1972):3638–3657, 2012.
 Roden et al. [2012] Jan Roden, Walter T Strunz, K Birgitta Whaley, and Alexander Eisfeld. Accounting for intramolecular vibrational modes in open quantum system description of molecular systems. The Journal of chemical physics, 137(20):204110, 2012.
 Roden et al. [2011] Jan Roden, Walter T Strunz, and Alexander Eisfeld. Nonmarkovian quantum state diffusion for absorption spectra of molecular aggregates. The Journal of chemical physics, 134(3):034902, 2011.