Simulation of Quantum Dynamics of Excitonic Systems at Finite Temperature: an efficient method based on Thermo Field Dynamics.

Quantum electron-vibrational dynamics in molecular systems at finite temperature is described using an approach based on Thermo Field Dynamics theory. This formulation treats temperature effects in the Hilbert space without introducing the Liouville space. The solution of Thermo Field Dynamics equations with a novel technique for the propagation of Tensor Trains (Matrix Product States) is implemented and discussed. The methodology is applied to the study of the exciton dynamics in the Fenna-Mathews-Olsen complex using a realistic structured spectral density to model the electron-phonon interaction. The results of the simulations highlight the effect of specific vibrational modes on the exciton dynamics and energy transfer process, as well as call for careful modeling of electron-phonon couplings.

Introduction

Unraveling the role of quantum effects in the time evolution of various molecular systems and assemblies under realistic conditions at ambient temperatures is a key problem of modern physical and biological chemistry[1]. Long range energy and charge transfer in natural as well as in artificial systems are among the most important processes in which quantum coherent motion can be of relevance[2-5]. However, a proper understanding of these processes is often hampered by the impossibility to properly simulate the evolution of quantum systems with many degrees of freedom.

The quasi-adiabatic path integral (QUAPI)[6] and the hierarchical equations of motion (HEOM)[7-9] are among the most successful numerically exact methods for density matrix propagation. However, both methods become numerically demanding at low temperature[10] and when the Hilbert space of the system is very large[8, 10–15]. A number of approximate methods based on density matrix formalism is also available, but their range of validity can be very limited and system dependent[16-25]. Numerically accurate evolution of large systems has also been described by the density matrix renormalization group (DMRG) methodology, and the associated time-evolution algorithms[26, 27].

Wave function propagation methods employing a basis set representation, such as the multiconfiguration time-dependent Hartree (MCTDH) method and its multilayer extension (ML-MCTDH)[28, 29], Gaussian based MCTDH and other basis set methods[30-33], are powerful tools at very low temperature[34], but become unhandy in high temperature cases, as their application requires a statistical sampling of the initial conditions and faces both theoretical and computational difficulties[35-37]. On the other hand, basis set methods are very versatile, and capable of handling a large variety of Hamiltonian operators[38, 39].

The development of alternative approaches to the simulation of many body quantum dynamics at ambient conditions is therefore indispensable for better understanding and proper exploitation of quantum effects in nanosystems. In this work we discuss a novel theoretical methodology based on Thermo Field Dynamics theory[40, 41] that combines an accurate Hamiltonian description of quantum dynamics at finite temperature with the flexibility of a basis set representation[42]. We then apply this methodology to the study of exciton dynamics in the Fenna-Matthews-Olsen (FMO) complex, which has nowadays become a “guinea pig” of exciton dynamics theory and quantum biology[1, 2, 43, 44]. The exciton dynamics at FMO has been simulated by all major numerically exact quantum methods, e.g. QUAPI[45], HEOM[46, 47], ML-MCTDH[48]. Several explicit parameterizations of the bath spectral density accounting for the impact of intra- and inter-molecular modes on the FMO dynamics have been developed[49, 50]. In the present work, we simulate exciton dynamics in FMO at ambient temperature modeling the electron-phonon interaction with a realistic spectral density obtained from experimental data[51].

Results

Quantum Dynamics at Finite Temperature

Problems formulated in the quantum mechanics language require the calculation of the expectation value of some dynamical variable, A where ρ(0) is the initial density matrix of the system, and A(t) = e Ae − is the Heisenberg representation of the operator A, H being the system Hamiltonian ( = 1). The trace operation implies a weighted sum over all the thermally accessible states. In most molecular systems the energies of the electronic degrees of freedom are usually much higher than the vibrational energies. The effect of a finite temperature is then to create a thermal population of excited vibrational states, while only one electronic state of the entire system, |e〉, is tangibly populated. Within the validity of this condition we can safely employ the approximationHere Z is the proper partition function and ρ vib is the equilibrium Boltzmann distribution of the vibrational degrees of freedom, which, in the present work, is described using harmonic approximation, and β = 1/k T, where T is the temperature of system and k the Boltzmann constant. Consequently, the trace operation involves only a summation over a thermal distribution of vibrational states, where are the creation (destruction) operators of the k-th bosonic degree of freedom with frequency ω , and the cyclic invariance of the trace operation has been used for the symmetrization. Following the Thermo Field Dynamics approach[40] the above trace can be evaluated by introducing a set of auxiliary boson operators and their corresponding occupation number states , and rewriting the summation as[52] The dummy tilde variables do not affect the expectation value since A(t) is independent of them, and the states form a complete orthonormal set. We notice that in the above summation the numerical values of {n } and are identical. Defining the so called thermal vacuum statewhere represent the vacuum state of the ensemble of physical and tilde bosons, andthe expectation value 3 becomes[53] The above equation is an extension of the fundamental result of Thermo Field Dynamics[40, 54, 55] and the transformation e − is often referred to as Bogoliubov thermal transformation[56]. Equation 6 can be equivalently written as (see Methods section, and ref. 42)where the wavefunction |φ(t)〉 satisfies the Schrödinger equationwith the thermal operators

