Tight-binding model of the photosystem II reaction center: application to two-dimensional electronic spectroscopy

We propose an optimized tight-binding electron–hole model of the photosystem II (PSII) reaction center (RC). Our model incorporates two charge separation pathways and spatial correlations of both static disorder and fast fluctuations of energy levels. It captures the main experimental features observed in time-resolved two-dimensional (2D) optical spectra at 77 K: peak pattern, lineshapes and time traces. Analysis of 2D spectra kinetics reveals that specific regions of the 2D spectra of the PSII RC are sensitive to the charge transfer states. We find that the energy disorder of two peripheral chlorophylls is four times larger than the other RC pigments.


Introduction
Photosystem II (PSII) is a unique biological system that is capable of oxidizing water to molecular oxygen [1]. PSII is a complicated pigment-protein structure, containing D1, D2, CP43 and CP47 subunits. Peripheral complexes, such as light harvesting complex II (LHCII) absorb sunlight and transfer energy to the D1-D2 proteins, where the initial charge separation occurs.
The structure of the PSII reaction center (RC) has been resolved by x-ray crystallography [2][3][4][5]. The D1 and D2 branches contain eight pigments responsible for the long-wavelength absorption: two primary chlorophylls P D1 and P D2 (from analogy with the RC of photosynthetic bacteria often called the special pair), two accessory chlorophylls Chl D1 and Chl D2 , two pheophytins Pheo D1 and Pheo D2 and two peripheral chlorophylls Chlz D1 and Chlz D2 . The simplified picture of the arrangement of pigment molecules in the PSII RC is presented in figure 1. Charge separation proceeds only through the D1 branch [6].
The structure of the PSII RC plays a very important role in determining its charge transfer (CT) properties. Compared to the bacterial RC, the special pair pigments are more separated and slightly tilted with respect to each other [7]. Due to these differences in the structural arrangement, the lowest energy excited state is no longer localized on the special pair as is the case in the bacterial RC and corresponds to a coherent superposition of several pigment molecules from the PSII RC, which is known as the multimer state [8]. Moreover, the contribution from all pigments to this multimer state is not uniform; it strongly depends on the realization of static disorder as determined by the protein conformation. These spectral properties also help to increase the redox potential, which allows PSII to perform water splitting.
Numerous optical spectroscopy measurements have been performed on the PSII RC [9][10][11]. Important information about the optical properties of different RC pigments was obtained from the spectra of mutant PSII complexes [12,13] and different preparations of RC complexes [14,15]. Photon echo experiments suggested that Chl D1 is the primary electron donor and that Pheo D1 is the primary electron acceptor [16]. The latter was supported by transient absorption (TA) measurements [17]. Stark spectroscopy experiments showed that the Stark spectra [33] suggested that the CT state associated with the special pair is strongly mixed with pigment states. It was concluded that the Chl D1 pathway is responsible for the fast ∼700 fs component of the kinetics and that the special pair pathway gives the slower ∼3 ps component [34]. Modeling with only one CT pathway was not able to reproduce the experimental kinetics. Regarding electron separation, these simulations [33,34] extended the Frenkel exciton model to include CT states on the same level as additional excitonic modes. However, the fermionic properties of electrons and holes are not accounted for in this approach, so the double excitations are not defined and the approach in principle cannot account for excited state absorption. Consequently, the recent simulations of the 2D spectra of the PSII RC [37], based on the Novoderezhkin model, were inherently limited-they did not treat double excitations.
The solution to this problem is to start from the electron-hole model [38]. In this approach, which was developed for semiconductors [39,40], the electrons and holes are core particles that comprise both molecular excitations and CT states. Double, triple, etc excitations are generated by the fermionic creation and annihilation operators. Excited state absorption relevant for pump-probe or 2D spectroscopy is thus properly accounted for.
In this paper, we employ the two-band tight-binding (TB) approach [38] to treat both molecular excitations and CT states on an equal footing and optimize the model parameters with respect to the 2D spectroscopy signals. The model additionally applies some restrictions on the parameters. For instance, in the Frenkel exciton model the energy fluctuations are usually assumed as independent for different molecules. In the TB model, additional spatial correlations between these fluctuations involving CT states become important. We also replace the nuclear fluctuation spectral density, which had a complicated fine structure in the simulations based on the Novoderezhkin model [32-34, 37, 38], with a simple broad function that accounts for high-frequency contributions in a single term. Simulations based on our model reproduce both absorption and 2D spectra of the PSII RC at 77 K. Also, we obtain a reasonable match between experimental and simulated peak kinetics.

