;

1 Introduction↩︎

Simulating dynamics of complex systems has been and remains a computational challenge that has given rise to a flowering of algorithms and software. Since measured observables often depend on a limited number of degrees of freedom (DoF), a customary starting point is the partition of the complete system into a main subsystem that gathers the active DoFs and an environment, that is, the system-bath model [1]–[3]. From the quantum mechanical method, the first step is the definition of a model Hamiltonian symbolically separated in three parts related to the relevant subsystem, the reservoir, and the system-bath interaction. Since the work of Caldeira-Leggett [4], the reservoir is usually an ensemble of harmonic oscillators, and the interaction term is bilinear in the system and the bath operators. In this pioneering model, two electronic states (denoted spin, i.e. spin \(1/2\)) are coupled to a harmonic reaction coordinate explicitly included in the system Hamiltonian and coupled to a set of secondary bath oscillators. A great alternative already discussed in a seminal paper on electron transfer [5] restricts the system to the two-electronic states at a given reference nuclear geometry and coupled to a new set of non interacting oscillators obtained by a unitary transformation. This generates the so-called « star model » with all the oscillators (phonon modes or all the molecular and solvent vibrational modes) coupled to the electronic subspace, while the first model is generally referred to as the reaction coordinate mapping method [6]–[11]. The star model has been used to simulate molecular excited state dynamics for instance excitonic states in light-harvesting molecular structures [12], quantum dots [13], or excited states in chromophores coupled by non-adiabatic interstate interactions [14]–[16].

We focus here on the dynamics of electronic excitation transfer via a conical intersection, which leads to strong non-adiabatic coupling [17], [18]. When the dynamics is treated by the wave function approach, the linear vibronic coupling (LVC) model is the basic tool for many simulations [19]–[25] with Multi Layer Multi Configuration Time Dependent Hartree (ML-MCTDH) [26] or Matrix Product State (MPS) computational methods [16], [27]–[29] involving some hundreds of vibrational DoFs. In the open quantum system simulations with the hierarchical equations of motion (HEOM) [30]–[36], the simulation of dynamics via a conical interaction is challenging since at least three baths are involved: two baths make the energy of the electronic states fluctuate and a third bath strongly alters the electronic coupling. A discussion between the two alternative models including or not active modes in the open system was presented in the reference [14]. The partition retaining a two-dimensional subspace of nuclear degrees of freedom where the degeneracy of the conical intersection is lifted coupled to a residual bath is adopted in ref. [10]. The star model where all the nuclear degrees of freedom are in the baths leads to stronger system-bath couplings and requires more computational resources [14], [15].

When the environmental effects are taken into account, the central tool is the spectral density matrix, whose diagonal terms indicate the strength of the system-bath coupling at each bath frequency and off-diagonal terms are related to the correlation of the effects of a mode at a given frequency on two different states. Indeed, the intramolecular modes often show strong correlations, in the sense that displacement of these modes either reduces or increases the energy gap between distinct electronic states. The spectral density matrix may be constructed from discrete \(\textit{ab}\) \(\textit{initio}\) data, broadened to generate a continuous distribution and possibly sampled again to generate fewer modes than in the initial model for wave packets simulations [37], [38]. An alternative is to directly obtain a continuous matrix by transforming correlation functions calculated by molecular dynamics, and then possibly proceed with sampling [28], [29], [39]–[41]. This approach allows the explicit inclusion of solvent effects, vibration damping, and even conformational changes to be accounted for via the effective (continuous) spectral function matrix that contains information about the energy fluctuations arising from all environmental influences.

Bath correlation is ubiquitous in the modelling of condensed-matter systems and has been considered mainly in the context of excitons in multi-site Hamiltonian simulating light-harvesting systems or biological environment [12], [42]–[50], in qubit physics [51], [52] and also condensed matter, i.e. electron-phonon interaction in solids/metals/semiconductors - the basis of pairing in superconductivity. These strong (spatial) correlations in bath fluctuations naturally arise from the extended (delocalized, plane wave) nature of lattice phonons [53], but can be still simulated in thermalized time-evolving density operator with orthogonal polynomials algorithm (TEDOPA) by including long-range interactions [54]. The problem arises both in sample-based wave packet simulations by MCTDH [17], [38] or MPS [28], [29] and in dissipative dynamics with memory kernels treated by HEOM [30]–[36], [55]–[58].

In this work, we address the central issue of bath de-correlation procedure by HEOM and T-TEDOPA formalisms. A first idea could be to propose a transformation of the collective bath modes to diagonalize the correlation matrix or the associated spectral density matrix. However, whether the transformation used to diagonalize the spectral density matrix is frequency dependent or frequency independent is a crucial point within the HEOM framework. We shall mention situations where this diagonalization is possible and propose different approximate strategies useful to reduce the computational cost by involving a simple set of independent environments compatible with the standard HEOM formulation. Interestingly, we will show that this frequency dependence could be taken into account in T-TEDOPA. We will compare this procedure directly based on the frequency dependent diagonalization with a recent alternative method using T-TEDOPA in a simulation of a conical intersection[59]–[61]. The goal is then to extract parameters of a system-bath coupling terms from the spectral density matrix calculated by molecular dynamics. This strategy leads to an operator that explicitly couples a given electronic state to different baths as it is the case for spatial correlations. This illustrate the flexibility to write the system-bath coupling with transformed modes and operators. This will be the heart for the de-correlation procedure in the HEOM or T-TEDOPA formalisms that have different numerical constraints. As illustration, we reconsider the LVC model of a conical intersection (sec.7) in a symmetric meta-substituted planar poly(phenylene-ethylene) (PPE) dimer already used in different applications with MCTDH [22] or HEOM [62], [63]. The dynamics of two excited states depend on three baths: the two tuning baths that make the energy gap with the ground state fluctuate and the coupling bath, which causes the variation of the electronic interstate coupling. In recent HEOM simulations [62], [63], the two correlated tuning baths were replaced by a single shared bath with a modified system operator by an approximate frequency independent transformation of the baths. Here, we intend to justify this approach.

We first simulate the excited state dynamics by HEOM (secs.3 and 4) including all the crossed terms which is rarely done in HEOM simulations [55], [56] where approximate treatments are used [12], [42], [46], [62], [63]. Then, we justify and discuss the validity of the approximate strategy, which is very useful to reduce the computational resources in HEOM simulations (sec.5). The efficiency of the T-TEDOPA method to apply the exact frequency dependent transformation leading to a single shared bath in favorable cases is presented in sec.6. This strategy in T-TEDOPA is then an alternative to both the exact correlated HEOM method and to another approximate de-correlation procedure recently proposed [28], [29], [64].

2 Fully (anti) correlated baths↩︎

We start with a simple example of a two-level system described by the LVC model in which each state \(\left| n \right\rangle\) is coupled to a bath of oscillators as schematized in fig.1. The partition of the full model leads to the generic Hamiltonian \[H={{H}_{S}}+{{H}_{B}}+{{H}_{SB}} \label{eq:Htotal}\tag{1}\] where \({{H}_{S}}=\sum\nolimits_{n=1}^{2}\epsilon_n{\left| n \right\rangle }\left\langle n \right|\) is the system Hamiltonian, \({{H}_{B}}\) is the environment Hamiltonian, i.e., the sum of the harmonic baths \({{h}_{n}}=\frac{1}{2}\sum\nolimits_{j}{\left( p_{n,j}^{2}+\omega _{n,j}^{2}q_{n,j}^{2} \right)}\) expressed in mass weighted coordinates. in the LVC model, the \(q_j\) are normal modes assumed to be the same in each electronic state. By discarding any electronic interstate coupling, the two tuning baths are coupled only diagonally to the system, the system-bath coupling may always take the form of a sum of bilinear terms \[{{H}_{SB}}=\sum\nolimits_{n=1}^{2}{{{S}_{n}}}{{B}_{n}} \label{HSBnew}\tag{2}\] where the system-bath coupling operators are projectors on the basis states \({{S}_{n}}=\left| n \right\rangle \left\langle n \right|\). The collective bath modes are linear combinations of the bath coordinates \[B_{n}=\sum\nolimits_{j}{\kappa_{j}^{(n)}q_{j}} \label{eq:Collectivebathmode}\tag{3}\] where \(\kappa _{j}^{(n)}\) is the gradient making the energy of the corresponding state fluctuate.

