Excitonic Energy Transfer of Cryptophyte Phycocyanin 645 Complex in Physiological Temperature by Reduced Hierarchical Equation of Motion

Abstract

Recently, many researches have shown that even photosynthetic light-harvesting pigment-protein complexes can have quantum coherence in their excitonic energy transfer at cryogenic and physiological temperatures. Because the protein supplies such noisy environment around pigments that conventional wisdom expects very short lived quantum coherence, elucidating the mechanism and searching for an applicability of the coherence have become an interesting topic in both experiment and theory. We have previously studied the quantum coherence of a phycocyanin 645 complex in a marine algae harvesting light system, using Poisson mapping bracket equation (PBME). PBME is one of the applicable methods for solving quantum-classical Liouville equation, for following the dynamics of such pigment-protein complexes. However, it may suffer from many defects mostly from mapping quantum degrees of freedom into classical ones. To make improvements against such defects, benchmarking targets with more accurately described dynamics is highly needed. Here, we fall back to reduced hierarchical equation of motion (HEOM), for such a purpose. Even though HEOM is known to applicable only to simplified system that is coupled to a set of harmonic oscillators, it can provide ultimate accuracy within the regime of quantum-classical description, thus providing perfect benchmark targets for certain systems. We compare the evolution of the density matrix of pigment excited states by HEOM against the PBME results at physiological temperature, and observe more sophisticated changes of density matrix elements from HEOM. In PBME, the population of states with intermediate energies display only monotonically increasing behaviors. Most importantly, PBME suffers a serious issue of wrong population in the long time limit, likely generated by the zero-point energy leaking problem. Future prospects for developments are briefly discussed as a concluding remark.

Introduction

Photosynthetic process begins with the transformation of light energy to electronic excitation energy, which is transferred through the pigment-protein complex. This excitonic energy transfer (EET) is an important topic in understanding natural photosynthesis and in constructing artificial photosynthetic materials. Quantum yield in the EET can be very high, even close to unity in some instances, and many researches have made continuous progresses in finding the reason. Recently, long-lived quantum coherence in EET of photosynthetic system is observed in some photosynthetic bacteria at cryogenic and physiological temperatures.1-4 It is a significant discovery because the pigment-protein complex is so complicated that a common sense will expect that the energy would be transferred through incoherent decay of the initially excited state.5 Accordingly, many studies followed to elucidate the reasons behind such unexpected results.6-9

Some of the pigment-protein complexes that display quantum coherence are from cryptophtyes. Cryptophtyes has not only chlorophyll but also phycobiliproteins that can absorb almost all regions of visible light.10 For example, Phycocyanin PC645 pigment-protein complex has a primary role of the EET in photosynthesis of Chroomonas CCMP270. It consists of two distinct monomer which has one dihydrobiliverdin (DBV), two phycocyanobilins (PCB), and one mesobiliverdin (MBV) together with four protein subunits enfolding pigment groups.11 These bilins have maximum absorption wavelengths at 585, 645, and 622 nm. Because the phycobiliprotein can absorbs light with longer wavelength than chlorophyll, the algae can perform photosynthesis efficiently even in deep sea environment.10 Experimentally, it turned out that the complex shows quantum coherence in room temperature (294 K).3 In this experiment, the central DBV dimer was initially excited, and the excitonic energy passed through MBVs and finally excited PCBs.3

Theoretically, such a pigment-protein complex system is usually treated with a subsystem-bath model. The finite electronic states of pigments can be treated as a finite subsystem (often call just a “system”), and then the protein complex affects the subsystem as a bath. Thus, the model is composed of a subspace that follows quantum mechanics and the remaining subspace that is described by classical mechanics. Of course, the two subspaces “ cross-talk” to each other and the total system follows the quantum-classical mechanics. The quantum-classical Liouville equation (QCLE) is the most widely adopted starting point for treating such a system.12 Unfortunately, even the most simple bath model has almost an infinite number of degrees of freedom and the computational cost can become easily unaffordable. Therefore, many approximate methods have been proposed as alternatives such as iterative linearized density matrix (ILDM)13 approach and Poisson bracket mapping equation (PBME)14,15 formalism. Reduced hierarchical equation of motion (HEOM)16-18 is another propagation method for the subsystem-bath model. Unlike ILDM or PBME, reduced HEOM has many restrictions and cannot be applied to general subsystem-bath models. However, it is numerically exact and can be used to supply a reference to other approximation methods.

