Photosystem II (PSII) is the only biological system capable of splitting water to molecular oxygen. Its reaction center (RC) is responsible for the primary charge separation that drives the water oxidation reaction. In this work, we revisit the spectroscopic properties of the PSII RC using the complex time-dependent Redfield (ctR) theory for optical lineshapes [A. Gelzinis *et al.*, J. Chem. Phys. **142**, 154107 (2015)]. We obtain the PSII RC model parameters (site energies, disorder, and reorganization energies) from the fits of several spectra and then further validate the model by calculating additional independent spectra. We obtain good to excellent agreement between theory and calculations. We find that overall our model is similar to some of the previous asymmetric exciton models of the PSII RC. On the other hand, our model displays differences from previous work based on the modified Redfield theory. We extend the ctR theory to describe the Stark spectrum and use its fit to obtain the parameters of a single charge transfer state included in our model. Our results suggest that $ChlD1+PheoD1\u2212$ is most likely the primary charge transfer state, but that the Stark spectrum of the PSII RC is probably also influenced by other states.

## I. INTRODUCTION

Photosynthesis begins with the absorption of solar light by specific pigments, (bacterio)chlorophylls and carotenoids, that are bound to light-harvesting proteins. Subsequently the energy is transferred to reaction centers, where the excitation energy triggers a fast charge separation, through which it is transduced into chemical potential. Oxygenic photosynthesis in plants, cyanobacteria, and algae begins in photosystem II (PSII), which is the only known biological complex with strong enough redox potential to split water. PSII is a dimeric supercomplex, consisting of major and minor light harvesting complexes and the reaction center (RC). Excitation within the light-harvesting antennae of PSII is transferred to the RC, leading to the oxidation of water, the reduction of plastoquinone, and the formation of a proton gradient across the thylakoid membrane. The PSII is thus the main producer of oxygen on Earth.^{1,2}

The PSII RC itself contains six chlorophyll (Chl) and two pheophytin (Pheo) molecules arranged into two almost symmetric branches, corresponding to the D1 and D2 proteins.^{3} The pigment arrangement in the PSII RC is shown in Fig. 1. Two closely situated Chls, P_{D1} and P_{D2}, are known as the “special pair” (from the analogous attribution in the bacterial RC). Each branch also contains an accessory chlorophyll Chl_{D1/D2} and a pheophytin Pheo_{D1/D2}. Two peripheral chlorophylls Chlz_{D1/D2} are about 25 Å away from the other six RC pigments, and they do not participate in the primary charge separation. *In vivo*, the electron transfer proceeds through the D1 branch to Pheo_{D1} and then further to plastoquinone Q_{A} and finally to Q_{B}.^{2} The most often studied PSII RC preparation (the so-called D1D2cyt*b*_{559} complex), however, lacks the oxygen evolving complex, the redox active tyrosine Tyr_{Z} and the quinone acceptors,^{4,5} thus allowing the electrons to accumulate on Pheo_{D1}. Meanwhile, the hole remains mostly localized on P_{D1}.

The energy and charge transfer processes in the isolated PSII RC, the PSII core (which includes the CP43 and CP47 antenna complexes), or the whole PSII have been extensively studied over the last few decades, as reviewed in Refs. 1 and 6–12. Nonetheless, the pathways of charge separation and the identities of the participating states are still under debate. Many studies claim that the primary electron donor is the Chl_{D1} pigment; thus, the primary charge transfer (CT) state is $ChlD1+PheoD1\u2212$.^{13–17} In addition, a series of studies^{18–20} have highlighted the importance of the special pair CT state $PD2+PD1\u2212$, which has been proposed to be strongly mixed with the exciton states, providing an alternative charge separation pathway. These suggestions have to be tested against an increasing amount of experimental data using appropriate theoretical modeling.

Most of the theoretical work on the PSII RC is based on the exciton concept.^{21} The first exciton model that captured some features of the PSII RC spectroscopic data was the so-called multimer model.^{22} It was formulated by taking the structural data of the bacterial reaction center into account and additionally assuming equal transition energies for all the core pigments. Later, as the PSII structural data became available with increasing resolution,^{3,23–25} the theoretical models were refined,^{16,18,20,26–30} and many spectroscopic properties of the PSII RC were successfully described. Some models have included one or more CT states phenomenologically, treating them as additional exciton states.^{18,20,26,28} Only a few studies have accounted for the relevant properties of electrons and holes.^{29,31}

Although there have been attempts to obtain the PSII RC site energies directly from the structural data,^{30,33,34} the majority of the theoretical models derive their parameters from fits to the optical spectra. Since exact calculations are almost always too computationally expensive, approximate methods have to be applied to calculate optical lineshapes and dynamics. We have recently analyzed various approaches to calculating absorption lineshapes^{35} and have found that the complex time dependent Redfield (ctR) theory gives results very close to the exact ones. On the other hand, the modified Redfield theory, used in Refs. 18, 20, 26, and 28 and in our previous work on the PSII RC^{29}, was found to perform significantly worse. The Markovian Redfield approach, used in Refs. 16, 27, and 30, showed improvements but did not perform as well as the ctR theory. These findings raise some doubts about the optimization of the parameters of the previous models of the PSII RC^{16,18,20,26–29} and motivate our use of ctR theory for describing its spectroscopic properties.

The charge transfer states, being optically dark, are insensitive to most spectroscopic measurements like absorption, linear, or circular dichroism. Nonetheless, due to their strong static dipole moments, they are sensitive to the presence of static electric fields and thus can influence the results of Stark spectroscopy. Stark experiments have been performed on the PSII RC,^{36} but the results were not easy to understand without a detailed theoretical description. Previously, the PSII RC Stark spectrum was fitted using the modified Redfield theory.^{18} Even though significant progress in understanding the observed features was achieved, a check of the obtained conclusions, especially the parameters of the CT states, using a more accurate lineshape theory is definitely called for.

In this work, we revisit the spectroscopic properties of the PSII RC using the ctR theory.^{35} We construct a model that includes a single CT state together with molecular excitations. We obtain model parameters by fitting a number of linear optical spectra. We also extend the ctR theory to calculate the Stark spectrum, which is then used to tune the parameters of the CT state. Due to the more advanced lineshape theory, in our current model, we improve upon the parameters obtained in our previous work based on the modified Redfield theory.^{29}

## II. THEORETICAL BACKGROUND

### A. Tight-binding Hamiltonian

To describe the molecular system, we use the tight-binding Hamiltonian,^{29,31} which is an extension of the Frenkel exciton model^{21} that enables the inclusion of CT states. This method is similar to other recent descriptions of Frenkel excitons and CT states in molecular systems.^{37–39} In Secs. III and IV, this model will be applied to the PSII RC, but its formulation is independent of any particular molecular system.

