B800-to-B850 relaxation of excitation energy in bacterial light harvesting: All-state, all-mode path integral simulations

Significant efforts in physical chemistry and chemical physics are currently focused on understanding excitation energy transfer (EET) in molecular aggregates. At the heart of this venture is the molecular apparatus of photosynthetic bacteria and plants specifically designed to harvest solar energy, known as light harvesting complex.¹ Most work has focused on the inner- and outer-antenna light harvesting complexes of purple photosynthetic bacteria, known as LH1 and LH2, the analogous complexes of photosystem I found in higher plants, and the Fenna–Matthews–Olson (FMO) complex. Ever since their discovery and characterization,^2–6 light harvesting complexes of pigments and biomolecules have been used as prototypes for understanding intermolecular energy transfer riddled with the complex interaction of electronic and vibrational degrees of freedom whose interplay is responsible for the observed EET. Through the intriguing architecture of many light harvesting complexes, energy appears to be “funneled” to the photosynthetic reaction center with unmatched efficiency. However, several questions still remain unanswered. Hidden within the remarkable energy transfer efficiency of nature’s machinery is the question about the role of vibrational motion in facilitating or impeding electronic energy transfer pathways. On the one hand, quantum mechanical energy transfer between nonadiabatically coupled electronic states exhibits rich effects of delocalization, superposition, and interference. At the same time, the molecular vibrations within large light harvesting complexes and their encapsulating biomolecular environments induce dynamic disorder, which leads to damping of quantum coherence. Moreover, vibronic modes resulting from strong exciton–vibration interaction can modify EET time scales in nontrivial ways. From a more practical perspective, clarifying the role of vibrations^7–10 in biological energy transfer is important for designing artificial systems for efficient solar energy harvest and storage.

Experimental work in this field has benefited from the advent and use of multidimensional spectroscopic techniques. Apart from the characterization of light harvesting complexes using linear absorption spectra, two-dimensional electronic spectroscopy developed in the 1990s has been used to investigate the explicit mechanisms of ultrafast energy transfer.^11–16 Recently, two-dimensional electronic–vibrational spectroscopy has emerged as a powerful tool for the selective characterization of vibronic timescales.^17–19 

Almost all available theoretical tools have been applied to the dynamics of EET in light harvesting complexes.^20–50 Perhaps not surprisingly, theoretical and computational investigations of EET in these large systems have encountered several challenges. Beyond the prevalence of nonadiabatic coupling, which precludes the possibility of simpler propagation on Born Oppenheimer potential surfaces, the scales of the electronic couplings, site asymmetries, static disorder, vibrational frequencies, reorganization energies, and thermal energy (at or around room temperature) fall within a narrow range in the typical EET regime. Thus, perturbative approaches, such as Redfield or Forster approximations, are rendered unreliable.^22,26,29 While simple one- and two-vibration models^{39,40,43,47,49} provide valuable insights, they lack the essential effects of decoherence from the nuclear motion on a broad range of time scales, which is essential to energy transfer. Moreover, due to the large size and complexity of the biomolecular apparatus, even within a harmonic treatment of all vibrational degrees of freedom, the size of the total Hilbert space grows too rapidly, rendering traditional wavefunction methods rather inefficient. Classical trajectory based simulations have been performed with surface hopping methods,⁴⁴ while other treatments exploit the Meyer–Miller–Stock–Thoss mapping^51,52 in combination with linearized semiclassical initial value representations^33,53 or partially linearized path integral approximations.³⁴ Fully quantum mechanical calculations have been performed using hierarchical equations of motion^54,55 to investigate the EET dynamics of the FMO complex^28,29,45 with the assumption of Drude bath models, and subsequent calculations^42,47 have added one or more vibrational or vibronic modes to the model bath.