The modified Hamiltonian operator is defined aswhere is any operator acting in the vibrational tilde space. Equations 7, 8, 9 and 10 are the main theoretical result of the work. In the methodology described above the evaluation of the thermal average 〈A(t)〉 can be reduced to the solution of the thermal Schrödinger equation 8 with the Hamiltonian specified by Eq. 9, followed by the computation of the desired expectation value.

Since the coupling with the tilde space doubles the number of nuclear degrees of freedom, and since a thermal environment can be realistically mimicked only using hundreds or thousands degrees of freedom, the solution of the time-dependent Schrödinger equation 8 requires efficient numerical methods, suitable to treat a large number of dynamical variables[38]. Here we follow our recently proposed methodology[42] and represent the full vibronic wavefunction using the so-called Tensor-Train (TT) format (Matrix Product States, MPS, in the physics literature)[57-59]. Equation 8 is then solved using a methodology based on the time-dependent variational principle (TDVP) recently developed by Lubich, Oseledets and Vandereycken[59]. The reader is referred to the original papers[58, 60] for a detailed analysis of the TT decomposition, and to the Methods section (see also ref. 42).

Exciton Dynamics in the FMO complex

In order to apply the above methodology to the study of exciton dynamics at finite temperature in the FMO complex we consider a widely employed model Hamiltonian in which a set {|n〉}, of coupled electronic states interacts linearly with a phonon bathHere |n〉 describes an electronic state with the excitation localized on the n-th pigment, ε is the electronic energy of state |n〉, J are electronic couplings, ω are the frequencies of the bath of harmonic oscillators, and the parameters g determine the strength of the electron-phonon coupling. In the above notation the index k labels all the vibrations of the system.

The modified thermal Hamiltonian which controls the finite temperature dynamics is readily obtained applying the Bogoliubov thermal transformation to the Hamiltonian operator 11[42] Here we have used to exploit the invariance properties of the thermal Bogoliubov transformation[56]. At T → 0 the mixing parameters θ become zero, sinh (θ ) → 0, the coupling to the tilde space disappears, and the standard Schrödinger equation is recovered as expected. For high-frequency modes, , sinh (θ ) ≈ 0 and cosh (θ ) ≈ 1 even at room temperature. As a rule of thumb high-frequency modes need not be incorporated into the tilde Hamiltonian. This leads to additional reduction of the active space and computational savings.

The exciton part of the FMO Hamiltonian, the site energies ε and electronic couplings J , have been retrieved from ref. 61. Vibronic coherences are essentially determined by the distribution of the bath vibrational frequencies and their coupling constants g . Here we consider the general case in which the excited states of each BChl are independently coupled to a bath of N phonons (uncorrelated baths). Consequently, in the seven site FMO model, we have 7N vibrations, characterized by the coupling parameters g . Since a single pigment is excited in each electronic state, only N components of g with are nonzero for a given n. These parameters are assumed to be the same for all BChls and are conveniently specified by the so called bath spectral density[49, 50] We also point out that in our methodology the values of the parameters g can be determined using any method of choice, and no particular benefit is derived from the above specific assumption. Often the Drude-Lorentz spectral density is used, J(ω) = 2γλω/(γ 2 + ω 2), which however deviates substantially from the measured spectral density at 4 K[51]. Recent theoretical analyses suggest that the use of a structured spectral density can lead to quite relevant changes in the vibronic dynamics of the system[47, 48, 62, 63].

Following the very recent work by Schulze et al.[48], we model the electron-phonon interaction by discretizing the experimental spectral density of ref. 51 with N = 74 vibrations uniformly distributed in the range [2,300 cm−1]. This way the times corresponding to the frequency ω  = 2 cm−1 and the line spacing Δω = (ω  − ω )/N are safely beyond the observed time evolution of the system. Accordingly, our model consists of 74 vibrations per molecular site, thus 518 overall vibrational degrees of freedom which are doubled to 1036 due to TFD methodology. The numerical approach to the solution of this problem is described in the Methods section.