We consider the ground and the single-excited state manifolds that are sufficient to calculate linear spectra. The global ground state is denoted by *g*, and its energy is set to zero. In the single-excited state manifold, we have the molecular (pigment) excitation states *n*^{*}, with excitation localized on a single site *n*, and the charge transfer states *m*^{+}*n*^{−}, which are not accounted for in the Frenkel exciton model. To describe the charge transfer states, we introduce separate creation and annihilation operators for electrons and holes. Thus $e^n\u2020$ $e^n$ creates (annihilates) an electron on site *n*, and $h^n\u2020$ $h^n$ creates (annihilates) a hole on site *n*. Then molecular excitation (ME) states are described by $e^n\u2020h^n\u2020|g=|n*$ and CT states by $e^n\u2020h^m\u2020|g=|m+n\u2212$. The electron and hole creation and annihilation operators satisfy the Fermi commutation relations.^{29,31}

The total Hamiltonian is composed of Hamiltonians of the system, the bath, and their interaction,

The electronic degrees of freedom compose the system that is described by the tight-binding Hamiltonian

Here $tmne$$tmnh$ is the electron (hole) hopping rate between LUMO (HOMO) orbitals of different pigments. Diagonal elements $tmme/h$ correspond to the orbital energies. *N* is the number of molecules, and the total number of states in the single excited state manifold is *N*_{tot} = *N*^{2}. *J*_{mn} is the resonance coupling between the molecular excitations on sites *m* and *n*. The term $Vmneh$ corresponds to the electron-hole Coulomb attraction energy between sites *m* and *n*. Note that electron-electron and hole-hole repulsion terms are not needed to describe the single excitation manifold. For this Hamiltonian, the energies of the excited states are

Here $\epsilon n*$ is the energy of the ME state *n*^{*} (site energy) and $\epsilon m+n\u2212$ is the energy of the CT state *m*^{+}*n*^{−}. Note that for calculations only the energies and not their constituent terms need to be specified. The couplings between different excited states are

Thus, couplings involving the CT states are related, which is not the case if the CT states are treated as additional molecular excitation states.^{18,20,26}

All vibrational degrees of freedom of the molecules themselves and their environment compose the bath and are represented by an infinite set of harmonic oscillators with the Hamiltonian,

Here $\omega j$ is the frequency, $p^j$ is the dimensionless momentum, and $x^j$ is the dimensionless coordinate of the *j*th bath mode. We set $\u210f=1$ throughout the paper. The system-bath coupling is assumed to be linear in the bath coordinates

Here the two terms are separated for future convenience. The parameter $dnn,jeh$$dmn,jeh$ is the dimensionless displacement of the equilibrium configuration of the *j*th bath mode between the ground state and the excited state of the *n*th molecule (*n*^{+}*m*^{−} CT state). This formulation differs from the previous applications of the tight-binding model, where the bath was coupled to electron and hole orbitals and not to the excited states themselves.^{29,31} The present approach is more natural when describing excitations in a molecular aggregate and not in a solid state system.

We introduce bath operators $F^n=\u2212\u2211j\omega jdnn,jehx^j$ and $F^mn=\u2212\u2211j\omega jdmn,jehx^j$, enabling us to rewrite the system-bath coupling Hamiltonian as

For simplicity, we subsequently denote both the ME and the CT states by letters with tildes $\xe3,b\u0303,\u2026$ Then summation over all single excited states is expressed as

Here $z\xe3$ denotes any quantity related to the state $\xe3$. The reorganization energies of the excited states are defined as $\lambda \xe3=\u2211j\omega jd\xe3,jeh2/2$.

For a harmonic bath, the key quantities that describe the bath effects are the correlation functions,^{40–42}

Here

defines the operator in the interaction representation. The equilibrium bath density operator is

with $\beta =kBT\u22121$ being the inverse temperature in energy units, and *k*_{B} is the Boltzmann constant. The spectral densities

describe the bath in the frequency domain. They show how strongly the system states are coupled to the specific frequency bath modes. From the spectral densities, the bath correlation functions can be calculated as

The reorganization energies are also related to the spectral densities

Diagonalization of the system Hamiltonian [Eq. (2)] leads to the introduction of the generalized exciton basis (denoted as such, because it includes not only Frenkel but also CT excitons),

with $\psi e,\xe3=\u27e8e|\xdb\u2020|\xe3\u27e9$ being the elements of transformation matrix $\xdb\u2020$, which satisfies

Here the superscript EB (SB) means the generalized exciton (site) basis. We assume that the matrix elements $\psi e,\xe3$ are real. In the generalized exciton basis, the bath correlation functions are

with $F^e1e2=\u27e8e1|\u0124SB|e2\u27e9$. The lineshape functions, required to describe the optical lineshapes (see Sec. II B), are defined as

When static disorder is present in the system, the energies of the states are different in each particular case of the ensemble. Thus in the system Hamiltonian [Eq. (2)], one needs to change $Vnneh\u2192Vnneh+\delta \epsilon n*$ and $Vmneh\u2192Vmneh+\delta \epsilon n+m\u2212$, with $\delta \epsilon n*$ and $\delta \epsilon n+m\u2212$ being random numbers from a Gaussian distribution with zero mean and some chosen standard deviations $\sigma n*$ and $\sigma n+m\u2212$. All observable quantities then have to be averaged over disorder, i.e., $X\xaf=Xdis$.

### B. Calculations of optical spectra

#### 1. Linear absorption

For calculations of the linear absorption, one has to include the interaction with a time-dependent electric field in the Hamiltonian, which induces the transitions between the ground and excited states,

Here $\mu \u2192^$ is the electric dipole operator, $\mu \u2192n$ is the transition dipole vector for pigment *n*, and $o\u2192$ is the polarization vector of the electric field. For now we disregarded the permanent dipole moments (the diagonal elements of the dipole operator). Note that the CT states are dark and carry no transition dipole strength. We also define the scalar dipole moment operator $\mu ^=\mu \u2192^\u22c5o\u2192$ for later convenience.

Assuming linearly polarized incoming light, the linear absorption spectrum is calculated as a one-sided Fourier transform of the linear response function,^{35,40,42–44}

The linear response function is formally a tensor