The path integral formulation of quantum dynamics^56,57 leads to efficient tools for performing fully quantum mechanical calculations on condensed phase processes. Perhaps, the most appealing feature of the path integral formalism is the ability to integrate out any number of harmonic degrees of freedom exactly,⁵⁸ at any temperature, without increasing the computational cost. The oldest real-time path integral method for system-bath Hamiltonians, the quasi-adiabatic propagator path integral (QuAPI) algorithm,^59,60 has been used to simulate the dynamics of the FMO complex.^36,37,61 The quantum–classical path integral^62,63 (QCPI) that (while applicable to anharmonic environments) yields exact, fully quantum mechanical results if the bath consists of harmonic degrees of freedom has been applied to bacteriochlorophyll (BChl) dimers⁶⁴ with highly accurate exciton–vibration coupling parameters obtained from fluorescence experiments.⁶⁵ More recently, the modular path integral (MPI) methodology^66–70 was used to simulate the EET dynamics in chains containing up to 19 BChl units as well as the B850 ring of the LH2 complex, treating explicitly the two electronic states and 50 normal mode intramolecular vibrations in each pigment.⁵⁰ Recent work developed a small matrix decomposition of the path integral^71–73 (SMatPI), which eliminates the storage needs of the QuAPI algorithm, allowing the simulation of larger systems and long-memory processes. The SMatPI algorithm offers an additional, very powerful tool for studying EET in molecular aggregates.¹⁰ 

In this paper, we report numerically exact simulations of EET dynamics between the B800 and B850 rings of the LH2 complex at room temperature, following photoexcitation of various eigenstates of the single excitation Hamiltonian. We treat the singly excited 8-unit B800 ring and the 16-unit B850 ring explicitly, including all state-to-state couplings and their interaction with all vibrational modes using the SMatPI methodology. Since the initial condition is given by an electronic eigenstate, the EET from the B800 to the B850 ring is entirely enabled by the motion of the nuclei; thus, the absence of dynamical approximations in the treatment of exciton–vibration coupling is crucial for obtaining reliable results.

In Sec. II we introduce the exciton–vibration Hamiltonian, along with its parameterization for the LH2 complex, and give a short overview of the SMatPI algorithm. In Sec. III, we present the results of our simulations by showing the populations of the exciton eigenstates that follow the excitation of four eigenstates of the exciton Hamiltonian with different outer vs inner ring compositions, along with significant coherences. When the initial state has density distributed primarily over the B800 ring, the relaxation process is found to involve an initial ultrafast (<0.1 ps) step of intra-B800 energy redistribution, followed by a considerably slower (∼1 ps) component associated with B800–B850 transfer. Interestingly, this diversity of vastly different exciton–vibration relaxation time scales is resolved to a smooth and perfectly exponential decay of the total B800 ring population, with an overall 1/e time of ∼0.5 ps, in good agreement with experimental findings. Only the very rapid component is observed if the main mixed character B850* state (an exciton with primarily B850 composition that lies within the B800 band) is initially excited, with almost no transfer of energy to the outer ring. We also obtain the equilibrium populations of the electronic eigenstates by evaluating the imaginary time path integral and discuss the anomalously large populations that the mixed character eigenstates exhibit. We conclude in Sec. IV with a summary of our findings and additional comments. In the  Appendix, we present detailed information about the convergence of the SMatPI calculations.

The B850 ring of Rhodopseudomonas molischianum comprises 16 chromophores. The B800 ring is made of half as many, i.e., n = 8 units. The crystal structure of these systems is well characterized, and the energies and couplings of the chromophores are available from electronic structure calculations.^5,74 Figure 1 shows the structure of the composite 16 + 8 chromophore system and the couplings among singly excited states.

We adopt the Frenkel exciton model⁷⁵ to construct the Hamiltonian and study the energy transfer dynamics. Each chromophore k = 1, …, n is treated as a pair of ground and excited electronic states, denoted $gk$ and $ek$⁠, with energies $εkg$ and $εke$⁠, which are coupled to vibrational degrees of freedom. We focus on the dynamics within the singly excited subspace of the aggregate⁷⁶ and designate as $k$ the state in which only monomer k is excited, i.e.,

etc. Within the normal mode approximation, the ground and excited electronic states of each BChl are described in terms of two diabatic harmonic potential surfaces,

and

where $q̂ik$ and $p̂ik$ are the positions and momenta of the normal modes in monomer k, respectively, and ω_ik and c_ik are the respective vibrational frequencies and exciton–vibration coupling strengths. The latter are related to the Huang–Rhys parameters S_ik and the dimensionless distance parameter s_e (which only rescales the exciton–vibration couplings based on the separation between the potential minima) through the expression

