On the performance of a photosystem II reaction centre-based photocell

We investigate the performance of a theoretical photosystem II reaction centre-inspired photocell device through the framework of electron counting statistics. In particular we look at the effect of a structured vibrational environment on the mean current and current noise.


Supplementary Note 1: PSIIRC parameters Electronic Parameters
The state energies in wavenumbers for the photocell model are given in Supplementary Table S1. The electronic PSIIRC Hamiltonian H el used for the excited state dynamics can be constructed using the energies of the single excitation states of the six core chromophores and two charge transfer (CT) states of the Chl D1 charge separation pathway 26 along with the couplings from Supplementary Table S2.

Spectral Densities
The parameters for the low energy phonon environment, described using a Drude spectral density J D (ω) (see Eq. 1 in the main paper) are λ D = 35cm −1 and Ω D = 40cm −1 , which are the reorganisation energy and cut-off frequency respectively 21,26 . The parameters for the high energy modes used in the full PSIIRC spectral density J(ω) are given in Supplementary Table S3 21,28 for which a damping constant γ j = 10cm −1 is used for all modes 28 and the reorganisation energy associated with each mode is given by λ j = ω j s j . The frequencies, Huang-Rhys factors and damping constants for the broadened modes in spectral densities J 1 (ω) and J 2 (ω), used to approximate the 48 underdamped modes in nonperturbative calculations, are ω 1 = 342cm −1 , s 1 = 0.4, γ 1 = 100cm −1 , and ω 2 = 742cm −1 , s 2 = 0.32, γ 2 = 100cm −1 . The reorganisation energies for the broadened modes are given by λ i = ω i s i . To quantify the stronger coupling of the CT states to the bath, the reorganisation energies and line broadening functions of the primary (I) and secondary (α) CT states were scaled using values of ν I = 3 and ν α = 4 respectively 26 .

Supplementary Note 2: Theoretical Methods for Electronic Dynamics
The different parameter regimes found within PSIIRC systems make it a challenging computational task to accurately describe the exciton and charge transfer dynamics in a self-consistent manner. Two frameworks for the photocell ) A hybrid model where the exciton dynamics is described using the non-perturbative hierarchical equations of motion (HEOM). (ii) A Pauli master equation where the exciton dynamics is described using modified Redfield theory. Both these frameworks share a common description of the incoherent rates for charge transfer, optical excitation and coupling to leads which can be obtained from a Lindblad dissipator. More detail is given on the theoretical basis of these frameworks below.

Hamiltonian Dynamics
The Hamiltonian of the system and bath is H = H el + H B + H I . The electronic Hamiltonian H el = i |i i| + ij T ij (|i j| + |j i|) is given in Supplementary Tables S1 and S2. The exciton states and associated energies are given in Supplementary Table S4. Additionally, H B = kh ω k b † k b k is the bath Hamiltonian describing a continuum of bosonic modes and H I = ki g ik |i i|(b † k + b k ) defines the interaction between the system and bath. The bath is characterised by a spectral density J(ω) = k g 2 k δ(ω − ω k ) for which several forms are investigated in the main paper (see Eqs. 1 and 2). The bath is assumed identical for each site apart from a scaling for the CT states taking into account their stronger coupling to the bath. Different open quantum systems approaches, depending on the framework considered, are used to determine the dynamics of the electronic system under the influence of the bath.

Charge Transfer: Förster Theory
Primary and secondary charge separation were modeled using Förster theory 42,43 (equivalent to Marcus rates in the excited state) which is perturbative with respect to the electronic coupling between states and describes incoherent electron transfer events. The Förster rate expression has the general form where V ab is the electronic coupling between states a and b and S ab is the spectral overlap between these states. This overlap can be expressed equivalently in both the time and frequency domain Supplementary Figure S1: Comparison of perturbative theories with the hierarchical equations of motion for downhill transfer rates in a dimer system 70 . The electronic system Hamiltonian is H = (∆ /2)σ z + V σ x which is interacting with a Drude spectral density with reorganisation energy λ = 100cm −1 and cutoff frequency Ω c = 53cm −1 . The rates are calculated for different electronic coupling strengths (a) The inset of (b) shows the dependence of transfer rate on reorganisation energy λ for an energy gap ∆ = 10cm −1 . where are the absorption and fluorescence lineshapes respectively. The line broadening functions are defined g(t) = g D (t)+ j g j (t) where g D (t) and g j (t) are the line broadening functions for the Drude mode and the jth underdamped mode respectively. For a Drude . The summations over Matsubara frequencies ν k = 2π β k can usually can truncated at a suitable k depending on the temperature. For primary charge transfer we have transfer from an exciton state (X) to the primary CT state (I). In this case the electronic coupling between X and I is V XI = n∈X c n (X)V nI where c n (X) is the amplitude of site n in exciton X and V nI is the electronic coupling between site n and CT state I. The exciton reorganisation energy has the form λ X = i |c i (X)| 4 λ and similarly the line broadening function g X (t) = i |c i (X)| 4 g(t) where c i (X) are the amplitudes of exciton X at site i and λ and g(t) are the site reorganisation energy and line broadening function. The primary CT reorganisation energy and line broadening function are λ I = ν I λ and g I (t) = ν I g(t).
Secondary charge transfer rates are calculated in a similar way though the reorganisation energies and line broadening functions of both the donor and acceptor are now scaled.

