Primary Charge Separation in the Photosystem II Reaction Center Revealed by a Global Analysis of the Two-dimensional Electronic Spectra

The transfer of electronic charge in the reaction center of Photosystem II is one of the key building blocks of the conversion of sunlight energy into chemical energy within the cascade of the photosynthetic reactions. Since the charge transfer dynamics is mixed with the energy transfer dynamics, an effective tool for the direct resolution of charge separation in the reaction center is still missing. Here, we use experimental two-dimensional optical photon echo spectroscopy in combination with the theoretical calculation to resolve its signature. A global fitting analysis allows us to clearly and directly identify a decay pathway associated to the primary charge separation. In particular, it can be distinguished from regular energy transfer and occurs on a time scale of 1.5 ps under ambient conditions. This technique provides a general tool to identify charge separation signatures from the energy transport in two-dimensional optical spectroscopy.


I. MODEL HAMILTONIAN
The molecular structure of the photosystem II (PSII) reaction center (RC) has been resolved by X-ray crystallography and is well known [1][2][3]. It has eight cofactors embeded in a protein matrix. For the modeling, we describe the optical transition in each cofactor in terms of individual quantum two-level systems. We use the matrix elements of the electronic coupling between the pigments and the initial values of the site energies of Ref. [4]. We add four different charge separated (CS) states, denoted by P + D1 P − D2 , Chl + D1 P heo − D1 , P + D1 Chl − D1 and P + D1 P heo − D1 . The coupling between the CS states were resolved by transient absorption spectroscopy [5,6]. We assume the electronic coupling to be constant and further refine the site energies of the Frenkel exciton (FE) and the CS states by fitting the linear absorption and circular dichroism spectrum at different temperatures to the experimentally determined results. The final Hamiltonian matrix in the optimized form and the single-excited state manifold reads All entries are given in units of cm −1 and have been shifted by 14700 cm −1 by the rotating-wave approximation. The associated index of Hamiltonian matrix are labeled as: P D1 , P D2 , Chl D1 , Chl D2 , P heo D1 , P heo D2 , Chlz D1 , Chlz D2 , Here, ǫ m * n * is the site energy of the doubly excited state formed by the FE-FE states. ǫ m * n + k − and ǫ m + n − k + l − are the site energies associated to the FE-charge separated states and the CS-CS states, respectively. J m * n * ,m ′ * n ′ * is the electronic coupling of the FE-FE states between the states m * n * and m ′ * n ′ * . Moreover, J m + n − k + l − ,m ′ * n ′ * is the coupling between the CS-CS state m + n − k + l − and the FE-FE state m ′ * n ′ * . In addition, J m * n + k − ,m ′ * n ′+ k ′− is the coupling between the FE-CS states m * n + k − and m ′ * n ′+ k ′− , and J m + n − k + l − ,m ′ * n ′+ k ′− is the coupling between the CS-CS state m + n − k + l − and the FE-CS m ′ * n ′+ k ′− . Finally, J m + n − k + l − ,m ′+ n ′− k ′+ l ′− is the coupling between the CS-CS states m + n − k + l − and m ′+ n ′− k ′+ l ′− . In this calculation, we assume the K mm,nn = 0, K mm,kn = 500 cm −1 and K nm,jk = 2000 cm −1 .