The total exciton–vibration Hamiltonian of the aggregate becomes⁴³ 

where

with

where $Jkk′$ is the electronic coupling between pigments k, k′. The electronic energy $Ek$ and the exciton–vibration term H_ev in Eq. (2.6) comprise the excited state terms for molecule k and the ground state terms for the rest of the chromophores, respectively. In Eq. (2.5), we also express the electronic Hamiltonian in terms of its eigenstates φ_i and eigenvalues E_i.

All electronic parameters were taken from the parameterization of the exciton Hamiltonian for the isolated LH2 complex in photosynthetic bacteria [shown in Fig. 1(c)] by Tretiak et al.⁷⁴ The magnitudes of couplings span a large range (3.2–363 cm⁻¹), where the strongest interactions exist within the B850 ring (363, 320, and 102 cm⁻¹) and the weakest are seen in the B800 (25–3.2 cm⁻¹). The strongest interactions between the chromophores of different rings (53 and 38 cm⁻¹) lie in between the two sets of intra-ring couplings. The larger B850 ring shows weak dimerization. The chromophores within each dimer are labeled as α and β, where chromophores labeled α have an identical excitation energy as the B800 molecules, while the excitation energy of those labeled β is lower by 82 cm⁻¹. The intradimer coupling (363 cm⁻¹) is larger than the coupling (320 cm⁻¹) between neighboring BChls that belong to two different dimers. Non-nearest neighbor interactions become progressively smaller with intermonomer separation, with the longest–range interactions spanning four monomers in the B850 ring and three monomers in the B800 ring. This distribution of electronic couplings and excitation energies are reflective of the subtle differences in local environments of pigments and create interplays of varied timescales. We discuss the eigenstates of the two-ring Hamiltonian in Sec. II C.

The vibrational parameters for each BChl are encoded in the spectral density function $Jω$⁠, which basically consists of 50 discrete vibrational modes that in a combined experimental and theoretical study⁶⁵ were found to couple strongly to the electronic states of an isolated BChl molecule. While the explicit inclusion of these discrete modes (shown in gray) using spectroscopic Huang–Rhys factors ensures a highly accurate treatment of the intramolecular vibrations of the pigment complexes, it does not capture the effects of the solvent or the encapsulating protein, which are responsible for irreversible dynamics. In this work, we treat the discrete molecular vibrations explicitly while adding a weakly coupled dissipative bath of the well-known Ohmic form,⁷⁷ 

with ω_c = 200 cm⁻¹ and ξ = 0.4 to capture the dissipative effects from the low frequency fluctuations of the environment. The background bath contributes a smaller reorganization energy (160 cm⁻¹), while the reorganization energy of the discrete modes is equal to 217 cm⁻¹. For consistency, we also broaden slightly the discrete components while maintaining the highly structured form and preserving the reorganization energy. The low frequency shape qualitatively resembles the vast diversity of environmental spectral densities^26,37,78,79 previously optimized for the LH2 and FMO complexes and has a minor effect on the subpicosecond dynamics of the BChl rings compared to the discrete vibrations. The overall spectral density, along with its components, is shown in Fig. 2(a). We also show in Fig. 2(b) the time correlation function of the total harmonic bath for each BChl at 300 K,

where β = 1/k_BT. The highly structured spectral density gives rise to persistent oscillations in the correlation function, which survive up to ∼1.5 ps.

The eigenstates of the exciton Hamiltonian and their energies reveal the range of the purely electronic time scales of EET. Figure 3 shows the energy eigenvalues, along with the spatial distribution of the eigenstate probability densities. Since by far the strongest couplings exist within the B850 ring, the lowest (φ₀ − φ₆) and the highest (φ₁₆ − φ₂₃) energy eigenstates are delocalized on the inner ring, forming the lower and upper B850 bands. Furthermore, since the dominant couplings are between nearest neighbors, the band splitting (1488 cm⁻¹) is not far from the 4J = 1452 cm⁻¹ value attained for a tight binding Hamiltonian in the n → ∞ limit. However, significant deviations from the tight-binding eigenvalue pattern appear due to the large range of coupling values and nonlocal interactions. For example, the two groups of B850 states have different spreads, with the higher band showing smaller spacings. In addition, because the B800 couplings are weak, the eigenstates (φ₇ − φ₁₅) that lie closest to the site energies correspond primarily to pigments on the B800 ring (forming the B800 band), with the exception of φ₁₃ (which has 78% of its density spread over the B850 ring) and φ₁₀ (which has 34% B850 character). The first of these states (and to some extent the second as well) is associated with the nearly dark B850* band.²⁵ The dynamical manifestations of the eigenstate manifolds and their interaction with vibrational modes will be investigated in Sec. III.