Previously, we have calculated the EET dynamics in PC645 by PBME.19 Here, we perform HEOM simulation to provide more accurate results of the same dynamics, with a purpose of providing benchmark references toward further developing PBME. We provide more detailed derivation of the equation, and compare the results of dynamics simulations of the same subsystem-bath complex with HEOM and PBME. Future prospects are also discussed in view of how to utilize the present results.

 

Mathematical Formulation

The formulation of reduced hierarchy equation of motion (HEOM) is already published.20,21 For completeness, the outline of the equation and an algorithm for solving this equation numerically will be explained here. With HEOM, the time-dependent solution is evaluated for the following Hamiltonian, which describes the excitation energy transfer:

The subsystem Hamiltonian, Hs represents the pigment part of photosynthetic pigment-protein complex, and the bath Hamiltonian, Hb describes the protein environment. Hs-b corresponds to the interaction of the subsystem and the bath. Hs is formulated by the Frenkel exciton model,

The basis of the subsystem is a set of single pigment excited states. Namely, | K〉 means a state with the excited k-th pigment with other pigments in their ground states. Ek is the excitation energy of the k-th pigment, and Jkl is excitonic coupling between k-th and l-th pigments. A generally accepted model for describing environment is a collection of phonons linearly coupled to excitation energies of the corresponding pigments:

Here, pξ and qξare dimensionless momentum and position of the ξ-th phonon mode. The interactions between the system and bath are described as

Vk is defined as Vk =|k〉〈k| , and uk = −∑ξcξkqξ, where cξk is the linear coupling coefficient between the subsystem and the bath. All the phonon modes of the bath are assumed to be independent from each other. The Liouvillian of the total system Ltotal, the subsystem Ls, the bath Lb, and the bathsubsystem coupling Ls-b are defined to match their corresponding Hamiltonian components. The time evolution of the total density matrix is described as

The formal solution of this equation can be written as the following:

By reducing the system into the subsystem, the evolution of the operators is simplified. For this, the bath part is ensembleaveraged, and only the subsystem part remains for detailed descriptions. The density matrix of the total system ρtotal(t) is reduced to ρ(t) by ρ(t) = Trb{ρtotal(t)}. For simplicity, we assume that the total density matrix at the initial state can be constructed as a tensor product of subsystem and bath components: ρtotal(0) = ρ(0)⊗exp(−βHb)/Z, where Z is a (canonical) partition function defined as Z = Trb exp(−βHb)). The interaction picture is considered here such that Hs and Hb evolve operators, and Hs-b propagates wavevectors. The tilde sign denotes this interaction picture nature.

With the assumption of the decomposition of initial total density matrix, the reduced composition of Liouvillian and total density matrix can be decomposed into the reduced time evolution operators and as

and

where the bath average of an operator is . When the exponential of the superoperator is explicitly expanded,

From , the n-th order term in the above equation is

when we adopt a vectorized time notation, ds = dsn···ds2ds1· There are a total of 2n individual terms in the commutator complex in the above equation. These can be formally enumerated as a sum as

where each αi can take either 0 or 1. Of course, the causality condition is taken care of by the time ordering operator T+ such that the set of is rearranged in time-increasing order. By adopting Eq. (4),

Now, it is easy to imagine that the last bath averaging term in the above equation will vanish when n is an odd integer due to the symmetry of the harmonic oscillator. In fact, from Wick's theorem, 22 it can be shown that this bath averaging term can be reduced into products of correlation terms from the Gaussian nature of the harmonic oscillators:

