Absorption and Circular Dichroism Spectra of Molecular Aggregates With the Full Cumulant Expansion.

The exciton Hamiltonian of multichromophoric aggregates can be probed by spectroscopic techniques such as linear absorption and circular dichroism. To compare calculated Hamiltonians to experiments, a lineshape theory is needed, which takes into account the coupling of the excitons with inter- and intramolecular vibrations. This coupling is normally introduced in a perturbative way through the cumulant expansion formalism and further approximated by assuming a Markovian exciton dynamics, for example with the modified Redfield theory. Here, we present the implementation of the full cumulant expansion (FCE) formalism ( J. Chem. Phys. 142, 2015, 094106) to efficiently compute absorption and circular dichroism spectra of molecular aggregates beyond the Markov approximation, without restrictions on the form of exciton-phonon coupling. By employing the LH2 system of purple bacteria as a challenging test case, we compare the FCE lineshapes with the Markovian lineshapes obtained with the modified Redfield theory, showing that the latter presents a less satisfying agreement with experiments. The FCE approach instead accurately describes the lineshapes, especially in the vibronic sideband of the B800 peak. We envision that the FCE approach will become a valuable tool for accurately comparing model exciton Hamiltonians with optical spectroscopy experiments.

Introduction

The optical and photophysical properties of chromophore aggregates such as light-harvesting complexes are determined by the interactions between the monomer excitations and by the coupling between these excitations and the vibrational degrees of freedom.[1,2] These interactions are quantified, respectively, by the exciton couplings between transitions localized on individual chromophores and by the spectral density of the exciton–phonon coupling. Together, these interactions determine the excited-state dynamics of light-harvesting systems.

Excitonic and vibronic couplings can be parameterized by fitting the linear and nonlinear spectra of the aggregates, guided by prior knowledge of their structure.[3−6] Detailed high-resolution structures of chromophoric aggregates have also enabled the use of quantum chemistry (QM) methods to compute a priori both excitonic Hamiltonians and vibronic couplings.[1,7−14] However, these quantities still need to be compared with experiments, by simulating the aggregate spectra.

While the spectrum of a single chromophore can be evaluated exactly starting from the corresponding spectral density,[15] the mixing of electronic states resulting from the exciton couplings prevents a simple calculation of the overall spectrum without further approximation. In the delocalized basis of electronic eigenstates, the exciton–phonon coupling is not diagonal anymore; the off-diagonal exciton–phonon coupling cannot be treated exactly. Numerically exact methods, such as the hierarchical equations of motion (HEOM)[16−19] and the stochastic path integral (sPI),[20] are computationally very expensive and limited to benchmark calculations on model systems. On the other hand, perturbative approaches based on Redfield equations have been extensively used,[12,13,21,22] but their applicability range is limited.[23] The main concern about the Redfield-like approaches arises from the Markovian and secular approximations; in particular, the Markovian approximation is known to break down when vibrational modes are in resonance with electronic energy gaps.[23] Several strategies have been devised to overcome the Markov/secular approximations.[23−32] Among these, the full second-order cumulant expansion (FCE) offers the most general form of non-Markov/secular approach within perturbative treatment of the off-diagonal exciton–phonon coupling.[28] For linear absorption properties, the FCE is equivalent to the second-order perturbative time-convolutionless (TCL2) quantum master equation.[28] Despite the approximation to the real non-Markovian dynamics, the TCL2 method performs well for linear absorption spectra, although it is qualitatively inaccurate for nonlinear spectra.[33]

In this contribution, we implement the FCE approach to compute absorption and circular dichroism (CD) lineshapes of molecular aggregates, without restrictions on the functional form or the site-dependence of spectral densities. We show an application of the method to the LH2 system of purple bacteria,[34,35] using excitonic parameters[10] and spectral densities[13] recently calculated by some of us. We demonstrate that a Redfield-based approach suffers some limitations in describing the absorption lineshape of LH2, while the FCE approach yields an excellent agreement with the experiments, offering a new method to compare QM exciton calculations with absorption and circular dichroism spectra.

Methods

Hamiltonian

The total Hamiltonian of the excitonic aggregate, , can be written as , where the electronic part is the exciton Hamiltonianwhere is the vertical excitation energy of chromophore i, and V is the electronic coupling between the excited states of chromophores i and j. The eigenstates μ of the exciton Hamiltonian are expressed in the site basis aswhere the excitonic coefficients cμ = ⟨j|μ⟩ are the eigenvectors of the exciton matrix.

The bath correlation function C of chromophore n is determined by the spectral density

Absorption and Circular Dichroism