The parameters g  cosh (θ ) and g  sinh (θ ) entering the thermal Hamiltonian govern the coupling of the electronic subsystem with physical and tilde bosonic degrees of freedom. Hence, it is tempting to introduce the spectral densitieswhich describe the electron-vibrational couplings in the physical (subscript p) and tilde (subscript t) subspace. As temperature goes to zero, J (ω) → J(ω) and J (ω) → 0. The two spectral densities are reported in Figure 1 at 77 K and 300 K. The comparison of lower and upper panels in Figure 1 reveals how effective electron-vibrational coupling increases with temperature, notably for lower-frequency modes.

Figure 1

Effective site spectral densities J (ω) and J (ω) describing the coupling of the physical and tilde bosonic degrees of freedom with the electronic subsystem at different temperatures. (a,b) 77 K, (c,d) 300 K.

Figure 2 shows the total time-dependent populations p (t) of seven (n = 1–7) BChl molecules of the FMO complex (standard numbering of the FMO cofactors is used). The populations are evaluated by Eqs 7 and 8 for A = A  = |n〉 〈n| so that p (t) = 〈A(t)〉. The initial excitation is assumed to be initially localized on site 1. In all panels, p 1(t) and p 2(t) exhibit pronounced oscillations, as expected[45, 46].

Figure 2

The time evolution of the electronic populations p (t) of seven (n = 1–7) BChl molecules of the FMO complex at different temperatures indicated in the panels. The initial excitation is localized on site 1.

At T = 0 K (upper panel) the populations are in perfect agreement with the results obtained by Schulze and coworkers using ML-MCTDH[48].

At T = 77 K (middle panel) p 3(t) drops to about 0.6 at t = 1 ps. On the other hand, no pronounced difference in the behaviors of p 1(t) and p 2(t) at T = 0 and 77 K is observed. In the language of spectral densities defined in Eq. 14, it means that the contributions of the lower-frequencies vibrational modes (which are strongly temperature dependent) are quite significant in the dynamics of p 3(t) already at 77 K, while are less pronounced in the dynamics of p 1(t) and p 2(t).

If temperature increases up to 300 K (lower panel) p 3(t) further decreases to 0.3 at t = 1 ps, and the oscillatory components of p 1(t) and p 2(t) are significantly reduced but still visible. As can be seen from Figure 3, which shows an enlargement of the lower panel of Figure 2 at longer times, small amplitude beatings of p 1(t) and p 2(t) are clearly observable even after 700 fs. Such long-lived beatings at ambient temperature have not been reported in models employing an approximate spectral density in the Drude-Lorentz form[46], Ohmic form[45] or Adolphs-Regner (single peak) form[45]. The beatings revealed in the present work at T = 300 K are due to the strongly structured spectral density and frequency-dependent coupling between the electronic subsystem and the vibrational degrees of freedom.

Figure 3

Enlarged section of the lower panel of Fig. 2.

To elucidate how the spectral densities of Figure 1 affect the fraction of BChls that are significantly occupied during the time evolution of the system, we compute the inverse participation ratio Π(t), defined as[64, 65] It is easy to show that Π(t) = 1 for a completely localized exciton wavefunction, while Π(t) = N site (7 in the present case) for a perfectly uniform state. Therefore, Π(t) can be considered as an effective length, measuring the spatial extent of the exciton wave function over the aggregate. Figure 4 shows the computed Π(t) for the FMO complex at different temperatures. At T = 0 K Π(t) has a strong quantum behavior showing an oscillatory increase for the first 400 fs which is followed by an oscillatory decrease to a value of 2 at t = 1 ps. This is an indication of the exciton self-trapping. Therefore, a small number of sites are accessible to the systems during its evolution at T = 0 K, as is also evident from the population dynamics in Figure 2. For T = 77 K, the qualitative behavior of Π(t) remains the same but the number of accessible sites increases to 3 at t = 1 ps. At room temperature the number of accessible sites increases significantly and the effective length of the exciton is about 5.5 at t = 1 ps. The effect of a finite temperature is thus not only to provide a decoherence mechanism but also to increase the number of sites simultaneously accessible for the energy transfer process and to destroy the exciton self-trapping (cf. ref. 65).

Figure 4

Inverse participation ratio Π(t) as a function of time; (−) 300 K, (− −) 77 K, (− ·) 0 K.

Summarizing, we have developed a new theoretical and numerical approach for the determination of time dependent properties in large molecular aggregates. The methodology is based on Thermo Field Dynamics theory and requires doubling of all the thermalised degrees of freedom. The evaluation of an observable 〈A(t)〉 is reduced to the calculation of a simpler expectation value 〈φ(t)|A |φ(t)〉 where |φ(t)〉 satisfies the Schrödinger equation 8. The methodology has been implemented in the framework of the Tensor Train/Matrix Product State representation of the wave function, and using a novel technique for the numerical integration of Tensor Trains based on the time-dependent variational principle[59]. We have successfully applied this methodology to the study of quantum dynamics of energy transfer in the FMO complex. The present results draw attention to the importance of an accurate modeling of the bath spectral density. A highly structured spectral density is required to observe long-lasting oscillations in 〈A(t)〉 at room temperature which are absent when a simple Drude-Lorentz model[46] or Ohmic model[45] are employed.