Here $\rho ^gg=|gg|\rho ^Beq$ is the ground state equilibrium density operator, $o\u21921$ and $o\u21922$ denote the polarizations of the incoming light and registered signal, and $\u2022or$ denotes orientational averaging, which has to be performed for randomly oriented systems. In this expression, we have assumed the rotating wave approximation (RWA). For isotropic systems, the tensorial nature of the response function can be neglected and scalar response function can be defined as $R1t=R1t;o\u2192,o\u2192$, with $o\u2192$ being an arbitrary polarization vector.

In the generalized exciton basis, the response function is expanded as

Assuming the secular approximation (considering only terms *e*_{1} = *e*_{2}, thus no coherence transfer) and using the cumulant expansion^{40} and the ctR theory for the off-diagonal fluctuations, we obtain^{35}

Here $\mu \u2192e1g=\u27e8e1|\mu \u2192^|g\u27e9$ is the excitonic transition dipole moment, and $ge1t=ge1e1,e1e1t$. The ctR theory term is

#### 2. Linear dichroism

Linear dichroism (LD) spectroscopy is performed on oriented samples (they are no longer isotropic). For simplicity we orient the molecular coordinate system in such a way that the membrane normal lies along the *z* axis. The linear dichroism spectra then corresponds to a difference between two absorption measurements: one where the incoming light is polarized parallel to the membrane plane, and the other, when the incoming light is polarized perpendicular to the membrane plane^{46}

The required absorption spectra are then defined by

The response functions are calculated using the same approximations as above. The relevant orientational average is

Finally, the LD spectrum can be calculated as

with

#### 3. Circular dichroism

The circular dichroism spectrum is defined as a difference between the absorption of left- and right-circularly polarized light^{21,46,47}

To properly describe the CD spectrum, one has to take into account the coupling with the magnetic field. Therefore, in this case, the system-field coupling Hamiltonian is

Here $m\u2192^$ is the magnetic dipole moment operator and $o\u2192B$ is the polarization vector of the magnetic field. Following Ref. 44, we use an effective molecular magnetic dipole moment defined as

Here $R\u2192n$ is the position vector corresponding to the center of the pigment *n*. The CD spectrum can also be calculated from the response function

The relevant response function is found to be^{44}

For an isotropic medium, we can again set $RCD1t=RCD1t;o\u2192,o\u2192$, though note that the same vector $o\u2192$ is taken as both arguments. Using the same level of theory as above, the response function is

Here $m\u2192e1g=\u27e8e1|m\u2192^|g\u27e9$ is the excitonic magnetic dipole moment.

We note that this level of theory neglects the intrinsic CD of the pigments. In addition, the CD spectra can be affected by couplings to the higher energy states,^{21,48} which we ignore in the present treatment.

#### 4. Stark spectroscopy

Stark spectroscopy relies on the effect of an applied static electric field on an absorption spectrum and is defined as the difference between the spectrum with the field on and the field off^{49}

Theoretically, the key issue is to calculate the absorption spectrum with the presence of an electric field. For this, the system-field Hamiltonian is extended to also include the static field

Here the static electric field is denoted by $E\u2192$, and we assume that the electric dipole operator $\mu \u2192^$ stays constant even in the presence of the electric field—therefore we are describing only the linear Stark effect. The term $\u2212\mu \u2192^\u22c5E\u2192$ can be added to the system Hamiltonian because it is constant. Here we can no longer neglect the diagonal elements of the dipole moment operator. In fact, they cause a shift in the excited state energies,

Here,

is the difference between the permanent dipole moments of the excited and the ground states and will be called the static dipole moment for convenience. We have to subtract the ground state energy because we have set it to zero. Thus, for simplicity, we redefine the dipole moment operator

with $\xce$ being the identity operator.

In this work, we follow Ref. 50 and use perturbation theory to calculate the Stark spectrum. To this end, we first write the expression for the linear response function in the presence of the static electric field

Here $o\u2192E$ is the polarization vector of the static electric field, and $E$ is its magnitude. Note that here the interaction picture is assumed with respect to $\u0124$ [see Eq. (1)], contrary to Eq. (13). If we expand this expression in powers of $E$ (this is done by expanding the time-ordered exponentials and collecting the terms), then the zeroth order term gives the usual absorption spectrum, the first order term vanishes for isotropic samples, and the second order term is responsible for the Stark signal.

Thus, the Stark spectrum response function is

Here we have set $o\u21921=o\u21922=o\u2192E$. Expansion of these three terms in the eigenstate basis of the system Hamiltonian produces many contributions to the Stark spectrum. Nonetheless, the majority of them are nonresonant, i.e., oscillating with factors other than $e\u2212i\omega \xafet$, $\omega \xafe$ being the band gap between the ground and the single excited state manifolds. In fact, only a single contribution survives the RWA (this limit was denoted as pure RWA in Ref. 18),

Here we have denoted $\mu abE/E=\u27e8a|\mu \u2192^\u22c5o\u2192E/E|b\u27e9$. Then, using the secular approximation and performing the second order cumulant expansion, we obtain the following expression, which includes the ctR theory terms to describe the off-diagonal fluctuations:

If we assume magic angle $cos\theta m.\u2009a.=1/3$ between the static and probing electric fields, then the orientational averaging results in18^{,}45

The Stark spectrum is finally calculated as

Expression (47) is one of the major results of the present work.

## III. MODEL OF THE PSII RC

Since the PSII RC has 8 pigments, in addition to 8 ME states, there could be 56 possible CT states. If we limit the CT states to the particular pigments (Chl_{D1}, Pheo_{D1}, P_{D1}, and P_{D2}) that have been proposed to be involved in charge separation and take into consideration the directionality of the charge transfer, this number is significantly reduced to 6. Restricting the charge transfer states to neighboring pigment pairs further reduces this number to 3. Given that the majority of the spectra that we fit (absorption, LD, CD) are not significantly affected by the dark CT states, we include a single CT state in our model. We note that only the Stark spectrum is sensitive to CT states, and previous attempts to fit the Stark spectrum of the PSII RC included a single CT state.^{18} We performed calculations that considered either the $PD2+PD1\u2212$ or $ChlD1+PheoD1\u2212$ CT states and from the fits of the Stark spectrum (see Sec. IV for details) we have chosen $ChlD1+PheoD1\u2212$ as the primary CT state. Thus, our model includes 9 states in the single excitation manifold.

We obtain the site energies of the pigments from the fitting procedure. The resonance couplings between the pigments were taken from Ref. 51, where they were calculated using the Poisson TrEsp method.^{52–54} The CT state energy and the coupling between CT state and molecular states were also obtained from the fits. Formally, there are two couplings involving a single CT state,

but we have chosen $tmne=tmnh$ for simplicity. The entire single excitation Hamiltonian matrix is given in Table I.