We seek the time evolution following photoexcitation of an eigenstate of the exciton Hamiltonian. Oscillator strengths for transitions to the various exciton states of the two-ring system are not available. However, it is instructive to investigate the EET process from several initial conditions. If photoexcitation produces eigenstate φ_i, the time evolution of the population of exciton eigenstate φ_f is given by the sum

where

is the reduced density matrix (RDM) in the “site” representation of singly excited BChl molecules. The trace in Eq. (2.10) is with respect to all degrees of freedom except the electronic states, $ρ̂b(0)$ denotes the initial density of the vibrational and background bath, and the subscripts m′, m″ and k′, k″ represent single BChl excited states. The off-diagonal elements of the RDM, or “coherences,” which develop from the evolution from an excited eigenstate, are given by

Within the Franck–Condon approximation, the vibrational degrees of freedom remain in the electronic ground state upon photoexcitation.

The RDM at the time t = NΔt, where Δt is the path integral time step, is calculated using the SMatPI decomposition of the QuAPI expression.⁸⁰ The latter has the form

where the path integral variables $ik±$ label the auxiliary states of the n-state Hamiltonian, which define paths in the space of electronic states, $Kik±,ik−1±$ are short-time forward-backward propagators for H_e, and

is the QuAPI-discretized⁸¹ influence functional,⁵⁸ which incorporates the effects of the vibrational degrees of freedom on the electronic system at the specified temperature. The influence functional contains nonlocal interactions that couple the path integral variables at different times. This temporal nonlocality or “memory” prevents the evaluation of Eq. (2.12) through step-by-step iteration,^82,83 commonly used to propagate a wavefunction or density matrix in the space of all degrees of freedom (which is practical only in very small systems). The QuAPI algorithm^59,60 relies on the decrease in nonlocal influence functional couplings (i.e., finite “memory”) to evaluate the RDM through iterative propagation of a tensor that spans the memory length. The SMatPI decomposition^71,72 eliminates the need for tensors, obtaining the RDM through a sum of small matrix products,

where L is the memory length in units of the path integral time step and M^(N,N−r) are n² × n² matrices. Equation (2.14) is an analytically derived, numerically exact decomposition of the RDM, which requires minimal storage that is comparable to that of the RDM of interest. Furthermore, the total cost of the SMatPI algorithm is only that of a single QuAPI propagation step. Finally, we note that the memory length L may exceed the entanglement length r_max of the path integral variables (which determines the cost of the SMatPI calculation),⁷³ further decreasing the cost of propagation. These features of the SMatPI methodology allow its application to the 24-state LH2 Hamiltonian with the highly structured spectral density of each BChl monomer. The RDM in the eigenstate representation is obtained from the computed values of $ρ̃$ and the eigenvector matrix of the electronic Hamiltonian, as shown in Eqs. (2.9) and (2.11).

The highly structured spectral density shown in Fig. 2 gives rise to persistent small-amplitude oscillations that decay extremely slowly. Despite this, the eigenstate populations shown in Sec. III are converged to within 0.01 or better. The convergence is discussed in detail in the  Appendix, which shows the variation of the populations with the memory length up to L = 200. Convergence is aided by the large size and a large number of vibrations of the LH2 complex, which lead to small fluctuations of the RDM elements with memory and a short entanglement length. Further verification of convergence was obtained by comparing to equilibrium values through imaginary-time path integral calculations on the full exciton–vibration Hamiltonian.

The transfer of energy among eigenstates of the exciton Hamiltonian is entirely mediated by exciton–vibration coupling. Even though the spectral density written in the site basis is identical for all BChl molecules, the interaction of exciton eigenstates with the vibrational modes is not uniform and, thus, is far from simple. As a consequence, the time scales and pathways of energy relaxation depend not only on the relation between electronic energy gaps and the BChl spectral density but also on the eigenstate components, giving rise to considerable complexity that is manifested in the highly nontrivial EET evolution obtained from the path integral calculations discussed below.