The methodology developed in the present work offers new qualitative insights into the dynamics of excitonic systems at finite temperatures. The time evolution of these systems is governed by the thermal Hamiltonian of Equation 12, which looks like a standard excitonic Hamiltonian H of Equation 11, but contains twice as many vibrational degrees of freedom: the physical ones and the auxiliary (tilde) ones. Hence, all nontrivial dynamic effects are governed by three sets of the parameters: the electronic couplings J as well as the temperature-dependent electron-vibrational couplings g  cosh (θ ) and g  sinh (θ ). The parameters g  sinh (θ ) determine the coupling of the electronic degrees of freedom with the vibrational tilde degrees of freedom, and regulate the effective number of vibrational modes of the system. If T = 0 K, g  sinh (θ ) = 0 and the tilde variables are totally decoupled. At higher temperature both g  cosh (θ ) and g  sinh (θ ) tend to contribute on an equal footing. Since the system dynamics is purely Hamiltonian, the time evolution of any observable 〈A(t)〉 is the result of pure dephasing of the wave packet formed by a large number of vibronic eigenstates of the Hamiltonian . All oscillatory behaviors in 〈A(t)〉 are therefore vibronic by definition, since they are caused by the combined effect of electronic and temperature-dependent electron-vibrational couplings. Hence, a proper simulation of the time evolution of excitonic systems at finite temperature requires a careful modeling of electron-phonon interaction.

Methods

Derivation of TFD Schrödinger equation

Equation 6 can be transformed into a convenient Schrödinger representation by first rewriting it aswhere is any operator acting in the tilde vibrational space. The choice of the gauge is dictated exclusively by computational convenience and does not affect the expectation value 〈A(t)〉. Hencewhereand the modified Hamiltonian operator is defined as

Equation 17 is clearly equivalent towhere the wavefunction |φ(t)〉 satisfies the equationWe point out that in order to obtain a numerical solution of the Schrödinger equation 21 the Hamiltonian must have an analytical representation or a form which is suitable for numerical treatment. This can be accomplished by expanding in power series of creation-annihilation operators (or position and momentum operators) and using the fundamental relations[40] The transformed Hamiltonian depends on temperature through the parameters θ .

Quantum Dynamics with Tensor-Trains

Let us consider a generic expression of a state of a d dimensional quantum system in the formwhere |i 〉 labels the basis states of the k-th dynamical variable, and the elements are complex numbers labeled by d indices. If we truncate the summation of each index i the elements represent a tensor of rank d. The evaluation of the summation 24 requires the computation (and storage) of n terms, where n is the average size of the one-dimensional basis set, which becomes prohibitive for large d. Using the TT format, the tensor C is approximated aswhere G (i ) is a r  × r complex matrix. In the explicit index notationThe matrices G are three-dimensional arrays, called cores of the TT decomposition. The ranks r are called compression ranks. Using the TT decomposition 25 it is possible, at least in principle, to overcome most of the difficulties caused by the dimensions of the problem. Indeed, the wave function is entirely defined by d arrays of dimensions r  × n  × r thus the required storage dimension is of the order dnr 2.

In a time-dependent theory the cores G (i ) are time dependent complex matrices whose equations of motion can be found by applying the time-dependent variational principle (TDVP) to the parametrized form of the wave functionThe resulting equations of motion can be written in the formand provide an approximate solution of the original equation on the manifold of TT tensors of fixed rank, . In equation 28, is the orthogonal projection into the tangent space of at |Ψ(G(t))〉. We refer the reader to refs 59 and 66, where the explicit differential equations are derived and their approximation properties are analyzed, and to ref. 67 for a discussion of time-dependent TT/MPS in the theoretical physics literature.

Several techiques exist to compute the time evolution of TT/MPS[59, 68–70]. Here we adopt a methodology recently developed by Lubich, Oseledets and Vandereycken, which combines an explicit expression for the projector and an extremely efficient second order split projector integrator specifically tailored to the TT format[59]. The computations presented in this paper have been performed using a code based on the software library developed by Oseledets and coworkers.

All results presented in this paper are numerically converged in the sense that the ranks of the TT cores are increased until no significant variations are observed in the solution. The dynamics at 300 K required an average value r  = 40 to obtain fully converged results.