. | P_{D1}
. | P_{D2}
. | Chl_{D1}
. | Chl_{D2}
. | Pheo_{D1}
. | Pheo_{D2}
. | Chlz_{D1}
. | Chlz_{D2}
. | $ChlD1+PheoD1\u2212$ . | $E\lambda $ . | $\sigma $ . |
---|---|---|---|---|---|---|---|---|---|---|---|

P_{D1} | 15 080 | 158 | −27.3 | −41.8 | −4 | 12.6 | 0.5 | 0.6 | 0 | 663 | 80 |

P_{D2} | 158 | 15 015 | −46.8 | −22 | 15.1 | −3 | 0.6 | 0.6 | 0 | 666 | 80 |

Chl_{D1} | −27.3 | −46.8 | 14 800 | 3.5 | 43.5 | −2.2 | 1.7 | −0.1 | 125 | 675.5 | 80 |

Chl_{D2} | −41.8 | −22 | 3.5 | 15 010 | −2.4 | 41.7 | −0.1 | 1.8 | 0 | 666 | 80 |

Pheo_{D1} | −4 | 15.1 | 43.5 | −2.4 | 14 950 | 1.5 | −2.5 | −0.2 | 125 | 669 | 80 |

Pheo_{D2} | 12.6 | −3 | −2.2 | 41.7 | 1.5 | 14 850 | −0.2 | −2.6 | 0 | 673.5 | 80 |

Chlz_{D1} | 0.5 | 0.6 | 1.7 | −0.1 | −2.5 | −0.2 | 14 985 | 0.1 | 0 | 667.5 | 80 |

Chlz_{D2} | 0.6 | 0.6 | −0.1 | 1.8 | −0.2 | −2.6 | 0.1 | 15 065 | 0 | 664 | 80 |

$ChlD1+PheoD1\u2212$ | 0 | 0 | 125 | 0 | 125 | 0 | 0 | 0 | 14 270 | 701 | 550 |

. | P_{D1}
. | P_{D2}
. | Chl_{D1}
. | Chl_{D2}
. | Pheo_{D1}
. | Pheo_{D2}
. | Chlz_{D1}
. | Chlz_{D2}
. | $ChlD1+PheoD1\u2212$ . | $E\lambda $ . | $\sigma $ . |
---|---|---|---|---|---|---|---|---|---|---|---|

P_{D1} | 15 080 | 158 | −27.3 | −41.8 | −4 | 12.6 | 0.5 | 0.6 | 0 | 663 | 80 |

P_{D2} | 158 | 15 015 | −46.8 | −22 | 15.1 | −3 | 0.6 | 0.6 | 0 | 666 | 80 |

Chl_{D1} | −27.3 | −46.8 | 14 800 | 3.5 | 43.5 | −2.2 | 1.7 | −0.1 | 125 | 675.5 | 80 |

Chl_{D2} | −41.8 | −22 | 3.5 | 15 010 | −2.4 | 41.7 | −0.1 | 1.8 | 0 | 666 | 80 |

Pheo_{D1} | −4 | 15.1 | 43.5 | −2.4 | 14 950 | 1.5 | −2.5 | −0.2 | 125 | 669 | 80 |

Pheo_{D2} | 12.6 | −3 | −2.2 | 41.7 | 1.5 | 14 850 | −0.2 | −2.6 | 0 | 673.5 | 80 |

Chlz_{D1} | 0.5 | 0.6 | 1.7 | −0.1 | −2.5 | −0.2 | 14 985 | 0.1 | 0 | 667.5 | 80 |

Chlz_{D2} | 0.6 | 0.6 | −0.1 | 1.8 | −0.2 | −2.6 | 0.1 | 15 065 | 0 | 664 | 80 |

$ChlD1+PheoD1\u2212$ | 0 | 0 | 125 | 0 | 125 | 0 | 0 | 0 | 14 270 | 701 | 550 |

Our model takes into account the disorder value of the energies of the pigments and CT states. We assume the disorder value for all pigments to be the same, and the CT state to have different disorder values. In both cases, the disorder value was described by a Gaussian distribution with standard deviations obtained from the fits (given in Table I).

The strengths of the transition dipole moments of Chl and Pheo pigments are taken to be 4.4 and 3.4 D, respectively.^{27,29} The transition dipole vectors are assumed to be directed from *N*_{B} to *N*_{D} atoms (for nomenclature see, e.g., Ref. 61) as obtained from the 3WU2 crystal structure.^{3} The CT state was assumed to have no transition dipole strength. The static dipole moments are assumed to be rotated by $15\xb0$ in the direction of the *N*_{C} atom from the *N*_{B} to *N*_{D} axis (this was based on the calculations of absorption difference spectra corresponding to the $PD1+PheoD1\u2212$ CT state, see Sec. IV for details). Their strengths are taken to be 1.55 and 1.3 D for Chls and Pheos, respectively. The static dipole moment of the CT state was calculated by assuming charges fully localized at the geometric centers of the four *N* atoms of the involved pigments. This gives a static dipole strength of 49.3 D for the $ChlD1+PheoD1\u2212$ state. The transition dipole moments, the static dipole moments, and the geometric centers of the pigments of the PSII RC are given in Table II. Here the *z* axis is perpendicular to the membrane plane, with data taken from the 3WU2 crystal structure^{3} obtained from the OPM database.^{55}

. | $\mu x$ . | $\mu y$ . | $\mu z$ . | $\Delta \mu x$ . | $\Delta \mu y$ . | $\Delta \mu z$ . | R_{x}
. | R_{y}
. | R_{z}
. |
---|---|---|---|---|---|---|---|---|---|

P_{D1} | 1.08 | −3.74 | 2.05 | 0.32 | −1.07 | 1.07 | 37.87 | 4.05 | 11.99 |

P_{D2} | −0.66 | 3.73 | 2.24 | −0.13 | 1.03 | 1.15 | 32.57 | 9.55 | 11.12 |

Chl_{D1} | −0.33 | 4.06 | 1.67 | 0.29 | 1.34 | 0.73 | 30.65 | −2.13 | 9.30 |

Chl_{D2} | 0.11 | −4.04 | 1.74 | −0.36 | −1.31 | 0.75 | 40.19 | 15.45 | 8.24 |

Pheo_{D1} | 1.65 | −0.55 | −2.92 | 0.84 | 0.10 | −0.99 | 33.73 | −4.13 | −0.31 |

Pheo_{D2} | −1.77 | 0.19 | −2.90 | −0.89 | −0.19 | −0.93 | 36.67 | 16.71 | −1.41 |

