How Markovian is exciton dynamics in purple bacteria? Journal of Chemical Physics

We investigate the extent to which the dynamics of excitons in the light-harvesting complex LH2 of purple bacteria can be described using a Markovian approximation. To analyse the degree of non-Markovianity in these systems, we introduce a measure based on ﬁtting Lindblad dynamics, as well as employing a recently introduced trace-distance measure. We apply these measures to a chromophore-dimer model of exciton dynamics and use the hierarchical equation-of-motion method to take into account the broad, low-frequency phonon bath. With a smooth phonon bath, small amounts of non-Markovianity are present according to the trace-distance measure, but the dynamics is poorly described by a Lindblad master equation unless the excitonic dimer coupling strength is modiﬁed. Inclusion of underdamped, high-frequency modes leads to signiﬁcant deviations from Markovian evolution in both measures. In particular, we ﬁnd that modes that are nearly resonant with gaps in the excitonic spectrum produce dynamics that deviate most strongly from the Lindblad approximation, despite the trace distance measuring larger amounts of non-Markovianity for higher frequency modes. Overall we ﬁnd that the detailed structure in the high-frequency region of the spectral density has a signiﬁcant impact on the nature of the dynamics of excitons. ' 2017 Author(s). All article content, except where otherwise noted, is licensed under a Creative Commons Attribution (CC BY) license (http://creativecommons


I. INTRODUCTION
Photosynthesis provides the energy for nearly all life on Earth. The conversion of light energy into biologically useful chemical energy is a complex process performed by a variety of different organisms. Energy is captured in the form of photons and transported to a reaction center, where it is converted into useful chemical energy. The capture and transport process, generically known as light harvesting, is of great interest due to its near unit quantum efficiency 1,2 and because of recent suggestions of long-lasting quantum coherent effects based on two-dimensional (2D) femtosecond spectroscopic experiments. [3][4][5][6][7][8] Understanding how light harvesting achieves such high efficiency, and the role of quantum coherence in this efficiency, has been the goal of much theoretical and experimental work in the past few years, with heavy focus on the interaction of the excited states with their surrounding environment.
Within the body of theoretical work, many interesting and general insights into the potential role of the environment have been gained through the use of the Lindblad master equation. [9][10][11][12][13][14][15][16][17][18][19][20][21] Ideas such as noise-assisted transport, 9,10,[12][13][14]16 quantum locking, and momentum rejuvenation 14 have been used to suggest how interactions with the environment can help speed up energy transport in networks of chromophores. a) n.linden@bristol.ac.uk b) fred.manby@bristol.ac.uk However, it is still not certain whether these phenomena play a role in the real biological context, where the assumptions inherent in the master-equation approach may not apply. It has been shown, for physically realistic values of the environment coupling, that Lindblad and the related Redfield equations fail to capture the correct exciton-transfer rates. 22 The question we address here is whether or not a priori Lindblad equations, parameterised in any manner, are capable of reproducing the dynamics of exciton transfer. Previous research into this question has resulted in differing conclusions. [23][24][25][26] In particular, within a dimer extracted from the Fenna-Matthews-Olson (FMO) complex, small amounts of non-Markovianity were observed for an overdamped Brownian spectral density; 23 however, in calculations using the quasiadiabatic path-integral approach (QUAPI) on the full FMO complex with realistic spectral densities, no non-Markovianity was found. 24 Here we look at the degree of non-Markovianity in a dimer parameterised to reflect the B850 ring of light harvesting complex two.
The role of underdamped intra-molecular modes has been the focus of much attention recently, and in particular modes that are resonant with energy gaps between excitonic eigenstates. 11,[27][28][29][30][31][32][33] In some pigment-protein complexes, these modes show significant coupling to the excitons. [34][35][36][37] This has led to theoretical studies that provide evidence for the important role they may play in assisting energy transport and in producing the long-lived coherences seen in 2D spectroscopy experiments. [27][28][29][30][31][32][33] The phenomenon whereby these underdamped molecular vibrations induce long-lived coherences and drive transitions between excitonic eigenstates has been demonstrated using a semi-classical treatment of the intra-molecular mode, 28 which it was claimed demonstrates that the fundamental mechanism is not dependent upon non-Markovian dynamics. However, in Ref. 32 it is argued that, in models of dimers relevant for photosynthetic systems, efficient transport can show, and benefit from, non-classical effects, so there may be a significant role for non-Markovian effects. The extent to which non-Markovianity is increased through coupling to resonant modes is investigated in this paper.
Much progress has been made in developing metrics for quantifying non-Markovianity in the dynamics of open quantum systems. For a thorough review of these measures, we refer the reader to Ref. 38. Broadly speaking these measures can be categorised into those that look at properties of the dynamical maps that induce the time evolution of the system; [39][40][41][42][43] monotonicity of distance measures between two time-evolved trajectories of the system; [44][45][46][47] and properties quantifying the correlations between the system of interest and an ancillary system. 40, [48][49][50] Other measures are based upon the sign of decay rates of a Lindblad master equation in its canonical form 51 and properties of the affine transform induced by the dynamical map. 52 Some of these measures are complete in the sense that they capture any deviation away from the Markovian regime, whilst others capture only particular aspects of non-Markovianity. The complete measures either require information on the dynamical map or the use of an ancillary system which incurs extra, often impractical, computational cost. To tackle this issue we have developed a measure of non-Markovianity based upon finding the closest fitting Lindblad master equation. This measure has the benefits of being bounded (between 0 and 1), and being related in a direct and straightforward way to the feasibility of using the Lindblad approximation. This paper is structured as follows: In Section II we discuss the measures we will use to quantify the degree of non-Markovianity, introducing our measure based upon complete parameterisation and fitting of a Lindblad operator. In Secs. III, IV, and V, we apply these measures to chromophore dimers parameterised to model those of the light-harvesting complex LH2. In Section III we look at the effect of a broad phonon bath on the degree of non-Markovianity, which we model with the hierarchical equations-of-motion (HEOM) method. We show that despite small amounts of non-Markovianity, the dynamics is well-captured by a Lindblad approach, provided the site coupling in the excitonic Hamiltonian is a free parameter. In Section IV we treat the same dimer system coupled to high frequency intra-molecular modes. We present our calculations of the couplings to these modes which were computed using density functional theory calculations of a bacteriochlorophyll-a (BChl-a) molecule. We investigate how the degree of non-Markovianity varies with the number of discrete modes as well as the frequency of the discrete modes. We find that the degree of non-Markovianity increases with the number of discrete modes for both measures of non-Markovianity. We also find that, according to the Lindblad fitting measure, vibrational modes resonant with the gap between the excitonic eigenstates induce the most non-Markovianity. Finally, in Section V we combine both models and investigate more realistic structured spectral densities using HEOM to model the phonon bath and include a discrete intra-molecular mode. Here we find that the discrete vibrational modes contribute significantly to the degree of non-Markovianity, calling into question the use of Lindblad evolution to model exciton dynamics.