Tight-binding Hamiltonian
In the TB electron-hole model [38], each pigment is represented by two electronic orbitals: the highest occupied molecular orbital (HOMO) and the lowest unoccupied molecular orbital (LUMO). The model is constructed starting from operatorsê † m (ê m ) that create (annihilate) an electron in the LUMO orbital of site m and onĥ † m (ĥ m ) that create (annihilate) a hole in the HOMO orbital. These operators satisfy Fermi commutation relations The Hamiltonian is given by [38] Here t e mn (t h mn ) is the electron (hole) hopping rate between LUMO (HOMO) orbitals of different pigments. W d mn is the dipole-dipole-type resonance interaction between excitons on sites m and n. V e mn = V (r m − r n ) is the electron-electron Coulomb repulsion between sites m and n, V h mn = V (r m − r n ) is the hole-hole Coulomb repulsion between sites m and n, and V eh mn = V (r m − r n ) is the electron-hole Coulomb attraction between sites m and n. Note that V e mn = V h mn = V eh mn . In this model, the last term represents K couplings [41], which are responsible for the shifting of double-excited state energies. These shifts were not taken into account in [38].
The single-excitation manifold consists of single electron-hole pair states, denoted The Hamiltonian matrix elements for these states are In the double-excitation manifold, there are two electrons and two holes in the state n |g with k > l and m > n. Three kinds of doubly excited states are possible: FE-FE states m * n * , FE-CT states m * n + k − and CT-CT states m + n − k + l − . Various types of single-and double-excitation Hamiltonian matrix elements are given in appendix A.
The single-excitation eigenstates are defined by the transformation while for the double-excitation eigenstates we have The transformation matrices ψ e,kl and f,klmn are calculated by diagonalizing the single-and double-excitation blocks of the system Hamiltonian. We assume that all elements of these matrices are real. Additional vibrational degrees of freedom constitute the thermal bath, which is modeled by an infinite set of harmonic oscillators, that are linearly coupled to electron and hole orbitals, Here ω j ,p j andx j are the frequency, dimensionless momentum and corresponding dimensionless coordinate of the jth bath mode, respectively. d e n j (d h n j ) is the dimensionless displacement of the equilibrium configuration of the jth bath mode between the ground and excited states of the electron (hole) at the nth site. Note that we seth = 1 throughout the paper. We define bath operatorsF e m = − j ω j d e m jx j andF h m = − j ω j d h m jx j . Then the system-bath interaction Hamiltonian is expressed aŝ In this picture the bath induces energy fluctuations of HOMO and LUMO levels. All bath-induced effects to system dynamics are determined by the following bath correlation functions [42]: HereF denotes operatorF x k in the interaction representation and x = e or h, whilê is the bath equilibrium density operator, β = (k B T ) −1 is the inverse temperature in energy units and k B is the Boltzmann constant. The spectral densities are often introduced to characterize the bath [43]: Since the spectral densities are independent of temperature, they are the most convenient bath characteristics. In practice, the spectral densities are modeled as continuous functions of Figure 2. Scheme of the 2D spectroscopy experiment.
frequency. The corresponding correlation functions and the spectral densities are related by the fluctuation dissipation theorem [43]: where x = e, h or eh.
Given this system, the external classical optical field creates or annihilates molecular excitations:Ĥ whereμ is the dipole operator and µ m is the transition dipole vector for site m. In this formulation the electric field does not couple with CT states. The total Hamiltonian is then a sum of system, system-bath, bath and system-field terms.