Chlz_{D1} | −1.91 | −0.64 | 3.91 | −1.01 | −0.32 | 1.13 | 31.92 | −27.24 | 8.07 |

Chlz_{D2} | 2.07 | 0.98 | 3.76 | 1.06 | 0.41 | 1.05 | 38.91 | 39.35 | 4.76 |

. | $\mu x$ . | $\mu y$ . | $\mu z$ . | $\Delta \mu x$ . | $\Delta \mu y$ . | $\Delta \mu z$ . | R_{x}
. | R_{y}
. | R_{z}
. |
---|---|---|---|---|---|---|---|---|---|

P_{D1} | 1.08 | −3.74 | 2.05 | 0.32 | −1.07 | 1.07 | 37.87 | 4.05 | 11.99 |

P_{D2} | −0.66 | 3.73 | 2.24 | −0.13 | 1.03 | 1.15 | 32.57 | 9.55 | 11.12 |

Chl_{D1} | −0.33 | 4.06 | 1.67 | 0.29 | 1.34 | 0.73 | 30.65 | −2.13 | 9.30 |

Chl_{D2} | 0.11 | −4.04 | 1.74 | −0.36 | −1.31 | 0.75 | 40.19 | 15.45 | 8.24 |

Pheo_{D1} | 1.65 | −0.55 | −2.92 | 0.84 | 0.10 | −0.99 | 33.73 | −4.13 | −0.31 |

Pheo_{D2} | −1.77 | 0.19 | −2.90 | −0.89 | −0.19 | −0.93 | 36.67 | 16.71 | −1.41 |

Chlz_{D1} | −1.91 | −0.64 | 3.91 | −1.01 | −0.32 | 1.13 | 31.92 | −27.24 | 8.07 |

Chlz_{D2} | 2.07 | 0.98 | 3.76 | 1.06 | 0.41 | 1.05 | 38.91 | 39.35 | 4.76 |

We assume that fluctuations of the energies are uncorrelated, thus $CAB\u2032\u2032\omega =CA\u2032\u2032\delta AB$. All pigments (both Chls and Pheos) are described by the same spectral density consisting of two parts, $Cm*\u2032\u2032\omega =CLF\u2032\u2032\omega +CHF\u2032\u2032\omega $. In our previous work,^{29} we assumed that the low frequency part was given by a Debye spectral density. However, this form is unsuitable for low temperature spectra calculations due to the divergence of the total Huang-Rhys factor.^{62} Instead we adopt the superohmic spectral density that is free from this deficiency, consisting of two terms,

The two cutoff frequencies were chosen to be $\omega 1c=20\u2009cm\u22121$ and $\omega 2c=80\u2009cm\u22121$. The values for reorganization energy were obtained from the fits of spectra and found to be $\lambda 1=15\u2009cm\u22121$ and $\lambda 2=35\u2009cm\u22121$. The high frequency part is represented by 8 underdamped vibrational modes,

The damping parameter was chosen to be the same for all modes, $\Gamma \u22121=1\u2009ps$. The frequencies *W*_{j} and the Huang-Rhys factors *s*_{j} were taken from $\Delta FLN$ measurements of Chl *a* in solution^{56} and are given in Table III. We included only the modes with center frequencies less than 750 $cm\u22121$, as higher frequency modes are not expected to influence system dynamics, being completely nonresonant. For the CT state, we keep the high frequency part of the spectral density the same but introduce an additional multiplicative factor for the low frequency part, as CT states are expected to couple more strongly to the environment. Thus, $CCT\u2032\u2032\omega =\eta CLF\u2032\u2032\omega +CHF\u2032\u2032\omega $. From the fits of the Stark spectrum, we obtained $\eta =12$ for the $ChlD1+PheoD1\u2212$ state.

Mode . | 1 . | 2 . | 3 . | 4 . | 5 . | 6 . | 7 . | 8 . |
---|---|---|---|---|---|---|---|---|

$Wj,\u2009cm\u22121$ | 93 | 189 | 262 | 352 | 389 | 519 | 573 | 745 |

s_{j} | 0.012 | 0.010 | 0.013 | 0.023 | 0.015 | 0.016 | 0.017 | 0.034 |

Mode . | 1 . | 2 . | 3 . | 4 . | 5 . | 6 . | 7 . | 8 . |
---|---|---|---|---|---|---|---|---|

$Wj,\u2009cm\u22121$ | 93 | 189 | 262 | 352 | 389 | 519 | 573 | 745 |

s_{j} | 0.012 | 0.010 | 0.013 | 0.023 | 0.015 | 0.016 | 0.017 | 0.034 |

## IV. RESULTS

### A. Fits of optical spectra

Most of the model parameters were obtained from a fit of nine linear spectra shown in Fig. 2. In particular, we fitted absorption at 6, 77, and 298 K, linear and circular dichroism at 6 and 298 K, and absorption difference spectra with exchanged pheophytins at 6 K. The latter spectra correspond to the so-called RC1x and RC2x complexes,^{57,58} where in the RC1x complexes, 80% of the Pheo_{D2} is estimated to be replaced by 13^{1}-OH-Pheo, and in the RC2x complexes, all Pheo_{D2} and 50% of Pheo_{D1} are replaced. When calculating RC1x/RC2x spectra, we took this into account accordingly. The energies of the exchanged Pheos were also free parameters of the fit. In the fitting procedure, we minimized the sum of the deviations from the experimental data at each wavelength for each spectra and added them with different weights (1.0 for Abs 6 K and 298 K, 2.0 for Abs 77 K, 5.0 for Abs RC-RC1x/RC2x 6 K, 0.05 for LD 6 K, 0.005 for CD 6 K and 298 K, and 0.01 for LD 298 K). Smaller weights were chosen for some spectra due to the poorer quality of the fit. Before calculating the differences, the calculated spectra were normalized in the following way: absorption spectra were normalized by setting the maximum amplitudes to unity, CD/LD spectra were normalized in such a way that the difference between the experimental and calculated spectra would be the smallest, and for absorption difference RC and RC1x/RC2x spectra were normalized to unity separately (this was done also to the experimental data) and then subtracted. Fitting was performed using R software^{63} with a differential evolution algorithm^{64} as implemented in Ref. 65.

We obtain excellent fits for all the absorption experiments, and good qualitative fits for circular dichroism (the CD spectrum of PSII RC is nonconservative, while our theoretical approach can give only conservative spectra in the *Q*_{y} region). There is some quantitative difference for the calculated linear dichroism spectra, with positive feature $\u223c680\u2009nm$ having less amplitude than in the experimental data. Overall, our results are of similar quality to previous studies on the PSII RC.^{16,18}