Figures 4 and 5 show the populations of the exciton eigenstates at T = 300 K, obtained with four different choices of initial preparation, which corresponds to photoexcitation to eigenstates φ₇, φ₁₀, φ₁₂, and φ₁₃. These states are chosen for their different compositions in terms of excited BChl participation. The first of these, state φ₇, is the lowest energy eigenstate that is delocalized over the entire B800 ring. Exciton state φ₁₂ is confined to the outer ring as well but has a density on only four excited pigments. The other two eigenstates, φ₁₀ and φ₁₃, are spread over both rings, with φ₁₃ composed mostly of excited B850 chromophores. As mentioned in Sec. II C, the mixed-character state φ₁₃ is the dark B850* state that is believed to play a special role in inter-ring energy transfer. With each initial condition, the populations of initially unexcited degenerate pairs are identical.

In the first three cases (initial conditions φ₇, φ₁₀, and φ₁₂), the population of the initially excited state shows two distinct time scales. A rapid decay is observed within ∼70 fs, which is followed by an abrupt transition to a much slower relaxation. This two-component behavior is in line with the observed decay pattern of the excited B800 state in various light harvesting complexes of purple bacteria.¹ Relaxation is achieved through the motion of the nuclei, which enables coupling among exciton states. Rapid transitions to other B800 states are mediated by the low-frequency components of the spectral density, which are resonant to intraband spacings. During the initial fast step, significant transient population is built on the other exciton eigenstates with B800 character (φ₇ − φ₉, φ₁₁, φ₁₂, φ₁₄, and φ₁₅), indicating some redistribution of excitation energy on the outer ring. A smaller population accumulates on state φ₁₀, and much less energy transfer to φ₁₃ is observed, consistent with the considerable and large, respectively, B850 components of these states. Therefore, a fraction of the excitation arrives to the B850 ring via the mixed character B850* excitons. The populated B800 states are seen to gradually decay over a longer time of ∼1 ps, transferring their energy to states of the lower B850 ring. We note that some uphill energy transfer (to eigenstates within the B800 cluster) also seems to be observed within this time interval. Transient bath-induced excitation (in addition to relaxation) has been seen in two-level systems.⁸⁴ 

A careful examination of the eigenstate populations within the B800 band reveals that the sequence of energy transfer does not always follow the order of energy eigenvalues. In particular, the two mixed states φ₁₀ and φ₁₃ are populated later than other close-lying states. As discussed earlier, the EET dynamics involves the complex interplay among exciton energy gaps, spectral density structure, and exciton–vibration coupling, which attains a very complex form in the basis of eigenstates. On the other hand, it is seen that the energy transfer patterns resulting from the initial excitation of eigenstates φ₇ and φ₁₂ are very similar, despite the rather different pigment composition of these states (delocalized vs spread over four monomers) and the fact that φ₁₂ belongs to a degenerate pair. Additional insights for these observations are offered in Sec. III C, where we examine the equilibrium values of the eigenstate populations dressed by exciton–vibration mixing.

Following the rapid redistribution of excitation energy among the eigenstates of the B800 band, inter-ring exciton relaxation to B850 eigenstates is observed. The energy cascade within the inner ring eigenstate manifold generates nearly monotonically increasing populations over a much longer time scale of 0.5–1 ps.

The EET dynamics is seen to be quite different when the initial excitation corresponds to the mixed B850* state φ₁₃. The second, slow relaxation component is practically absent in this case, as expected for a state that is composed primarily of inner ring BChl excitations. The excitation energy is now directly transferred to the lower B850 states with minor transient population accumulation. Furthermore, the exciton relaxation is considerably faster with this initial condition, as the slow ring-to-ring transfer step is largely avoided. The relaxation of the other mixed character state φ₁₀ exhibits some of these dynamical features as well, but to a smaller extent, because of its smaller B850 component. We revisit these patterns in Sec. III C, with a discussion from the perspective of energetics.