Justification of Perturbative Treatment for Charge Transfer Dynamics
The parameter regime for primary charge separation in the Chl D1 pathway is such that a perturbative approach to the dynamics using Förster theory is a good approximation to the exact dynamics. In the site basis the electronic coupling between Chl D1 and Phe D1 and the primary CT state is 70cm −1 while the CT reorganisation energy ranges from 105cm −1 for a Drude spectral density J D (ω) to 1769cm −1 for the full spectral density J(ω). A comparison of transfer rates calculated using HEOM and Förster theory is given in Supplementary Fig. S1 70 , which shows downhill transfer rates for a dimer system in various parameter regimes. The relevant regime comparable to primary CT in PSIIRC is the right-hand side of plot a, where the electronic coupling is less than both the reorganisation energy and the energy gap. In this region, Förster theory and HEOM give transfer rates in close agreement. The small electronic coupling and large energy gap mean transfer must occur through transient resonant states of the system rather than through coherent evolution.

Incoherent Excitation and Coupling to Leads
For the incoherent rates of electron hopping between the source lead of the system a value of Γ L = 201.6cm −1 was used. The rate to hop from the system into the drain lead (Γ R ) was varied for the current and Fano factor calculations as a function of voltage.

Hybrid Model: Hierarchical Equations of Motion
The exciton relaxation dynamics in photosynthetic systems can be described using the non-perturbative hierarchical equations of motion (HEOM) 34,35,36 . This treats the exciton-vibrational dynamics exactly but is computationally expensive and requires the hierarchy to be truncated at some suitable level. In order to include additional incoherent processes in the photocell dynamics a hybrid model must be used. In this case a Lindblad dissipator is included to capture the coupling of the excitons to the rest of the photocell 68 . In particular the dissipator describes the excitation of the system from the ground to exciton manifold, primary and secondary charge transfer, and coupling to the leads.
Due to the numerical difficulty of describing the system-bath dynamics the full PSIIRC spectral density cannot be included in this framework. Instead, a Drude mode and a broadened underdamped mode (see Fig. 2 from the main paper) overlapping the exciton energy gaps is used to approximate the full spectral density. In the main paper the photocell current and power output for this approximate spectral density are compared with a photocell having just a Drude bath.
The infinite hierarchy of coupled differential equations including the Lindblad dissipator has the form: where σ (n,m − ,m + ) are auxiliary density matrices indexed by the vectors n = {n 1 , n 2 , ..., n N } to account for the Drude part of the spectral density and m − = {m − 1 , m − 2 , ..., m − N }, m + = {m + 1 , m + 2 , ..., m + N } for the underdamped mode. The system density matrix ρ s (t) is given by the matrix indexed by zero vectors, σ (0,0,0) . n i±1 indicates the vector n with element i → n i ± 1 with similar definitions for m ± i±1 . The coherent system dynamics and incoherent rates for charge transfer, excitation and coupling to the leads are accounted for within the Liouvillian L where H e l is the electronic Hamiltonian and D(σ) = p γ p [a p σa † p − 1 2 {a † p a p σ}] is a Lindblad dissipator. γ p are the incoherent rate for coupling between the exciton states and the rest of the photocell with a † p (a p ) the associated jump operators. The vectors ν D , ν 0 and ζ 0 contain the Drude cutoff frequency Ω D , the mode damping constant γ 0 and the quantity ζ 0 = ω 2 0 − 1 4 γ 2 0 respectively, for each site of the system. The hierarchy above is given in the rescaled form which allows for efficient numerical computation of the time dynamics and steady state. The coefficients c D i , c − i and c + i are the leading coefficients of the exponential expansion of the bath correlation functions for Drude (D) and underdamped Brownian oscillator (±) spectral densities. Matsubara terms were not explicitly considered as they had little effect on the dynamics at 300K and increased the cost of computation. However a temperature correction term which approximates the effect of the Matsubara frequencies was The operators for coupling to the lower tiers are •} the superoperators representing commutator and anticommutator relations involving site i system-bath coupling operator V i = |i i|.