The absorption spectrum is given by the half-sided Fourier transform of the transition dipole correlation function.[26] In the site basiswhere ρg is the ground-state density matrix of the bath and the i, j indices run on the chromophores of the system. Applying the Condon approximation, the transition dipoles can be taken out of the bath tracewhere I(t) is the absorption tensor of the aggregate[28]and M = · is the dipole strength matrix.

In the same way, the CD spectrum is related to the electric dipole–magnetic dipole correlation function[26]where denotes the imaginary part. Neglecting the intrinsic magnetic moment of each chromophore, m can be written aswhere and R are the excitation energy and the position, respectively, of chromophore j. One can then define the rotatory strength matrix[26]and express the CD spectrum asBoth the absorption and the CD lineshapes are determined by the absorption tensor I(t), which will be evaluated using FCE. The lineshape tensor can be transformed into the exciton basis aswhere cμ are the coefficients of the exciton states. In the same way, the dipole strength matrix and rotatory strength matrix can be transformed into an exciton basis. However, the matrix r given in eq is nonsymmetric and gauge-dependent. A gauge-invariant expression can be obtained by substituting the excitation energy with the geometric average and symmetrizing the matrix

Full Cumulant Expansion for Arbitrary Spectral Densities

The absorption tensor I(t) can be written in the exciton basis via the full second-order cumulant expansion given in ref (28) aswhere the FCE lineshape matrix K is defined asHere, n denotes a site; μ, ν, and α refer to exciton states; Xμα = (cμ)*cα; ℏωμα = Eμ – Eα; and C(t) is the autocorrelation function of the energy gap.

The numerical integration in eq can be made easier by separating the case ωμν = 0. First, we evaluate Fμαν(t) by integrating by partsIn the case where ωμν = 0, we cannot separate the integral as done below in eq , but noting thatwe reformulate eq as (μ = ν)where we have defined the auxiliary tensors Gνα(t) = ∫0 eC(t2)dt2 and Hμα(t) = ∫0t2eC(t2) dt2 Next, we consider the case of ωμν ≠ 0. In this case, the integral can be separated into two terms

Numerical Implementation

The calculation of Kμν(t) is the computationally most intensive part of the FCE evaluation, especially when the number of distinct spectral densities increases. At the beginning of the calculation, exciton weights Wμα = XμαXαν are computed and summed over the chromophores that have the same autocorrelation function C(t)At each time step, only the auxiliary tensors Gμα and Hμα are propagated from their previous values. In this implementation, we use the trapezoid ruleThe two tensors are allocated as Nsite2 × NSD, where NSD is the number of unique spectral densities in the calculation. The tensor Fμαν is computed with eq or 17 when ωμν is smaller than a threshold.

Finally, the lineshape matrix is computed, by summing over the NSD different spectral densities, givingThe numerical calculations described in this Section were implemented in a Fortran 90 code available for download under the LGPL license agreement.[36]

Results

The LH2 system of purple bacteria[34,35] provides an ideal test case for the FCE approach. LH2 comprises 27 bacteriochlorophyll a (BChl) pigments, which give rise to a typical absorption spectrum in the near infrared. The BChls in LH2 are organized into two rings, responsible for the absorption peaks at 800 nm (B800 ring) and at ∼850 nm (B850 ring), respectively. The B850 ring contains nine αβdimers of tightly packed BChls, whose strong nearest-neighbor couplings delocalize the excitation and shift the peak frequency to the red. Conversely, the BChls in the B800 ring are more separated and weakly coupled, though the effect of this coupling on the band shape is not negligible.[3,13] LH2 has been extensively investigated theoretically, with both modeling studies[3,6,37−41] and quantum chemical calculations,[10,13,42−44] and experimentally, by a variety of spectroscopic techniques.[45−50]

To accurately compare the FCE lineshapes to experiments on LH2, we need reliable parameters for the exciton and exciton–phonon couplings. Some of us have employed a polarizable multiscale quantum chemical model to accurately compute the site energies and excitonic couplings for LH2.[10] We have shown that excitonic parameters computed on the crystal structure (@CRY) represent the exciton structure of LH2 at very low temperature, whereas the parameters computed along a molecular dynamics simulation (@MD) represent LH2 at room temperature. In this work, we will employ the two sets of parameters from ref (10) to describe the spectral lineshapes at different temperatures.