These behaviors help elucidate the role of the B850* state φ₁₃ (and to some extent of φ₁₀) in the EET pathways. Starting from a B800 excitation, the energy is rapidly transferred to the other close lying exciton states, including φ₁₃, and subsequently transferred to the lower B850 states. The EET to the B850 band out of outer ring states is slow, but it is very fast starting from the mixed B850* state. In the latter case, the rapid redistribution of the energy on the B850 ring depletes the B850* state(s), causing only minor population accumulation above the long-time equilibrium value. Thus, the B850* states provide an additional pathway for exciton relaxation and, in turn, for B800-to-B850 EET.

Although the time axis in the panels that show the individual state populations extends to 1.5 ps for visual clarity, we note that we have computed the populations up to 5 ps, a time interval sufficiently long for all populations to have settled to their equilibrium values. The propagated RDM preserves the norm to four decimal places and the calculated long-time populations are the same for all initial conditions (and also consistent with expected room temperature Boltzmann occupancies), further verifying the accuracy of the SMatPI results. A systematic analysis of convergence is given in the  Appendix.

Also shown in Figs. 4 and 5 are the total populations of the two rings, obtained by adding the populations of the corresponding excited monomer states. It is seen that 90% of the energy of eigenstates with primarily B800 composition is eventually transferred to the B850 ring and that the transfer is almost complete within 1 ps. Even though the population of the initially excited state is characterized by an ultrafast and a much slower component, the total ring populations exhibit exponential dynamics. This is seen in Fig. 6, which shows the population of the B800 ring for the two initial conditions shown in Fig. 4 (eigenstates φ₇ and φ₁₂) relative to its equilibrium value on a logarithmic scale. In both cases, the 1/e time of 0.5 ps obtained from the path integral results is in very good agreement with experimental observations of 0.6–1 ps at room temperature for various LH2 complexes.^1,22,85–89 Given the highly accurate BChl Huang–Rhys factors and the numerically exact nature of the path integral calculations, this agreement also provides strong evidence for the accuracy of the electronic coupling parameters reported in Ref. 74. When the initially prepared eigenstate has a substantial B850 character, equilibration is achieved faster.

Additional insights are offered by studying the off-diagonal elements of the RDM, which are often referred to as “coherences.” The bulk of these have values close to zero during the entire time evolution, while others show interesting time dependence. Figure 7 shows some of the coherences $Qijm(t)$ with the largest values for the initial states considered above, and they also exhibit the two time scales identified in the population dynamics. While the largest values and rapid, oscillatory evolution are observed within the initial 0.1 ps, Fig. 7 shows that the dynamics of these coherences lasts for the entire time of population transfer, i.e., ∼1 ps. This suggests that the evolution creates long-lived coherent superpositions of some eigenstates, primarily of eigenstates that are delocalized over one of the two rings.

We point out that the real part of the coherences does not approach zero at equilibrium. This is because the eigenstates of the exciton–vibration Hamiltonian are linear combinations of direct product states. These superpositions are reflected in the long-time values shown in Fig. 7.

The strong mixing of exciton and vibrational states in the equilibrium density matrix of the LH2 is further seen by comparing the equilibrium eigenstate populations of the LH2 complex $pieq$ predicted from the bare electronic Hamiltonian,

to those obtained from the composite exciton–vibration Hamiltonian,

which we evaluated using the imaginary time path integral formulation of the Boltzmann operator.^90–92 We also calculate vibration-dressed energy eigenvalues $Ẽi$ by writing

Figures 8 and 9 show the bare and vibration-dressed equilibrium populations and eigenvalues of the 24 electronic eigenstates. The value for the lowest state, $P0eq=0.32$⁠, is in excellent agreement with the long-time values of the converged dynamical results, while the population obtained from the bare electronic Hamiltonian gives $p0eq=0.41$⁠, which is in error by 28%. This large change is a consequence of the cumulative vibrational reorganization of 24 BChl pigments in this complex. Figure 8 shows that while $Pieq<pieq$ for i = 0, 1, 2, the trend reverses for higher lying states (as required to satisfy normalization). This implies that the spread of the vibration-dressed energies is smaller than the spread of electronic energies. This result is not unexpected, based on the well-known two-level system behavior where high-frequency bath modes reduce the tunneling splitting, leading to slower system-bath dynamics.⁹³ However, while the eigenstate populations that correspond to the bare electronic Hamiltonian decrease with an increasing quantum number, Fig. 8 shows that the mixed B850* states have higher population values than those preceding them in sequence, violating the expected trend. In particular, as seen in Fig. 9, the vibrationally dressed level $Ẽ10$ is shifted to the lower edge of the B800 band, and the level $Ẽ13$ is shifted even lower, within the range of the B850 states.

