Single Molecule Spectroscopy of Monomeric LHCII: Experiment and Theory

We derive approximate equations of motion for excited state dynamics of a multilevel open quantum system weakly interacting with light to describe fluorescence-detected single molecule spectra. Based on the Frenkel exciton theory, we construct a model for the chlorophyll part of the LHCII complex of higher plants and its interaction with previously proposed excitation quencher in the form of the lutein molecule Lut 1. The resulting description is valid over a broad range of timescales relevant for single molecule spectroscopy, i.e. from ps to minutes. Validity of these equations is demonstrated by comparing simulations of ensemble and single-molecule spectra of monomeric LHCII with experiments. Using a conformational change of the LHCII protein as a switching mechanism, the intensity and spectral time traces of individual LHCII complexes are simulated, and the experimental statistical distributions are reproduced. Based on our model, it is shown that with reasonable assumptions about its interaction with chlorophylls, Lut 1 can act as an efficient fluorescence quencher in LHCII.


I. EXCITON MODEL
The interaction of the pigments, i. e. the system, with the vibrations, i. e. its environment, is treated by second order perturbation theory. There the bath is completely described by its correlation function or, equivalently, in the frequency domain, by the spectral density of the bath vibrations. We use the spectral density from Ref. [1]. It is constructed from one overdamped oscillator representing the slow protein vibrations and 48 high-frequency intrapigment modes (for more context for the spectral density and excitonic model see Ref. [2]): (1) The (temperature dependent) correlation function given by this spectral density is assumed to be uncorrelated between individual sites and differs only between Chl a and Chl b, while the difference is in the coupling strength ν n = ν a/b : The time-dependent correlation function is then obtained by Fourier transform of (2). We use the same values of the spectral density parameters as in Ref. [1]. The difference between the vertical, Franck-Condon transition of the pigments, which are called site energies in this text, and their 0-0 transitions is given by the reorganization energy due to the interaction with the bath: Because the pigments are strongly coupled, the preferred basis of calculations is the excitonic basis in which the system Hamiltonian is diagonal with eigenvalues, exciton energies, ω i0 .
All quantities, including correlation functions C(t), reorganization energies λ and transition dipole moments µ have to be transformed into the excitonic basis: The position of the zero phonon lines of the excitonic transitions is then The spectral lineshapes are calculated by 2nd order cumulant expansion employing the socalled lineshape functions: The lineshape function is conveniently expressed in terms of the spectral densities (1 − cos(ωt)) + i (sin(ωt) − ωt) .
The absorption spectrum is calculated as where the absorption lineshape is Here Γ i is the population relaxation rate from state i. The fluorescence is similarly given as where P i is the steady-state population of state i and the fluorescence lineshape is The population transfer rates in the Redfield theory are obtained as The population relaxation rates of chlorophylls are a result of energy transfer and radiative and non-radiative decay, i.e. Γ i = j k ji +Γ i , whereΓ i = n |c n i | 2Γ n ,Γ n = 1/τ n is inverse lifetime τ n of site n excited state, which is taken to be 3 ns for all n.

II. BULK SPECTRA MEASUREMENT
To check for sample degradation a control bulk measurement of absorption and fluorescence was performed. In Fig. 1 the measured spectra are given together with the spectra of LHCII trimers taken from [3]. Their excellent agreement confirms absence of degradation in course of the experiment and also justifies the usage of the trimer spectra for the bulk spectra modeling in the main text.