A. Time non-local quantum master equation
For the numerical simulations of the linear and 2D spectra, we have applied the time non-local (TNL) quantum master equation [7,8]. The time evolution of the total density matrix ρ(t) at time t which includes the system and the bath is governed by the Liouville-von Neumann equation with the Liouville superoperator L, according to (h = 1) The total Hamiltonian H tot = H s + H b + H sb + H ren includes the system, the bath, and the interaction and renormalization terms. For the system with a single degree of freedom x and the bath consisting of an ensemble of harmonic oscillators H b = N j=1 p 2 j /(2m j ) + m j ω 2 j x 2 j /2 , the coupling between them is assumed to be of the form H sb = f (x) N j=1 c j x j with some real function f (·). The projection scheme of Nakajima and Zwanzig [9] separates the dynamics of the system and the bath. The bath is assumed in a thermal state represented by the canonical density operator ρ eq b = exp(−βH b ) with a given temperature T = (k B β) −1 . Applying the projector P = ρ eq b tr b with tr b ρ eq b = 1 and its orthogonal complement Q = 1 − P yields a formally exact quantum master equation for the time evolution of the reduced system density operator ρ s (t) = tr b ρ(t) [7] in the form ofρ Here, ρ tot (0) is the total density operator of the system and bath at initial time. The Liouville superoperators L s , L sb and L ren are associated with corresponding Hamiltonian operators. Moreover, L eff s · = −i[H s + H ren , ·] and T is the time ordering operator [10]. Next, we expand the correlated thermal equilibrium state up to the first order in the system-bath coupling and obtain with the respective partition functions Z tot = tr exp(−βH tot ), Z b = tr b exp(−βH b ), and Z s = tr s exp(−βH s ). Next, we take the trace over the system degrees of freedom on both sides of Eq. (5) and get Here, χ = (1/Z s )tr s f (x)e −βHs with the coupling function f (x) defined below Eq. (3). As well established, the bath and coupling parameters can be collected in the spectral density This enters in the bath correlation function with the real part a(t) and imaginary part b(t). After inserting Eqs. (5) and (6) into Eq. (4), we express the last three terms of Eq. (4) by a(t) and b(t) in the form The potential renormalization is given in terms of the spectral density as µ = ∞ −∞ dω 2π J(ω)/ω. In order to obtain an analytic form of the bath correlation function, any given spectral density (in our particular case, we use the standard Ohmic form) can be approximated by a sum of Lorentzian-like spectral terms [11,12] according to The spectral amplitude p k , the frequency Ω k and the width Γ k follow from the expansion of the original function in terms of the Lorentzian shapes. Inserting the expanded form of J(ω) into Eq. (8) results in with the Matsubara frequencies ν k = 2πk/β. Moreover, n ′ is the number of Matsubara frequencies used. Next, we rewrite the correlation functions as a(t) = nr k=1 α r k e γ r k t and b(t) = ni k=1 α i k e γ i k t with n i = 2n, n r = 2n+n ′ . Then, we define new auxiliary "density matrices" which incorporate both memory effects and initial correlations according to The time-retarded Eq. (4) (first term) can be then deconvoluted into a set of coupled first-order equations aṡ ρ r k (t) = (L s (t) + γ r k )ρ r k (t) + L − ρ s (t), k = 1, . . . , n r , ρ i k (t) = (L s (t) + γ i k )ρ i k (t) + L + ρ s (t), k = 1, . . . , n i .
This set of coupled time local quantum master equations was used for the calculations of the quantum dynamics and the resulting spectra. It generates an effective time-nonlocal quantum dynamics due to the coupled differential equations.

B. Linear absorption and circular dichroism spectra
We have used the first-order transition dipole moment correlation function to calculate the absorption and CD spectra defined by Here, µ is the transition dipole moment and the subscript g refers to performing the trace over a thermally equilibrated bath distribution. R m,n is the distance vector between the monomers m and n. The correlation functions can be calculated as µ(t)µ(0) g = tr s {µtr b [e −iHt µρ g e iHt ]}, and r m,n ·µ m ×µ n (t) g = tr s {R m,n ·µ m ×tr b [e −iHt µ n ρ g e iHt ]}.