A. Measuring non-Markovianity
Before exploring measures of non-Markovianity, it is worth considering a definition for Markovian dynamics introduced in Ref. 38. Let ρ(t) be the density matrix describing the dynamical evolution of the state of some quantum system. The dynamical evolution from a state ρ(t 1 ) to a state ρ(t 2 ) at times t 2 ≥ t 1 can be described by the dynamical map that induces the evolution, E t 2 ,t 1 ρ(t 1 ) = ρ(t 2 ), where E is a completely positive and trace-preserving (CPTP) map and does not depend on the state upon which it acts. A family of such linear quantum maps E t 2 ,t 1 induces Markovian dynamics if for all t 2 ≥ t 1 ≥ t 0 it fulfills the composition law Breuer et al. 47 introduced a convenient measure for non-Markovianity based on the trace-distance metric. The trace distance quantifies the difference between two states, ρ 1 and ρ 2 , and is defined as It satisfies 0 ≤ D ≤ 1, and D = 0 if and only if ρ 1 = ρ 2 , and D = 1 if and only if ρ 1 is orthogonal to ρ 2 (defined as Tr ρ † 1 ρ 2 = 0). This metric has the property that it is a contraction under a CPTP map E, From our definition of Markovian dynamics above, a Markovian evolution is induced by a CPTP map and due to the composition law in Eq. (1), any possible division of that evolution into smaller time steps is also a CPTP map. Thus the trace distance between two initial states undergoing Markovian dynamics will monotonically decrease and any increase in the trace distance indicates non-Markovian effects. This monotonic decrease can be intuitively understood as arising from the entirely dissipative nature of the bath in the Markovian regime. Both initial states continuously lose information to the environment as they tend towards the same steady state. Only when there is a back-flow of information from the environment can this process be temporarily reversed and the distance between the states be increased. We would like to construct a measure for this property and quantify an answer to the question: given some dynamical process that takes ρ(0) → ρ(t), how non-Markovian is it? From the above property, Breuer, Laine, and Piilo (BLP) 47 concluded that the degree of non-Markovianity can be quantified by calculating the extent to which the trace distance increases between two density matrices starting at different initial states, where θ, the Heaviside function, ensures that the derivative only contributes to the integral when it is positive. This quantification is dependent upon the pair of initial states chosen. To remove this dependence, a search over all possible pairs of initial states can be performed and the maximum defined to represent the degree of non-Markovianity for the system, The computational difficulty of performing such a maximisation can be reduced by taking into account recently proved conditions for initial-state pairs that maximise N BLP , 53 which state that the maximum occurs for initial states that are orthogonal and lie on the boundary of the state space (which is equivalent to having a zero eigenvalue). An alternative approach is to consider typical initial states of the system. For the case of light harvesting, the true nature of the initial state is still debated; 26 however, for simplicity it is often assumed that the exciton is initially localised at a single chromophore. Thus, as well as looking at N BLP , we also investigate the typical amount of non-Markovianity, based on the fully localised initial states for a dimer system. We note that in such cases the degree of non-Markovianity may be an underestimate.
We will refer to this measure as the BLP measure of non-Markovianity and, in order to distinguish it from other uses of the trace distance in this paper, we will refer to the trace distance calculation associated with it as the BLP trace distance.
Although very useful, the disadvantages of the BLP measure are two-fold. First, there are some forms of non-Markovianity that the BLP measure will not capture. 38,51,54 In particular, we write a state ρ of dimension d in the form where α i are the elements of a real d 2 1 dimensional vector α known as the Bloch vector and F i are a set of orthogonal traceless hermitian operators (e.g., the Pauli matrices for d = 2). The action of any dynamical map on ρ can now be written in the form of an affine transform of the Bloch vector, It may be shown that the trace distance between two evolving states is independent of the vector t. 54 For this reason the BLP metric is sometimes called a witness rather than a measure of non-Markovianity. 38 Second, since the measure is not bounded above, it is not clear for what values of the measure it can be concluded that Lindblad-type master equations are no longer a good approximation. It is this second point that motivates the introduction of a Lindblad fitting scheme.
Any quantum map L inducing a Markovian time evolution can be written in Gorini-Kossakowski-Sudarshan-Lindblad form, 38 d where the second term describes the dissipative effect induced upon the dynamics of the system by a Markovian bath, and where the L j form a complete set of orthogonal traceless operators and h is a positive matrix. For a two-dimensional system, one can represent the density matrix ρ in its Bloch vector form (see Eq. (6)), The action of a Lindblad operator can now be written as the affine transformation in Eq. (7).
It can be shown that for a two-dimensional system that M is a real symmetric matrix 55 and can therefore always be diagonalised. It is convenient to parameterise the eigenvalues of M with three real parameters q i such that the eigenvalues are (q 2 + q 3 )/2, (q 1 + q 3 )/2, and (q 1 + q 2 )/2. The conditions on the parameters can be derived by ensuring that the corresponding matrix h in Eq. (9) is always positive. 56,57 These result in the conditions, We then explore the complete set of Lindblad dissipators by varying q i , t i , and the 3 × 3 orthogonal matrices (parameterised by three real parameters β 1 , β 2 , β 3 corresponding to angles) that diagonalise a general 3 × 3 symmetric matrix M.
It should be noted that this parameterisation characterises fully only those Lindblad evolutions with constant rates. This simple parameterisation, consisting of nine real parameters, makes it possible to perform an optimisation to find the Lindblad master equation that best fits a given quantum evolution. Given some quantum evolution of a two-state system ρ(t), we can find the closest-fitting Lindblad evolution ρ L (t) by minimising the function, with respect to the parameters t i , q i , and β i and subject to the constraints in Eq. (10). In Eq. (11) we have employed the trace distance to quantify the difference between the exact evolution and the Lindblad evolution at each time, t. The integral quantifies the total difference during the trajectory, resulting in Q being the mean trace distance between the two trajectories. This measure acts as a complementary measure of non-Markovianity to the BLP measure presented earlier. It has the advantages that it is bounded between 0 and 1 and that it directly quantifies the extent to which the non-Markovianity present in the system renders the Lindblad master equation ineffective at capturing the dynamics. Examination of the properties of the best-fitting Lindblad also provides some insight into the effects that non-Markovianity has on the dynamics of the system.
As for the BLP measure, we can show that the largest values of non-Markovianity from this fitting measure will be for pure initial states. Suppose our initial state is a mixed state, where σ i are some set of pure-state density matrices and p i the associated probabilities. Let the dynamical maps that induce a HEOM and Lindblad evolution from time 0 to t be given by E H t and E L t , respectively. The trace distance between these two trajectories at time t is then given by By the convexity property of the trace distance, we can show that this is always less than or equal to the weighted sum of the trace distance between the time evolved pure state density matrices, Thus, the total trace distance over the evolution of the mixed state must always be less than or equal to the trace distance over the evolution of one or more of its pure-state components.
In order to find the Lindblad equation that best captures the underlying processes of a physical system, it is not sufficient to minimise the function Q for a single initial state as the resulting Lindblad may not capture the correct dynamics for other initial states. Thus we compute the best fitting Lindblad, averaged over initial states, However, in practice we found that whether one fits to many initial states or just to one initial state, the resulting Q value is almost identical in both cases. In practice we used 105 pure initial states selected to uniformly cover the surface of the Bloch sphere.