III. RANDOM WALK MODEL
In the original description by Valkunas [4] and, in more detail, in [5], the protein diffusive motion was described as a continuous-time random walk (CTRW) on a two-dimensional potential energy surface. Here we somewhat simplify this description in the following way.
The two generalized coordinates represent fast and slow degrees of freedom and thus, when following the slow dynamics, the fast fluctuations can be adiabatically eliminated. It can be shown that when we consider the potential dependence on the fast coordinate the same in both the on and off states for simplicity, i.e. setting x 0 = 0, λ 2 /λ 1 = 1 in [4], the fast coordinate can be completely removed. This leaves us with effectively one-dimensional problem.
After transforming into dimensionless coordinates y γ 1 k B T → y, the parabolically approximated minima of the potential for the on (1) and off (2) state are Here γ 1,2 determine the steepness of the potential and the second minimum is shifted by y 0 along the slow coordinate and by U 0 along the potential energy.
The protein then performs random walk on the respective surface, with probabilities to tunnel to the other surface k i is a rate of falling from the i−th potential, the first exponential term reflects the energy gap law and the min term ensures detailed balance condition. The coefficient α can be treated as a constant and ω 0 is a characteristic frequency of the protein environment vibrations responsible for the tunneling.
The protein diffusion under this conditions can be described either by the CTRW or by a discrete random walk (RW). The former approach was employed by Valkunas et al. [4,5].
However, we believe that for our purpose it is better to solve this problem as a discrete RW for two reasons. First, the coupled Smoluchovski equations for the CTRW on the two potentials can be decoupled only for conditional probabilities, i.e. assuming that the system was in the opposite state in the previous interval, and, in the same time, employing the same, equilibrium initial condition for each dwell time. When continuously modeling the trajectory of a single protein, we do not have to include the resetting after switching and also the conditioning will be inherent, as the system is observed being in a particular state. And second, when we want to simulate the individual intensity time traces, it is more natural to really follow the trajectories of the individual proteins on their PES.
If we want to follow the time trace of every single molecule, we should follow its particular trajectory. The blinking statistics will then be recovered by averaging over a large number of molecules, exactly as in the experiment. To the purpose of following trajectories of individual proteins, we need to describe its discrete RW (DRW) in the potential. We will denote probability of going right (left) as p (q). In the symmetrical RW we have p = q = 1 2 . If the protein is at a coordinate y, in the next step it will move with probability p to y +a and with probability q to y − a, where a is the length of the step. Inspired by classical approach by van Kampen [6], we augment the position dependent probabilities in the presence of the potential U(y) as We note that defined probabilities defined in Eq. (17) reflect the detailed balance condition . Considering a small step a, we can use Taylor expansion in y Now, considering that p(y) + q(y) = 1, we get Using the potential form (15), we get for the probabilities The length of the step a on respective surface can be related to the transformed diffusion coefficient D 1,2 : where ∆t is the time duration of the step. During ∆t the protein walks either to the left or right with the respective probability and with the probability κ i→j (y)∆t switches to the other surface.
Similarly to [4] the position near the potential intersection can be chosen as an initial condition. However, as we do not include resetting after switching in our approach, this parameter Valkunas [4, pH 8] This work  [4]. Considering the differences -LHCII monomers vs trimers, discrete vs continuous RW, no resetting vs resetting -the agreement is satisfactory.

We start with Liouville-von Neumann equation (LvNE) with a Hamiltonian
where is the matter Hamiltonian including part, H S , describing the relevant system of interest (electronic degrees of freedom (DOF)), thermodynamic bath H B , in which the system is embedded and system-bath interaction term H I . Interaction between the light and matter is described by semiclassical light-matter Hamiltonian in dipole approximation. We denote the full density matrix of the system by W (t). The LvNE reads as Our aim is to discuss approximations necessary to make a transition towards equations of motion for the so-called reduced density matrix (RDM) and ultimately to the exciton populations where |a are eigenstates of the Hamiltonian H S .