The fact that \(H_{SB}\) is a simple sum of terms does not mean that the baths are not correlated. Indeed, even if the normal modes \(q_k\) are independent harmonic oscillators, the operators \(B_i\) and \(B_j\) can still be correlated whenever they share common modes. The two baths emerging from the LVC Hamiltonian are naturally correlated: they are not distinct physical environments, but different projections of the same vibrational environment onto different electronic states.

The Gaussian statistical properties of the harmonic baths are characterized by the two-time correlation matrix \[{{C}_{lm}}(t-\tau )=\left\langle {{B}_{l}}(t){{B}_{m}}(\tau ) \right\rangle \label{eq:C40t41new}\tag{4}\] where \(l,m\in \left[\!\left[ 1,2 \right]\!\right]\) and \(\left\langle \centerdot \right\rangle\) is the average over the equilibrated baths at temperature \(T\) described by \[\rho _{B}^{eq}(T)=\exp \left( -\beta {{H}_{B}} \right)/T{{r}_{B}}\left[ \exp \left( -\beta {{H}_{B}} \right) \right] \label{eq:rhoeq}\tag{5}\] with \(\beta =1/{{k}_{B}}T\) and \({{k}_{B}}\) is the Boltzmann constant.

In the frequency space, the correlation matrix is linked to the temperature-dependent spectral density \[{{C}_{lm}}(t-\tau )=\frac{\hbar }{\pi }\int_{-\infty }^{\infty }{d\omega J_{lm}^{\beta }(\omega )}{{e}^{-i\omega (t-\tau )}} \label{eq:TFC40t41new}\tag{6}\] with \[J_{lm}^{\beta }(\omega )={{J}_{lm}}(\omega ){{n}_{\beta }}(\omega ) \label{eq:Jbeta}\tag{7}\] where \({{n}_{\beta }}(\omega )\) is the Bose function and \(J_{lm}^{{}}(\omega )\) is the temperature-independent spectral density that is first defined in a discrete representation: \[{{J}_{lm}}(\omega )=\frac{\pi }{2}\sum\nolimits_{j}{{g_{j}^{(l)}g_{j}^{(m)}}}\,\delta (\omega -{{\omega }_{j}}) \label{eq:spectraldensity}\tag{8}\] where \(g_{j}^{(l)}\) are the system-bath coupling given by \(\kappa _{j}^{(l)}\omega_j^{-1/2}\) for \(l=1,2\). The spectral density matrix contains information about the energy fluctuations arising from all environmental influences, including intramolecular or solvents modes. Intramolecular modes often show strong correlations, although some fluctuations affect the electronic states in effectively random (uncorrelated) ways (especially true of the solvent contributions). When constructing the spectral density matrix from the LVC model coefficients, the off-diagonal element \(J_{12}(\omega)\) becomes non-zero because the two electronic states couple to overlapping subsets of vibrational coordinates. However, when the spectral densities are smoothed by a broadening procedure, the continuous off-diagonal terms may become small due to cancellations between positive and negative mode-specific contributions. As a result, the effective bath correlation perceived by the dynamical model can become significantly weaker than the underlying LVC model.

The interesting inverse procedure followed in refs.[28], [29] starts from the spectral density matrix known from molecular dynamics calculations, including both the diagonal terms and the cross term \(J_{12}(\omega)\). The goal is then to reconstruct an explicit microscopic bath Hamiltonian that reproduces these correlations. To do this, one introduces a set of harmonic modes and assigns coupling coefficients such that each mode can couple simultaneously to several system operators. The \(H_{SB}\) operator is then similar to the one used in non-local correlations when a given electronic state is coupled to oscillators on different sites. [12], [42]. In practice, this means finding a set of couplings whose reproduce the spectral density matrix. Conceptually, this is a factorization problem that may be realized by different ways. The shared modes generated by this factorization are what produce the correlated baths.

One may consider two extreme situations where the baths are fully correlated, this means the collective modes are always in phase or fully anti correlated if they have opposite phases. The gradients are equal in each bath and have the same sign in the first case or opposite signs in the second one (see fig.1). Assuming fully correlated or fully anti-correlated baths means \(g_{j}^{(1)}=\pm g_{j}^{(2)}\;\forall j\). The \({{J}_{lm}}(\omega )\) matrix is then of rank 1 for each frequency and it may be diagonalized by a single frequency-independent matrix \[{{\mathbf{J}}^{\beta }}(\omega )=\mathbf{U}\mathbf{\Lambda} (\omega ){{\mathbf{U}^T}} \label{eq:DiagJbeta}\tag{9}\] where \(\mathbf{\Lambda (\omega )}\) is the diagonal matrix of the eigenvalues, which contains only one non-zero eigenvalue in this case. By this orthogonal transformation one may build new collective bath operators \({{\tilde{B}}_{j}}(t)=\sum\limits_{n}{U_{nj}}{{B}_{n}}(t)\) that remain uncorrelated at all times. Moreover, only the new bath associated with the non-zero eigenvalue is coupled to the system. The two initial baths are compressed in a single active shared bath, the other one being decoupled from the system. In the fully correlated or fully anti-correlated case, the shared bath is \[{{{B}}_{\pm }}(t)=\frac{1}{\sqrt{2}}\left( {{B}_{1}}(t)\pm {{B}_{2}}(t) \right) \label{eq:B43-}\tag{10}\] with the \(+\) or \(-\) sign respectively. The \(H_{SB}=\boldsymbol{s}^T\boldsymbol{b}\) operator (with \(\boldsymbol{s}\) and \(\boldsymbol{b}\) the column vectors of the operators and of the collective bath modes) may be written \(H_{SB}= \boldsymbol{s}^T\mathbf{U}\mathbf{U}^T\boldsymbol{b}\). The transformed system operator is then \[{{{S}}_{\pm }}=\frac{1}{\sqrt{2}}\left( {{S}_{1}}\pm {{S}_{2}} \right) \label{eq:Stilde}\tag{11}\] (see fig.1). Note that the standard spin boson model generally uses a single shared bath coupled to the system by the \({{\sigma }_{z}}\) Pauli matrix, which corresponds to a fully anti-correlated case.

3 Cross-correlation in master equation↩︎

