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.

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][3][4][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.
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][31][32][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][36][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

Results
Quantum Dynamics at Finite Temperature. Problems formulated in the quantum mechanics language require the calculation of the expectation value of some dynamical variable, where ρ(0) is the initial density matrix of the system, and A(t) = e iHt Ae −iHt 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 approximation Here 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 B T, where T is the temperature of system and k B the Boltzmann constant. Consequently, the trace operation involves only a summation over a thermal distribution of vibrational states, … n n where † a a , ( ) k k are the creation (destruction) operators of the k-th bosonic degree of freedom with frequency ω k , 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   † a a , k k and their corresponding occupation number states …   n n 1 2 , 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  n form a complete orthonormal set. We notice that in the above summation the numerical values of {n k } and  n { } k are identical. Defining the so called thermal vacuum state The above equation is an extension of the fundamental result of Thermo Field Dynamics 40, 54, 55 and the transformation e −iG is often referred to as Bogoliubov thermal transformation 56 where the wavefunction |ϕ(t)〉 satisfies the Schrödinger equation The modified Hamiltonian operator H is defined as where ∼ H vib 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 θ H 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  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 bath Here |n〉 describes an electronic state with the excitation localized on the n-th pigment, ε n is the electronic energy of state |n〉, J nm are electronic couplings, ω k are the frequencies of the bath of harmonic oscillators, and the parameters g nk 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 θ H which controls the finite temperature dynamics is readily obtained applying the Bogoliubov thermal transformation to the Hamiltonian operator 11 42 to exploit the invariance properties of the thermal Bogoliubov transformation 56 . At T → 0 the mixing parameters θ k become zero, sinh (θ k ) → 0, the coupling to the tilde space disappears, and the standard Schrödinger equation is recovered as expected. For high-frequency modes, θ  1 k , sinh (θ k ) ≈ 0 and cosh (θ k ) ≈ 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 ε n and electronic couplings J nm , have been retrieved from ref. 61. Vibronic coherences are essentially determined by the distribution of the bath vibrational frequencies and their coupling constants g nk . 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 nk . Since a single pigment is excited in each electronic state, only N components of g nk with = + − … k n N nN 1 ( 1) , , 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 nk 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 ω min = 2 cm −1 and the line spacing Δω = (ω max − ω min )/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 kn cosh (θ k ) and g kn sinh (θ k ) entering the thermal Hamiltonian θ H govern the coupling of the electronic subsystem with physical and tilde bosonic degrees of freedom. Hence, it is tempting to introduce the spectral densities which describe the electron-vibrational couplings in the physical (subscript p) and tilde (subscript t) subspace. As temperature goes to zero, J p (ω) → J(ω) and J t (ω) → 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 2 shows the total time-dependent populations p n (t) of seven (n = 1-7) BChl a 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 n (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 .
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.
To elucidate how the spectral densities of Figure 1 affect the fraction of BChl a s 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). 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 θ H 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 nm as well as the temperature-dependent electron-vibrational couplings g kn cosh (θ k ) and g kn sinh (θ k ). The parameters g kn sinh (θ k ) 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 kn sinh (θ k ) = 0 and the tilde variables are totally decoupled. At higher temperature both g kn cosh (θ k ) and g kn sinh (θ k ) 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 θ H . 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 (20) where the wavefunction |ϕ(t)〉 satisfies the equation We point out that in order to obtain a numerical solution of the Schrödinger equation 21 the Hamiltonian θ H must have an analytical representation or a form which is suitable for numerical treatment. This can be accomplished by expanding H in power series of creation-annihilation operators (or position and momentum operators) and using the fundamental relations 40 Quantum Dynamics with Tensor-Trains. Let us consider a generic expression of a state of a d dimensional quantum system in the form