A. Expansion in Orders of Electric Field
In order to make approximations only with respect to the system-bath interaction Hamiltonian, we first switch to the interaction picture with respect to the field free evolution of the system. We define Here, we used the field free evolution operator which propagates the density matrix from its initial condition at time t 0 to time t. Eq. (28) will be formally integrated, and we return to Schrödinger picture. This yields We obtained a recursive formula for the total density matrix W (t). By iterative reinsertion of Eq. (30) into itself, we obtain a series of terms in orders of the electric field E(t). We get where and for n > 0.
For any order of perturbation, we can write a differential equation instead of the integral expressions, Eq. (33). By differentiating Eq. (33) we get The second term on the right hand side of Eq. (34) is a source term which is a known function of time, if the electric field is known and the lower order density matrix is known. does not lead directly to a closed set of equations for the RDM ρ(t). Rather, one has to derive a master equation which includes terms describing phenomena emerging due to the reduction from a full set of DOF to a defined subset of the system's DOF.
We will formally separate the matter Hamiltonian H M into the relevant electronic DOF described by Hamiltonian H S , bath DOF described by the Hamiltonian H B and some interaction Hamiltonian H I . We define We switch into interaction picture with respect to the interaction free part of the matter Hamiltonian by definingW were This leads to a new equation of motion in a form where we defined a Liouville superoperator L I (t) corresponding to the interaction Hamiltonian in interaction picture (denoted by H I (t)): Here, A is an arbitrary operator, and we defined a source term in interaction picturē with We will now proceed along the lines of the standard projection operator technique (see e.g. [7,8]) by defining two projection superoperators P and Q (acting on the operators W (n) (t)) such that In particular we will use the well-known Argyres-Kelley projection superoperator which acts where w eq is the bath equilibrium density matrix. Applying these operators to Eq. (38) and eliminating QW (n) (t) we obtain a variant of the well-known Nakajima-Zwanzig identity: Here, The first term on the right hand side (r.h.s.) in Eq. (44) can be eliminated by an appropriate choice of the form of interaction operator H I . We assume that the system-bath interaction operator is of a pure dephasing form in the basis of the electronic states localized on the molecules of the complex (the site basis) We can always redefine this interaction operator in such a way that The reorganization energies can be always made part of the system Hamiltonian. We assume further on that Eq. (47) holds.
The second term on the r. h. s. of Eq. (44) is a source term by which the operator of n th order of perturbation is connected to the lower order operator. The third term is the so-called initial correlation term. Since the initial condition is always W (n) (t 0 ) = 0 (and t 0 is a time before any interaction with the field), this term also disappears. The fourth term has a very similar structure as the initial correlation term. It underlines the fact that the true initial condition for the W (n) is set by the source term, and, in the case that some correlation between the bath and the system is created during the process of excitation, it contributes to the excited state evolution.
The last term of Eq. (44) is responsible for the processes of energy relaxation and transfer and for dephasing of coherences. We will treat this term in line with Ref. [7] to turn it into a second order time local differential equation. We define formally a superoperator U(t, t 0 ) which describes an exact evolution of the projected statistical operator PW (n) (t).
In Eq. (44) we substitute PW (n) (τ ) =Ū(τ, t)PW (n) (t). Now PW (n) (t) can be taken out of the integral and the differential equation is local in time. Next we will approximate the integral by a second order expression in terms of L I (t). BothŪ(τ, t) and U Q (t, τ ) contain all orders of L I (t), however, in order to keep the integral in second order, we only need to use their zero order approximationsŪ(τ, t) ≈ 1 and U Q (t, τ ) ≈ 1 . Using the second order approximation in the last term of Eq. (44) we get a precursor to a derivation of the standard Redfield relaxation equations. We notice that all terms start with P and correspondingly, we perform the trace in all terms and divide out the equilibrium bath density matrix w eq . The result reads as where we defined a reduced source terms can be turned back to Schrödinger picture in a straightforward way. In the following section we will apply Eq. (48) to evaluate ρ (1) (t) and ρ (2) (t).