We now consider a generic d-level system coupling to \(N_{bath}\) harmonic baths \[{{H}_{SB}}=\sum\nolimits_{\alpha }^{{{N}_{bath}}}{{{S}_{\alpha }}}{{B}_{\alpha }}. \label{eq:HSB}\tag{12}\] The subsystem is described by the reduced density matrix \({{\rho }_{S}}(t)=T{{r}_{B}}\left[ \rho (t) \right]\) obtained by tracing the bath DoFs. In the interaction representation, the evolution of the system density matrix \({{\tilde{\rho }}_{S}}(t)\) is driven by a dynamical map, which involves a positive time ordered exponential denoted by \({{\mathcal{T}}_{+}}\) \[{{\tilde{\rho }}_{S}}(t)={{\mathcal{T}}_{+}}{{e}^{\int\limits_{0}^{t}{d\tau \left\langle L(\tau ) \right\rangle }}}{{\tilde{\rho }}_{S}}(0) \label{eq:rhoS}\tag{13}\] where \(\left\langle \centerdot \right\rangle\) means trace over the bath DoFs. In this representation, one has \(L(t)\centerdot =\frac{-i}{\hbar }\left[ {{H}_{SB}}(t),\centerdot \right]\) with \({{H}_{SB}}(t)=\sum\nolimits_{j}^{{{N}_{bath}}}{{{S}_{j}}(t)}{{B}_{j}}(t)\), \({{S}_{j}}(t)={{e}^{i{{H}_{S}}t}}{{S}_{j}}{{e}^{-i{{H}_{S}}t}}\) and \({{B}_{j}}(t)={{e}^{i{{h}_jt}}}{{B}_{j}}{{e}^{-i{{h}_jt}}}\). According to the Wick theorem, the cumulant expansion of Eq.(13 ) is exactly truncated at the second order if the statistics is Gaussian as it is the case for harmonic baths. Therefore, Eq.(13 ) becomes \[{{\tilde{\rho }}_{S}}(t)={{\mathcal{T}}_{+}}{{e}^{{{F}^{(2)}}(t,{{t}_{0}})}}{{\tilde{\rho }}_{S}}(0) \label{eq:rhocumul}\tag{14}\] with the second order cumulant \[{{F}^{(2)}}(t,{{t}_{0}})=\int_{{{t}_{0}}}^{t}{d\tau }\int_{{{t}_{0}}}^{\tau }{d\sigma }{{K}^{(2)}}(\tau ,\sigma )\] where the memory kernel is defined involving only two-time correlation functions \[\begin{align} & {{K}^{(2)}}(t,\tau )\centerdot =T{{r}_{B}}[L(t)L(\tau )\rho _{B}^{eq}\centerdot ] \nonumber \\ & =-\frac{1}{{{\hbar }^{2}}}\left\langle \left[ {{H}_{SB}}(t),\left[ {{H}_{SB}}(\tau ),\centerdot \right] \right] \right\rangle \label{eq:kernel} \end{align}\tag{15}\] and the trace is made over each bath. The master equation derived from the mapping of Eq.(14 ) reads: \[{{\dot{\tilde{\rho }}}_{S}}(t)={{\mathcal{T}}_{+}}\int\limits_{0}^{t}{d}\tau {{K}^{(2)}}(t,\tau ){{\tilde{\rho }}_{S}}(t).\] By inserting Eq.(15 ) in Eq.(14 ), one gets : \[\begin{align} & {{\hbar }^{2}}{{{\dot{\tilde{\rho }}}}_{S}}(t) \\ & =-\sum\limits_{l,m=1}^{{{N}_{bath}}}{\left[ {{S}_{l}}(t), \right.}\int\limits_{0}^{t}{d}\tau \left( {{C}_{lm}}(t,\tau ){{S}_{m}}(\tau ){{{\tilde{\rho }}}_{S}}(t) \right. \nonumber \\ & \left. -{{C}_{ml}}(\tau ,t){{{\tilde{\rho }}}_{S}}(t){{S}_{m}}\left. (\tau ) \right) \right]. \label{eq:mastereq} \end{align}\tag{16}\] Even if the system-bath coupling \({{H}_{SB}}\) (Eq.(12 )) is the sum of bath contributions that seem independent, the master equation (Eq.(16 )) rigorously involves the full correlation matrix of the bath collective modes \({{C}_{lm}}(t,\tau )\) (Eq.(4 )). The statistical properties are well known [55], [56] (they are summarized in the supplementary material). For two particular modes \(m,l\) the main relations are : \[{{C}_{ml}}(\tau ,t)={{C}_{lm}}(t-i\hbar\beta,\tau )=C_{lm}^{*}(t,\tau ).\]

4 HEOM for correlated baths↩︎

The crucial point of the HEOM method is the expansion of the correlation functions in terms of oscillatory decaying complex exponential functions [33], [65]–[68], each term being associated with an artificial decaying mode \[{{C}_{lm}}(t,\tau )=\sum\nolimits_{k}^{{{K}^{(lm)}}}{\alpha _{k}^{(lm)}}{{e}^{i\gamma _{k}^{(lm)}(t-\tau )}} \label{CdeTexpansion}\tag{17}\] and \[\begin{align} & {{C}_{ml}}(\tau ,t)=C_{lm}^{*}(t,\tau ) \nonumber \\ & =\sum\nolimits_{k}^{{{K}^{(lm)}}}{\bar{\alpha }_{k}^{(lm)}}{{e}^{i\gamma _{k}^{(lm)}(t-\tau )}} \label{CdeTconjg} \end{align}\tag{18}\] (The possibility of keeping the same exponent in expressions (17 ) and (18 ) is illustrated in the supplementary material [69]). The time non-local master equation (Eq.(16 )) is replaced by time local coupled equations among auxiliary density operators (ADOs) having the dimension of the reduced system matrix. The initial condition is assumed to be factorized \({{\rho }_{tot}(t=0)}={{\rho }_{S}}(t=0)\rho _{h_0}^{eq}(T)\) where \(\rho _{h_0}^{eq}(T)\) is the equilibrated bath interacting with the ground state with harmonic bath \(h_0\). This is justified for an ultra-short preparation by a Franck-Condon (FC) process that projects the ground equilibrium bath in each excited electronic state [12]. The ADOs will describe the behavior of the Gaussian environment. With the FC preparation, they are initially set to zero. Each auxiliary matrix \({{\tilde{\rho }}_{\mathbf{m}}}(t)\) is denoted by a vector of indices giving the occupation number in each artificial mode of the diagonal and off-diagonal correlation functions: \[\mathbf{m}=\left( n_{1}^{(11)}..n_{K}^{(11)},..,n_{1}^{(ll)}..n_{K}^{(ll)},..n_{1}^{(lm)}..n_{K}^{(lm)} \right)\] with \(l,m\in \left[\!\left[ 1,{{N}_{bath}} \right]\!\right]\) where \({{N}_{bath}}\) is the number of baths in the LVC model. Everything happens as if one considers an extended ensemble of baths. When the uncorrelated case involves \({{N}_{bath}}\) baths with spectral densities \({{J}_{ll}}(\omega )\) and \({{J}_{mm}}(\omega )\), one now adds the baths \({{J}_{lm}}(\omega )\) and \({{J}_{ml}}(\omega )\). Note that some uncorrelated baths may not contribute to the off-diagonal part, as it is assumed here in the LVC case for the interstate coupling bath.