We have adjusted the CT state parameters (its energy, coupling with molecular states, disorder values and multiplicative factor describing coupling with low frequency vibrations) from the fits of the Stark spectrum measured at 77 K,^{36} with results given in Fig. 3. During the fitting procedure, both the experimental and the simulated Stark spectra were normalized with respect to the absorption spectrum. We tried including either $PD2+PD1\u2212$ or $ChlD1+PheoD1\u2212$ CT states in our model. We found that a better fit was obtained with the $ChlD1+PheoD1\u2212$ state. Results with the $PD2+PD1\u2212$ state are also reasonable but are somewhat worse—the amplitude of the $\u223c680\u2009nm$ is underestimated, and a negative feature $\u223c700\u2009nm$ appears, which is absent in the experimental data. Therefore, we hold that $ChlD1+PheoD1\u2212$ should be the primary CT state, contrary to the conclusions obtained previously.^{18} We also note that calculations without any CT states roughly reproduce the amplitude of the $\u223c670\u2009nm$ feature but severely underestimate the amplitude of the $\u223c680\u2009nm$ feature. This is a clear indication that CT states significantly influence the Stark spectrum of the PSII RC, much more than in previous calculations.^{18}

### B. Calculations of independent spectra

Having performed the fits and obtained the optimal values of the parameters, we checked our model by calculating other linear spectra, for which experimental results are available in the literature.

We calculated the CD spectra of different PSII RC preparations (normal RC, RC1x, and RC2x) at 6 K. The comparison of calculations with experimental data (taken from Ref. 58) is presented in Fig. 4. Our simulations qualitatively reproduce the changes observed in the CD spectra upon the increase of exchanged Pheo (going from RC to RC1x to RC2x). First, the main positive feature at $\u223c680\u2009nm$ loses a significant amount of amplitude. Second, a positive feature at $\u223c675\u2009nm$ appears. Third, the negative feature at $\u223c670\u2009nm$ loses its amplitude. Fourth, a negative feature at $\u223c655\u2009nm$ appears.

We also calculated LD spectra at 6 K of RC, RC1x, and RC2x preparations. Comparison with experimental data is presented in Fig. 5. Again, our simulations qualitatively reproduce the main changes upon the increase of exchanged Pheo—appearance of negative feature at $\u223c655\u2009nm$, lack of changes in negative feature at $\u223c670\u2009nm$, and increase of amplitude of positive feature at $\u223c680\u2009nm$. Nonetheless, we have to note that the overall amplitudes of both the negative feature at $\u223c655\u2009nm$ and the positive feature at $\u223c680\u2009nm$ are considerably smaller in our simulations.

In Fig. 6, we compare the calculated triplet minus singlet absorption difference spectrum of the PSII RC with the experimental data.^{58} Calculations were performed assuming triplet localized at Chl_{D1} and ignoring its couplings and dipole strength. We obtain good agreement with experimental data, similar to other studies.^{16,18,30}

We calculated absorption difference spectra of normal RC preparations (here called RC6) and the so-called RC5 preparations, which lack one of the peripheral chlorophylls.^{66,67} Comparison with the experimental data is presented in Fig. 7. Experimental spectra were calculated by normalizing both RC6 and RC5 absorption at the $\u223c678\u2009nm$ peak and subtracting the spectra, then normalizing with respect to the RC6 spectra. Simulated spectra were obtained by calculating normal RC6 spectra, then calculating the spectra without the Chlz_{D1/D2} pigment, subtracting, and normalizing with respect to the RC6 spectra. Our calculations are consistent with Chlz_{D2} being lost in the RC5 preparations, with both peak position at $\u223c666\u2009nm$ and its half-width agreeing with the experimental data. The amplitude of the $\u223c666\u2009nm$ peak is higher in the calculations, which could be related to differences in the low temperature experimental absorption spectra, as the amplitude of $\u223c672\u2009nm$ peak varies in different reports.^{57–59,66}

In Fig. 8, we show a comparison of experimental and calculated RC5-$RC5PheoD1\u2212$ absorption difference spectra (here $RC5PheoD1\u2212$ means RC5 with reduced Pheo_{D1}), with experimental data taken from Ref. 68. Calculations were done assuming that in RC5 preparations the Chlz_{D2} is lost, consistent with our RC6-RC5 absorption difference calculations. When calculating spectra with reduced Pheo_{D1}, we neglected the couplings and dipole strength of this pigment. In addition, as suggested in Ref. 16, we took the electrochromic shifts into account. This was done by assuming a negative charge localized at the geometric center of the four nitrogen atoms of Pheo_{D1} and then calculating the interaction of that charge with the static dipole moments of the other pigments, that is,

Here $\mathit{\epsilon}r$ is the relative dielectric constant, which was assumed to be equal to 2.5 in our calculations. The CT state was neglected in the calculations for the reduced Pheo_{D1} since we expect that due to its huge energy shift, it would effectively decouple from the ME states. As seen from Fig. 8, we obtain excellent agreement for 77 K data, with both the main negative feature at $\u223c682\u2009nm$ and the positive feature at $\u223c673\u2009nm$ reproduced. In calculations for 277 K, we still reproduce the main negative feature but obtain a larger amplitude for the positive one compared to the experimental data.

We also calculated the absorption difference spectrum corresponding to the CT state $PD1+PheoD1\u2212$ taken from pump–probe measurements.^{19,60,69} This is presented in Fig. 9 for 10, 77, and 298 K. In calculations of the absorption with the CT state, we neglected the couplings and dipole strength of the involved pigments and accounted for the electrochromic shifts. The effective dielectric constant $\mathit{\epsilon}r$ was assumed to be equal to 2.5, as before. We obtain good agreement between experiment and calculations for all temperatures. It is important to note that agreement was obtained only by assuming an angle of $15\xb0$ between *N*_{B} and *N*_{D} axes towards the *N*_{C} atom for the static dipole moments of both Chls and Pheos.

## V. DISCUSSION

### A. Site energies

We obtained our site energies from the fit of nine linear spectra and then verified them by calculating more independent spectra. In this subsection, we will compare our results to other work. We have to note, however, that direct comparison of the site energies can sometimes be misleading, due to the type of spectral density employed and the amount of reorganization energy included in different models. Therefore, comparing the site energies relatively, e.g., between the corresponding pigments of different branches of the RC, is a better approach.

Our model finds that the lowest site energy is of the Chl_{D1} pigment, in agreement with many previous studies.^{16,18,29,30} This is consistent with the triplet minus singlet absorption difference spectrum calculations, which can be well described assuming that at low temperatures, the triplet is localized at the Chl_{D1} pigment. On the other hand, some studies assign a considerably higher energy for this pigment.^{70,71} Our fitting shows that the latter assignment is not consistent with the experimental data.

