Finite temperature quantum dynamics of complex systems: Integrating thermo‐field theories and tensor‐train methods

This review provides the fundamental theoretical tools for the development of a complete wave‐function formalism for the study of time‐evolution of chemico‐physical systems at finite temperature. The methodology is based on the non‐equilibrium thermo‐field dynamics (NE‐TFD) representation of quantum mechanics, which is alternative to the commonly used density matrix representation. TFD concepts are extended and integrated with the tensor‐train (TT) numerical tools leading to a novel and powerful theoretical and computational framework for the study of complex quantum dynamical problems. In addition, NE‐TFD techniques are extended to enable the study of dissipative open systems via a new formulation of the hierarchical equations of motion (HEOM) fully integrated with TT methodologies. We demonstrate that the combination of the TFD machinery with computational advantages of TTs results in a powerful theoretical and computational framework for scrutinizing dynamics of complex multidimensional electron‐vibrational systems. We illustrate the validity and the computational advantages of the developed methodologies by applying them to the study of quantum coherence effects in the energy‐transfer processes in antenna systems, to the analysis of fingerprints of vibrational modes in electron‐transfer and charge‐transfer processes in various model and realistic multidimensional molecular systems, as well as to simulation of other fundamental models of physical chemistry.

If mñi f g is an orthonormal basis of ℒ then mñjm 0ñ0 h i¼ δ mm 0 δññ0 X mn j mñihmñ j¼ 1: Further, the identity vector jIi is defined as This special vector allows to define a mapping between the dual space of H (i.e., the bra space) and the tilde space, indeed we have Using these relations it is possible to associate a vector of the ℒ space to each operator A acting in the space H Similarly, we can define a state vector jρ(t)i = ρ(t) j Ii, where ρ(t) is the density matrix of the system. Accordingly, the expectation value of A is defined as the scalar product The meaning of the above notation can be easily understood using the closure relation whence it is clear that the vector jAi is a linear combination of a basis of ℒ with coefficients given by the matrix elements A mn . Together with operators acting in the space H it is possible to define a set of operators acting in the space H . In particular, following Suzuki, 39 two operators A and B are weakly equivalent if and we write For each Hermitian operator A acting in the space H it is possible to define a tilde operatorÃ that is weakly equivalent to A as where the dag operator implies the Hermitian conjugation. Consequently, for Hermitian operators The tilde operator can be obtained from the original operators by the so-called tilde conjugation rules 40 and If A, B are two operators of the space H andÂ ¼ A ÀÃ † then proving the fundamental property of the twin-space formalism Finally, we notice that the identity vector jIi is invariant under any unitary transformation of the original Hilbert space basis set jmi. 40 Equation (12) allows to rewrite quantum statistical mechanics in a commutator-free way. Indeed, it is easy to demonstrate that the evolution of the vector jρ(t)i is given by the equation (ℏ = 1) whereH is a Hamiltonian operator identical to the physical Hamiltonian H but acting in the tilde spaceH. The super-operatorĤ ¼ H ÀH acts in the entire ℒ space and is the Liouville operator rewritten as a superoperator in the ℒ space. We refer the reader to the demonstration reported in the original articles by Schmutz 41 and by Suzuki. 39 When the initial condition of the system can be described by a Boltzmann distribution with a certain zero-order Hamiltonian operator H 0 we have where Z 0 is the partition function and β = 1/k B T, k B being Boltzmann's constant and T is the temperature. It is now worth comparing the above formulation with the more common double-space theory 42 for clarifying similarities and differences of the two approaches. In double-space theory the notation jmnii = jmihnj is often used to identify a state of a space that is the direct product of H and its dual. 42,43 Once we let the tilde space and the dual space coincide this approach and the twin formalism seem identical. 44 However, in standard double-space formulation what is often referred to as Liouville superoperator is a mere symbol for a commutator, L = [H, Á], which induces a Lie algebra in the Liouville space. In the twin-formulation outlined above the introduction of the tilde space and of its operators allows to eliminate all the commutators, overcoming the difficulties inherent in their evaluation, and replace them with the action of an operator in a Hilbert space. From a mathematical point of view the double vector jIi allows to map the dual of H into the tilde spaceH. Furthermore, as we shall see in the next section, the twin formulation permits a certain freedom (up to the gauge transformation) of the definition of the twin-space counterpart of the Liouville operator, which can be used to make mathematical formulations more compact and numerical simulations more efficient.