The coupled equations in the interaction representation are written by using compact operators: \[\begin{align} & {{{\dot{\tilde{\rho }}}}_{\mathbf{m}}}(t)=i\sum\limits_{l,m=1}^{{{N}_{bath}}}{\sum\limits_{k=1}^{{{K}^{(lm)}}}{n_{k}^{(lm)}}}\gamma _{k}^{(lm)}{{{\tilde{\rho }}}_{\mathbf{m}}}(t) \nonumber \\ & -\sum\limits_{l=1}^{{{N}_{bath}}}{{{\Phi }_{l}}(t)\sum\limits_{m=1}^{{{N}_{bath}}}{\sum\limits_{k=1}^{{{K}^{(lm)}}}{\tilde{\rho }_{\mathbf{m}_{k}^{+}}^{{}}(t)}}} \nonumber \\ & -\frac{1}{{{\hbar }^{2}}}\sum\limits_{l,m=1}^{{{N}_{bath}}}{\sum\limits_{k=1}^{{{K}^{(lm)}}}{n_{k}^{(lm)}}}\Theta _{k}^{(lm)}\tilde{\rho }_{\mathbf{m}_{k}^{-}}^{{}}(t) \label{eq:HEOMeq} \end{align}\tag{19}\] where \(\boldsymbol{m}_k^{\pm}={\left\{ ...n^{(lm)}_{k}\pm 1.... \right\}}\) represents the index vector where the occupation number of the \(k^{th}\) mode of bath \((lm)\) increases or decreases by one unit while all the other numbers remain unchanged. \({{\Phi }_{l}}(t)\centerdot =-\left[ {{S}_{l}}(t),\centerdot \right]\) with \(l=1,{{N}_{bath}}\) and \(\Theta _{k}^{(lm)}(t)\centerdot =\alpha _{k}^{(lm)}{{S}_{m}}(t)\centerdot -\bar{\alpha }_{k}^{(lm)}\centerdot {{S}_{m}}(t)\) with \(k=1,{{K}^{(lm)}}\). As illustrated in fig.2, in the uncorrelated case, the terms \({{\Phi }_{l}}(t)\centerdot\) involves a commutator of the system operator \({{S}_{l}}(t)\) only with raising ADOs of the corresponding bath with correlation function \({{C}_{ll}}(t,\tau )\) while in the cross-correlated case one must add commutators of this \({{S}_{l}}(t)\) with raising terms of the additional baths described by \({{C}_{lm}}(t,\tau )\). Similarly, without cross-correlation, only diagonal terms \(\Theta _{k}^{(ll)}(t)\centerdot\) arise. They involve the parameters of \({{C}_{ll}}(t,\tau )\) and the system operator \({{S}_{l}}(t)\) that acts on the lowering ADOs of this bath \(ll\). With cross-correlation, the supplementary term \(\Theta _{k}^{(lm)}(t)\centerdot\) involves the parameters of \({{C}_{lm}}(t,\tau )\) with action on lowering ADOs of the bath \(lm\) but with the system operator \({{S}_{m}}(t)\). Similarly, \(\Theta _{k}^{(ml)}(t)\centerdot\) implies bath \(ml\) and the system operator \({{S}_{l}}(t)\).

It is obvious that taking into account the cross-correlation involves \(N_{bath}^2\) so that the number of ADOs rapidly grows and computational resources become prohibitive when the hierarchy level is high. We now revisit and discuss a strategy to reduce the number of baths by introducing approximate uncorrelated shared baths.

5 Approximate de-correlation for HEOM↩︎

When the baths are not fully or anti-correlated as assumed in sec.2, the central issue in bath de-correlation is whether the transformation used to diagonalize the spectral density matrix is frequency dependent or frequency independent within the HEOM framework. We will show that this frequency dependence could be taken into account in T-TEDOPA. In general, the spectral density matrix \(J_{ij}(\omega)\) can always be diagonalized at each frequency, but if the corresponding transformation varies with \(\omega\), the bath operators become frequency dependent mixtures of the original operators. In that case, the transformed baths cannot generally be represented by a simple set of independent harmonic environments compatible with the standard HEOM formulation.

As a consequence, the practical possibility of bath de-correlation strongly depends on the spectral representation itself and we will distinguish the discrete or continuous representation.

(a) Discrete case. For discrete LVC spectral densities on a single site, the correlation matrix is \(J_{lm}^{\beta }(\omega_j )\), which is linked to the discrete spectral density \(J_{lm}^{{}}(\omega_j )\) (Eq.(7 )). When the baths are at the same temperature, we may discuss \(J_{lm}^{{}}(\omega_j )\) only. At each frequency \(\omega_j\), the latter is a dyadic matrix \[\mathbf{J}({{\omega }_{j}})={{\mathbf{g}}_{j}}\otimes {{\mathbf{g}}_{j}}\] built from the vectors \({{\mathbf{g}}_{j}^T}=\left( g_{j}^{(1)},..,g_{j}^{({{N}_{bath}})} \right)\) (Eq.(8 )). Accordingly, \(J_{lm}^{\beta }(\omega_j )\) has a single non-zero eigenvalue \[\begin{align} & \mathbf{\bar{J}}_{{}}^{\beta }(\omega_j )={{\mathbf{U}^T}}(\omega_j )\mathbf{J}_{{}}^{\beta }(\omega_j )\mathbf{U}(\omega_j ) \nonumber \\ & =diag({{\Lambda }_{1}}(\omega_j ),0..0) \end{align}\] with \({{\Lambda }_{i}}({{\omega }_{j}})=0,\forall i\in \left[\!\left[ 2,{{N}_{bath}} \right]\!\right]\). In the discrete case, there is only one non-zero eigenvalue \(\Lambda_1(\omega_j)\) at each frequency.

(b) Continuous case. After spectral broadening and possible fit by analytical expressions, however, positive and negative contributions may cancel over finite frequency windows, leading to a much smaller continuous cross spectral density \(J_{12}(\omega)\). The resulting diagonalization can therefore differ substantially between the discrete and continuous representations, potentially giving more non-zero eigenvalues or making the baths appear effectively uncorrelated only after smoothing. A deviation from this ideal dyadic behavior may also occur when the spectral densities are Fourier transforms of correlations functions computed by molecular dynamics [28], [29], [64].

De-correlating the baths is then in principle possible by working with the transformed collective modes and the transformed system operators as done in Sec.2. This would lead to a single effective shared bath associated to the \(\Lambda_1\) eigenvalue in the ideal case that preserves the dyadic structure of the discrete matrix. Working with uncorrelated baths and, even better, with a single shared bath, would substantially reduce the computational resources in the HEOM simulations. However, the frequency dependence of the system operators \(\bar{\boldsymbol{s}}=\boldsymbol{U}^T(\omega)\boldsymbol{s}\) and of the combination of the bath modes \(\bar{\boldsymbol{b}}=\boldsymbol{U}^T(\omega)\boldsymbol{b}\) is not suitable for HEOM. The goal is to find the optimal frequency independent transformation matrix \(\boldsymbol{U}_{\text{opt}}\) to minimize the correlations while ensuring a good approximation of the dynamics \[{\bar{\mathbf{J}}^{\beta}_{{{U}_{\text{opt}}}}(\omega)}=\mathbf{U}^T_{\text{opt}}\mathbf{J^{\beta}(\omega )}{{\mathbf{U}}_{\text{opt}}}. \label{eq:Jtransfoopt}\tag{20}\] The quality of \(\boldsymbol{U}_{\text{opt}}\) may be estimated by the remainder of the off-diagonal terms \[\varepsilon ({{\mathbf{U}}_{\text{opt}}})=\int{d\omega }{{\sum\limits_{i\ne j}{\left| {{\left( \mathbf{U}_{\text{opt}}^T{{\mathbf{J}}^{\beta }}(\omega ){{\mathbf{U}}_{\text{opt}}} \right)}_{ij}} \right|}}^{2}} \label{eq:error}\tag{21}\] Different possibilities may be examined:

(a) The approximate constant \(\boldsymbol{U}_{\text{opt}}\) matrix is denoted with an index \(\boldsymbol{U}_{\eta}\) when it is chosen at a representative frequency, for instance that of the dominant peak in the spectral density. This assumes that the spectral densities keep the same shape with a constant ratio among the couplings. For instance, in the two-bath model \({{g}^{(2)}}(\omega )=\pm \sqrt{\eta }{{g}^{(1)}}(\omega ) \space \,\forall \omega\) with \(\eta > 0\) and a constant sign. The cases where \(\sqrt{\eta }= 1\) are the fully correlated (if \({{g}^{(2)}}(\omega )= + {{g}^{(1)}}(\omega )\)) or fully anti-correlated (if \({{g}^{(2)}}(\omega )= - {{g}^{(1)}}(\omega )\)) baths, as discussed in Sec.2.

(b) Another choice, referred to as \(\boldsymbol{U}_\text{aver}\) does not postulate the constraint of a constant ratio between the couplings. The matrix is built from the average value of the eigenvectors of \(\mathbf{{J}}_{{}}^{\beta }(\omega )\).

(c) A third possibility is to take \(\boldsymbol{U}_{C(0)}\) given by the eigenvectors of the initial value of the correlation matrix \(\mathbf{C}(t=0)\). Indeed, each element is an average of the corresponding matrix element of \(\mathbf{J}^{\beta }(\omega )\) since it is its integral as a result of the Fourier transform relation (Eq.(6 )) .

(d) Finally \(\boldsymbol{U}_{\text{opt}}\) may be estimated by \(\boldsymbol{U}_\text{PCA}\) obtained by the principal component analysis (PCA) by making the singular value decomposition (SVD) of the \(N_{bath} \times N_{bath} \times N_{\omega}\) tensor where \(N_{\omega}\) is a large sampling of the frequency domain.

The exact and approximate strategies for a HEOM simulation are schematized in fig.3 for a model discussed below where two correlated tuning baths are diagonally coupled and one uncorrelated coupling bath is off-diagonally coupled to the excited states. The exact simulation involves the HEOM equations with correlated baths (Eqs.(19 )), illustrated in fig.3 (a). In some favorable cases, it can be followed with an approximate diagonalization of \(\mathbf{J}^{\beta }(\omega )\) by a frequency independent \(\mathbf{U}_{\text{opt}}\) matrix that provides two decoupled shared baths described by the new spectral functions \(\Lambda_1(\omega)\) and \(\Lambda_2(\omega)\) with in the ideal case \(\Lambda_1(\omega) >> \Lambda_2(\omega)\) so that \(\Lambda_2(\omega)\) may be discarded (fig.3 (b)).

In Eq.(21 ), we have estimated the validity of the approximate frequency independent transformation by the integrated Frobenius projection error. For a deeper estimation of the error, it could be interesting to generalize the error bound proposed in ref.[70] to quantify the impact on the expectation value of observables due to a variation of a single spectral density \(\Delta J(\omega)\). The error is expected to grow exponentially with time and is quantified by parameters that are related the correlation kernel (Eq.(6 )) by transforming \(\Delta J(\omega)\) in place of \(J(\omega)\). In the case of multiple correlated baths, the residual correlation matrix \(\mathbf{J^{\beta}(\omega )-{\bar{\mathbf{J}}^{\beta}_{{{U}_{\text{opt}}}}(\omega)}}\) (Eq.(20 )) is the source of error and the integrated strength of this residual bath correlation kernel derived from Eq.(6 ) could be involved in the error estimation as in the single-bath result.

6 T-TEDOPA for correlated baths↩︎

By the efficient T-TEDOPA mapping of Tamascelli et al.[59], [71], the initial thermalized mixed state is replaced by a pure state of the environment at 0 K. The initial factorized state becomes: \[\rho \left( t=0 \right)\text{ }={{\rho }_{S}}\left( t=0 \right)\otimes \ket{0 \ldots 0} \bra{0 \ldots 0}\] where \(\ket{0 \ldots 0}\) is the state with all the environmental modes in their ground state. The modes are defined by sampling the temperature-dependent spectral density \(J^{\beta}(\omega)\) on the extended frequency domain comprising negative frequencies. In the basic case where a given electronic state \(n\) is coupled to a unique bath by the coupling operator \({{S}_{n}}=\left| n \right\rangle \left\langle n \right|\), the system bath interaction in the continuous frequency domain may be written: \[{{H}_{SB}}=\int{d\omega\sqrt{{{J}^{\beta}_{nn}}(\omega )}}{{S}_{n}}q_{B_n}(\omega ).\] The infinite bath is mapped onto a discrete set of chain modes by an unitary transformation built with polynomials orthogonal with respect to the measure of the spectral density \(J_{nn}^{\beta}(\omega)d\omega\) [72]. The system interacts only with the first mode of the chain. All the coupling coefficients connecting the neighboring modes are obtained by a procedure based on the three terms recurrence relation of the polynomials [60], [72], [73]. The length of the chain must be adapted in order to avoid the return of the excitation to the system. The possible excitation in each mode is fixed by the Fock space of dimension \(d\). This configuration leads to the MPS representation formed by a chain of tensor \(d\)-dimensional cores bounded by the bond dimension \(r\).

When the system interacts with two correlated tuning baths, a transformation of the \(H_{SB}\) operator has been proposed recently [28], [29], [64]. This transformation uses a particular matrix \(\boldsymbol{A}(\omega)\) that factorizes the spectral density as followed: \[{{\mathbf{J}}^{\beta }}(\omega )=\mathbf{A}(\omega ){{\mathbf{A}}^{T}}(\omega ).\] The \(\boldsymbol{A}(\omega)\) matrix transforms the coupling operators \(\bar{\boldsymbol{s}}_A(\omega)=\boldsymbol{A}(\omega)\boldsymbol{s}\) and the bath modes \(\bar{\boldsymbol{b}}_A=(\boldsymbol{A}^{T})^{-1}\boldsymbol{b}\). The choice of the factorization matrix \(\boldsymbol{A}\) is not unique. In this proposal, the \(\boldsymbol{A}\) matrix is (by dropping the writing of \(\omega\)): \[\mathbf{A}=\left( \begin{matrix} \sqrt{{{G}_{1}}} & \sqrt{C{{G}_{2}}} \\ \sqrt{C{{G}_{1}}} & \sqrt{{{G}_{2}}} \\ \end{matrix} \right) \label{eq:matA}\tag{22}\] with \[{{G}_{1}}={{\left( J_{12}^{\beta } \right)}^{2}}/\left[ 2J_{22}^{\beta }(1-R) \right],\] \[{{G}_{2}}={{\left( J_{12}^{\beta } \right)}^{2}}/\left[ 2J_{11}^{\beta }(1-R) \right],\] \[C=(2J_{11}^{\beta }J_{22}^{\beta }/{{(J_{12}^{\beta })}^{2}})(1- R)-1,\] and \(R={{\left( 1-{{(J_{12}^{\beta })}^{2}}/(J_{11}^{\beta }J_{22}^{\beta }) \right)}^{1/2}}\). The two baths are independent and each bath couples to the two states by different spectral densities as illustrated in fig.4 \[\begin{align} & {{H}_{SB}}=\int{d\omega }\left( {{A}_{11}}(\omega ){{S}_{1}}+{{A}_{12}}(\omega ){{S}_{2}} \right){{q}_{{{B}_{A1}}}}(\omega ) \nonumber \\ & +\int{d\omega }\left( {{A}_{21}}(\omega ){{S}_{1}}+{{A}_{22}}(\omega ){{S}_{2}} \right){{q}_{{{B}_{A2}}}}(\omega ). \end{align}\] For each bath, the chain mapping is done with respect to \(G_1(\omega)\) and \(G_2(\omega)\) respectively. As a result, interactions between the other state and the chain modes may become long-ranged [28], [64].