Another somewhat controversial issue is the assignment of the Pheo_{D1/D2} energies. We find that Pheo_{D2} has the second lowest energy, while that of Pheo_{D1} is considerably higher. This is in general agreement with the findings of Ref. 30, although contrary to other reports in the literature.^{72} Since we can accurately reproduce not only absorption difference spectra RC-RC1x/RC2x but also absorption difference corresponding to reduced Pheo_{D1} and the CT state $PD1+PheoD1\u2212$, we hold that our assignment accurately reflects the system.

In previous work, the energies of the special pair pigments have mostly been found to be the highest of the pigments in the RC.^{16,18,20,26,29} Our findings suggest high, but not necessarily the highest energies for these pigments, with Chl_{D2} having a similar energy to P_{D1} and Chlz_{D2} to P_{D2}. Due to a relatively large resonance coupling between the special pair pigments, it is rather difficult to pinpoint their energies from a specific experiment. Another issue to consider is that due to very close proximity between the special pair pigments,^{3} their orbitals are likely overlapping. Thus, even the molecular excited states of P_{D1} and P_{D2} pigments might actually have partially delocalized or even CT character. Their excited state energies are likely to be only effective. Nonetheless, having simulated the absorption difference spectra corresponding to the $PD1+PheoD1\u2212$ CT state and obtained good agreement with experimental data, our energy for the P_{D1} pigment should be accurate within the limits of our model.

Regarding the peripheral chlorophylls, our results are in general agreement with other studies.^{16,18,29,30} Nonetheless, previous models often suggested very similar energies between the two Chlz molecules, while in the present work, we found a somewhat larger energy for Chlz_{D2}, which suggests that it is the pigment lost in the RC5 preparations.^{66,67} This is in agreement with the assignment from Ref. 30.

A convenient way to compare the site energies is to present the excited state energy distributions that show the energies of the eigenstates in which each pigment or CT state participates. These distributions are defined as

Here $\lambda ee,ee=\u2211\xe3\psi e,\xe34\lambda \xe3$ is the eigenstate reorganization energy. For our model, these distributions are presented in Fig. 10. In addition, we also present distributions from our previous model^{29} and from Ref. 30, where the energies of the sites were obtained by taking into account electrostatic interactions as calculated from the structural data.

It can be seen that there is good agreement between our current model and that of Ref. 30. We note that this agreement is because Eq. (55) accounts for differences in reorganization energy, thus even though our model has higher state energies, it also has larger reorganization energies and the combined result is very similar to Ref. 30. The only noticeable difference is in the special pigment energy distributions that are somewhat blueshifted in our model.

On the other hand, the differences between the excited state distributions of the present model and our previous work^{29} are larger, though the overall pattern is somewhat similar. The largest differences are for Chl_{D2}, Pheo_{D1}, and both Chlzs. One main difference is that the energy distributions for our previous work show more structure. This is because of the differences in disorder values between both models. Our previous work,^{29} following previous models based on the modified Redfield theory,^{18,20,26} assumed about three times smaller disorder than in the current model. Taking into account that presently we use a more advanced lineshape theory and a spectral density suitable for low temperature calculations and that the disorder value that we obtain is very similar to Refs. 16, 27, and 30, we believe that our current values are more accurate and that disorder in the previous studies^{18,20,26,28,29} might have been severely underestimated. These differences in disorder are expected to have significant influence on the time-resolved spectroscopic data and will have to be investigated more closely in the future.

### B. Stark spectrum and CT states

In this work, we extended the ctR theory to describe the Stark spectrum. This was done using the perturbative expansion in terms of the strength of the static electric field. Our expressions thus significantly improve upon Ref. 18, where the Markovian approximation and the modified Redfield lineshapes were used. Even though here we presented only the term corresponding to the RWA, other terms can be easily formulated. Nonetheless, our approach, though straightforward to compute, is computationally expensive, involving a double time integral over a triple sum. This means that it scales rather badly for larger systems. Additionally, the Stark spectrum is sensitive to the CT states and our results indicate that they possess very large disorder. Thus, Stark spectrum calculations need a large number of realizations to converge. We performed our fits with 25 000 realizations, and the obtained spectra were still somewhat noisy. The Stark spectra presented in this paper (Fig. 3) were calculated using 600 000 realizations. For comparison, 5000 realizations are sufficient to achieve good convergence for the absorption spectrum at 77 K. Although this suggests that our determined CT state parameters might possess some uncertainty, they still should be reasonably accurate.

Several conclusions can be made from our calculations of the Stark spectrum of the PSII RC. First, our results show that CT states influence the Stark spectrum more than was obtained previously.^{18} Without the CT states, we cannot obtain a reasonable agreement with experimental data. This implies that the coupling between the ME states and the CT states is somewhat larger than was estimated before.^{18} If so, this should affect calculations of the time-resolved spectroscopic data.

Second, our results indicate that *if* the Stark spectrum is mostly influenced by a single CT state, it is probably $ChlD1+PheoD1\u2212$ rather than $PD2+PD1\u2212$, contrary to what was previously concluded.^{18} Although the agreement between theory and experiment is poorer with the $PD2+PD1\u2212$ state than with $ChlD1+PheoD1\u2212$, the inclusion of the $PD2+PD1\u2212$ state does provide better agreement than when CT states are neglected entirely. The agreement that we obtain with the single $ChlD1+PheoD1\u2212$ state [Fig. 3(a)] is good but not excellent. It is likely that instead of a single state, more CT states are coupled to the ME states strongly enough to influence the Stark spectrum. We did not try to fit the experimental Stark spectrum with 2 CT states because the fit would have been very poorly constrained (fitting a single spectrum with eight free parameters). Instead we chose to make separate optimizations with the individual $ChlD1+PheoD1\u2212$ and $PD2+PD1\u2212$ states. Using the parameters obtained from the individual optimizations, if we include both CT states [see Fig. 3(c)] we obtain improved agreement of the $\u223c680\u2009nm$ feature [if compared to Fig. 3(a)] but at a cost of an appearance of a spurious negative feature at $\u223c700\u2009nm$.

Third, we have obtained very large disorder of the CT state ($\sigma ChlD1+PheoD1\u2212=550\u2009cm\u22121$), which is clearly visible in its energy distribution in Fig. 10. The obtained value of $\sigma ChlD1+PheoD1\u2212$ is about seven times larger than that of the ME states. This should be related to the fact that CT states have large static dipole moments, thus they are more susceptible to the changes in the environment. Also we note that our value of disorder is considerably larger than previous estimates.^{18,20,26,29} This could be due to the shortcomings of the modified Redfield theory, which was used in those studies. This result is in line with previous suggestions that disorder might control different charge separation pathways in the PSII RC.