C. Numerical calculations of two-dimensional spectra by highly parallelized message-passing-interface
For the calculation of two-dimensional spectra, we have applied the equation of motion-phase matching approach (EOM-PMA) established in Ref. [13]. In the EOM-PMA, the induced polarization in the direction of the photon-echo signal is calculated by the simultaneous propagation of three auxiliary density matrices ρ i=1,2,3 (t), each of which obeys a modified effective equation of motion according tȯ where V α (t, t α ) = XAe −(t−tα) 2 /(2Γ 2 ) e iωt , X is the transition dipole operator, Γ is the pulse duration, and ℜ is a relaxation superoperator. All three above master equations are here calculated by adopting the TNL method of Eq. (13) to the auxiliary density operators with the corresponding different time-dependent Hamiltonians. Then, the third-order induced polarization is obtained as where the brackets . . . indicate the evaluation of the trace. The total 2D Fourier-transformed spectrum is then given by the double Fourier transform of the photon-echo polarization signal with respect to the delay time τ = t 2 − t 1 and t according to Here, ω τ is the "coherence" frequency, ω t is the "detection frequency", and T = t 3 − t 2 is the "waiting" time.
For the concrete simulation of the 2D electronic spectra of the PSII reaction center, a message-passing interface (MPI) has been used to minimize the computation time. The coherence time window [−300 fs, 300 fs] was separated into time slices of length dτ = 10 fs. For further simplification, the inhomogeneous broadening by the static disorder was simulated with more than 100 realizations and each realization was calculated by one CPU. Thus, the total number of CPUs we need for one 2D electronic spectrum at a given waiting time is 61 × 100 = 6100 CPUs.

III. CORRELATED DISTRIBUTION GENERATED BY THE CHOLESKY DECOMPOSITION
The Cholesky decomposition is used [14] to generate cross-correlated static disorder and bath fluctuations for the PSII reaction center. To do so, we first define a correlation matrix A which determines the degree and the type of correlations between two elements. In our current model of the PSII reaction center, the correlation coefficients are given by From this, the L-matrix can be obtained by a Cholesky decomposition according to A = LL † . In order to generate spatially correlated static disorder, we define an array (∆ 1 , ∆ 2 , . . .) with the dimension being given by the number of exciton sites. Each ∆ i indicates the disorder for the ith site and is independently distributed according to a Gaussian. Similarly, to generate correlated bath fluctuations, we just need to define an array of fully independent baths, ( i x i , j x j , . . .). Then, the combined correlated static disorder and bath fluctuations are obtained as the product of the L-matrix and the two independent arrays. This approach has been successfully applied to generate correlated noise in a system-bath model [14].

IV. MULTIPLE PATHWAYS OF THE CHARGE DYNAMICS IN A TETRAMER
In this section, we report the results of the calculated 2D electronic spectra up to 12 ps for a tetramer which includes two CT states. Again, the TNL method is used. In order to describe the multiple pathways of the CT dynamics, we include two CT states with the site energies 0 cm −1 and −500 cm −1 . Then, the system Hamiltonian is given by As before, we have analyzed the calculated 2D spectra for different waiting times by the global fitting approach. Four different decay-associated spectra (DAS) can be resolved, with the associated decay time constants 39 fs, 422 fs, 3.04 ps and ∞, see Fig. 1. The fastest process with the time constant 39 fs is the fast electronic dephasing. The two other components with 422 fs and 3.04 ps reflect the evidence of the charge separation process which correspond to the CT from the excitonic states to the CS state with site energy 0 cm −1 , and, the secondary charge separation process from the site with 0 cm −1 to −500 cm −1 . This proves that the multiple pathways of the charge separation process in a tetramer can be fully revealed by the 2D electronic spectroscopy upon using the 2DDAS.

V. LONG-LIVED VIBRATIONAL COHERENCE IN THE PSII REACTION CENTER
Finally, we address the dynamics observed in our measurement of the PSII reaction center. In Fig. 2, we present the extracted time traces of the four selected peaks A, B, C and D, as marked in the 2D electronic spectrum at the waiting time T = 0fs. We clearly observe long-lived oscillatory signals which last over the full available window of the measurment time of 2 ps at the room temperature. This observation is in agreement with the previous observations reported in Refs. [15,16]. Based on our previous study [17], these small-amplitude oscillations are induced by ordinary vibrational coherence. In order to present the traces in one plot, the A-trace has been shifted downwards by -0.62. All the traces start at T = 50fs.

VI. DECAY-ASSOCIATED SPECTRUM OF 13.9 PS
We have experimentally resolved one decay component with a time scale of 13.9 ps, which can not be reproduced by the theoretical calculation, which is due to the required large, but unavailable computational resources. The resolved decay-associated spectrum is shown in Fig. 3, it only contains one diagonal peak at ω τ = ω t = 14800 cm −1 .