B. System bath model
To model the dynamics of excitons interacting with an environment, we adopt the standard system-bath Hamiltonian where the system Hamiltonian H S describes the dynamics of the exciton, H B describes the bath or environment, and H I describes the interaction between system and bath. The system is a chromophore dimer described by the commonly used Frenkel exciton Hamiltonian, where a ( †) i are annihilation (creation) operators of an excitation at bacteriochlorophyll (BChl) i, i is the energy of a single excitation at site i, and V ij is the electronic coupling between BChls i and j. Here we restrict ourselves to the one-exciton manifold, an approximation that is well-justified for certain photosynthetic organisms, such as purple bacteria, where the low intensity light conditions of their natural environments are such that the time scale between the arrival of photons is large enough that the probability of two excitations encountering each other during transport is very low. 58 The environment is modelled as a collection of quantum harmonic oscillators, with each BChl interacting with an independent set of bath modes of identical frequency and coupling strength. These vibrational modes are used to represent the influence of discrete intramolecular vibrational modes as well as the effect of solvent modes. This bath Hamiltonian is given by where b ( †) s,i are the annihilation (creation) operators of quanta for bath mode s of frequency ω s that is coupled to BChl site i.
The interaction Hamiltonian describes the coupling of the bath modes to the diagonal elements of the exciton Hamiltonian, Here g s is the coupling strength associated with bath mode s and is related to the reorganisation energy λ s = g 2 s / ω s and Huang-Rhys factor Λ s = g 2 s /( ω s ) 2 . The operator b † s,i + b s,i is proportional to the displacement operator of bath mode s of BChl i. Note that the summation over the index s applies in the case of discrete vibrational modes, but in the continuous case should be thought of as an integral.
The coupling parameters can be captured in a single function of frequency known as the spectral density J(ω), whose magnitude determines the coupling strength, This is particularly useful when the number of bath modes is infinite, for example, when modelling interaction with a solvent environment or a damped vibrational mode. For clarity we have adopted the notation for spectral densities as defined in Ref. 59. From now on in this paper, when referring to the spectral density, we will be referring to the antisymmetric component of the Fourie-Laplace transform of the energy gap correlation function C (ω), which is related to the definition above by