Pauli Master Equation: Modified Redfield Theory
An alternative framework used to model the PSIIRC photocell involves a Pauli master equation where modified Redfield theory is used to calculate the exciton relaxation dynamics 26 . The other incoherent rates for excitation, charge transfer and coupling to leads are the same as in the hybrid model framework. The reduced computational cost of the Pauli master equation allowed an investigation of the effect of the full PSIIRC spectral density on the behaviour of the photocell and also allowed the electron counting statistics of the photocell to be studied. Modified Redfield theory 43,69 describes population-to-population transfer among strongly coupled groups of chromophores. The transfer rates take into account the strong coupling to the phonon-bath but treat as a perturbation the off-diagonal coupling in the exciton basis i.e. H = k =k |k H el−ph kk k |, where |k are exciton states. The electron-phonon coupling matrix elements are given by H el−ph kk = N n=1 a kk (n)u n where a kk (n) is the spatial overlap of excitons |k and |k at site n. It is clear that modified Redfield should work well for systems in which we have quasi-localised excitons (ie. a kk (n) is small), which is the situation for the Chl D1 pathway in PSIIRC. The Markovian transfer rate between excitons a and b is given by the expression where the reorganisation energies and line broadening functions are given by is the amplitude of exciton p on site i, λ i is the reorganisation energy of site i and g i (t) is the line broadening function for site i. The rates for coupling to the leads, excitation and charge transfer are treated equivalently in both frameworks we consider. The Pauli master equation can be obtained from a Lindblad dissipator having the form D The jump operators a p (a † p ) causing population transfer between two states of the system give contributions to the time evolution of both the populations and coherences. Since the populations and coherences are decoupled from each other in the Lindblad formalism the coherences can be neglected and a master equation for the populations is obtained. This is given in the form of the Liouville space operator in Supplementary Eq. S8. The diagonal elements of the Pauli matrix have been omitted for compactness, but each column must sum to zero in order to satisfy a trace preserving time evolution. The rates k XY where X, Y = 1, .., 6 are the modified Redfield rates for transfer between excitons. The additional rates have been discussed in the previous sections where k XI(IX) and k Iα(αI) are Förster rates for primary and secondary charge separation and rates γn, γ(n + 1), Γ L and Γ R are for optical excitation and coupling to leads.

Supplementary Note 3: Comparison of Theoretical Frameworks for PSI-IRC Photocell Dynamics and Steady States
The two frameworks we consider in the main paper are compared here to justify the use of modified Redfield theory to describe exciton relaxation in the Pauli master equation framework. In Supplementary Fig. S2 the transient dynamics of the full photocell for the different spectral densities is presented, with the exciton populations summed to simplify the plot. For the Drude (J D (ω)) and single underdamped Browian oscillator (J 1 (ω)) spectral densities the hybrid model dynamics is well approximated by the Pauli master equation at both the initial stages and at longer times where the systems are relaxing towards a steady state. The Pauli dynamics with the full spectral density (J(ω)) is shown for completeness. The dynamics of exciton coherences as calculated with the hybrid model are shown in Supplementary Fig. S3. This shows the exciton coherences are small in amplitude and decay quickly towards a steady state close to zero. The use of modified Redfield which does not carry information about coherences between exciton states is justified on these grounds.
Supplementary Figure S2: Comparison of PSIIRC photocell population dynamics using the hybrid model and Pauli master equation. The exciton populations (X i ) have been summed to simplify the plots and the populations of the secondary CT state (α) are shown in the insets. The population of state β is not shown.
Additionally the steady states calculated by the two theories are compared in Supplementary Fig. S4 which shows the long time behaviour of the photocell. The steady state plays an important role in the counting statistics calculations. There is a close agreement between the two frameworks which suggests that calculations using the Pauli master equation should give results representative of the hybrid model behaviour.
Supplementary Figure S3: Exciton coherence dynamics within the hybrid model for spectral densities J D (ω) and J 1 (ω). Only coherences between exciton X 1 and the other excitons are shown.
Supplementary Figure S4 (S12) Introducing the pseudo-inverse R = QM −1 Q which is well defined since it is performed in the subspace orthogonal to the steady state. We then have which on substituting in Eq. S11 we get Equations S10 and S14 are the starting points for the recursive scheme. Expanding the zero eigenvalue λ 0 (χ), steady state |0(χ) and the perturbation ∆M(χ) about χ = 0 gives From this it is clear that the zero eigenvalue is equivalent to the cumulant generating function of the statistics. Substituting these expansions into Eqs. S10 and S14 give us the following expressions to calculate the cumulant to any order (n = 1, 2, 3, ...) The first and second order cumulant which are used to obtain the numerical results later on are explicitly shown below It is useful to note that M (n) = M J in the above since ∂ n ∂(iχ) n ∆M(χ)| χ=0 = ∂ n ∂(iχ) n M J (e iχ − 1)| χ=0 .

Supplementary Note 6: Counting Statistics of a Single Resonant Level
Equation 6 from the main text is plotted in Supplementary Fig. S6 to show the voltage dependence of the zerofrequency Fano factor for electrons passing through a single resonant level. The derivation of Eq. 6 is outlined in the following. The Liouvillian for the single resonant level system has the form 29,63 from which the steady state populations of the empty and occupied states can be calculated as ρ empty = Γ R Γ L +Γ R and ρ occupied = Γ L Γ L +Γ R . Substituting the expressions for the steady state populations into eV = E 0 + k B T ln( ρ occupied ρempty ) leads to Γ R = Γ L exp(−(eV − E 0 )/k B T ).
This can be used with the expression for the Fano factor in terms of the rates Γ L and Γ R 29 , to give the Fano factor as a function of V,