These results explain the relaxation cascade anomalies identified in Figs. 4 and 5. Starting with the B800 states φ₇ or φ₁₂, ultrafast vibrational redistribution transforms these to the dressed states, which are densely packed and give rise to rapid intraband electronic relaxation to close lying states, which lie either below or above. However, the dressed mixed states $φ̃10$ and $φ̃13$ that lie lower are separated by larger energy gaps and, thus, are populated later. On the other hand, when the initial excitation is placed on $φ̃13$⁠, the higher lying B800 states are minimally populated; thus, the slow relaxation process from the B800 to the B850 states is largely circumvented.

To summarize this section, we note again that when the composite system is excited to an eigenstate of the electronic Hamiltonian, no evolution (thus no energy transfer) would occur in the absence of vibrational motion. Exciton–vibration coupling with a variety of energy and time scales gives rise to relaxation and the transfer of energy to the B850 ring. The specifics of the mechanistic pathways for exciton relaxation and EET within the LH2 complex are complex and intricate, exhibiting rich facets of molecular exciton–vibration dynamics. We found that the energy relaxation cascade presents apparent anomalies that correspond to the mixed B850* states. However, these anomalies are explained by realizing that exciton–vibration coupling shifts these states to lower energies, close to the top of the lower B850 band.

In this paper, we reported fully quantum mechanical path integral simulation results for photosynthetic energy transfer between the B800 and B850 rings of the bacterial LH2 complex. Using a numerically exact path integral method, we avoided dynamical approximations and evaluated the real-time dynamics at room temperature within a normal mode treatment of all 50 coupled vibrational modes for each of the 24 BChl molecules, accounting for all coupling terms among chromophore excited states. The effects of the biomolecular environment were included using a dissipative bath of a form commonly used in open quantum systems. To the best of our knowledge, this simulation is the first to report fully quantum mechanical results in the two-ring LH2 system. The simulations allow the full characterization of the pathways, time scales, and mechanistic details of inter-ring energy transfer and present excellent agreement with available experimental results.

The results of our path integral simulations are summarized in Figure 10, which shows the distribution of excitation energy among the electronic eigenstates as a function of time. We found that the relaxation of excitation energy in the light harvesting complex from the B800 to the lower B850 band, occurs over approximately 1ps, and is characterized by different time scales and complex patterns that are not reducible to a simple picture constructed in terms of excitonic energy gaps. The exciton–vibration coupling, although quite simple in mathematical form for each BChl unit, is manifested in the state-to-state coupling in intricate ways, leading to intriguing relaxation pathways. When the initial excitation populates eigenstates localized within the outer ring, two sets of cascading relaxation time scales are identified. The faster components (∼70 fs) involve transient population of B800 band states, while a subsequent, considerably slower progression of transfer events is observed on a timescale of 1 ps, during which eigenstates of the inner ring are populated. One of the pathways for fast inter-ring energy transfer relies on the B850* states (of mixed B800–B850 character), whose electronic energies lie within the B800 band, although exciton–vibration coupling shifts these states to lower energies. A B850* excitation leads to even faster (<70 fs) energy redistribution within the B800 ring, and the slow B800-to-B850 energy-transferring component is not observed in this case. Despite the different times pertaining to the two cascading progressions, the overall B800-to-B850 energy transfer dynamics is perfectly exponential, characterized by a 0.5 ps 1/e time.

While it is evident that transitions between exciton eigenstates are possible only through the motion of the nuclei, the highly complex interplay among 24 coupled electronic states and 50 intramolecular modes in each BChl cannot be understood without resorting to computational treatments. This is further substantiated by the large deviation of long-time thermal populations of exciton states from their Boltzmann populations predicted by electronic eigenvalues and the presence of anomalies, which suggest a reordering of electronic eigenvalues due to significant exciton–vibration mixing. The fully quantum mechanical real-time path integral simulations shed light on the exciton–vibration dynamics in the light harvesting process of bacterial photosynthesis with unprecedented accuracy and clarity.