III. BROAD SPECTRAL DENSITY
The interaction of a BChl molecule with a protein and solvent environment can be modelled by coupling the system to a broad continuous spectral density. Here we explore the effect this interaction has on the degree of non-Markovianity in the dynamics of an exciton. To model this interaction, we use the HEOM method, 60 which permits an essentially exact treatment of quantum dynamics at finite temperature. The main disadvantage is that few spectral densities are computationally tractable with the method: here we employ the Brownian oscillator spectral density 61 where E λ is the total reorganisation energy, and γ is a damping constant, inversely related to the bath correlation time scale by γ = 1/τ c . This function determines the coupling parameter g s in the interaction Hamiltonian (Eq. (19)) through To obtain physically realistic values for these parameters, we compared this functional form to the spectral density fitted to experiments on the B777 monomer pigment of the LH1 complex by Renger and Marcus. 62 The B777 complex consists of a single BChl-a molecule bound to an α-helix protein, similar to the chromophore environment in the B850 ring of LH2. This consists of 18 BChl-a molecules closely arranged in a ring structure, each bound to α-helices. 63 Due to the similarity of the environment, it is reasonable to assume that the same shape function will apply well to the B850 ring of LH2, and this is what we have done. We configured E λ to equal the reorganisation energy of the B777 spectral density and fitted γ such that the peak frequencies matched.
The experimentally fitted and Brownian spectral densities are shown in Figure 1: the imperfect agreement between them has been well-documented. 59 It has been shown for this spectral density that decreasing γ (increasing the correlation time τ c ) results in greater amounts of non-Markovianity. 23 The effect of this parameter change on the shape of the spectral density is decreased coupling at higher frequencies and a sharper peak in the lower frequency regime. Similarly, here we see that the B777 spectral density has decreased coupling in the high frequency regime and a sharper peak in the lower frequency regime. Thus, by analogy, we suspect that the B777 spectral density would exhibit greater non-Markovianity than the Brownian spectral density used here.
Using the HEOM approach 60 as implemented in the Phi software package, 64 we modelled the effect of the Brownian spectral density on a symmetric dimer parameterised to reflect a dimer of the LH2 complex of purple bacteria. The site energies for the BChls were taken to be equal, and the coupling strength between them was set to V = V 12 = V 21 = 245 cm −1 (4.61 × 10 13 Hz). (Later on in Sections IV and V we also look at asymmetric dimers where the site energies differ by 1 − 2 = 200 cm −1 (3.77 × 10 13 Hz).) The initial state of the bath was taken to be a thermal initial state at temperature T = 288 K.
Here we look at the dynamics of this dimer and assess the degree of non-Markovianity using the Lindblad fitting measure. The fitting of the Lindblad was taken by minimising Q over 105 initial states chosen from a Lebedev grid 65 over a hemi-sphere of the Bloch sphere. Figure 2 (left column) shows the result for the initial state ρ(0) = |0 0| which is a typical result. It shows the dynamics of the site population ρ 11 for the HEOM calculation and the best fitting Lindblad. The upper row is for the case of fixed coupling strength V. The lower row is for the case where we allow this coupling strength to be optimised, along with the parameters that define the Lindblad operator. For the fixed Hamiltonian, we see significant differences between the two evolutions. For the Lindblad case, the frequency of oscillations in the site population is determined entirely by the system Hamiltonian. For the HEOM calculation however, the frequency of oscillations in the dynamics is determined by the system Hamiltonian as well as the shape of the spectral density and the strength of interaction with it.
The interaction with the Brownian spectral density, as parameterised here, produces faster oscillations in the dynamics compared to the isolated system, effectively increasing the coupling strength between the two chromophores. This is shown in the lower row of Figure 2 where an increase of the coupling strength of around 10% from V = 245 cm −1 (4.61 × 10 13 Hz) to V = 270 cm −1 (4.61 × 10 13 Hz) greatly improves the fit. The site populations show very similar values beyond 200 fs for this fit. Up to this point there is a small deviation between the two, indicating non-Markovian effects. In particular, the HEOM evolution shows slightly smaller-amplitude oscillations than the Lindblad evolution, corresponding to a smaller transfer of exciton density between the sites in this region. This result is seen clearly in the lower right panel of Figure 2, which shows the trace distance between the two evolutions. After around 200 fs, the density matrices are nearly identical between the two calculations, and even in the earlier dynamics the deviation is small.
The full range of Q values for all initial states used to find the optimal Lindblad is shown in Figure 3. The left figure is for the case of a fitted coupling strength V and the right figure shows the fixed result. Here again we can clearly see the large difference in values between the cases where the coupling strength is and is not optimised. In particular we see a broader spread for the fixed Hamiltonian case. This arises because for some initial states there are minimal oscillations in the dynamics, for example, the |+ +| or |− −| states. When the oscillations in the dynamics are of very low amplitude, it does not matter if the frequencies are mismatched; hence, these initial states fit well. This illustrates the need to allow the coupling strength to be a free parameter if we are to interpret the fitting values as a measure of non-Markovianity. Somewhat counterintuitively, when we do allow the coupling strength to vary it is the |− −| initial state that shows the greatest amount of non-Markovianity, despite being an initial state for which the site populations are initially at their equilibrium value. In this case most of the non-Markovianity arises in the dynamics of the off-diagonal elements of the density matrix. Allowing V to be fitted greatly reduces the value of Q and results in an average value of around 0.009.
The need to vary the system Hamiltonian shows the potential difficulties in finding the correct Lindblad master equation to model the dynamics of an exciton system so strongly coupled to the environment. However we can conclude that allowing the coupling strength between BChls to be optimised is necessary if we wish to assess the general applicability of Lindblad master equations and to quantify the degree of non-Markovianity. Thus in all subsequent Lindblad fitting optimisations, we have included the BChl coupling V in the set of optimised parameters. Allowing the site energies to vary does not appear to significantly impact the results.
We can examine the resulting form of the optimal Lindblad for this spectral density. If we write the Lindblad equation in the form of Equation (9)  . (24) We can see that the strongest contribution is the diagonal term for h 33 , which corresponds to pure dephasing in the site basis. This is the form typically used when modelling the interaction of the environment with excitonic systems. [9][10][11][12][13][14][15][16][17][18][19][20][21] It is clear however that there are significant contributions from other Lindblad terms and that a pure dephasing Lindblad is not nearly the optimal choice. Diagonalising h we find that the optimal Lindblad is dominated by two processes, The Lindblad operators corresponding to these processes are By construction, L 1 and L 2 are traceless and orthogonal (Tr(L i L † j ) = 2δ ij ). But neither is hermitian, and so they cannot simply be interpreted as dephasing in any basis. Figure 4 shows the trace distance used to calculate N BLP (ρ 1 , ρ 2 ), where ρ 1 (0) = |1 1| and ρ 2 (0) = |2 2|. Increases in the trace distance are seen throughout the evolution up until it reaches its equilibrium state, showing that non-Markovian dynamics persist beyond the 200 fs indicated in the Lindblad fitting result. However, it is clear that the most  (2)) as a function of time for the HEOM calculation in Figure 2 where ρ 1 (0) = |1 1| and ρ 2 (0) = |2 2 |. Regions of positive derivative correspond to non-Markovian dynamics. The total non-Markovianity N BLP (ρ 1 , ρ 2 ) = 0.14. significant increases occur during the first 200 fs; thus, the two measures show a reasonably good agreement. The total increase in the trace distance over the 1 ps of dynamics gives a total non-Markovianity of N BLP (ρ 1 , ρ 2 ) = 0.14, which is of a similar size to other work quantifying non-Markovianity for the Brownian spectral density. 23 Overall, we find that for these parameters of the Brownian spectral density, the exciton dynamics shows small amounts of non-Markovianity in both the BLP measure and the Lindblad fitting measure. The majority of the non-Markovianity occurs in the first 200 fs of the dynamics, where we see through the Lindblad fitting measure that it results in reducedamplitude oscillations in the site population, corresponding to the reduced initial transfer of exciton density to the neighbouring chromophore. However, the Lindblad fit is still capable of capturing the key features of the dynamics, suggesting that for these small values of the BLP metric, it may still be acceptable to adopt the Markov approximation provided the coupling strength in the Hamiltonian is treated as a free parameter.