Calculation of spectra
The linear absorption spectrum is simulated by summing up all possible transitions from the electronic ground state to all excited states. It is given by the absorptive part of the Fouriertransformed linear response function [43,44] The two-dimensional optical spectroscopy is a four-pulse four-wave-mixing experiment [43,45]. The sample is excited by three time-ordered laser pulses with wavevectors k 1 (first pulse), k 2 (second pulse) and k 3 (third pulse). These pulses are separated by intervals t 1 and t 2 , producing signals that are emitted along directions k S = ±k 1 ± k 2 ± k 3 in the time interval t 3 (see figure 2). The rephasing part is, however, detected in direction k I = −k 1 + k 2 + k 3 , while the non-rephasing part in direction k II = +k 1 − k 2 + k 3 . The 2D spectrum is obtained by performing a 2D Fourier transform t 1 → ω 1 and t 3 → ω 3 . The second time interval, t 2 , is left as an evolution parameter of the system.
In this work we model the so-called total 2D spectrum, which is the sum of rephasing W k I and non-rephasing W k II contributions. It is directly obtained in the pump-probe configuration. The total 2D signal is Using the response function formalism [43][44][45] and with the assumption of δ-shaped pulses, these terms are expressed as [46] Here S k I/II (t 1 , t 2 , t 3 ) is the contribution to the system response function corresponding to either the k I or the k II signal. In the following the electronic energy levels are included explicitly, while the nuclear part is divided into two parts-the diagonal fluctuations in the eigenstate basis are treated exactly via the second-order cumulant expansion, while the off-diagonal fluctuations are included perturbatively to second order. The electronic response functions are then modified by the eigenstate representation spectral lineshape functions, which encapsulate the effect of diagonal nuclear fluctuations [43]: Here α, β, γ and δ represent arbitrary single-or double-excitation eigenstates. C αβ,γ δ (ω) is the spectral density of the eigenstate basis. The relevant spectral densities are C e 1 e 2 ,e 3 e 4 (ω), which describes fluctuations in the single-excitation manifold, C e 1 e 2 , f 1 f 2 (ω), which describes correlations between fluctuations in single-and double-excitation manifolds, and C f 1 f 2 , f 3 f 4 (ω), which describes fluctuations in the double-excitation manifold. These spectral densities are given by a linear transformation of the real-space spectral densities defined by (15)- (17), as shown in appendix B. Population transport is significant for absorption due to lifetime broadening and is directly detected in the 2D spectrum. The transport is essentially described by the Pauli master equation approach. Its rates are given by the modified Redfield theory [47], which reduces to Förster and traditional Redfield theories in their respective ranges of validity [48]. The electron transfer theory of Marcus [49] is also reproduced by the modified Redfield theory [38]. It gives the following population transfer rate from state e to state e: Here the dot denotes the time derivative and is a type of reorganization energy of cross-correlation between the (Ĥ SB ) αβ and (Ĥ SB ) γ δ elements.
The absorption spectrum is calculated as [50] where e labels the single-excited eigenstate, µ eg = e|μ|g is the transition dipole moment between the ground state and the state e, • or denotes orientational averaging [44], ω eg = ε e − ε g is the energy difference of state e and the ground state g (usually ε g = 0) and ξ e = 1 2 e =e e R e e is the lifetime broadening. The third-order response functions for the 2D signals are calculated by summing up the relevant Liouville space pathways [43] as described in [44] (we used equations (141)-(153) with (166) and (179) and orientational averaging was performed using expressions from appendix L).