We model the exciton–phonon couplings of the BChls in LH2 with a different spectral density for each chlorophyll type. We take the spectral densities calculated in ref (13) for the three different BChls in LH2. These are comprised of a high-frequency part, directly calculated through a normal-mode analysis of the BChl in the protein, and a low-frequency part, modeled as a Drude oscillatorwhere j denotes the chromophore, and k denotes the vibrational modes of the chromophore, λc is the low-frequency contribution to the reorganization energy, and γc is the damping parameter of the overdamped Brownian oscillator. In ref (13), parameters λc and γc and the width σ of the static disorder were empirically adapted. The effect of static disorder was modeled by averaging the spectra over N realizations of the disordered exciton Hamiltonian, obtained by randomly choosing the site energies from a Gaussian distribution with standard deviation σ. The number N of realizations was deemed sufficient when the average spectrum of N/2 realizations was indistinguishable from the average of the remaining N/2. We thus used 5000 realizations for the 77 K spectra and 500 realizations for the 300 K spectra. The details of spectral densities and static disorder are reported in the Supporting Information.

Absorption Spectra

We first compare the LH2 absorption spectra with those obtained in ref (13) employing the modified Redfield (mR) equation. To do this, we compare in Figure the spectra obtained at 77 K from the same exciton Hamiltonian (@CRY) and disorder parameters. A noticeable shift is present, for both absorption bands, between the two predicted spectra due to the neglect of the off-diagonal reorganization energy in the mR theory. In particular, the mR equation neglects the imaginary part of the off-diagonal exciton–phonon coupling.[23] As already found in ref (23) for a model system, this off-diagonal reorganization shift can be different for different exciton states. In our case, the reorganization energy is different in the two bands because (i) we explicitly used different spectral densities for the B850 and B800 BChls, and (ii) the excitation is much more delocalized in the B850 ring than in the B800 ring. The difference in reorganization energies is reflected in the apparent energy gap between the B800 and B850 peaks, which is ∼60 cm–1 narrower in the FCE spectrum.

Figure 1

Absorption spectrum of LH2 at T = 77 K, obtained using FCE (red) and modified Redfield theory (blue), with the @CRY calculated exciton parameters (see text). The arrows indicate the energy gap in cm–1 between the two absorption maxima.

In Figure , we separately show the B850 and B800 band shapes, comparing the mR theory and FCE with the experiment.[5] We shifted and normalized all spectra on the experimental maxima to facilitate the shape comparison. Apart from the reorganization energy shift, the shape of the B850 band as predicted by mR theory is very similar to the one predicted by FCE, pointing to a small influence of non-Markovian effects in this band. In comparison, mR theory predicts a substantially symmetric shape for the B800 band, whereas the FCE band shape is strongly asymmetric. The origin of this asymmetry has already been debated and assigned to dipole strength redistribution due to exciton interactions and static disorder.[3,37] Here, using a more realistic description of exciton and exciton–phonon couplings, we find that non-Markovian effects, which are absent in the mR description but included in the FCE, enhance the asymmetry of the B800 band.

Figure 2

Absorption spectrum of LH2 at T = 77 K, obtained using FCE (red) and modified Redfield theory (blue) with the @CRY calculated exciton parameters (see text) and compared with the experiment at 77 K (dashed black line): (a) B850 band and (b) B800 band. Computed spectra are shifted to the experiment absorption maximum and normalized.

To better show the non-Markovian effects in the B850 and B800 bands, we reported in Figure S1 the absorption spectra from two sample realizations of the static disorder, computed with FCE and mR theory. The stronger static disorder of the B850 band gives rise to lineshapes with several narrow peaks, due to the splitting and mixing of the lowest exciton levels. These narrow peaks are similarly described by mR theory and FCE. On the other hand, the B800 band is still more asymmetrical with FCE than with mR theory, confirming the conclusions drawn from the averaged spectra.

As for room-temperature spectra, we find an even larger difference between the mR and FCE lineshapes (Figure ). In addition to the frequency shift noticed above, the FCE spectra present a narrower B800 band, with a pronounced vibronic structure. By contrast, the mR lineshapes show noticeable Lorentzian broadening at both the red and blue tails of the spectrum, due to the interexciton relaxation rates that are overestimated by the mR theory.

Figure 3

Absorption spectrum of LH2 at T = 300 K, obtained with the @MD calculated exciton parameters (see text) using FCE (red) and modified Redfield theory (blue). (a) Disorder parameters as in ref (13) and (b) disorder parameters slightly modified (σ = 270 cm–1, λc = 140 cm–1 for B850, and λc = 40 cm–1 for B800).

In the right panel of Figure , we present the same comparison, changing the width σ of static disorder for the B850 BChls, and the magnitude of the low-frequency component of the spectral density, λc (the B800 static disorder was kept constant). Now, the B850 lineshapes are more similar, but the mR spectrum still presents Lorentzian tails. Overall, the non-Markovian effects seem much stronger for room-temperature spectra, and the mR results strongly deviate from the FCE ones.