IV. DISCRETE SPECTRAL DENSITY
Here we investigate the effect that high-frequency, discrete modes of intra-molecular origin have on the degree of non-Markovianity in the exciton dynamics. We look at two key aspects. First, we investigate how the degree of non-Markovianity scales with the number of intra-molecular modes. Second, we study the effect of modes resonant with the gap in the excitonic spectrum by looking at the frequency dependence of the non-Markovianity for a single mode.
In order to get physically realistic parameters for the coupling strength of intra-molecular modes, we performed quantum chemistry calculations on an isolated bacteriochlorophyll-a molecule. The long phytyl chain was replaced with a methyl group to reduce computational expense; we anticipate that this truncation will not affect our conclusions. Using Gaussian 09 66 we used density functional theory (DFT) to optimise the geometry and compute normal modes at the B3LYP/TZVP level of theory. To determine the coupling strength between these vibrational modes and the electronic excitation, we used time dependent DFT (TDDFT) at the B3LYP/SVP level of theory to compute the transition energy from the ground to first-excited state at discrete points along each normal mode and fit this function to a straight line. The gradient d s of the resulting line gives the rate of change of excitation energy with displacement along the normal mode s, which the linear approximation is related to the coupling parameters that appear in Eq. (19) by Here m s is the effective mass of the normal mode s and ω s is the angular frequency. Thus for each vibrational mode of the molecule, we obtained the corresponding angular frequency and coupling strength. Figure 5 shows the corresponding Huang-Rhys factors for these modes.
To model these vibrational modes, we performed exact diagonalisation (with truncated boson number) of the systembath Hamiltonian introduced previously. For the dimer we consider here, the dimensionality of the bath can be reduced by defining a set of collective bath operators that act on both sites. These are defined by where modes described by c ( †) s influence both sites in phase, and r ( †) s modes act on each site in an equal and opposite manner. We can see this by substituting the inverse relations into the bath and interaction Hamiltonians to obtain (33) Since we have restricted ourselves to the one exciton basis, a † 1 a 1 + a † 2 a 2 is a constant of the dynamics. Thus displacement of the modes c ( †) s results in a uniform shift in the energy of both sites which has no effect on the dynamics. They are therefore The advantage of this approach is that for a small number of high-frequency modes ω k b T , the boson number can be truncated (without significant loss of accuracy) at a level at which exact diagonalization can still be performed.
The modes we are modelling here represent delta peaks in the spectral density. In the physical system, these delta peaks will be broadened due to interactions with the solvent environment, which acts to dampen the non-equilibrium oscillations of the vibrational mode. The time scale for this damping process depends upon the shape of the broadened peak. For Lorentzian broadening, the peak is defined according to the function where λ is the reorganisation energy, ω 0 is the frequency of the vibrational mode, and γ l the broadening factor such that C (ω) → δ(ω − ω 0 ) as γ l → 0. For this shape function nonequilibrium oscillations decay with e −γ l t . 27 Thus the results that follow are accurate for modes where the time scale for decay is long compared to the simulation time scale.
Here we look at how increasing the number of modes affects the degree of non-Markovianity in the excitonic subsystem. To do this, we perform exact diagonalisation on the full system-bath Hamiltonian and trace out the bath modes to obtain the excitonic reduced density matrix. For this calculation, due to the computational expense of the exact diagonalisation, we look only at the typical amount of non-Markovianity for initial states where the exciton is fully localised, i.e., we calculate N BLP (|1 1| , |2 2|) as well as the corresponding Q values for these initial states. Again the bath modes are taken to be in a thermal initial state at temperature T = 288 K. Figure 6 shows the amount of non-Markovianity present for a symmetric dimer coupled to N m discrete modes from our DFT calculations with frequencies ω > 10 14 Hz. The modes were ordered by strength of coupling, with the strongest coupling modes added first. The Huang-Rhys factors range from 0.012 to 0.004. The Fock space of the bath is capped to have at most N b = 15 quanta in total, which due to the small Huang- Rhys   FIG. 7. The degree of non-Markovianity as a function of the frequency of a single mode coupled to symmetric (left panels) and asymmetric (right panels) chromophore dimers. In each case metrics were extracted from 1 ps of dynamics. Each point in the upper plots represents the non-Markovianity for a particular pair of orthogonal initial states, whereas for the lower plots each point is for a single initial state. Red and arrows indicate the frequency degenerate with the gap between eigenvalues of the excitonic Hamiltonian. factors and high frequency of the modes is sufficient to capture essentially the exact dynamics. This was checked by increasing the boson number until convergence to a single result was observed. Each calculation is run for a total of 1 ps.
On the left of Figure 6, we can see that these discrete modes introduce a significant increase in the BLP trace distance corresponding to non-Markovian dynamics. We see similar results on the right hand side for our Lindblad-fitting metric. As the number of modes increases, the degree of non-Markovianity also increases. This is due to the increased total coupling between the bath and the system as each new mode is included. The decreasing gradient is due to the ordering of the modes that are added, starting from strongest to weakest coupling. These results suggest that the underdamped modes that exhibit prolonged non-equilibrium oscillations in the high-frequency region of the spectral density are likely to contribute in a combined way to increase the non-Markovianity of the system dynamics. In other words, increasing the amount of structure in the higher frequency regime of the spectral density will correspond with increased non-Markovianity.