C. First Order in Field
The first order reduced density matrix is driven by the first order source term s (1) which contains the zero order density matrix where W (0) (t) = w eq |g g|.
Here we already returned to Schrödinger picture. The equilibrium density matrix of the system does not evolve in time, and, correspondingly, the only time evolution in the source term originates from the electric field. As a consequence one can easily perform the trace to obtain s (1) (t) in a form where h.c. stands for Hermite conjugate, and µ kg = k|µ|g is the transition dipole moment between the ground state |g and an excited state |k . One can easily show that QS (1) (t) = 0, and thus the first order equations of motion reads as where is a relaxation tensor of a second order in system-bath interaction term.
For an infinite bath with a unstructured spectral density coupled weakly We could always choose t 0 sufficiently small to make the relaxation tensor constant. For pulses longer than the bath correlation time, even if we set the start of the dynamics, t 0 , before the center of the pulse so that E(t 0 ) < ǫ where ǫ is some chosen small value, the relaxation tensor will integrate into constant already during the action of the pulse. This is a consequence of the decoupling of the bath from the influence by external field which is otherwise indirectly mediated by the electronic DOF. This decoupling is a consequence of reducing our equations of motion to electronic DOF only. The only point in which the bath state is influenced in our description is the initial condition. Eq. (54) can therefore describe the dynamics correctly only if the excitation field is in a form of the Dirac delta function centered at some τ ′ . Than the start of the dynamics can be set to t 0 = τ ′ , and the bath can correctly react to the excitation of the system. One can demonstrate this idea by using Eq.
(33) to obtain ρ (1) (t) -this procedure involves no additional approximation. We obtain where is an operator on the same Hilbert space as ρ (1) (t). Eq. (57) can be written as where ρ 0 = |g g|, and U(t, τ ) is a formally exact evolution superoperator for the reduced density matrix. Eq. (56) corresponds to a sum (integral) of field free evolutions starting from an initial condition [µ, ρ 0 ] − at times τ weighted by the field amplitude E(τ ). The exact field free evolution has to involve an exact time dependent relaxation tensor R(t, τ ).
Notice that the second argument coincides with the start of the dynamics for a given field free propagation.
The discussion above underlines the fact that the RDM master equation, Eq. (54), can only be correct if the source term is of a Dirac delta form, or if the transient features of the relaxation tensor are unimportant. It turns out that for the first order density matrix ρ (1) (t) and its subsequent application to calculation of the second order RDM ρ (2) , the early dynamics is important. If we used Eq. (54) with constant relaxation (dephasing) tensor, we would limit the lineshape of the transitions to Lorentzian type. For more general and realistic lineshapes we need to use Eq. (56).

D. Secular Approximation in First Order Equations
For the first order term in the field, the transient features of the relaxation tensor turn out to be of crucial importance. They are responsible for the line shape of linear absorption.
The absorption spectrum of a multilevel molecular system can be calculated as (see e.g. [9]) Here the initial condition µρ 0 which consists entirely of optical coherences is propagated in time. In general the evolutions of different optical coherences are coupled. Nevertheless, the so-called secular approximation in which decoupling of individual coherences is postulated is very often valid. In particular, it can be shown that secular approximation is strictly valid in second order for a homodimer [10]. For molecules with similar reorganization energies, non-secular effects are very small, and only for systems were the participating states have very different reorganization energies one can find pronounced effects of coupled evolution of different optical coherences [11]. Because our interest in this work is with photosynthetic systems composed of pigments of similar kind (chlorophylls and bacteriochlorophylls) we can expect the secular approximation to work reasonably well. We can therefore assume that only the elements J ag (t, τ ) = a|J(t, τ )|g = U agag (t, τ )µ ag (60) of the operator, Eq. (57) are non-zero. In secular approximation, Eq. (54) splits into a set of independent equations for coherences. Green's function corresponding to the secular version of Eq. (54) gives a second order approximation to U agag (t, τ ). It is possible to show that for a bath composed of harmonic oscillators, second order master equation for an optical coherence gives an exact result [12]. It can be further shown that and representing a lineshape function of the transition |g → |a . For details see e.g. [8].
To conclude this section, in secular approximation we can obtain the reduced density matrix elements of the first order in field to all orders of the perturbation in the systembath coupling. This is a very good starting point for the derivation of the source term for the second order.