Due to the large difference between FCE and mR spectra, for simplicity, we compare only the FCE spectra with the experiments. For this comparison, we use the disorder parameters adopted in this work. In Figure , we show the B850 and B800 FCE lineshapes along with the experiment at 300 K.[10] While the overall broadening of the B850 band is reproduced, the asymmetry is not. The source of this asymmetry could be traced to the static disorder distribution and to the coupling of BChl states to charge-transfer states.[51] On the other hand, the shape of the B800 band is almost perfectly reproduced by the FCE approach. In particular, the vibronic tail, which could not be reproduced by the mR theory, fits the experimental spectrum almost perfectly.

Figure 4

Absorption spectrum of LH2 at T = 300 K, obtained with the @MD calculated exciton parameters (see text) using FCE (red), compared with the experiment at 300 K (dashed black line). Left: B850 band; Right: B800 band. The disorder parameters are σ = 270 cm–1, λc = 140 cm–1 for B850; σ = 40 cm–1 and λc = 40 cm–1 for B800. Computed spectra are shifted to the experiment absorption maximum and normalized.

In Figure S4 (Supporting Information), we show the dissection of the LH2 lineshape into B850 and B800 contributions. The vibronic sidebands are clearly visible in the B800 spectrum and reflect in the total LH2 lineshape. In addition, there is a non-negligible contribution of B850 pigments also in the B800 region, which slightly alters the B800 lineshape in the full spectrum. To characterize the effect of the weak exciton coupling on the B800 lineshape, we compare in Figure S5 the B800 lineshape with the monomer lineshape obtained with the same spectral density and static disorder. The exciton coupling gives rise to a small but noticeable red shift and to a more asymmetrical shape of the main absorption band. In addition, a small intensity redistribution from the vibronic tail to the main peak can be noticed.

Circular Dichroism Spectra

The CD spectra of LH2 computed at 77 K and at 300 K with FCE are reported in Figure , along with their experimental counterparts. Despite the excellent agreement in the absorption lineshapes with the experiment at both temperatures, the computed CD spectrum at 300 K deviates considerably from the experiments, especially in the B850 couplet. It was already observed in ref (10) that the intensity of the negative B850 band was underestimated by the calculations. Indeed, the sign and magnitude of the B850 couplet were shown to be extremely sensitive to small changes in the orientation of the transition dipole moments of the B850 BChls.

Figure 5

Circular dichroism spectra of LH2 computed using FCE (red) and compared with the experimental spectra (dashed black line): (a) spectrum at 77 K from the @CRY data; the computed spectrum was shifted by 1403 cm–1 to match the position of the B800 absorption and (b) spectrum at 300 K from the @MD data; the computed spectrum was shifted by −343 cm–1 to match the position of the B800 absorption. All spectra were normalized to one at the negative B800 peak.

Overall, the qualitative shape of the CD spectrum at 77 K is similar to that obtained in ref (10) using only Lorentzian broadening. On the contrary, in the 300 K spectrum computed here, the B800 couplet is somewhat narrower than the prediction from MD calculations in ref (10). It should be noted, however, that here we are employing a frozen arrangement of the chromophores, as calculated on the crystal structure, thus neglecting the fluctuations of couplings and transition dipoles, which affected the spectra predicted in ref (10).

Discussion and Conclusions

We have shown that the inclusion of non-Markovian effects and reorganization energy shifts through FCE theory has a considerable impact on the two exciton bands of LH2. The most striking differences are seen in the B800 band, for which the Redfield-based theories are known to have limitations arising from the near-degeneracy of exciton states.[3] This effect is less pronounced in the B850 band due to the larger couplings and static disorder, which increase the separation between exciton energies. The reorganization energy shifts, on the other hand, can be easily taken into account also in Redfield-like theories, by including the imaginary part of the relaxation tensor.[4,27]

The FCE approach yields an excellent agreement with the experiment for the B800 absorption, including its vibronic tail. This tail arises from the discrete part of the spectral density in eq , which in our case was calculated by means of multiscale QM methods.[13] The spectral densities employed in this work present several peaks in the region at around 750 and 1200 cm–1 and over 1500 cm–1, which contribute to the vibronic tail. The strong underestimation of this tail by mR theory can be traced back to the Markov treatment of all off-diagonal exciton–phonon interactions.[20,23] In fact, the mR theory treats only the diagonal part of the exciton–phonon coupling in a time-dependent fashion, whereas the off-diagonal part is reduced to a simple lifetime symmetric broadening.