Inspired by the procedure relevant in HEOM, one could propose to use the factorization with the orthogonal matrix built with the spectral decomposition of \(\mathbf{J}(\omega)\) and the diagonal matrix \(\mathbf{\Lambda }(\omega )\) \[\begin{align} & \mathbf{J}_{{}}^{\beta }(\omega )=\mathbf{U}(\omega )\mathbf{\Lambda }(\omega ){{\mathbf{U}^T}}(\omega ) \nonumber \\ & ={{\mathbf{A}}_{D}}(\omega )\mathbf{A}_{D}^{T}(\omega) \end{align}\] with \[\boldsymbol{A}_D(\omega)=\boldsymbol{U}(\omega)\sqrt{{\mathbf{\Lambda}}(\omega)}. \label{eq:matAD}\tag{23}\]

By using \(\boldsymbol{A}_D(\omega)\) for the transformation, the system operator becomes \(\bar{\boldsymbol{s}}_{A_D}(\omega)=\boldsymbol{A}_D(\omega)\boldsymbol{s}\). One finds a similar scheme with two uncorrelated baths, each of which is coupled to the two states with different spectral densities. However, in a favorable case when the dyadic structure of the spectral density matrix is nearly conserved, a single eigenvalue \(\Lambda_1(\omega)\) dominates (\(\Lambda_1(\omega)>>\Lambda_2(\omega)\)) and a single shared bath might be sufficient with \[{{H}_{SB}}=\int{d\omega }\left( {{c}_{1}}(\omega ){{S}_{1}}+{{c}_{2}}(\omega ){{S}_{2}} \right){{q}_{{{B}_{A_{D1}}}}}(\omega ) \label{eq:HSBTEDOPA}\tag{24}\] where \({{c}_{n}}(\omega )=\sqrt{\Lambda_{1}(\omega )}U_{n1} (\omega )\) (\(n=1,2\)) and \({q}_{{{B}_{A_{D1}}}}\) represents the new bath coordinate. The eigenvalue is included in the system operator while it is in the transformed correlation function in HEOM. It is worth noting that the T-TEDOPA procedure allows to account for the full frequency dependence of the transformation \(\boldsymbol{U}(\omega)\) at the cost of the long-range interactions. It is an interesting alternative to exact correlated HEOM simulations mainly when a single eigenvalue of the spectral density matrix dominates. If \(\Lambda_2(\omega)\) is not negligible, a second term must be added in Eq.(24 ) and the procedure becomes similar to the other ones using \(\boldsymbol{A}(\omega)\) with two baths interacting with the two states. However, it is not obvious which of the transformations \(\boldsymbol{A}\) or \(\boldsymbol{A}_D\) is the best. In principle, the spectral transformation with \(\boldsymbol{A}_D\) de-correlates the baths. The \(\boldsymbol{A}\) procedure also uses independent baths but with a rotation of the spectral modes. The \(\boldsymbol{A}\) matrix may be written \(\boldsymbol{A}=\mathbf{U}\sqrt{\mathbf{\Lambda}}\mathbf{O}\) where \(\mathbf{O}\) is an orthogonal matrix. This combination of spectral modes may induce some dynamical correlation. The length of the chains depends on the shape of the coefficients of the \(S_n\) operators provided by the different transformations. It seems that slowly varying functions are more likely to avoid long chains. A kind of compromise must be found between the complete de-correlation of the baths and the efficiency of the chains.

We will illustrate that, at least in favorable cases, the factorization with \(\boldsymbol{A}_D(\omega)\) giving a single shared bath is possible. In HEOM, the only way to avoid the full correlated simulation is the choice of an optimal frequency independent transformation \(\boldsymbol{U}_{opt}\). T-TEDOPA could in principle treat the full correlated case by using the factorization with \(\boldsymbol{U}(\omega)\sqrt{\mathbf{\Lambda}(\omega)}\). A simplification may also be used by choosing an optimal frequency independent \(\boldsymbol{U}_{opt}\) matrix. This approximation has already been shown in ref.[63].

7 Correlation in a LVC model of a conical intersection↩︎

We consider a molecular model with a ground and two excited states coupled by a conical intersection, the ground state being decoupled from the excited manifold. In a diabatic representation, the generic LVC model Hamiltonian centered at the Franck-Condon Point (FCP) geometry reads :

\[\begin{align} & H=\left( \begin{matrix} {{\varepsilon}_{0}} & 0 & 0 \\ 0 & {{\varepsilon}_{1}} & 0 \\ 0 & 0 & {{\varepsilon}_{2}} \\ \end{matrix} \right)+\left( \begin{matrix} {{h}_{0}} & 0 & 0 \\ 0 & {{h}_{1}} & 0 \\ 0 & 0 & {{h}_{2}} \\ \end{matrix} \right) \nonumber \\ & +\left( \begin{matrix} 0 & 0 & 0 \\ 0 & \sum\nolimits_{j}{\kappa _{j}^{(1)}{{q}_{j}}} & \sum\nolimits_{j}{{{\kappa}'_{j}}{{q}_{j}}} \\ 0 & \sum\nolimits_{j}{{{\kappa}'_{j}}{{q}_{j}}} & \sum\nolimits_{j}{\kappa _{j}^{(2)}{{q}_{j}}} \\ \end{matrix} \right) \label{eq:LVCmodel} \end{align}\tag{25}\] where \({{\varepsilon}_{n}}\) (\(n=0,1,2\)) are the diabatic potential energies at the FCP for each normal mode. The vibrational Hamiltonian of electronic state is supposed to be harmonic. By neglecting any Duschinsky rotation, the normal modes are assumed to be the same in each electronic state. Only their equilibrium position \(q_{j,eq}^{(n)}\) is different. In this harmonic approximation, the excited electronic energies at FCP \({{\varepsilon}_{n}}={{\varepsilon }_{n}^0}+{{\lambda }_{n}}\) (\(n=1,2\)) are the sum of the energy at the minimum for each mode \(q_{j,eq}^{(n)}\) and the renormalization energy \({{\lambda }_{n}}=\frac{1}{2}\sum\nolimits_{j}{\omega _{j}^{2}}q_{j,eq}^{(n)2}\). \(\kappa _{j}^{(n)}=\omega _{j}^{2}\left( q_{j,eq}^{(0)}-q_{j,eq}^{(n)} \right)\) are the gradients of the diabatic potential energy curves at the FCP and the sum runs over all the tuning modes. \({{\kappa}'_{j}}\) are the gradients of the electronic interstate coupling. The diagonal and off-diagonal sums define collective modes having the dimension of an energy: the collective tuning modes \(B_{n}^{{}}=\sum\nolimits_{j}^{{}}{\kappa _{j}^{(n)}q_{j}^{{}}}\) (\(n=1,2\)) that are coupled to the system by \(S_n=\ket{n} \bra{n}\) and the interstate coupling collective mode \({{B}_{3}}=\sum\nolimits_{j}^{{}}{{{\kappa}'_{j}}}q_{j}^{{}}\) coupled by \(S_3=\ket{1}\bra{2} + \ket{2}\bra{1}\). The two baths \(B_1\) and \(B_2\) are correlated as illustrated in Sec.2 for the fully correlated case. We assume that the coupling bath \(B_3\) is completely uncorrelated. It is here a simple spectator bath for the correlation problem but it plays a significant role in the non-adiabatic dynamics.

8 Illustrative Application↩︎