This work was supported by the National Science Foundation (Award No. CHE-1955302). Some of this research was part of the Blue Waters sustained-petascale computing project, which is supported by the National Science Foundation (Award Nos. OCI-0725070 and ACI-1238993) and the state of Illinois. Blue Waters is a joint effort of the University of Illinois at Urbana-Champaign and its National Center for Supercomputing Applications.

The authors have no conflicts to disclose.

The data that support the findings of this study are available from the corresponding author upon reasonable request.

Here, we show the detailed convergence of the SMatPI calculations. The structured BChl spectral density gives rise to a correlation function with persisting oscillations [Fig. 2(b)], leading to very long memory. Because of the continuous spectral background shown in Eq. (2.7), the oscillations in the correlation function eventually (at 1.5 ps) decay to zero.

We have investigated the convergence in several hypothetical two-ring complexes of varying sizes with identical couplings between BChl units. The “8 + 4” complex, comprising 8 pigments in the inner ring and 4 in the outer, is the smallest structure that contains all interactions present in the natural “16 + 8” LH2 complex, yet contains half as many excited states, allowing more efficient and systematic convergence testing.

Figure 11 shows SMatPI results for the populations of eigenstates φ₀ and φ₁ when the initially excited state is φ₆ (which has mixed character, thus is analogous to φ₁₀ and φ₁₃ of the 16 + 8 complex). In each panel, we show results obtained by varying the memory parameter L for fixed entanglement length r_max. With r_max = 3, increasing the memory length at first leads to a significant monotonic change in the populations. However, as the memory length is further increased, the eigenstate populations begin to vary in a nonmonotonic fashion, remaining confined within a narrow band that shrinks as L is increased. As shown in Fig. 11, the population of φ₀ oscillates with a maximum (long-time) spread of 0.02 when L varies in the range of 10–20, which decreases to 0.005 for 100 < L < 110. Upon increasing L to 200, the population curves become identical to four decimal places. This behavior is a consequence of the highly structured form of the BChl spectral density, which creates persisting recurrences in the correlation function of slowly decreasing amplitude. Similar observations can be made for the population of eigenstate φ₁. Because of the lower population of this state, the fluctuations as a function of memory length are even smaller. The populations of other eigenstates are very small and their convergence is even better at all times.

Similar trends are obtained with r_max = 4 and 5. In both cases, the populations obtained with sufficiently large L lie within narrow bands of comparable width. As shown in Fig. 11, SMatPI results with r_max = L = 4 and 5 lie within these bands as well. The large L correction to the r_max = 5 population is smaller than 0.01. This shows that while the memory of the BChl vibrational bath is very long, the entanglement is quite short and, thus, is well within the capabilities of the SMatPI algorithm. We note that these favorable trends, which allowed us to obtain converged results in the presence of the persisting oscillatory character of the spectral density, can be attributed to the large size of this system.

As a further confirmation of convergence, we have computed equilibrium populations by evaluating the Boltzmann matrix elements for the composite electronic–vibrational Hamiltonian using the imaginary-time path integral formulation.^90–92 The equilibrium populations, shown as dashed lines in Fig. 11, are seen to be in excellent agreement with the converged results. For comparison, we also show equilibrium populations obtained in the absence of vibrational degrees of freedom. In line with the observations discussed in Sec. III C, it is seen that the inclusion of exciton–vibration coupling modifies the equilibrium values of the state populations in a profound way, altering the value of the lowest eigenstate by ∼28% here as well. The two mixed B800–B850 states in this smaller cluster are φ₆ and φ₇. Again, these two states have anomalously large populations through exciton–vibration coupling. While their equilibrium populations are small (<0.03), the dynamical SMatPI results capture the population ordering correctly.

The convergence of the SMatPI results is even better in the larger 16 + 8 complex. In this case, the largest population (that of φ₀) has a considerably smaller value because there are twice as many eigenstates within approximately the same energy band. Figure 12 shows that the spread of results obtained with L = 10 −20 is significantly smaller in this case.