It has been suggested that nonsecular/non-Markov theories take into account localization effects, as they tend to reduce the exciton delocalization effects predicted by Markov/secular Redfield theory.[26] In particular, these localization effects are expected to reduce the intensity of the CD bands, as suggested in ref (26) for a chlorophyll dimer. To verify this hypothesis, we compared the CD spectrum calculated at room temperature by FCE and mR theories (Figure S2 in the Supporting Information). On the one hand, the B800 couplet is more intense for FCE than for mR theory, contrary to the expectations. We attribute this effect to the narrower overall shape predicted by FCE for the B800 band. On the other hand, a larger asymmetry in CD intensities is found for the mR lineshape between the B850 and B800 couplets. The B850 band has an overall more positive CD, and the B800 band has a net negative CD. In the FCE lineshape, this asymmetry is reduced, which points to reduced delocalization between the B850 and B800 rings by including non-Markov and nonsecular terms.

We have computed the spectra of LH2 at different temperatures using distinct exciton Hamiltonians, which have been derived in ref (10) from crystal structure calculations (@CRY) and molecular dynamics calculations (@MD). Indeed, there is spectroscopic evidence that the nearest-neighbor couplings in the B850 ring reduce when increasing the temperature.[47,52] We also assessed how the LH2 spectra change between 77 and 300 K without changing the exciton Hamiltonian. To this end, we computed the LH2 spectra at 77 and 300 K using only the @CRY exciton Hamiltonian (Figure S3a,b in the Supporting Information). By keeping the exciton Hamiltonian fixed, it is not possible to reproduce the blue shift of the B850 band, which is, however, apparent in the experiments (Figure S3c,d). Therefore, we conclude that two different exciton Hamiltonians are appropriate for low and room temperature.

Calculation of lineshapes in light-harvesting complexes has been an active research area. Theoretical methods are primarily based on second-order quantum master equations, including various versions of Redfield equations. In this context, the FCE formalism is the most basic second-order form, from which various Redfield equations can be derived with further approximations. For example, one subtle issue is the time order of the system–bath interaction operator, which is strictly imposed in the cumulant form adopted in FCE but approximately obeyed in the differential form adopted in TCL2 or TC2 quantum master equations.[53] In addition to accuracy, FCE also provides a convenient starting point to calculate emission spectra and the Forster energy transfer rate in molecular aggregates. However, we note that calculating emission spectra presents additional challenges related to the entanglement between the system and the bath created when the former relaxes in the excited state, assuming that a factorized initial condition leads to a qualitative difference from the exact results.[28] A comparison of several second-order approximations and FCE has been reported in ref (28) both analytically and numerically. This study further establishes its validity using the important light-harvesting system of LH2 and extends the FCE method to the calculation of CD spectra.

Correct treatment of the absorption lineshape is fundamental to assess the quality of a computed exciton Hamiltonian by comparing optical spectra with experiments. We have shown that, for LH2, the apparent frequency gap between the B800 and B850 bands is strongly influenced by the approximations introduced by the mR theory. The same effect might also appear in other exciton systems, where different exciton states have different degrees of localization. In the major and minor antenna complexes of higher plants, for example, there are pigment clusters in which excitation is delocalized over two or three chromophores and other more isolated pigments where the excitation is virtually localized.[54−56] We envision that a more accurate lineshape theory such as the FCE will help refine the exciton Hamiltonian models of these antenna complexes.

From our results, it appears that the modified Redfield approach introduces a bias in the estimation of lineshape, which might reflect in an incorrect assessment of exciton Hamiltonian parameters when comparing optical spectra with experiments. Gelzinis et al. have reached similar conclusions by studying a model dimer with different lineshape theories.[23] They noticed that Redfield-like approaches have been used to extract Hamiltonian parameters from experimental spectra and argued that these fitted parameters should be reassessed with care. Indeed, lineshapes have been customarily modeled with Redfield-like approaches, both for directly fitting site energies and couplings[6,22,47,56] and for refining exciton Hamiltonians calculated by atomistic methods.[54]

In conclusion, we have presented the implementation of the full cumulant expansion approach to computing the absorption and CD lineshapes of molecular aggregates, overcoming the usual Markovian and secular approximations made in the widely employed modified Redfield theory. Our implementation does not impose any restriction on the functional form of the spectral density, allowing efficient computation of linear spectra for systems with more than 20 states. As an example, we have shown the application of the FCE method to the LH2 antenna of purple bacteria, demonstrating that the FCE overcomes the limitations of the modified Redfield theory and better reproduces absorption lineshapes, including vibronic sidebands. The FCE approach will represent a valuable tool to achieve a quantitative comparison between calculated exciton Hamiltonians and experimental optical spectra.