| Finite-temperature Schrödinger equation
TFD stems from the simple observation that an expectation value in the twin-space formalism can be cast into a wavefunction-like form using the identity where we have used the property of Equation (7), and the state vector with the initial condition Therefore the use of the twin space allows to recast all quantum statistical mechanics into a wave function theory in a double Hilbert space ℒ ¼ H H À Á . We notice that the density vector jρ(t)i and the new wave function jψ(t)i obey the same evolution equation but with different initial conditions.
Let us now consider a system that is prepared at thermal equilibrium at t = 0, that is, Z being the partition function. We further observe that in most molecular systems the characteristic energy of electronic states are much higher than those of nuclear (vibrational) states, therefore we can assume that composite electronvibrational system at ambient temperature resides in the ground state of the electronic subsystem (A), j0 A i, and possesses thermal distribution over the DoFs of the vibrational subsystem (B), that is H B being the vibrational Hamiltonian, in the ground electronic state, and Z B the corresponding partition function. The assumption (20) is rooted into Born-Oppenheimer approximation, it has a huge number of applications in molecular physics and spectroscopy, 1 and is the key to the description of molecular aggregates, 46-48 electron-vibrational 49 , and vibration-rotational 50 systems. If the above assumption holds true, it is possible to show that Equation (17) can be rewritten in the equivalent form 51 in which the new Hamiltonian operator, H, is defined as whereh B is any operator acting in the tilde space of the vibrational subsystem, which can be chosen as an arbitrary gauge to help solve the mathematical problem. Equation (21) should be solved with the initial condition where, according to Equation (18) in which jki labels the vibrational states of the original Hilbert space and jki the states of the corresponding fictitious tilde space. Therefore, the solution of the original Liouville -von Neumann equation for ρ(t) with the initial condition (20) is reduced to the solution of the Schrödinger equation (21) in the extended vector space spanned by the basis ni ki ki É È , where jni labels the electronic basis set of system A, and { j ki, jkig the physical and tilde vibrational basis sets, respectively. In this special formulation of TFD theory we have explicitly taken advantage of the separation of the energy ranges of the electronic and vibrational subsystems (Equation (20)), avoiding the use of tilde variables for the entire system and employing them only for the thermalized DoFs of subsystem B (cf. Refs. 52 and 53).
The thermal vacuum vector j0 B (β)i possesses multiple nonzero components, and its use in practical simulations can be rather difficult, requiring imaginary time propagation. [54][55][56] However, the key advantage of the present method lies in the possibility to have a compact analytical representation of the thermal vacuum j0 B (β)i. 57 Indeed, instead of the solution of the Schrödinger Equation (21) with the initial condition (23) it is much more practical to introduce a unitary transformation e ÀiG in the ki ki É È subspace obeying the identity where j 0 B0B i is the ground state in the ki ki É È subspace. Equation (25) defines what can be called generalized thermal Bogoliubov transformation. 58 Its application to Equation (21) yields the transformed TFD Schrödinger equation where and the initial condition corresponds to the global vacuum state. The expectation value of any operator A (acting in the physical, {j nij ki}, vector space) can now be rewritten as We underline that the expectation values evaluated via Equations (15) and (30) are identical. The TFD machinery with thermal Bogoliubov transformation gives an alternative representation of quantum mechanics which is fully equivalent to the traditional density-matrix representation. The last step required to complete our theory is to provide an explicit analytical form of the thermal Bogoliubov transformation and a set of rules to obtain the transformed Hamiltonian operator H θ . The explicit form of the transformation 25 is well known for two cases of notable interest, that is, for ensembles of free bosons and free fermions. [21][22][23]39,59 Extensions to more general systems can also be obtained although their practical use can be rather cumbersome. 51 The problems we wish to tackle in this work require, in most cases, that the physical observable of interest be averaged over a thermalized ensemble of nuclear vibrations which are modeled as harmonic oscillator states, that is H B of Equations (20) and (24) is where a † j (a j ) are the creation (annihilation) Bose operators ( a j , a † j 0 h i ¼ δ jj 0 ) and ω j are vibrational frequencies. The collection of vibrational DoFs that defines the above Hamiltonian operator includes both high frequency modes typical of molecular systems, and low frequency "environment" modes. After the seminal work of Caldeira and Leggett on quantum dissipative systems, 60 the importance of this model of the environment in various applications is hardly overestimated.
The operator of the thermal Bogoliubov transformation corresponding to the Hamiltonian of Equation (31) reads 22 where Hence, thermal Bogoliubov transformation introduces thermal noise into the physical system by coupling it to the fictitious tilde system through the temperature-dependent mixing parameters θ j . If the Hamiltonian H and the operatorh B are polynomials in creation-annihilation operators (which is the case in many real world applications) then the explicit form of H θ can be constructed by using the fundamental relations 22,53,61 e iG a j e ÀiG ¼ a j cosh θ j À Á þã † j sinh θ j À Á ð34Þ Therefore, the transformed Hamiltonian H θ depends on temperature through the parameters θ j . Equations (26)-(35) are the working tools of our new approach to finite temperature quantum dynamics in molecular systems.

| Finite-temperature electron-vibrational Hamiltonian
Let us now focus on a special type of the Hamiltonian operator describing a set of coupled electronic states interacting with a phonon bath Here jni labels an electronic state of the system with the energy ε n , V nm are electronic couplings, ω k are the frequencies of the bath harmonic oscillators, and the parameters g nk determine the strength of the electron-phonon coupling. The Hamiltonian (36) has a large number of applications, such as the fundamental description of molecular aggregates properties, 46-48 the analysis of molecular processes in the linear vibronic-coupling theory, 49,62 electron-transfer and charge-transfer properties in molecular systems. 11,63,64 In what follows we choose the gauge of Equation (22) as Applying the Bogoliubov transformation to the Hamiltonian operator (36) we obtain the temperature-dependent operator H θ ¼ e iG He ÀiG ¼ In deriving the above expression we used the invariance property 22 e iG a † n a n Àã † nã n À Á e ÀiG ¼ a † n a n Àã † nã n : The operator H θ of Equation (38) consists of two parts: a modified physical Hamiltonian in which the linear coupling terms are multiplied by cosh(θ k ) factors, and the vibrational tilde Hamiltonian. The excitation of the tilde vibrations is caused by the linear coupling terms g kn sinh(θ k ). Sinceh B enters Equation (38) with a negative sign, vibrational excitations in the tilde space correspond to a flow of energy from the physical system to the fictitious tilde system. It is this type of coupling that accounts for thermal noise. Once the explicit structure of H θ has been determined the evaluation of the expectation value of any observable hA(t)i can be obtained from Equation (30) after the solution of the TFD Schrödinger Equation (26).
The explicit dependence of the Hamiltonian operator on the temperature through the parameters θ k enables a smooth transition from the low temperature to the high temperature case. 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, θ k ( 1, 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 (see Section 4.3). This a priori selection of the DoFs which need doubling cannot be used within the standard formulation of the Liouville-von Neumann equation for the density matrix: in that case the bra and ket spaces have by construction the same dimensionality.

| Basic theory
The solution of the TFD Schrödinger Equation (26) with the Hamiltonian (38) requires efficient numerical methods, suitable to accurately treat a large number of dynamical variables. Several techniques have been developed which can, at least in principle, overcome what has been termed the curse of dimensionality. 3,15 In our approach the TT decomposition, the simplest form of Tensor Network, has been adopted. Below we sketch the basic principles of the TT decomposition, and show how it can be applied to solve the thermal Schrödinger equation in twin-formulation. The reader is referred to the original articles 3,35,65 for a detailed analysis of the TT theory.
Let us consider a generic state jψi of a N dimensional quantum system having the form where ji k i labels the basis states of the k-th dynamical variable, and the elements C(i 1 , …, i N ) are complex numbers labeled by N indices. If we truncate the summation of each index i k the elements C(i 1 , …, i N ) represent a tensor of order N, where the word "order" means nothing but the number of DoFs. The evaluation of Equation (40) requires the computation (and storage) of p N terms, where p is the average size of the one-dimensional basis set that is usually much larger than 2, becoming therefore prohibitive for large N. Using the TT format, each element C(i 1 , …, i N ) of the tensor C is approximated as where each C k (i k ) is a r k À 1 Â r k complex matrix. In the explicit index notation The matrices C k are three dimensional arrays, called cores of the TT decomposition. The dimensions r k are called compression ranks. Since the product of Equation (41) must be a scalar, the constraint α 0 = α N = 1 must be imposed. In the matrix-product state (MPS) language the sizes of the ranks are referred to as bond dimensions. Using the TT decomposition 41 it is possible, at least in principle, to overcome most of the difficulties caused by the dimensionality of the problem. Indeed, the wave function is entirely defined by N arrays of dimensions r k À 1 Â n k Â r k thus requiring a storage dimension of the order Npr 2 . The meaning of the term Tensor Train can be easily understood by looking at the graph of Figure 1.
We notice that if all the ranks of the cores are 1, the TT format is equivalent to a Kronecker product. From this point of view it constitutes the simplest extension of a direct product to entangled systems. The higher the entanglement, that is the correlation between two or more DoFs, the larger the ranks of the TT cores connecting their respective indices. Indeed, a fundamental problem of the TT representation is to define the sequence of indices i k , k = 1, …, N with the minimum rank decomposition of the state vector, and of the operators. We do not aim to tackle this complex mathematical problem, however we can exploit the twin-space formalism to optimize the TT representation of the vector jψ(t)i.
The algebra of tensor trains has been thoroughly discussed elsewhere and efficient algorithm exists to perform basic operations such as scalar and matrix-vector products. For completeness here we briefly describe some of the basic operations, and refer the reader to more technical articles for details. 35 The addition of two vectors A and B in TT format can be easily implemented by merging the cores as in a direct sum. The cores of the new tensor C = A + B will be It follows that if the matrices A k (i k ) and B k (i k ) have ranks r k 0 and r k 0 0 , respectively, the rank of the core C k (i k ) will be r k ¼ r k 0 þ r k 0 0 . This increase in the ranks upon addition is a fundamental problem which is usually addressed by rounding the new TT with a prescribed accuracy (see Section 3.2). Matrix-by-vector multiplication is certainly one of the most important operations in linear algebra. In the manydimensional case this operation requires the contraction over a set of indices in the form The tensor X is said to be in TT-format if it is represented as Graphical representation of a tensor train. Each square node represents a core of the TT, and each vertical line represents an index i k of the tensor. Connecting lines represent the contractions over the indices α k X i 1 , ::, where X k (i k , j k ) are r k À 1 Â r k matrices, that is, in the explicit index notation If the vector B is in TT format then we have In this operation each core X k acts separately on the cores B k of the vector B. We notice that the TT-ranks of the resulting vector A are the product of ranks of the matrix and of the vector. This latter result severely limits the numbers of consecutive matrix-by-vector multiplications that can be performed and implies that after a multiplication a TTrounding must be performed to avoid rank growth.

| Tensor-train rounding
Suppose that the tensor A is in the TT-format with core ranks r k . This tensor can be for example the result of the addition of several tensors and therefore the ranks r k can be particularly large. To avoid waste of computer memory and time it is desirable to approximate A with a new tensor B in TT format such that Such a procedure is called rounding. The new tensor B will have ranks r k 0 ≤ r k while maintaining the prescribed accuracy ϵ. This problem is well known in the theory of TT decomposition and can be elegantly solved performing a sequence of singular value decompositions (SVD) of each core of the original tensor A retaining only singular values above a predefined threshold. The reader is referred to the original article for details about the implementation. 35

| Vibronic Hamiltonian operator in tensor-train format
In order to apply TT representation to solve the time-dependent Schrödinger Equation (21) both the state vector and the Hamiltonian operator have to be represented in TT format. The construction and storage of the Hamiltonian operator of Equation (36) in TT format can be performed by a sequence of additions of TT matrices followed by a rounding to keep the ranks to a reasonable value. If M is the number of electronic states and p k is the size of the truncated basis set for the k-th vibrational DoF the Hamiltonian operator (38) has the matrix form where H e is the electronic Hamiltonian matrix Q k and W k are the coordinate and energy matrices of mode q k in its energy eigenbasis representation I e is a unit matrix of size M Â M, S en is the matrix representation of the jnihnj electronic operator, and I k are identity matrices of size (p k Â p k ), k = 1, …, d.We now wish to provide a TT approximation, H ε of the Hamiltonian matrix (52) of the type such that To this end we observe that each direct product in Equation 52 is a tensor train with all ranks r k = 1, therefore their sum can be easily obtained using the rules of Equation 44. Since after addition the ranks of TT cores increase it is necessary to perform a TT rounding after each operation. Algorithm 1 describes the pseudocode for obtaining the TT approximation of the vibronic Hamiltonian matrix (52) with a prescribed accuracy ϵ.

| Application of tensor trains to time-dependent problems
In a time-dependent theory the cores are time-dependent complex matrices which must be determined from the numerical solution of the equations of motion of the system. 36,66,67 The simplest possible approach to tackle this problem is to use standard integrators of systems of ordinary differential equations (ODEs) combined with the algebra of tensor trains. Let us assume that the wave function can be represented as and that H is the Hamiltonian operator in TT format that controls its evolution. We could use a forward Euler scheme and obtain the solution of Equation (21) in the form As described in the preceding section the ranks of the cores of the vector y = HC are obtained by multiplying the ranks of H and those of C, and the sum C À iτHC results in a further increase of the ranks of the final vector. This scheme therefore requires an efficient rounding procedure after each time step. Almost always the truncation scheme allows for an a priori upper bound of the ranks r k . This truncation, however, results in artificial changes of the total energy and of the norm of wave function under unitary evolution. The so-called Time-Evolving Block Decimation (TEBD) technique, the Krylov subspace methods [68][69][70] belong to this class of methods. A recent implementation of the second-order Cranck-Nicholson integrator combined with the Alternating-Minimal-energy solver of linear systems in TT format has been proposed to control the increase in the ranks of the cores. 67 An alternative approach to tackle this problem is to apply the time-dependent variational principle (TDVP) to the parameterized form of the state given by Equation (57). 71 Since TTs of fixed rank form a closed manifold ℳ TT the resulting equations of motion can be written in the form where C labels all the cores of the TT representation (57), andP T ψ C t ð Þ ð Þ ð Þ is the orthogonal projection into the tangent space of ℳ TT at jψ(C(t))i. Equation (59) provides an approximate solution of the original equation on the manifold of TT tensors of fixed rank, ℳ TT . 65 This means that the projected evolution gives the best solution with a prescribed upper bound for the ranks of cores. This strategy is extremely appealing because it avoids, by construction, the growth of the ranks of the TT solution. The drawback of the methodology is that the accuracy of the solution must be verified a posteriori. This requires several calculations, with increasing the values of the ranks r k , to be performed, until a global convergence on a desired observable is reached.
We refer the reader to references, 36,72 where the explicit differential equations are derived and their approximation properties are analyzed, and to Ref. 73 for a discussion of time-dependent TT/MPS approximations in the theoretical physics literature.
In Sections 4 and 5.3 we employ the TDVP algorithm, and show how to control the behavior of the solution as a function of the average value of the ranks.

| Entanglement growth and structure of the tensor train
The TT approximation is effective only if the ranks r k of the cores are small. The structure of the train, that is the order of the indices of the tensor, can have a deep impact on the growth of the entanglement. A key aspect is that if two DoFs are highly entangled their indices should be as close as possible in the sequence that defines the TT representation of the wave function. The most desirable situation would correspond to a systems showing only nearest neighbor interactions and a structure of the tensor train that reflects the natural form of the Hamiltonian operator. This latter is dictated by the structure of the problem under examination, but several techniques exist to map the Hamiltonian of Equations (36) and (38) to chain-like models, thus easing the application of the TT format. [74][75][76] Yet, it has been recently demonstrated that this mapping might not be necessary at all, and sometime can even be counterproductive. 77 In all the applications presented in Section 4 we have found that ordering the electronic and vibrational DoFs of the Hamiltonian in Equation (27) with increasing values of their frequency ω i provides very good convergence properties of the TT representation, and the mapping to a chain form is not required.

| APPLICATIONS
The TFD-TT approach can be applied to study the time evolution of quantum observables of a wide range of chemicophysical systems. Below we describe some applications both to the study of model systems as well as to the characterization of energy and electron transfer problems in realistic many-dimensional molecular systems.
The Hamiltonian operators of all systems considered in this section are particular cases of the Hamiltonian of Equation (36) and of the thermal Hamiltonian of Equation (38).

| Spin-boson model
Starting from the general Hamiltonian of Equation 36 we have studied a prototypical spin-boson system 9 in which two degenerate electronic states, j1i and j2i, are coupled to low frequency vibrational modes where, as usual, σ z = ( j 1ih1j À j 2ih2 j ) and σ x = ( j 1ih2j + j 2ih1 j ), 2ϵ is the electronic energy difference between the two states, V is the electronic coupling, and the constants g k are responsible for the linear coupling between the ensemble of bosons of frequencies ω k and the spin. These parameters can be modeled starting from the so-called spectral density function, defined as In our model the Ohmic spectral density is employed with a cut-off frequency ω c = 53 cm À1 . The solution of Equation (26) using a basis set representation requires the discretization of the spectral density over a finite set of frequencies. We have adopted a non-uniform discretization procedure that ensures fast convergence with respect to the number of sampling points in which the frequencies and the coupling parameters are given by: 78 In our calculations the spectral density is discretized with d = 200 DoFs in the range (0,5ω c ]. Figure 2 shows the population P(t) of the initial upper electronic state at 300 K for the set of parameters reported in the caption. The TFD-TT results correspond to the full lines, while the blue dots correspond to the numerically exact populations computed in Ref. 79 via the HEOM methodologies (HEOM is a standard reference for benchmarking hightemperature simulations). Clearly, the TFD-TT and HEOM populations are in excellent agreement and virtually indistinguishable. The comparison unequivocally demonstrates the validity of the TFD-TT approach for this type of quantum dynamical problems.
The convergence properties of the numerical methodology are illustrated by Figure 3 which shows the population as a function of time for different values of the TT compression ranks. At a low temperature, T = 30 K, the convergence is achieved with very low compression ranks, while at higher temperature higher ranks are required. Since the required TT storage scales quadratically with the average TT rank, this amounts to an increased computational cost of the calculation.

| Reaction mode spin-boson
One of the key advantages of using a basis set approach is the possibility to describe a large variety of potential energy surfaces beyond linear electron-phonon couplings. This is important for the description of chemical reactions in which the surfaces can be highly anharmonic. In the following example the TFD-TT technique is used to study a system with bilinear couplings between electron and phonons. 9 We consider a system of two electronic states coupled to a single harmonic mode with the frequency Ω which is in turn coupled to a set of harmonic oscillators having Ohmic spectral density. The corresponding Hamiltonian, often referred to as reaction mode spin-boson, 80 can be written as where A, a k (A † , a † k ) denote the annihilation (creation) operators for the reaction and bath oscillators respectively, and the couplings coefficients λ k satisfy the relation The value of the parameter g determines the strength of the coupling between the high frequency mode and the electronic subsystem. In this model the phonons "drain" energy from the reaction coordinate and not directly from the spin system. As in the previous case the continuous spectral density has been discretized using 200 vibrations; and the frequencies ω k and couplings λ k are obtained from Equation (63).
In Figure 4 the electronic population P(t) at 30 K and 300 K is shown for two different values of the Kondo parameter α of Equation (65) (see figure caption for the Hamiltonian parameters).
At 30 K and small α the boson bath is not very effective in dissipating energy from the reaction coordinate. In this regime coherent oscillations of the electronic population persist at long times. When temperature increases to 300 K the damping is more evident but the oscillations remain underdamped. For larger α, the population at 30 K exhibits underdamped oscillations, while for 300 K the bath quenches the beatings at around 600 fs. The trend is natural, since the Kondo parameter controls the coupling of the reaction mode to the harmonic bath. The high-frequency modulation of the population dynamics is due to coherent vibrations of the reaction mode with a period 2π/Ω ≈ 22 fs.

| Energy transfer in supra-molecular assemblies
In the two preceding examples the vibronic couplings are derived from analytical functions which were devised to mimic mostly the behavior of low frequency vibrations of condensed phase systems. However, in most molecular systems this description does not hold, indeed, quite often the spectral densities are extremely structured, that is, absolutely not uniform. Providing an accurate description of finite temperature quantum dynamics in these systems can be a very challenging problem for most methodologies. In recent works 10,81 we have shown how the TFD-TT methodology can be applied to the study exciton dynamics in the Fenna-Matthews-Olson (FMO) complex, consisting of seven bacteriochlorophylls (BChl) embedded into a protein matrix. This type of problem can indeed be described by using the general model Hamiltonian of Equation (36).
The exciton part of the FMO Hamiltonian, the site energies ε n and electronic couplings V nm , have been retrieved from Ref. 82. Vibronic coherences are essentially determined by the distribution of the bath vibrational frequencies and their coupling constants g kn . For numerical convenience these parameters are assumed to be the same for all BChls (g kn = g k ). Recent theoretical analysis suggest that the use of a "structured" spectral density in the energy transfer process can lead to quite pronounced vibronic effects. 5,83-85 Therefore, following Schulze et al., 84 we have modeled the electron-phonon interaction by discretizing the experimental spectral density of Ref. 86 with d = 74 vibrations uniformly distributed in the range (2300 cm À1 ). This way the times corresponding to the frequency ω min = 2 cm À1 and the line spacing Δω = (ω max À ω min )/d are safely beyond the observed time evolution of the system. We notice that using this spectral density the overall number of physical and tilde vibrational DoFs is 1038 (14d). We are not aware of any simulation at finite temperature with such a large number of DoFs.
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 DoFs. 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 5 at 77 K and 300 K. The comparison of lower and upper panels of the figure reveals how effective electron-vibrational coupling increases with temperature, notably for lower-frequency modes. Figure 6(a) shows the total time-dependent populations p n (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 Equations (30) and (26) for A = A θ = jnihnj so that p n (t) = hA(t)i. The initial excitation is assumed to be localized on site 1. In all panels, p 1 (t) and p 2 (t) exhibit pronounced oscillations, as expected. 87,88 At T = 0 K (upper panel) the populations are in perfect agreement with the results obtained by Schulze et al. using ML-MCTDH. 84 At T = 77K (middle panel) p 3 (t) drops to about 0.6 at t = 1 ps. On the other hand, there is no pronounced difference in the behaviors of p 1 (t) and p 2 (t) at T = 0 and 77 K. In the language of spectral densities of Equation 66, it means that the contributions of the lower-frequencies vibrational modes (which are strongly 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. Such long-lived beatings at ambient temperature have not been reported in models employing an approximate spectral density in the Drude-Lorentz form, 87 Ohmic form or Adolphs-Regner (single peak) form. 88 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 DoFs.
To elucidate how the spectral densities of Figure 5 affects the fraction of BChls that are significantly occupied during the time evolution of the system, we have computed the inverse participation ratio Π(t), defined as 46,89 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 6(b) 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 6(a). 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 ambient 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. 89).

| Electron-transfer in photosynthetic reaction centers
One of the key advantages of the TFD-TT methodology developed in Section 2 is that it can combine the best of wavefunction and density-matrix approaches into a single theoretical framework. Indeed, in our formulation the tilde space comprises only thermalized low frequency DoFs, therefore whenever the thermal mixing parameter θ of Equation (33) falls below a certain threshold the corresponding tilde variable is totally disentangled from the evolution of the physical system and can be omitted from the dynamical problem (see also Equation (29)). The application of TFD-TT to the study of the dynamics of the electron-transfer (ET) between the accessory bacteriochlorophyll (B A ) and the bacteriopheophytin (H A ) in bacterial reactions centers, 11 clearly shows this special feature of the methodology.
Previous numerical studies have shown that this process can be modeled as a radiationless transition which involves mainly intramolecular vibrations which carry most of the reorganization energy. 19,64 However, quantum dynamical results were obtained only for models including a reduced space of nuclear vibrational coordinates, and the effect of finite temperature was not taken into account. 19 Following a common approach, 90,91 the ET process is described by the Hamiltonian of Equation (36) with only two states, jAi, j Fi to describe the electronic subsystem. These latter are represented by the direct product of the neutral and anionic forms of the two isolated molecules, that is, j Ai ¼j B À A i j H A i and j Fi ¼j B A i j H À A i. The Hamiltonian (36) can then be rewritten in the form where ϵ A , ϵ F are the electronic energies of the states jAi, jFi respectively, V is their electronic coupling. The linear vibronic couplings g k can be obtained from the computation of the equilibrium geometry, and of the matrices of normal modes of vibrations of the two electronic states via well-known relations. [92][93][94][95] Our model comprises 267 vibrational modes, and the parameters adopted in here are taken from our previous work. 96 Figure 7 shows the temperature-dependent spectral densities defined in Equation (66). As already said, the coupling with the tilde space is negligible for high frequency modes (ω > 2000 cm À1 ) and not reported. It is immediate to see that as temperature increases the coupling with the low frequency part of the spectrum increases, while the high frequency region is left almost unaffected. This observation enables to analyze the relevance of temperature effects for each single degree of freedom and to reduce the computational costs by a priori removing some of the tilde DoFs from the Hamiltonian.
In the present case the physical number of DoFs is 267, which should be doubled to 534 upon inclusion of the tilde space. However, a large fraction of high frequency modes has a negligible vibronic coupling, g k sinh(θ k ), thus it is possible to reduce the overall number of nuclear DoFs to 400 without any loss in the accuracy of the model. Figure 8(a) shows the convergence behavior of the TT methodology for different ranks of the cores. For sake of simplicity all cores have the same ranks, although different values are, in principles, allowed. As can be seen, the small rank approximation provides a good description of the dynamics only for short times. More specifically for r = 10 the dynamics is accurate up to 60 fs, while for r = 20 the dynamics is in almost quantitative agreement with the exact result up to 250 fs. Increasing r to 60 provides only slight modifications in the long-time tail of the decay. Indeed, after 800 fs the discrepancies between the population decay curves with r = 60 and r = 50 have an average relative deviation of about 5%. In all the calculations a basis set of harmonic oscillator eigenfunctions has been used with maximum quantum number 20 for all the DoFs, which guarantee converged results. 11 Figure 8(b) shows the electronic population of the initial state j B À A H A i as a function of time at different temperatures. Temperature effects are not dramatic, as expected, since most of the vibronic activity is associated with high frequency vibrations. Increasing the temperature from 10 K to 77 K results in a very small decrease of the decay rate. This F I G U R E 7 Effective site spectral densities J p (ω) and J t (ω) describing the coupling of the physical and tilde bosonic DoFs with the electronic subsystem at different temperatures. (a,b) 77 K, (c,d) 300 K effect cannot be described in the framework of the classical ET theory and is very likely due to the highly quantized nature of the vibrational density of states at very low temperatures. A modest increase in the population decay rate upon increasing T from 77 K to 298 K is observed. This effect can be attributed to the increased number of accessible vibronic states. The overall behavior is similar to what is found in other ET processes where large vibronic couplings are associated to highly quantized modes. 97-99

| Theoretical formulation
As discussed in the preceding sections the twin-space formulation of the Liouville-von Neuman equation enables the development of the TFD Schrödinger equation. Furthermore, using this formalism it is possible to define a reduced state vector analogous to the familiar reduced density matrix. This was first discussed by Arimitsu and Umezawa 44 using projection operator techniques and leads to a set of equations which share the same formal structure as the standard reduced density matrix approaches, but make explicit use of tilde operators. We have recently demonstrated that this procedure can be used to derive an equation of motion for the reduced density matrix based on the hierarchical solver technique developed by Tanimura and Kubo. 100,101 Here we briefly sketch the main steps of the derivation and refer the reader to the original article 12 for further details. Let us consider a system described by a Hamiltonian operator where A is the subsystem of interest, B is the "bath," that is, the irrelevant subsystem, V is their coupling, and H 0 = H A + H B . In twin-space formulation the Liouville super-operator is thus given bŷ whereĤ A ¼ H A ÀH A and so on. The HEOM theory by Kubo and Tanimura 100 is certainly one of the most important theoretical developments for the study of open systems. It was originally formulated as a methodology to describe a system interacting with a non-Markovian Gaussian environment, and has later been extended to treat other types of system-bath interactions. 100,102 In its current formulation it has been applied to study the dynamics of a large variety of chemico-physical processes including energy and electron transfer 83,87,103-105 as well as heat transport 106 and electron transport. 7 The reader can refer to the recent work by Tanimura 4 for a more detailed analysis of the HEOM methodology.
A fundamental assumption of HEOM theory is that the system-bath interaction can be factorized as where S k and Q k are system and bath operators respectively. Furthermore, the bath operators are described as a linear combination of position operators q j of harmonic oscillators The coupling super-operator is thus given byV Under this conditions it is possible to demonstrate that 107,108 where jρ A (t)i I is the reduced density matrix of subsystem A obtained by tracing over all the possible states of the "irrelevant" subsystem B, andK V t ð Þ ¼ e iĤ 0 tV e ÀiĤ 0 t being the coupling between subsystems A and B in the interaction picture. Differentiating Equation (74) one obtains 109 It is fundamental to note that the variable τ in the integral above ranges over all times, so that, for τ < s the time ordering operator mixesV τ ð Þ with all the terms of the expansion of the exponential operator exp À R t 0K 2 ð Þ I s ð Þds , making it impossible to obtain an explicit equation for jρ A (t)i I . The HEOM method provides a way to disentangle the above equation in the special case of a Gaussian bath.
Indeed, it is possible to demonstrate that after some easy manipulations the second order cumulant can be written as (see Ref. 12 where S t ð Þ,S t ð Þ, and Q(t) are operators in the interaction picture. In the limit of a continuous distribution of bath modes the correlation function hQ k (t 1 )Q k (t)i can be written as where J k (ω) is the k-th bath spectral density. Furthermore, if Equation (78) can be represented as series of the form c kj e Àγ kj jtÀt 1 j ð79Þ where c kj and γ kj are complex coefficients, the second order cumulant reduces tô in which we have introduced the super-operatorsŜ As a last step of our derivation we define a set of auxiliary state vectors 100,110 where m = {m kj } is a set of non-negative integers. Here, the index k labels the number of independent bath spectral densities J k (ω), and the index j labels the number of terms in the expansion of Equation (79). It is readily verified that the vector jρ A (t)i I , describing the physical state of our system, corresponds to the auxiliary state vector having all indices m kj = 0, that is, j ρ A t ð Þi I ¼j ρ 0 A t ð Þi. The above definition takes into account the scaling factors originally proposed by Shi and coworkers which improves the numerical stability of the final system of equations. 110 HEOM are readily derived upon repeated differentiation of the jρ m i with respect to time. Moving to the Schrödinger representation the set of equations is obtained, where m ± 1 kj = (m 10 , …, m kj ± 1, …), and the explicit time dependence of the auxiliary vectors has been dropped. The price to pay for disentangling the time ordering operation of Equation (74) is that HEOM constitutes an infinite set of first-order ordinary differential equations. Fortunately, using the hierarchy it is possible to devise very efficient truncation schemes which enable obtaining highly accurate results with a finite system. The reader is referred to the original articles for the derivation of an optimal truncation scheme. 100,111 In the above derivation we have not considered low-temperature corrections which can be included straightforwardly from a direct application of the original approach suggested by Ishizaki and Tanimura. 111 To further simplify the structure of the HEOMs we follow Tanimura 102 and introduce a set of vectors jmi = j m 10 m 11 …m 1K m 20 …m MK i, and their corresponding boson-like creation-annihilation operators and the auxiliary density vector and rewrite the hierarchical equations of motion in the compact form with the initial condition given by jη(0)i = j ρ A (0)i j 0i.

| Tensor-train representation of the auxiliary reduced density vector
Tensor train representation can be applied to the density vector jηi in a manner that is similar to that used for the treatment of the thermal Schrödinger equation. If d is the dimensionality of the original Hilbert space of our system, that is the number of DoFs of the H A Hamiltonian operator, and if the dissipative environment is described using M uncorrelated spectral densities J k (ω) each expanded into K Matstubara terms, the vector jη(t)i of Equation (86) can be considered as a tensor with N = 2d + KM indices. Therefore, one possible way to represent jη(t)i in TT format is to employ a product of 2d + MK low rank matrices. For sake of simplicity, in the following only one Matsubara term is considered for each spectral density, K = 1. The generalization to K > 1 requires a slightly more involved notation. If we label with μ k (m k ) the TT core matrices associated with the ith spectral density, and with ρ 2f À 1 (i f ), (ρ 2fjf ) the TT core matrices associated with the fth physical (tilde) DoFs, jη(t)i can be written in TT form as In the above expression only the component with {m k = 0, k = 1, …, M} is required for the computation of physical observables, that is The alternation of the indices i k ,j k provides a very convenient scheme for the computation of expectation values of observables. Indeed, the identity vector can easily be written as a direct product which can be directly translated into the TT format, since, as already mentioned in Section 3.1, the Kronecker product can be considered a TT in which all the ranks of the cores are 1. From the above expression it is immediate to derive TT representations of the vectors corresponding to observables using Equation (3). We notice that a different order of the indices would make the summations not separable increasing the overall computational cost. Finally, the use of tilde operators allows to easily perform grouping of variables that are strongly correlated. 13

| Application of twin-space TT-HEOM to model systems
The methodology described above has already been successfully applied to the study of spin-boson and exciton-polaron systems. Here we briefly report on an illustrative application to a charge transfer problem between two identical molecular sites. Both sites are linearly coupled to a set of seven nuclear vibrations. The parameters of the vibronic model are reported elsewhere 12 and have been used to describe the charge-transfer (CT) process in a pentacene dimer. 112 In Figure 9(a) the population dynamics of the two site system is reported for two different values of the bath reorganization energy λ = 300 and 90 cm À1 . In both cases the characteristic bath frequency is γ = 53 cm À1 . Due to the complexity of the model and to the large number of vibronically active DoFs it is not easy to disentangle the different contributions to the population dynamics. The initial fast decay of the populations is very likely due to pure excitonic couplings 113 while the small oscillations with periods of about 20 fs, clearly evident at longer times, are caused by the vibronic activity of several high frequency modes. For λ = 90 cm À1 the CT dynamics is underdamped while an overall overdamped behavior can be observed for λ = 300 cm À1 .
Finally, the convergence properties of the numerical methodology are illustrated in Figure 9(b) where the populations of the two states of the dimer model as a function of time for different values of the TT compression ranks are compared. Considering the relatively small number of DoFs the ranks necessary to reach a converged dynamics are considerably large. Figure 9(c) also shows the norm of the state vector as a function of time. As can be readily seen for very small ranks the norm drastically decreases, by about 20%, after a very short time and even for very large ranks (r = 115) there is a 1% loss after 800 fs. This behavior can be traced back to the combination of the reduced Liouville equation with the TDVP solver. 114,115 Indeed, this methodology preserves the norm hΨj Ψi and the mean value hΨj Xj Ψi during the evolution governed by Equation (87). Both quantities however have no direct physical meaning. In our formalism the true norm to be preserved during the evolution is As can be expected this problem introduces artifacts which can be alleviated by increasing the ranks of the TT cores. However, this easily becomes a limiting factor since the required TT storage scales quadratically with the TT ranks. Very recently Shi et al., 116 following the original idea of Heller 115 have suggested to correct the norm conservation problem by adding a constraint via an additional Lagrange multiplier.

| FURTHER DEVELOPMENTS
The TFD-TT framework offers ready-to-use methodology for the simulation of quantum dynamics of a broad class of electron-vibrational systems with many DoFs at finite temperature. The methodology can be applied to various polyatomic species and molecular aggregates of practical interest. Apart from that, two other directions should be highlighted. (a) TFD-TT simulations, if converged, provide numerically accurate quantum dynamics. Such simulations can therefore be employed for benchmarking evolutions of selected quantum systems and testing different approximate quantum/quasiclassical/semiclassical simulation protocols. (b) The TFD-TT methodology is valid, without any alterations, for time-dependent Hamiltonians. It is thus tempting to apply it to driven systems, to problems of quantum control, and to simulations of nonlinear femtosecond spectroscopic signals. For example, Beguši c and Vaníček have recently combined TFD machinery with the single-trajectory semiclassical thawed Gaussian approximation, which allowed them to simulate third-order response functions and spectroscopic signals on the fly. 117 The TFD-TT methodology is not only an efficient suite of algorithms and codes for the simulation of quantum dynamics. This methodology, in which temperature is incorporated into the Hamiltonian, offers a new paradigm of the description of molecular species. Any electron-vibrational system specified by the Hamiltonian H, the density matrix ρ(t), and the initial thermal vibrational distribution ρ(0) is equivalently described by the thermal Hamiltonian H θ with twice higher number of vibrational modes possessing temperature-dependent electron-vibrational couplings, but characterized by the thermal wave function jψ θ (t)i prepared initially in the global electron-vibrational vacuum state j ψ θ 0 ð Þi ¼j 0 A i j 0 B0B i. Putting differently, the TFD methodology maps the Schrödinger equation driven by the Hamiltonian H at zero temperature to the TFD Schrödinger equation driven by the thermal Hamiltonian H θ at finite temperature. This analogy is useful for the development of TFD picture of various aspects of quantum dynamics of molecular systems. For example, one can study semiclassical evolution generated by the TFD-version of the van Vleck-Gutzwiller and Herman-Kluk propagators, or introduce temperature-dependent Franck-Condon and Huang-Rhys factors for electron-vibrational transitions. Note also that quantum dynamics driven by the TFD Schrödinger equation can directly be simulated by the MCTDH methods. This observation uncovers deep interconnection between the MCTDH and TFD-TT methods, as well as calls for their application to the same systems for comparing numerical efficiency of the two methods.
As has been noted above, the TFD-TT methodology can readily be generalized to a broader class of Hamiltonians in which site energies ε n , electronic couplings V nm , and electron-vibrational coupling constants g kn are polynomials in vibrational creation-annihilation operators. For example, a conical intersection between the electronic states n and m driven by the vibrational coupling mode j corresponds to V nm $ a j þ a † j . Hence finite-temperature quantum dynamics of a larger number of important for applications molecular systems can be scrutinized by the TFD-TT machinery. The TFD-TT approach can also be extended to rovibrational systems, where vibrational subsystem should be considered as fast, while rotational subsystem should be thermalized and treated as slow.
The present methodology is based on the Born-Oppenheimer approximation of Equation (20), which implies that the electron-vibrational system is initially prepared in the ground state of the electronic subsystem with thermal distribution over the DoFs of the vibrational subsystem. The extension beyond the Born-Oppenheimer approximation, which is valid for initially correlated thermal electron-vibrational systems, has been theoretically developed in Ref. 51 and needs to be tested in explicit simulations.
Finally, it looks promising to apply the TFD-TT methodology to study inter-molecular and intra-molecular vibrational energy redistribution and transfer. These problems, which are conceptually rooted into venerable Fermi-Pasta-Ulam model, require accurate simulation of quantum dynamics of nonlinearly coupled oscillators in the electronic ground state. Very recently, the first step in this direction was made in Ref. 118, where the TT methodology was applied to investigate energy transfer between two localized vibrational modes nonlinearly (via Fermi-resonance) coupled to the chain of harmonic oscillators. Since high-frequency modes were considered only, thermal effects were insignificant in that study. The inclusion of low-frequency modes into similar models and simulation of the ensuing nonlinear multiscale quantum dynamics via TFD-TT methods is currently under way.
As for the twin-space TT-HEOM formulation the field of application can be potentially very similar to that of the TFD-TT approach.

| CONCLUSION
In this work, we have presented a comprehensive overview of recent theoretical development which enables accurate finite-temperature quantum dynamical simulations of systems with multiple electronic and nuclear DoFs. The approach is based on TFD theory and TT decomposition. TFD describes temperature effects by the coupling of the system to a fictitious bosonic bath, so that the number of nuclear DoFs is doubled. In a basis-set representation, the timedependent TFD wave function is an array of the dimension M Â d 2 , where M and d are the number of electronic and vibrational basis DoFs respectively. The density matrix describing the same problem has the dimension of M 2 Â d 2 . The TFD wave function offers, therefore, a more compact way of information storage in comparison with the density matrix, notably for systems with multiple electronic states. Furthermore, high-frequency modes need not be coupled to the fictitious bosonic bath. This leads to additional reduction of the active space and computational savings. A new timedependent TFD Schrödinger equation has been proposed and its numerical solution via TT/MPS representation of the vibronic wave function has been discussed in detail. Our approach to the application of TFD has provided a new view of quantum dynamics at finite temperature, by introducing temperature-renormalized electron-vibrational couplings and spectral densities. This picture shows explicitly which DoFs are thermalized, and which are not.
The use of TT decomposition enables to handle a large number of variables: the storage of order N tensors in TT/ MPS format scales linearly with N. The results of our numerical simulations of energy and charge transfer processes clearly show that the methodology is very accurate and robust, and it can be applied to problems with many DoFs at any temperature. The present approach is based on a basis set representation of the wave function and can in principle be applied to realistic chemical dynamics problems using computed potential energy surfaces.
In standard wave function methodologies, finite temperature effects are taken into account by averaging the quantity of interest over different initial conditions. Therefore, the distribution of initial conditions must be sampled via Monte Carlo methods, and for each initial state a separate dynamical problem must be solved. Both these aspects limit the applications of such methods to large systems at room temperature. Here we have shown that using TFD theory it is possible to describe finite temperature effects without the need to solve a large number of independent dynamical problems with different initial conditions. The increased computational cost, due to the doubled number of nuclear DoFs, can be kept under control by using a TT/MPS representation of the wave function, which is very accurate for the cases discussed in this article.
TT have revealed to be a very promising approximation of the wave function of vibronic systems, allowing to describe molecular processes with multiple coupled electronic states. The very simple structure of the tensor network makes the TT approximation suitable for the development of automatic procedures to optimize the sequence of indices which provide an optimal entanglement growth, although this latter point is still a subject of active research.
Finally, we have shown that by employing the twin formulation of quantum statistical mechanics it is possible to extend the use of numerical methods suitable for wave-packet propagation directly to the study of reduced density matrix evolution in open systems. The proposed approach is based on the twin-space formulation of quantum statistical mechanics and introduces a significant benefit to the actual numerical calculations of expectation values of dynamical operators of system variables. A key advantage of the mathematical description of HEOM reported here is that it allows to fully exploit the TT formalism in the double Hilbert space.

CONFLICT OF INTEREST
The authors have declared no conflicts of interest for this article.

DATA AVAILABILITY STATEMENT
Data sharing is not applicable to this article as no new data were created or analyzed in this study.