Model of the photosystem II (PSII) reaction center (RC)
The PSII RC consists of eight pigment states, so we can have up to 8 × 8 = 64 singly excited states (8 molecular excitations and 56 CT states). However, in our simulations we include only a few CT states-those that have been suggested to participate in the charge separation process. We thus include six CT states that describe two available electron transfer pathways [19,34]: We denote the first pathway as the Chl D1 pathway and the second as the special pair pathway. Now we consider the elements of the single-excitation Hamiltonian. First, we note that complete knowledge of all the parameters defined in the TB Hamiltonian (3) is not necessary for calculations of optical spectra: only their specific combinations resulting in state energies and couplings (see (A.1)-(A.5)) need to be specified. The energies of the single-excited states in our model were obtained as follows. The energies of the core pigments and P + D2 P − D1 state were taken from [34], where they were obtained from a fit of various low-temperature spectra. The energies of the peripheral chlorophylls were based on [33], while those of the remaining CT states were based on [34,38]. The energies of the pigments and of the CT states were fine tuned to better fit the absorption and the 2D spectra at zero t 2 delay.
Considering the couplings of the single-excited states, those between the Frenkel excitations were taken from [51], where they were calculated using the TrEsp method [52] from the 1.9 Å resolution structural data [5]. In order to obtain the couplings involving CT states, we have assumed that the electron and hole hopping rates are equal, that is t e mn = t h mn . Then, we have taken t e P D1 P D2 = 45 cm −1 based on J P * D1 ,P + D2 P − D1 from [34], where it was extracted from a fit of various low-temperature spectra, and we have also taken t e Chl D1 Pheo D1 = 70 cm −1 and t e P D1 Chl D1 = 40 cm −1 from couplings J Chl * given in [34], where they were estimated from the fit of 77 K TA kinetics. We also assume that t e P D2 Chl D1 = t e P D1 Pheo D1 = 0 due to the relatively larger distance between the pigments. Then all the couplings involving CT states were calculated using (A.4) and (A.5). The single-excitation Hamiltonian used in our model is presented in table 1.  Now we turn to the double-excitation block of the electronic Hamiltonian. It involves Frenkel-type couplings, electron and hole hopping terms, electron-electron and hole-hole repulsion and electron-hole attraction terms and K coupling terms. Frenkel-type couplings and electron and hole hopping terms define the off-diagonal elements of the double-excitation Hamiltonian block. These are directly related to the single-excitation Hamiltonian elements as shown in appendix A. The electron-electron repulsion, hole-hole repulsion, electron-hole attraction and K couplings determine the shifts of double-excitation state energies, i.e. they appear on the diagonal. The repulsion and attraction terms can be evaluated as was done in [38]. The values of K couplings need to be additionally specified. Since the K couplings depend on the product of static dipoles of the states [41], and CT states have much bigger static dipoles than FE states [33,34], the K couplings may be very different depending on the spatial configurations of electrons and holes in specific double-excited states. We can thus identify the FE-FE-type state where electron and hole pairs are two Frenkel excitations. We assume that the K couplings in such FE-FE states are relatively small and set K mm,nn = 0 [41]. The next type of double excitations are those composed of two CT states. Due to the large dipole value of a CT state, the K coupling strengths should be huge and they induce large shifts of the energies of CT-CT states, effectively decoupling them from the rest of the double-excited states. As the optical fields are resonant with the FE transitions, the CT-CT states are off-resonant and are ignored in our model. Finally, we can have the double-excited state composed of one CT state and one Frenkel exciton. We assume that the energies of all CT-FE-type states are shifted by the same value K , which can be determined from the fit of the simulated 2D spectra. K = 100 cm −1 was obtained from our simulations. Given these assumptions, all the double-excitation Hamiltonian elements were calculated using (A.7)-(A.15).
Dipole moments for pigments have been obtained from crystallographic data [5], assuming that the transition dipole is pointing from N B to N D atoms for both chlorophylls and pheophytins. The strengths of the transition dipoles are 4.4 and 3.4 Debye for chlorophylls and pheophytins, respectively [53]. The transition dipole moments used in the calculations are presented in table 2.
We next describe the model of the bath. The bath-induced fluctuations are separated into fast and slow modes. The slow modes are slower than the experiment timescale, so they can be Here the tilde denotes random variables, ε m * = ε m * , ε m + n − = ε m + n − andx e/h m =t e/h mm − t e/h mm is the random shift of electron or hole orbital energy from the mean. We next assume that electron and hole orbital energy shifts at the same site have the same standard deviations, and can be correlated, Here σ m is the standard deviation of orbital energies on site m and ν m ∈ [−1, 1] is the correlation coefficient between LUMO and HOMO orbital energies. Orbital energies on different sites are uncorrelated; thus Note that if shifts of electron and hole orbital energies are completely anticorrelated, the Frenkel excitation energy has no fluctuation. In that case an increase of the energy of the HOMO orbital would be accompanied by a decrease of the LUMO orbital energy (or vice versa), resulting in no fluctuation of the FE state. Standard deviations of single-excited state energies were optimized to fit both absorption and 2D spectra and are presented in table 2. We have taken 15 16 . This choice of ν values means that P + D2 P − D1 and Chl + D1 Pheo − D1 have, respectively, two and four times wider disorder than the core pigments, as was suggested in [34]. Since other pigments do not participate in the CT states present in our model, the rest of the σ m and ν m do not have to be specified explicitly. In our formulation, CT states have correlated disorder with other states involving the same pigments. These correlations are described in appendix C. Note that in our model, peripheral chlorophylls have significantly larger disorder than the rest of the pigments.
The same formulation holds for fast fluctuations, which are responsible for the homogeneous broadening of the spectra, however, their correlation coefficients should be different. We assume that the electron and hole orbital energy fluctuations at the same site are described by the same spectral densities and may be correlated: The correlation parameter is η m ∈ [−1, 1]. The fast fluctuations of electron and hole orbital energies at different sites are uncorrelated. We assume that all FE states have the same spectral density Then the spectral density of CT states is The values η P D1 = η P D2 = − 1 3 and η Chl D1 = η Pheo D1 = − 2 3 reproduce χ P + D2 P − D1 = 1.5 and χ Chl + D1 Pheo − D1 = 3.0 as in [34]. While these coupling values are not as large as in some other studies [54,55], they were inferred in earlier simulations in the framework of modified Redfield theory [33]. Application of a more advanced theory, capable of reproducing polaronic effects due to strong system-bath coupling [56], especially for the CT states, might require revision of these values. This specifies the spectral density of all the CT states of our model in terms of C (ω). The values of χ mn are given in table 2. Note that for FE states χ mm = 1. As is the case with static disorder, fast fluctuations of CT states are correlated with other states involving the same pigments-these correlations are given in appendix C.
The fast fluctuations are now described by a single spectral density C (ω). We assume it to consist of two terms: The first term describes coupling with low-frequency overdamped vibrations modeled by the Debye spectral density [57] with λ L = 40 cm −1 and γ L = 40 cm −1 . These values were obtained from simulations of absorption and 2D spectra, as the low-frequency part of the spectral density gives the homogeneous broadening. The second, damped, term represents coupling with highfrequency bath modes and was obtained by reducing the spectral density of Novoderezhkin et al [34]. Their function includes 48 vibrational modes, obtained from fluorescence line narrowing experiments [58], thus giving experiment-based detailed vibrational sidebands in the spectrum. Our parameters λ H = 500 cm −1 , γ H = 500 cm −1 and ω H = 1250 cm −1 lead to the same total reorganization energy for FE states λ = λ L + λ H = 540 cm −1 as in [34]. This spectral density is presented in figure 3(a) alongside the spectral density used in [34]. Additionally, for comparison we also plot the spectral density of Raszewski et al [31,53], which was used for the PSII RC, as well as the one by Jang and Silbey [59]. These two last models are similar in the range 50-1000 cm −1 but have different low-frequency behavior (see the inset in the figure). Additionally they both lack the high-frequency component of Novoderezhkin's model, responsible for the large reorganization energy. Single monomer absorption spectra at 77 K without static disorder calculated using the spectral densities shown in figure 3(a) are presented in figure 3(b). Our spectral density gives very similar vibrational sidebands at this temperature compared to Novoderezhkin's model. The homogeneous absorption linewidth of Jang's and Raszewski's models is similar to ours, while Raszewski's model additionally has a sharp δ-function-like peak, resembling the zero-phononline (ZPL) indicated by an arrow in figure 3(b) (this is related to ω 5 behavior of this model at ω → 0).
Our formulation allows us to describe the ensemble-averaged eigenstate properties. The excited state energy distributions of our model, which depend on contributions of these states to the eigenstates of the system and their energy spectrum, are defined by [7,31,53]   where m denotes either the FE or CT state in the site basis, e labels the single-excited eigenstate and • dis denotes averaging over disorder. We include reorganization energies λ ee,ee because the bath induces relatively large reorganization shifts. These distributions calculated using our Hamiltonian and disorder parameters are presented in figure 4. They incorporate not only energetic disorder of the states but also excitonic effects, as is most clearly seen from the fact that both special pair pigments contribute to two dimeric states. A few things should be noted. First, the Chl * D1 state has a lower energy than the lower-energy state of the special pair. Therefore, it should be energetically favorable to start charge separation from this pigment, especially at lower temperatures. Another thing to consider is that since CT states have quite wide energy distributions, it is clear that any electron transfer timescales in the RC are very strongly dependent on particular realizations of disorder. Also we note that despite having considerably larger disorder than the rest of the pigments, the peripheral chlorophylls do not show wide energy distributions, which can be attributed to their electronic states being mostly localized.