where the summation runs over all possible ways of forming pairs out from 2n operators for the bath averaging. In the classical case, this is equivalent to Isserlis' theorem23 of multivariate normal distributions, composed by the Gaussian harmonic oscillator position distributions. This fact also brings an important recursive characteristic between consecutive non-vanishing terms A2n and A2n+2:

This leads to the time evolution equation for the reduced density

where is involved in the evolution by bath coupled to j-th site,

The superoperator notations are defined as (commutator) and (anticommutator). The derivations of the recursion relationship and the formulation of can be found in the Appendix. The symmetrized correlation function of uk(t), namely , and the response function, describe the fluctuation of the excitation energy of the k-th site and the dissipation of the bath coupled to the k-th site. Then, the spectral density of the k-th site, Jk(ω), is defined by the imaginary part of the inverse Fourier transformation of the response function

where .Since the response function is odd by definition, the Ƒ−1 {χ(t)}(𝜔) is a pure imaginary function, and Eq. tells us that the spectral density is also an odd function. Therefore,

The quantum fluctuation-dissipation theorem24 constructs the relationship between the symmetrized correlation function and the response function χk(t) and the spectral density Jk(𝜔) as

The reorganization energy λk can be reproduced with spectral density by using a relationship . We construct the bath model with Drude-Lorentz spectral density, which approximately follows Markov process as . This spectral density describes overdamping Brownian motion and was widely used in experimental and theoretical studies.15,25-27 Then, the response function has the forms as in t > 0, and the correlation function is approximated provided that the temperate is high enough. Namely, as . With this condition, the time evolution of is given by

The operators in this equation are defined as

As () does not depend on s, Eq. (20) can be represented in a SchrÖdinger picture as

Here, σ is introduced as an auxiliary operator, transformed from the reduced density matrix by Θ, defined as

with vector index n. From its definition, σ(0,0,…,0)(t) is identical to ρ(t). The time evolution of the auxiliary operator is related to the other auxiliary operators with higher and lower index as

When the index (n1, n2,…, nN) of the auxiliary operator is 0, the third term in the above equation vanishes. Because the upper bound of the index is infinite in theory, the increasing index of the auxiliary operators must be truncated to calculate the time evolution equation in finite time. If the characteristic value of the subsystem part, 𝜔e, is smaller than (), the decay by the bath is faster than the evolution of the subsystem Hamiltonian. This approximation simplifies Eq. (25) as

if

The characteristic value can actually be estimated as the spectral range of the electronic Hamiltonian. Although σ is generated from ρ , it can be evaluated simultaneously with it using their time evolution equations. Therefore, by solving a set of linear differential equations with ρ and σ(1,0,…,0), σ(0, 1, …, 0), …, σnf , we can get the time evolution of the density matrix. The number of auxiliary operators needed to perform HEOM is

 

Dynamics of PC645

We assumed the identical bath to different sites, τ1(= γ1 −1) = τ2 = … = τN = τ and λ1 = λ2 = … = λN. Previous studies used the sum of two Drude-type spectral densities with different relaxation time.1319 However, for simplicity, we used only one density in this study, with shorter relaxation time τ = γ−1 = 50 fs together with the same reorganization energy as before as λ = 260 cm−1. We also adopted the excitation energy of each bilin and their excitonic coupling from the data in our previous study.19 For comparison, we calculated the dynamics by PBME with the same spectral density as in HEOM simulation. At least in the PBME case, we could observed that switching from the composite spectral density to a single Drude-type density has only a negligible effect on the dynamics.

Figure 1 shows the population dynamics of PC645 from HEOM simulations in a sub-picosecond scale. We can see a fast decay of the initially excited DBVc with strong coherence with DBVd. It is expected because the coupling between DBVs are the strongest in all couplings between bilins. The population decay of the dimer component is compensated by a nearly linear increase in the populations of other bilins. The coherent oscillation stays up to almost 300 fs. This dynamics resembles the results from previous PBME simulation at least in a qualitative sense.

Figure 1.The population dynamics of all bilins in PC645 at room temperature (300 K) by HEOM.