A. Frequency dependence
Here we investigate how the frequency of a discrete vibrational mode affects the degree of non-Markovianity. For completeness and since we are looking for resonances with the excitonic system, we study an asymmetric dimer ( 1 − 2 = 200 cm −1 (3.77 × 10 13 Hz)) as well as a symmetric one ( 1 = 2 ). In each case we couple the dimer to a single mode with a Huang-Rhys factor, S = 0.1. We perform a maximisation over 105 initial states chosen from a Lebedev grid 65 over a hemi-sphere of the Bloch sphere, with the antipodes of these points forming their orthogonal counterparts. The bath mode is taken to be a thermal state at T = 288 K. The BLP trace distance is then calculated between orthogonal pairs and the Markovianity measure N BLP (ρ 1 , ρ 2 ) (Eq. (4)) is calculated. We also apply our Lindblad fitting metric Q , which is applied to the all trajectories resulting from the initial states used for the BLP measure. The results are shown in Figure 7. Figure 7 shows the results for both the symmetric and asymmetric dimers using the BLP measure N BLP and Lindblad fitting measure Q. For the BLP metric in both the symmetric and asymmetric cases, higher-frequency modes are found to lead to more non-Markovian dynamics. The faster time scale of dynamics for the high-frequency modes results in a more rapid exchange of information between the system and the environment and therefore a greater total exchange over a picosecond of dynamics. Using this measure alone one might conclude that there is nothing significant about the non-Markovianity induced by resonant modes. However, for the Lindblad fitting measure Q, we find a maximum in the degree of non-Markovianity at the resonant frequencies. This implies that although there is a greater total increase in the trace distance for the case of higher frequency modes, the contribution of that increase to deviations from the Lindblad model is smaller.
To understand why this is, we examined the time evolution of the trace distance used to calculate the BLP measure. Figure 8 shows three typical examples for a low-frequency mode (top panel), a resonant mode (middle panel), and a highfrequency mode (bottom panel). It can be seen that although the total increase in the trace distance is greatest for the high frequency modes, the amplitude of the trace distance oscillations is greatest for the near resonant modes. These large amplitude oscillations in the BLP trace distance have a significant effect on the dynamics of the exciton, resulting in irregular oscillations in the site populations that the Lindblad equations are unable to capture. The essential role resonant modes could play in assisting energy transport, which is to drive transitions between excitonic eigenstates and is observed for an essentially Markovian semi-classical treatment of the resonant mode. 28 What we see here is that for a full quantum treatment, these resonant vibrational modes induce large-amplitude non-Markovian oscillations that correspond to a maximal deviation from the Lindblad approximation. Whether this strong non-Markovianity contributes to assisting energy transport is an open question. We have also shown that when trying to determine the applicability of Lindblad master equations using the BLP metric, the most important feature is the amplitude of the tracedistance oscillations, rather than the total increase in the trace distance. This is essentially an extension of the typical justifications for a Markovian approximation but extended to apply to strict non-Markovianity measures. For small oscillations in the BLP trace distance, it is easy for the monotonically decreasing trace distance of a Lindblad evolution to provide a good fit. Similarly if the time scale for the oscillations is small compared to the time scale of the phenomena of interest, then the dynamics can be well modelled by Lindblad dynamics.