Simulations of spectra of the PSII RC
Our model of the PSII RC is able to reproduce the experimental 77 K absorption spectrum as shown in figure 5. Two peaks at ∼670 and ∼680 nm and their lineshapes are well reproduced, especially at the lower-energy shoulder. The deviation on the high-energy part of the spectrum, starting at ∼665 nm, may be related to vibrational sidebands, which we include approximately.  [14].
The simulations allow to exclude one of the two CT pathways from the model. In simulations without the Chl D1 pathway we removed the Chl + D1 Pheo − D1 state, and in simulations without the special pair pathway we removed the P + D2 P − D1 , P + D2 Chl − D1 and P + D1 Chl − D1 states. The absorption spectrum without either the Chl D1 or the special pair CT pathway is shown by points in figure 5. It can be seen that the absorption spectrum of the PSII RC is not sensitive to CT states, as can be expected. Removing CT pathways only reduces the amplitude of the ∼670 nm peak.
The main results of this work are simulations of the 2D spectra of the PSII RC at 77 K. A comparison of experimental and simulated 2D spectra for various t 2 values is presented in figure 6. In both theoretical and experimental data, the strongest peak is at ∼680 nm. The peak at ∼670 nm is much weaker. That is correctly reproduced in our model. The lower crosspeak is present in both simulations and experiment. Its amplitude and overall shape are captured by our model. Our simulations show the upper crosspeak, which, however, is absent from the experimental data. Instead, experimental spectra show a negative feature that is due to excited state absorption. That is related to double-excited states. Since our model treats the K couplings as one effective parameter, it is not able to precisely reproduce double-excited state features. This is of no consequence, however, when talking about single-excited state dynamics. Another potential explanation for lineshape distortion at early times, particularly for the crosspeaks, may come from directional filtering or signal reabsorption effects [60] in our experimental data. While our data were taken in the pump-probe geometry with a 4 • crossing angle, there might still be a directional effect. We also note that our model correctly captures both homogeneous and inhomogeneous broadening.
The time evolution of the 2D spectra is also correctly described by our model. Both simulations and experiment show decay of the ∼670 nm diagonal peak. The amplitude of this peak decays quickly within 1 ps and then continues to decay slowly up to 100 ps. The ∼680 nm peak also shows gradual decay with increasing t 2 . The lower crosspeak shows an increase in   amplitude with increased waiting time, signaling the presence of energy transfer. The upper crosspeak is absent from experimental data during the entire spectral evolution, while it remains present in the simulations. At long times (t 2 = 100 ps), the presence of this peak in the simulated 2D spectrum is the only difference from the experimental data. Simulations without either of the CT pathways show that only some regions of the 2D spectra are sensitive to CT states included in the model. Overall, all sets of simulated 2D spectra are quite similar.
To demonstrate the time evolution of the 2D spectra more explicitly, we present the kinetics of selected points of the spectra (see figure 7). We have also calculated kinetics without the Chl D1 or the special pair CT pathway. Points for the Diag1 peak were selected at maximum amplitude at the diagonal, the Diag3 peak represents 0.4 times that maximum value and the Diag2 peak is in the middle of the other two diagonal peaks. Crosspeaks CP1 and CP2 were Kinetics of selected peaks from 2D spectra of the PSII RC. Squares denote experimental data, circles-simulations, triangles-simulations without the Chl D1 CT pathway and stars-simulations without the special pair CT pathway. Kinetics for all points were normalized for both experiment and simulations by assuming the value at t 2 = 50 fs to be equal to 1.
chosen as crossing points of lines going through the diagonal peaks. Note that because of these definitions, the wavelengths of the peaks from experimental and simulated spectra do not match precisely-see figure 6 for the exact position of the selected peaks. For simulations without either of the CT pathways, the positions of selected peaks were chosen to be the same as in full simulations.
Regarding the Diag1 peak kinetics, our model does not match the initial experimental time evolution, yet gives the correct final amplitude. Simulations without either of the CT pathways show larger deviations from the experimental data. The Diag2 kinetics are well reproduced by our model for times up to 500 fs. Simulations with one CT pathway removed do not show any noticeable difference. The kinetics of the Diag3 peak are reproduced very well by our model. For this peak, simulations without the Chl D1 CT pathway show poorer agreement with experimental data, if compared to full simulations or simulations without the special pair pathway.
Turning to the crosspeak kinetics, we can see that the evolution of the CP1 peak is reproduced well by our model. Simulations without the Chl D1 pathway show poorer agreement with experiment while those without the special pair pathway are even poorer. For the CP2 kinetics, our simulations miss the initial fast component, yet the rest of the kinetics are very similar to the experiment. Simulations without either of the CT pathways are not much different than the full simulations, except for the final amplitude.