Figure 2 displays the comparison between the dynamics from HEOM and PBME. Comparing (a) and (b), it is apperant that the decay rate is faster in HEOM than that in PBME and the lifetime of the coherence is rather shorter. The most promenent difference appears in two MBVs of (c) and (d). In PBME, the populations of MBVs almost linearly increase as diplayed by PCBc158. However, in HEOM, the population of MBVs increases quite fast in a short time regime (~200 fs) and the population of PCBc158 starts to be larger than MBV populations after ~500 fs. Considering the Boltzmann factor, the populations of DBVs will be the least, followed by those of MBVs and PCBs. The energies of the same class of chromophores are of course similar. Namely, each group of DBVs, MBVs, and PCBs are close in energy. However, the energy gaps between different classes are much larger. For example, the DBV-MBV gap is 864 cm−1 and the MBV-PCB gap is 484 cm−1 respectively on average, with the MBV energy in the middle. Thus, we speculate that ~500 fs time scale is too short to induce a thermal equilibrium among these different classes of chromophores. However, this timescale is still long enough to display equibrium- like behavior inside the same class of chromophores (within DBVs, within MBVs, and within PCBs). The largest absolute coupling between any DBV and any MBV (43.9 cm−1) is somewhat larger than the maximum of absolute coupling between DBVs and PCBs (−46.8 cm−1). We may attribute the phenomenon that the populations of DBVs are transferred more to MBVs than to PCBs in the short time regime to this coupling strength difference. However, DBV and MBV energies are closer than in DBV/PCB pairs, and this closer match will likely induce better condition for more feasible energy transfer. The excitations of MBVs do not show oscillations and their transfer to PCBs is quite monotonic.

Figure 2.Comparison of population dynamics by HEOM (left) and PBME (right).

From an operational point of view, the computational cost of HEOM can be estimated from the ratio between the spectral range 𝜔e of Hamiltonian and γ. In PC645, 𝜔e = 1834 cm−1 is approximately 17 times of γ. If we consider ║n f║ ∞ ~1.5𝜔e/γ, the number of auxiliary operators needed to perform HEOM is almost 18 million. In the PC645 case, considering that the size of a single auxiliary operator is 8¥8, the size of the data needed at one time step of integration is almost 109 complex numbers. This translates to 16 Gigabytes of memory usage when the double precision representation is adopted for complex numbers. The computational time and the storage requirement increase rapidly with the number of elements in the auxiliary operators. Although utilizing parallel computation platforms with distributed memory will alleviate the computational cost dramatically as HEOM algorithm can be trivially parallelized, the high order scaling with respect to the size of the auxiliary operator (or the dimensionality of the subsystem) will be too severe for applying HEOM to any systems much larger than the presently studied PC645.

 

Conclusions

HEOM is known to give the accurate solution to harmonic oscillator bath models with linear coupling to subsystems. We have redisplayed the detailed working formula of HEOM together with a model Hamiltonian of PC645 that has been adjusted for application to HEOM. The excitonic state populations evolved qualitatively similar to the previously reported evolution patterns from PBME simulations. The coherent oscillation between populations on monomers of DBVs was apparent and their transfer to PCBs was confirmed. Distinct transfer characteristics in HEOM simulations compared to PBME results were that excitonic energy passed through MBVs before it was transferred to PCBs. The density elements on MBVs increased faster in short time regime, after which they decreased with increasing density elements on PCBs. This aspect was not captured in PBME simulations and populations on MBVs and PCBs increased together.

The PBME method is an approximation and possesses quantitative error. The error also is known to become severer as the time scale increases. Most importantly, the method cannot attain thermal equilibrium condition in the long time limit and the final populations do not follow Boltzmann distributions. The HEOM method does achieve the thermal equilibrium (confirmed also with the PC645 case) and shows quantitatively accurate results in both short-time and long-time regimes. The improvement of the original PBME method should reproduce these aspects. For example, the dynamics of PC645 should display an increase followed by a decrease in the populations of MBVs, instead of monotonic increases. The behaviors observed from HEOM simulations can definitely be used as references for improving PBME during future studies.