E. Photo-induced Initial Correlation Term in Second Order Equations
Second order field perturbation of the density matrix enables us to calculate populations of excited states relevant for SMS. In the experiment, we observe fluorescence from excited states populated due to relaxation dynamics subsequent to the excitation by external field.
According to Eq. (41), the second order source term has a form This result is general, i.e. independent of our previous discussion of secular approximation.
We have discussed the validity of secular approximation for first order equations in preceding section. The validity of the secular approximation in the second order term, however, has to be discussed separately. In the excited state band, non-secular terms can fine tune the final, quasi-equilibrium (electronically excited) state of the system reached after the propagation (see e.g. [9]). If secular approximation is not valid, the coherences do not decay to zero.
If, however, the working basis is chosen well for the particular problem, e.g. delocalized excitonic basis is chosen in case of weak system-bath interaction, or localized basis is chosen for strong system-bath coupling, the need for non-secular theory can be minimized. Most importantly, large body of theoretical understanding of photosynthesis is based on the secular approximation. We will therefore apply secular approximation in second order equations to simplify our treatment.
Unlike in the case of the first order RDM, second order contribution QS (2) (τ ), which enters the photo-induced initial correlation term (see Eq. (50)), is in general non-zero.
Before we derive the second order reduced equations of motion in Section IV F of this SI, we will treat below this additional term.
We are interested in a contribution of the photo-induced initial correlation term to populations of excited states, and we will therefore express this contribution in the basis of the eigenstates of the electronic Hamiltonian. For populations we can keep working in interac- where we made the same approximation as in the relaxation tensor, namely, that U Q (t, τ ) ≈ 1. The term PS (n) (τ ) which is a part of the QS (n) (τ ) = (1 − P)S (n) (τ ) term is by definition of P proportional to w eq and tr bath H I ab (t)w eq = 0 because of Eq. (47). The Q projector is therefore effectively equal to unity, here. Using the system-bath interaction Hamiltonian elements in eigenstate basis ∆V ab (t, t 0 ) = n a|n n|b ∆V n (t, t 0 ), we can write the pulse correlation term as A straightforward but tedious derivation using second order cumulant expansion and rotating wave approximation leads to an expression in terms of bath correlation functions C ab,cd (t) = 1 2 tr bath {∆V ab (t)∆V cd (0)w eq } , and the corresponding line shape functions The expression reads as where It is easy to show that K aa (t, τ ) = 0 and therefore all contributions to the population from the pulse induced correlation term are of the non-secular form. The light excites coherence between states |a and |b and the term of the first order in system bath interaction (∆V ab (t)) is responsible for converting this contribution to a population of state |a . All terms K ab for a = b oscillate. Therefore, we can ignore the contribution ofĪ aa to the populations in our work for the same reason we ignore other non-secular contributions. For slow field envelops E(t) and resonant pulses (the case in SMS) the termĪ aa (t, t 0 ) contains an integral over oscillating functions and it is itself oscillating. Contributions from different times therefore cancel and the overall contribution to populations is small.
Besides the smallness of the contribution due to oscillation, one can find situations wherē I aa yields zero for other reasons. Because K ab is proportional to correlation functions we can see that K ab is zero if the delocalized states |a and |b do not share the same sites. For instance, when the excitonic basis coincides with the site basis (vanishing coupling between sites) the termĪ aa is automatically zero. Another case yielding zero is a case of homodimer.
There we can assume that g aa,aa (t) and g bb,bb (t) are the same, and we can also assume that C ab,bb (t) = C(t) n a|n n|b | b|n | 2 = ξ ab,bb C(t), where the bath correlation function C(t) ≡ C n (t) is the same for all (both) sites. The second equals sign in Eq. (73) defines the coefficients ξ ab,bb . The term K ab then consists of two terms proportional to (ξ ab,aa − ξ ab,bb ) and (ξ ab,aa + ξ ab,aa ), respectively. These terms are both zero.

F. Second Order Equations of Motion
Similar to Section IV C, we can use Eqs. Nevertheless, because populations change on a much slower time scale than coherences, and because we are in fact only interested in populations at long times after excitation, a constant relaxation tensor obtained by the limit R eff (t) = lim t 0 →−∞ R(t, t 0 ) is expected to be a reasonable approximation for the description of our SMS experiment. The effective relaxation tensor does not depend on the initial condition time, but it still weakly depends on the running time of the experiment (as discussed in the main text). Neglecting the photoinduced initial correlation term, we obtain the following equations of motion for the second order RDM: Applying secular approximation, writing out explicit form of the source term s (2) (t), and selecting only equations for eigenstate populations, we arrive at Eq. (8) of the main text.