Discussion
Our optimized electron-hole TB model properly defines the single-and the double-excitation manifolds and allows simulations of the linear absorption and the 2D spectra. The simulated results show very good correspondence with the experiment. Yet there are some disagreements, which may be related to several issues. One of them is the application of the modified Redfield theory. The theory implies the existence of excitons and describes all properties in the exciton basis. In our opinion, this is the most crude approximation in our simulations since we have demonstrated that polaronic effects may be very relevant when exciton reorganization energies are significant [56].
Another important issue in the modeling of the PSII RC is the construction of the doubleexcitation manifold. Application of the TB model allowed us to correctly incorporate CT states in our model. However, since we treated K couplings as a single shift parameter, we have obtained an approximate picture of the excited state absorption. This resulted in deviations of the upper crosspeak and its time evolutions. A more detailed modeling is necessary. This could be achieved by calculating the K couplings explicitly, as was done for the FMO aggregate in [41].
Even though we use a novel theoretical formulation-the TB model-we still use many parameters determined in earlier work on the PSII RC. Notice that the TB model reduces unambiguously to the Frenkel exciton model when CT states are excluded. Hence, we must ensure that our model remains consistent with earlier modeling of various experiments. The parameters of the CT states have been optimized by Novoderezhkin et al [33,34,61]; our simulations should be consistent with these as well. Hence, CT state parameters were adapted to the theoretical framework of the TB model. Since the TB model is more detailed, some parameters remain undefined, e.g. K couplings in the double-exciton manifold. These have been adjusted to match the present experiment. As the microscopic structure of the whole protein is known, the complete construction of the TB Hamiltonian parameters should be possible as was done for the Frenkel exciton model for FMO [62], PSI [63], as well as for LH2 [64] and LH3 aggregates [65].
An important result of our simulations is that we set the disorder of the peripheral chlorophylls to be four times wider than for the rest of the pigments. Earlier PSII RC models had the same disorder for all pigments [31][32][33]53]. Peripheral chlorophylls are a considerable distance away from the core pigments (see figure 1) and their surrounding protein environment might be different, thus inducing stronger disorder. To check the validity of our assignment, we also calculated the absorption difference spectrum of full RC complexes with those lacking one of the peripheral chlorophylls, denoted in the literature as the RC5 complexes [14]. A comparison of simulated and experimental spectrum at 77 K is given in figure 8. It shows that while our simulations show wider bleaching at ∼667 nm, the difference is not that big. This implies that our assignment of larger disorder to the peripheral chlorophylls is consistent with experimental data.
In simulations, we calculate the 2D spectra without including laser pulse properties, e.g. finite pulse duration and finite bandwidth. In experiments these two properties may affect the peak profiles due to two effects. Firstly, the finite pulse duration creates the conditions in experiments where different pulses overlap and so other resonant contributions to the signal may squeeze in, apart from the rephasing and non-rephasing ones. However, this is a small contribution as shown in [46]. Secondly, the amplitudes of the peaks, which are shifted away from the pulse central frequency, become reduced to some degree by weaker laser field intensity in this specific spectral region. However, we only show experimental data outside of the pulse overlap region. We note that our simulations were done by assuming δ-shaped pulses; hence, our simulated spectra are independent of the laser pulses. The simulations, thus, show pure system properties and reveal all possible details of the system, which could be available in an ideal experiment. Here we do not attempt to reproduce spectral distortions that result from finite pulses but rather explore the predictions of our model free from such effects.
In recent years there has been a lot of discussion about the observed beatings in 2D spectra. Although usually attributed to electronic coherences [21,24,[66][67][68], there has recently been an influx of work to analyze beatings that are of vibrational origin [28,29]. The amplitude of observed beatings can decrease dramatically due to the presence of static disorder, as has been clearly shown in simulations of 2D spectra of J-aggregates [69]. Recently, it was shown that 2D spectroscopy can distinguish between decoherence and inhomogeneous dephasing [70]. In experimental 2D spectra of the PSII RC, no beatings were clearly resolved [30]. However, the use of shorter pulses and polarization-dependent studies may reveal such features. Our modeling does not show any beatings either. Due to the presence of disorder, this cannot rule out the possible existence of noticeable electronic coherences. However, we also performed simulations of 2D spectra of the PSII RC for a single realization of disorder, which showed that all beatings decay within a few tens of femtoseconds. Therefore, we can conclude that due to large reorganization energies in the PSII RC we do not have any noticeable coherent dynamics and that observed energy and CT timescales and kinetics are not affected by it. Of course, this conclusion must be supported by explicitly including vibrational modes.
One of the biggest advantages of 2D spectroscopy is its ability to estimate homogeneous broadening, coming from system-bath interaction, and inhomogeneous broadening, coming from disorder. As seen from figure 6, our model captures both of these aspects very well. The spectra at different temperatures are shifted vertically for better visibility. Simulated results were averaged over 6000 realizations of disorder. Experimental data at 77 K were taken from [14] and at 150 and 277 K from [72].
The most important parameter of the system-bath interaction is the bath spectral density. Its low-frequency part characterizes the homogeneous broadening of the single-pigment absorption spectrum (see figure 3(b)). In our model for a decaying bath correlation function in a finite time, e.g. in τ c , we get the lineshape function g (t) → t, when t > τ c . The homogeneous broadening is then described by Re = k B T · C (ω)/ω, in the limit when ω → 0. Neglecting the second term in equation (39), which is a small contribution at ω → 0, our spectral density gives Re = 2λ L k B T /γ L ∼ 100 cm −1 at 77 K temperature. Radiative and non-radiative relaxation as well as energy transfer should be taken into consideration when describing the homogeneous linewidth in a molecular aggregate. In this case each excitonic transition has its own homogeneous broadening. The homogeneous broadening and the spectral density behavior at ω → 0 become non-trivial at low (4-10 K) temperatures, where the ZPL becomes important [71]; however, this region is beyond the scope of our present study. We also note that our model is able to explain the temperature dependence of the absorption spectrum of the PSII RC, as shown in figure 9, where we obtain a good match of both the peak pattern and spectral lineshapes. This gives further credence to our choice of spectral density form and parameters.
The observed time evolution of the 2D spectra of the PSII RC contains information about the CT processes in the system. Recent publications propose the existence of two distinct electron transfer pathways [19,34]. We developed the improved model by using the TB Hamiltonian to represent both molecular excitations and CT states on an equal footing. Our model gives reasonable agreement with the experimental kinetics of 2D spectra. Therefore, we can make observations regarding the CT processes and pathways. Our simulations show that there is little difference in the simulated kinetics of Diag1, Diag2 and CP2 peaks when we remove either the Chl D1 or the special pair pathway. Only the Diag3 and CP1 peak simulated kinetics are different when different CT pathways are removed. This means that only specific regions of the 2D spectra of the PSII RC are sensitive to the CT states, which are included in the model. This can also be seen from figure 6, where we present simulations of 2D spectra with either the Chl D1 or the special pair pathway removed. We note that our simulations show that inclusion of both CT pathways provides better agreement with experimental data.

Conclusions
We have constructed a TB model Hamiltonian of the PSII RC and used it for simulating the 2D spectra. Our model properly includes CT states, and describes double excitations as well as correlated energy disorder and fluctuations. We find that the peripheral chlorophylls have four times larger disorder than the other RC pigments.
We showed that our model is capable of reproducing the majority of features present in the experimental 2D spectra. Excellent agreement of homogeneous and inhomogeneous broadening in experimental and simulated spectra gives credibility to our choice of bath spectral density and disorder parameters. Analysis of 2D spectra kinetics showed that inclusion of the two electron transfer pathways helps us obtain better agreement with experiment.
In the double-excitation manifold, the Hamiltonian matrix elements are e k e l h m h n |Ĥ S |e k e l h m h n = t e kk δ ll + t e kl δ lk + t e lk δ kl + t e ll δ kk δ mm δ nn + t h mm δ nn + t h mn δ nm + t h nm δ mn + t h nn δ mm δ kk δ ll + V e kl + V h mn − V eh km − V eh lm − V eh kn − V eh ln + K km,ln δ kk δ ll δ mm δ nn