We present an illustrative example based on the \(\textit{ab}\) \(\textit{initio}\) LVC model of the 1,3-bis (phenylethynyl)benzene dimer where a conical intersection couples the excited states \(\text{S}_1\) and \(\text{S}_2\) (see fig.5). The corresponding potential energy surfaces and inter-state coupling have been computed using DFT (ground state) and TD-DFT (excited states) at the CAM-B3LYP/6-31\(+\)G\(*\) level of theory calculated by B. Lasorne and coworkers [22], [74]. The minimum-energy conical intersection (MECI) of the seam was characterized with the branching-space vectors associated to the halved gradient difference vector (GD) and the derivative coupling vector (DC). This is referred to as the m22 dimer in the following. All LVC data (gradients along the normal modes at the FC geometry) have already been used in previous simulations [22], [62] and are given in the supplementary material of these works.

This is an example where the tuning baths obtained from Eq.(25 ) are strongly correlated and for which the approximate treatment of a shared bath has been used in previous HEOM applications [62], [63] without really justifying it. The LVC Hamiltonian is built by retaining only the planar geometry. At the FCP, the two \(\text{S}_1\) and \(\text{S}_2\) states are very close, at 4.43 and 4.47 eV. The tuning and coupling modes have the A\(_1\) or B\(_2\) symmetry of the C\(_{2v}\) point group, respectively. The two states may be excited via orthogonal transition dipoles, \(\mu_y = 3.96\) \(ea_0\) for \(\text{S}_1\) and \(\mu_z = -1.87\) \(ea_0\) for \(\text{S}_2\) (the axes are represented in fig.3).

The continuous spectral densities of the tuning (\(J_{11}\), \(J_{22}\)) and coupling (\(J_{33}\)) baths and the crossed spectral density (\(J_{12}\)) are given in fig.6. In the absence of information provided by molecular dynamics in a solvent, we use a smoothing strategy [37], [38] by broadening each discrete data by a Lorentzian function \[\delta (\omega -{{\omega }_{j}})\sim\frac{1}{\pi }\frac{\Gamma}{{{\left( \omega -{{\omega }_{j}} \right)}^{2}}+{{\Gamma }^{2}}}. \label{eq:delta}\tag{26}\] The phenomenological width \(\Gamma={80}\) cm\(^{-1}\) is chosen to match the spectral enlargement introduced in a previous MCTDH simulation of linear spectra [22]. For the HEOM simulations, the spectral densities are fitted here by Tannor-Meier Lorentzian functions[63], [68], [75], [76] that provide an analytical expression of the parameters in the \(C_{lm}(t)\) expansion (Eqs.(17 ) and (18 )). The main relations are summarized in the Supplemental Material. The dominant peaks arise at the high frequency symmetrical or antisymmetric stretching modes of the acetylenic bond. At this frequency, the gradients have the same sign, which explains the large positive value of the crossed term. After smoothing, the crossed spectral density \(J_{12}(\omega)\) has no negative value in this case. Note that the convolution procedure effectively provides negative values when most gradients have different signs. The positively correlated tuning baths are related to what is generally referred to as a sloped conical intersection.

We will compare the two exact simulations at T \(=\) 298 K driven by the correlated HEOM equations (Eq.(19 )) and by the de-correlated T-TEDOPA method using the factorization by the \(\boldsymbol{A}_D(\omega)\) matrix (Eq.(23 )). For HEOM, the main constraints is the very fast increase of the number of ADOs with the number of decay modes of all the correlation functions. It is given by \({{N}_{\text{ADO}}}=\frac{\left( L+K \right)!}{L!K!}\). Five correlation functions are involved in this application, four for the two correlated tuning bath and one for the tuning bath, \(K=\sum\nolimits_{l,m=1}^{2}{{{K}^{(lm)}}}+{{K}^{(33)}}\) (see fig.3). Each spectral density is fitted by three Lorentzian functions, each of which leads to two decay modes in the correlation function. The contribution of the Matsubara terms coming from the Bose function may be fitted by two exponential functions for each correlation function, as done in refs. [63], [77] (the parameters of this fit at 298 K are given in the SI). These terms are correctly taken into account in the approximate treatment using \(\mathbf{U}_{opt}\) but not in the correlated case (Eqs.(19 )). In our simulation of the correlated case, the number of artificial modes is thus \(K^{(lm)}=3 \times 2= 6\) (number of Lorentzian functions per bath times number of decay modes per Lorentzian). At level \(L = 7\) of the hierarchy, the number of ADOs generated by the total number of artificial modes \(K=5 \times 6 =30\) (number of baths times number of modes per bath) is already \(N_{\text{ADO}} = 10,295,472\). Note that the approximate method with a single tuning shared bath requires significantly less computational effort. The number of correlation functions is then only 2, one for the shared bath and one for the tuning bath. With 6 decay modes per bath, a simulation at the same level \(L = 7\) demands only \(N_{\text{ADO}} = 50,388\) This reduction of the computational resource allows us to take into account the Matsubara terms by the fitting procedure. With two supplementary functions per bath, the number of modes is now 8 and at level \(L = 7\), which involves \(N_{\text{ADO}} = 116,280\). Note that a filtering procedure could reduce the ADO numbers during the dynamics [62], [65]. T-TEDOPA dynamics is done with a rank \(r=8\), a Fock space \(d=10\), 100 modes for the coupling bath and 90 modes for the tuning bath(s). Dynamics with the \(A_D\) procedure involves a single shared tuning bath. The storage then requires \(N_{TEDOPA}= (100+90) \times 10 \times 8^2= 121,600\) (total number of modes times the dimension of the Fock space times the square of the rank) complex elements. In the \(A\) approach, there are two tuning baths and this involves \(N_{TEDOPA}= (100+180) \times 10 \times 8^2= 179,200\) elements. In terms of storage resources, the T-TEDOPA method is clearly more effective than HEOM (discussed here in the standard implementation that could be notably improved by the tensor train method). A comparison among different implementations was already made in ref. [63]. The \(A_D\) strategy is also more advantageous, at least when the density matrix has only one non-zero eigenvalue. Even if the comparison of computational time strongly depends on the hardware used, in our simulations, T-TEDOPA was always more efficient in computational time.

The demanding resources for solving the correlated HEOM at high level stimulate the use of the approximate shared bath when this is possible. Obviously, the ratio between the two spectral densities \(J_{22}(\omega)/J_{11}(\omega)\) depends on the frequency and one has to choose the best frequency independent transformation \(\boldsymbol{U}_{opt}\) for the approximate de-correlation. The eigenvector components \(U_{n1}(\omega)\) (\(n = 1,2\)) of \(\boldsymbol{J}^{\beta}(\omega)\) associated with the dominant eigenvalue \(\Lambda_1(\omega)\) are shown in the upper part of fig.7 (\(\Lambda_2(\omega)\) is negligible). The lower panel of fig.7 shows the elements of the \(\boldsymbol{A}_D(\omega)\) (Eq.(23 )) and \(\boldsymbol{A}(\omega)\) (Eq.(22 )) matrices used in the T-TEDOPA applications.

Figure 8 gives the dominant eigenvalue \(\Lambda_1(\omega)\) of \(\boldsymbol{J}^{\beta}(\omega)\). The second eigenvalue \(\Lambda_2(\omega)\) is negligible. We compare this dominant eigenvalue with the approximate values \(\bar{\Lambda}_1(\omega)\) obtained by the different strategies searching for a constant transformation matrix \(\boldsymbol{U}_{opt}\) presented in sec.5. In the strategy based on the the projection of the principal vector of the PCA, we use a tensor of 1000 vectors \((J^{\beta}_{11}(\omega_k),J^{\beta}_{22}(\omega_k),2J^{\beta}_{12}(\omega_k))\) with a regular sampling and the main PCA vector maximizes the variance. The errors (Eq.(21 )),