V. STRUCTURED SPECTRAL DENSITY
Having looked at the non-Markovianity in the presence of a broad spectral density to model a low-frequency phonon bath, and a discrete spectral density to model intra-molecular modes, we now turn our attention to a more realistic model that includes both. To do this, we include two environments: a discrete vibrational mode that is treated quantum mechanically on the same footing as the system Hamiltonian, and a continuous bath as before.
The discrete bath mode is described as where r ( †) is the annihilation (creation) operator for the collective vibrational mode defined in Section IV. The continuous bath mode, modelled by H B 2 and H I 2 , is treated through the HEOM formalism; as before coupling strengths g s are defined by the Brownian spectral density (Eq. (22)).
To model the dynamics of the total Hamiltonian H = H S + H B 1 + H I 1 + H B 2 + H I 2 , we used H S + H B 1 + H I 1 as the system Hamiltonian in the Phi software package. This effectively creates a system of n sites, where n is the dimension of the Hilbert space of the Hamiltonian after truncation of the boson number for mode r. Each site represents a state of the excitonic dimer and a state of the bath mode. We then ensured the HEOM bath coupled only to the excitonic system by correlating coupling terms on sites with identical exciton states. We then traced out states associated with mode r to recover the excitonic density matrix. FIG. 10. The BLP trace distance as a function of time for the combined HEOM and discrete mode calculation in Figure 9 where the initial states were ρ 1 (0) = |1 1| and ρ 2 (0) = |2 2 |. The non-Markovianity measures N BLP = 0.9.
Initially we chose the strongest coupling mode from our DFT calculations whose frequency is ω > 10 14 Hz. Figure 9 shows the result for a mode frequency ω = 2.26 × 10 14 Hz, Huang-Rhys factor S = 0.012, and for the |1 1| initial state which is a representative example. The optimal Lindblad was found by minimising Q over 105 pure initial states as in Sections III and IV A. Here we see clearly the non-equilibrium oscillations that the undamped vibrational modes induce on the system dynamics. In the HEOM calculations where only the phonon bath is included, the system reaches equilibrium at around 700 fs. Here we see that beyond this point nonequilibrium oscillations persist due to the undamped vibrational mode. This can also be seen in the BLP trace-distance evolution in Figure 10, where oscillations continue throughout the evolution.
The BLP measure results in a non-Markovianity of N BLP (ρ 1 , ρ 2 ) = 0.9 which is around seven times larger than for the system coupled only to the phonon bath. The Lindblad fitting on the other hand shows a smaller increase in non-Markovianity of approximately a factor of two, with Q = 0.013. An examination of Figure 9 illustrates the disparity, as although the total increase in the trace distance is greater, the effect on the dynamics has been minimal. This again illustrates that the amplitude of increases in the trace distance is the decisive factor in determining the applicability of Lindblad master equations.
In order to interpolate between the regime of a single mode and the regime where many intra-molecular modes are coupled to the system, we looked at the effect of a single effective vibrational mode which combines the coupling strength of neighbouring vibrational modes into one. 13 modes were combined between frequencies 1.31 × 10 14 Hz and 1.47 × 10 14 Hz to form the single bath mode modelled by the Hamiltonians in Eqs. (37) and (38). The frequency of the effective mode is given by the mean frequency and the Huang-Rhys factor is given by the sum of the Huang-Rhys factors of the selected modes (this commonly used approach 34,37,67 conserves the average number of quanta excited in the bath). The resulting mode has frequency ω = 1.4 × 10 14 Hz and Huang-Rhys factor S = 0.026. Again we performed a minimisation of Q over 105 different pure initial states and showed the |1 1| initial state as a representative example.
The results are shown in Figures 11 and 12. It can be seen that for a mode of this coupling strength and frequency we start to see a more significant deviation from the Lindblad models. The site populations no longer show the characteristic oscillatory exponential decay but include irregularities throughout the evolution. The amplitude of the long lived non-equilibrium oscillations is now much greater than before. The average trace distance between the Lindblad evolution and the HEOM calculation is Q = 0.04, and the BLP measure of non-Markovianity gives a value of N BLP (ρ 1 , ρ 2 ) = 1.2. Taking into account the differences between the exact calculation and the Lindblad seen here and our results earlier in Figure 6, where we showed that non-Markovianity increases as a more detailed structure is included in the spectral density, and considering our DFT calculations which show that there are many more intra-molecular modes that have been excluded from this particular calculation, the evidence suggests that a Lindblad equation would no longer be able to reproduce the dynamics for a realistic spectral density.