Overall, we would like to suggest that calculations of the Stark spectrum of the PSII RC should be made not to extract the values of the system parameters (due to the computational expense) but to constrain them. Thus, models which include a number of CT states needed to calculate the charge separation dynamics^{20,29} should be tested against the experimental Stark spectrum using the level of theory developed in this paper. Since each CT state necessarily adds a number of parameters to the model, their values could be severely constrained by their sensitivity to the Stark spectrum.

We would also like to draw attention to our calculations of the absorption difference spectra corresponding to the CT state $PD1+PheoD1\u2212$ at different temperatures (Fig. 9). The good agreement between theory and experiment implies that even at room temperature the relaxed pump-probe spectra correspond to this specific state. Therefore, the energy gap between this state and all the other CT states should be relatively large, a few times more than the thermal energy at room temperature. This information should be useful in simulating time-resolved spectroscopic data.

### C. Fluorescence

In previous theoretical studies on the PSII RC, fitting of relaxed fluorescence was employed to obtain the model parameters.^{16,18,26} In our view, this approach is problematic and could lead to erroneous results. Since a normal RC performs charge separation with high efficiency, it should barely fluoresce. Thus, the majority of the fluorescence may originate from those reaction centers that cannot perform charge separation, due to their particular realization of static disorder that provides an unfavorable energy landscape. Fits of the relaxed fluorescence would then depend not on the normal RC, but only on a specific subensemble. A better approach to calculations of the fluorescence of the RC may be to simulate the time-resolved fluorescence and then integrate over all times. Unfortunately this requires including all the relevant CT states in the model and thus is out of scope of the present study.

### D. Adequacy of the theoretical approach

In the present work, we used the ctR theory to describe the optical lineshapes of the PSII RC. We previously tested the ctR theory on a molecular dimer^{35} and found that it provided excellent agreement with an exact approach. In addition, it was found to be considerably more accurate than the modified Redfield approach, which has often been used for calculations of the PSII RC optical spectra.^{18,20,26,28,29} The ctR theory was also somewhat more accurate than the Markovian Redfield approach (denoted by us as cR theory in Ref. 35), which was also applied to the PSII RC.^{16,27,30,73} As the cR and the ctR theories give reasonably close results, it is not surprising that our present work is in general agreement with the latest work where the cR theory was used.^{30}

We note that the ctR theory is formulated using the secular approximation. In usual circumstances, non-secular terms are expected to make small contributions to spectral lineshapes.^{35,44} On the other hand, when reorganization energies are very large^{74} or differing between the states, like between ME and CT states,^{75} the effects of non-secular terms might increase considerably. On the other hand, they should be strongly reduced, when mixing between the ME and CT states is relatively small. This is the case for our model, due to the large energy gap between the ME and CT states and huge disorder of the latter. Thus we expect that non-secular terms would play only a minor role in our model. Another issue is that the inclusion of non-secular terms is non-trivial and would increase the computational costs substantially, especially for the Stark spectrum.

Here we assumed the same form of the spectral density for CT as for ME states, with the only difference being a multiplicative factor for the low frequency term. Even though the actual shape of the spectral density of the CT state might be considerably different, there is no way to estimate it. Thus, our assumption is the simplest one currently possible. Note that we used the same spectral density for simulations of data ranging from 4 K to room temperature. This is based on the harmonic bath approximation. We believe it is valid, as our simulations match experimental data at all temperatures. Additionally, we note that the degree of anharmonicity was found to be small for a photosynthetic Fenna-Matthews-Olson complex.^{76}

We also assumed independent fluctuations for all states. Naturally, fluctuations of ME and CT states sharing the same molecules may be correlated. Nonetheless, currently there seems to be no reliable way of estimating these correlations for the PSII RC. On the other hand, parametrization in terms of correlation coefficients would increase the number of free parameters in our model. Thus, we chose not to do this, for simplicity.

Recently a focus has shifted to the question and relevance of quantum coherence, mostly related to the vibrational degrees of freedom, in the PSII RC. This interest was stimulated by simultaneous publications demonstrating coherent oscillations in two-dimensional electronic spectra.^{77,78} The possibility of coherence enhanced charge separation was either inferred^{78} or demonstrated in a proof of concept manner.^{77} To describe the influence of vibrations, one must either include them explicitly in the Hamiltonian^{79} or in the spectral density. These descriptions should be equivalent if accurate theories are used. The former approach, however, better captures the influence of vibrations if perturbative approaches are employed. However, the computational complexity then increases significantly. In the present work, we include several vibrational modes in the spectral density. The ctR approach allows us to qualitatively capture their influence,^{35} which is not the case for the cR or modified Redfield approaches. Of course, an exact approach, like the HEOM theory^{80} or stochastic calculations of path integrals,^{74} should be used if possible. Nonetheless, we believe that the ctR theory provides a very good balance between accuracy and computational cost. The latter is very important to enable extensive fitting of experimental spectra, as was done in this work.

## VI. CONCLUSIONS

We have revisited the spectroscopic properties of the PSII RC using the ctR theory. We constructed a model based on the tight-binding Hamiltonian, including a single CT state. The parameters of the model were obtained by fitting nine linear spectra and then further tested by calculations of other independent spectra. Our simulations are in good to excellent agreement with all the calculated spectra. The parameter values of our model differ from previous models based on the modified Redfield theory and are mostly in agreement with recent structure-based estimations that accounted for electrostatic effects. Additionally, we extended the ctR theory to calculate the Stark spectrum and applied it to the PSII RC. Our results showed that no agreement with experimental data could be obtained without an inclusion of CT states. We concluded that $ChlD1+PheoD1\u2212$ is the primary CT state, but the Stark spectrum is probably also influenced by other states. Finally, we suggest that calculations of the Stark spectrum should be used to constrain the CT state parameters obtained by fitting the time-resolved data. Our future work will be directed towards simulating the 2D spectrum of the PSII RC using the ctR theory.

## ACKNOWLEDGMENTS

This work was supported by the Research Council of Lithuania (LMT Grant No. MIP-080/2015). J.P.O. acknowledges support from the Office of Basic Energy Sciences, the US Department of Energy under Grant No. DE-SC0016384. Computations were performed using the resources of the High Performance Computing Center “HPC Sauletekis” at Faculty of Physics, Vilnius University.