which should vanish with the exact \(\boldsymbol{U}(\omega)\) matrix, are very small \(5.07 \times 10^{-5}\) for \(\boldsymbol{U}_\text{aver}\) (strategy (b)), \(3.28 \times 10^{-5}\) for \(\boldsymbol{U}_{C(0)}\) (strategy (c)) and \(4.51 \times 10^{-5}\) for \(\boldsymbol{U}_\text{PCA}\) (strategy (d)).

In \(\boldsymbol{U}_{\eta =1}\) the coefficients \(U_{11}\) and \(U_{21}\) are equal to 0.707106 while in \(\boldsymbol{U}_{C(0)}\), \(U_{11}=0.70707\) and \(U_{21} = 0.70714\), which are roughly the same. This approximation with \(\boldsymbol{U}_{\eta =1}\) is an acceptable approximation in this example. It corresponds to the well known transformation of the modes allowing to describe the conical intersection in terms of three effective modes. The first with linear couplings \(\kappa _{j}^{+}=(\kappa _{j}^{(1)}+\kappa _{j}^{(2)})/2\) is associated with the variation of the average energy. The second one is the tuning mode related to the energy gap between the two excited states with couplings \(\kappa _{j}^{-}=(\kappa _{j}^{(1)}-\kappa _{j}^{(2)})/2\) and the third one is the inter-state coupling mode [17], [78]–[82].

We now present the dynamics driven by the correlated baths in HEOM and by T-TEDOPA with the shared bath resulting from the \(\boldsymbol{A}_D(\omega)\) factorization. They are compared with some approximate methods. The population \(P_2\) in the upper excited state is plotted in fig.9 for an initial superposed state with equal weights in \(\text{S}_1\) and \(\text{S}_2\) at \(T=298~\)K (\(P_1=1-P_2)\). The agreement between the correlated HEOM and the new T-TEDOPA procedure with \(\boldsymbol{A}_D(\omega)\) is very satisfactory given the absence of the Matsubara terms in HEOM. We first compare with the approximate fully correlated case, this means the shared bath generated by the \(\boldsymbol{U}_{\eta=1}\). The early dynamics is similar. The difference from the exact result is more significant with respect to the relaxed values (the effect is still larger with \(\boldsymbol{U}_{\eta=0.8}\) given in the supplementary material). This difference comes from a modification of the renormalization energy \(\lambda_2\) when modifying the effective spectral density \(J_{22}\). Indeed, in the harmonic approximation, the energy at FCP is related to the minimum energy by \(\epsilon_m = \epsilon_m^0 + \lambda_m\). The variation of the difference \(\lambda_1-\lambda_2\) involves a modification of \(\epsilon_1^0-\epsilon_2^0\) and therefore of the equilibrated populations. Finally, we consider the T-TEDOPA simulation with the \(\boldsymbol{A}(\omega)\) factorization. As expected, the \(\boldsymbol{A}_D(\omega)\) method based on the exact diagonalization of \(\boldsymbol{J}^{\beta}(\omega)\) provides results closer to the HEOM results. The \(\boldsymbol{A}_D(\omega)\) strategy gives a behavior similar to the fully correlated case.

The decay of the coherence of a superposed initial state with equal weights is shown in fig.10. We present the same simulations as in fig.9 and also the uncorrelated case that assumes two independent baths. The decoherence decay is much slower when taking into account the bath correlation. In a non-adiabatic wave packet approach, the electronic coherence is linked to the overlap of the nuclear wave packets evolving in each excited state. In the case of a sloped conical, the gradients have the same sign, the wave packet components evolve in the same direction, which promotes long-term overlap. This effect is well described by the HEOM simulation with correlated baths or by T-TEDOPA with the \(\boldsymbol{A}_D(\omega)\) formalism. The approximate fully correlated case with \(\boldsymbol{U}_{\eta = 1}\) is also acceptable in this example at least during the first fifty femtoseconds. As for the populations, the MPS method with the \(\boldsymbol{A}(\omega)\) provides results similar to the fully correlated case. On the contrary, neglecting the correlation is a very bad approximation predicting a too fast decay.

9 Concluding remarks↩︎

In this work, we revisit the crossed correlation of tuning baths in a LVC model to tackle the problem of continuous spectral densities with HEOM and T-TEDOPA simulations. Indeed, in a dense vibrational environment, dynamics of the excited electronic states is driven by correlated fluctuations, both in the electrons-nuclei partition [46] and in the main vibrational coordinate mapping [48].

The correlated baths have been addressed early in the HEOM community [55]–[57]. In the applications on excitation energy transfer in photosynthetic aggregates, the correlation concerns baths on different sites. It is generally treated by assuming a simplified expression for the cross-correlation function [12], [31], [42]. Here, we consider a local LVC model with a conical intersection for which an approximate correlation matrix was used in our previous works [24], [62], [63].

For the HEOM simulations, we compare the exact correlated propagation with an approximate method that de-correlates the baths in favorable cases by using an optimal frequency independent transformation of the tuning bath modes. The latter is chosen from an average of the eigenvectors of the \(\mathbf{J}^{\beta}(\omega)\) matrix or from the eigenvectors of the initial bath correlation matrix \(\mathbf{C}(t=0)\). When the dyadic structure of the discrete \(\mathbf{J}^{\beta}(\omega)\) matrix is preserved during the smoothing procedure, a single eigenvector dominates that is associated with a single shared bath. This approximation is then very interesting since it requires far fewer computational resources than the correlated HEOM when the hierarchy level is high.

Methods using the discretization of the spectral densities for MCTDH applications also face the correlation problem [38]. The difference with T-TEDOPA mainly lies in the sampling of the temperature dependent spectral density and in the technical way of building the chain of modes. Here, we have shown that the T-TEDOPA method allows in principle to use the exact frequency dependent transformation of the modes diagonalizing the \(\mathbf{J}^{\beta}(\omega)\) matrix. However, HEOM and T-TEDOPA have very different numerical constraints. T-TEDOPA can manage the de-correlating frequency dependent transformation but the constraint lies in the length of the chains and in the long range interaction. The transformed system-bath Hamiltonian may be obtained by different factorization matrices using the uncorrelated spectral modes as suggested here or a combination of these modes as done in a recent application where the spectral densities are extracted from molecular dynamics [28], [29], [64]. When the factorization matrix is not the spectral one, some dynamical correlation remains but the numerical price of the chains may be reduced. In favorable cases, T-TEDOPA using the exact spectral factorization is a very efficient alternative to the exact correlated HEOM simulation.

Supplementary Material↩︎

The Supplementary Material gathers the parameters of the Tannor-Meier Lorentzian functions fitting the different spectral densities. We summarize the main characteristics of the bath correlation functions, specifically the case of the expansion of the complex conjugate function (Eq.(17 )). We also give a supplementary example of the dynamics by different approaches.

10 Acknowledgments↩︎

We warmly acknowledge Dr B. Lasorne and Dr J. Galiana for the \(\textit{ab}\) \(\textit{initio}\) data used to calibrate the two-bath model. We also thank Dr A. Jaouadi for helpful discussions. B.LD acknowledges the support of iSiM Sorbonne. EM acknowledges support from ANR project NQESim (Grant No. ANR-23-CE29-0024-01). AWC acknowledges support from the ANR project RADPOLIMER ( ANR-22-CE30-0033 ).

Data availability↩︎

The data that support this study are available on request from the authors. All the T-TEDOPA simulations were done with the MPSDynamics Julia package [83] that can be found along with the documentation at https://github.com/shareloqs/MPSDynamics.

References↩︎