A. Frequency dependence
Finally we look at how the non-Markovianity depends upon the frequency of the discrete vibrational mode in the FIG. 13. Variation in the degree of non-Markovianity for the Brownian spectral density and a discrete mode as a function of the frequency of the discrete vibrational mode. The vibrational mode couples with Huang-Rhys factor S = 0.1. In the upper plots, each point represents the non-Markovianity for a different pair of orthogonal initial states. For the lower plots each point is for a single initial state. Red represents the frequency which is on resonance with the difference in energy between the eigenstates of the excitonic Hamiltonian. presence of a broad spectral density. The calculation performed is equivalent to that in Section IV A, but with the additional presence of the Brownian spectral density. Figure 13 shows the result for a symmetric and asymmetric dimer. The degree of non-Markovianity overall is significantly reduced, by approximately an order of magnitude, compared to the calculations without the broad spectral density. The overall trends however are similar. The BLP non-Markovianity measure again shows that high frequency modes induce the greatest total trace-distance increase, due to the faster dynamics and therefore faster exchange of information. Again the Lindblad fitting measure shows that there is a maximum in the non-Markovianity. Interestingly, however, this maximum is no longer at the frequency resonant with the system Hamiltonian (represented by the red points) but instead occurs at slightly higher frequencies. The shift seems to correspond to the shift in the system coupling strength V required to fit a Lindblad to the broad spectral density in Section III. If we take the resonant frequency of this fitted Hamiltonian (represented by the arrows), we find that it matches well with the new maximum in the non-Markovianity. This shows that the effective Hamiltonian of the system is really changed by the interaction with a broad spectral density. It is possible that if light harvesting systems have evolved to take advantage of resonance with strongly coupling discrete modes, they will have evolved such that the discrete mode is resonant with this effective system Hamiltonian.

VI. CONCLUSIONS
We have studied the excitonic dynamics of a chromophore dimer interacting with a vibrational environment in order to assess the applicability of Lindblad-type master equations for modelling photosynthetic energy transfer. The environment was modelled using both HEOM, to model a broad phonon bath, and exact diagonalisation, to treat discrete vibrational modes of intra-molecular origin. To assess the validity of Lindblad master equations, we employed a recently developed method, based upon the trace distance, to measure the degree of non-Markovianity by quantifying increases in the trace distance between two evolutions. We also introduced and applied a measure of non-Markovianity which is based upon a general parameterisation of a Lindblad master equation for the twosite case. This general parameterisation allows one to quantify the distance between a general two-state quantum evolution and the best-fitting Lindblad evolution.
For a Brownian spectral density, parameterised to fit the spectral density of the B777 complex, small amounts of non-Markovianity are observed using both metrics. In particular, the Lindblad fitting metric shows that this small amount of non-Markovianity results in slightly reduced-amplitude oscillations in the site populations; however, overall the main features of the dynamics are well-captured by a Lindblad equation.
We found that to find the optimal Lindblad evolution, the system Hamiltonian had to be adjusted in such a way that the coupling between the chromophores was increased. This is direct evidence that interaction with an environment can speed up the transfer process between chromophores by introducing faster exciton transfer between the two sites.
However, whether similar adjustments to the Hamiltonian would be required for a chromophore system of more than two sites is unknown. What this does show, however, is that it is not possible to accurately model exciton transfer by simply employing a Lindblad formulation with a system Hamiltonian parameterised separately.
We also investigated the significance of the frequency of vibrational modes and whether any particular frequencies introduced greater non-Markovianity. We investigated this with and without the presence of a broad spectral density. We found that for the trace-distance metric of non-Markovianity, higher frequency modes introduce the most non-Markovianity. This is because higher frequency modes introduce faster oscillations in the trace distance resulting in a greater total increase over a given time period.
For our Lindblad fitting measure of non-Markovianity, however, we found that, for the case excluding the broad spectral density, modes resonant with the difference between eigenvalues of the system Hamiltonian showed the greatest non-Markovianity. The difference between the two measures can be understood in terms of the amplitude of oscillations in the trace distance. The modes resonant with the excitonic system induce the largest amplitude oscillations in the trace distance. Thus the amplitude of oscillations in the trace distance is the decisive feature by which to determine whether or not the Lindblad approximation is a good fit, as opposed to simply observing the total trace-distance increase. This result applies generally beyond applications to photosynthetic systems. Whether the observed large amplitude oscillations in the trace distance for resonant modes contributes to the recently observed increase in energy transfer rates [27][28][29][30][31][32][33] is left as an open question.
For the case where the broad spectral density was included, we found that the maximum non-Markovianity occurred at slightly higher frequencies than the resonant frequency of the system Hamiltonian. We showed that this shift is most likely due to the change in the effective Hamiltonian of the system caused by the interaction with the broad spectral density. The Hamiltonian resulting from the Lindblad fit to the broad spectral density, which contained a parameterised coupling strength V, was shown to have an eigenstate energy difference closer to the frequency of maximum non-Markovianity.
Finally, we looked at whether a Lindblad master equation would be capable of modelling the dynamics of a system coupled to a more realistic spectral density. In particular, we studied the effect of including some of the discrete vibrational modes that are responsible for the structure in the highfrequency region of the spectral density. For the dimer model interacting with only discrete vibrational modes, we found that as the number of discrete vibrational modes was increased, the degree of non-Markovianity also increased. For a model that included a single discrete intra-molecular mode and a broad spectral density, the impact on the degree of non-Markovianity was not very significant. However there are many such discrete intra-molecular modes, so in order to bridge the gap we modelled a single effective mode which incorporated the coupling strengths of 13 modes of a similar frequency. Here we found that the evolution deviates from the Lindblad approximation,