E Mangaud

B Lasorne

) O Atabek

M Desouter-Lecomte

Statistical distributions of the tuning and coupling collective modes at a conical intersection using the hierarchical equations of motion

)) whereby the information contained in the auxiliary HEOM matrices is exploited in order to derive the nuclear dissipative wave packet, i.e. the statistical distribution of the displacement of the two tuning and coupling collective coordinates in each electronic state and the coherence. This allows us to visualize the distribution, all along the non-adiabatic decay. We explore a large parameter space for a symmetrical conical intersection model and a symmetrical initial Franck-Condon preparation. Some parameters could be controlled by external fields while others are molecule dependent and could be designed by molecular engineering. We illustrate the relation between the strongly coupled electronic and bath dynamics together with a geometric measure of non-Markovianity.

I. INTRODUCTION

A conical intersection (CI) of two electronic potential energy surfaces is one of the key regions where the Born-Oppenheimer approximation strongly breaks down [START_REF] Herzberg | Intersection of potential energy surfaces in polyatomic molecules[END_REF][START_REF] Michl | Electronic Aspects Of Organic Photochemistry[END_REF][START_REF] Klessinger | Excited States and Photochemistry of Organic Molecules[END_REF][START_REF] Yarkony | Diabolical conical intersections[END_REF][START_REF] Domcke | Conical Intersection, Electronic Structure, Dynamics and Spectroscopy[END_REF][START_REF] Worth | Beyond Born-Oppenheimer: Molecular dynamics through a conical intersection[END_REF][START_REF] Bersuker | The Jahn-Teller Effect[END_REF][START_REF] Baer | Beyond Born-Oppenheimer: Electronic Nonadiabatic Couling Terms and Conical Intersections[END_REF][START_REF] Matsika | Nonadiabatic events and conical intersections[END_REF][START_REF] Domcke | Conical Intersection: Theory, Computation and Experiment[END_REF] . It occurs in N-atom systems with N ≥ 3. According to the noncrossing rule [START_REF] Neuman | Uber merkwürdige diskrete Eigenwerte. Uber das Verhalten von Eigenwerten bei adiabatischen Prozessen[END_REF] , two adiabatic electronic states may be energydegenerate along a (3N -8) dimensional seam when they have the same spin multiplicity or a (3N -7) one when they are of different spin multiplicity. In the former case, the other two degrees of freedom form a subspace, called branching space, in which the adiabatic potential energy surfaces appear as two cones touching at the CI (degeneracy is lifted to first order over the plane tangent to the branching space). This "funnel" is likely to induce ultrafast decay since the non-adiabatic coupling becomes infinite at the CI. The latter is thus at the core of photochemical processes in excited states and is giving rise to a wide research area, both in advanced time-dependent spectroscopy [START_REF] Kowalewski | Simulating coherent multidimensional spectroscopy of nonadiabatic molecular processes: From the infrared to the x-ray regime[END_REF][START_REF] Schuurman | Dynamics at conical intersections[END_REF][START_REF] Wu | Two-dimensional electronicvibrational spectroscopic study of conical intersection dynamics: An experimental and electronic structure study[END_REF][START_REF] Palacino-González | Analysis of transient-absorption pump-probe signals of nonadiabatic dissipative systems: "Ideal" and "real" spectra[END_REF] and in numerical simulations [START_REF] Garavelli | Potential-energy surfaces for ultrafast photochemistry Static and dynamic aspects[END_REF][START_REF] Robb | A computational strategy for organic photochemistry[END_REF][START_REF] Olivucci | Computational Potochemistry[END_REF][START_REF] Worth | Using the MCTDH wavepacket propagation method to describe multimode non-adiabatic dynamics[END_REF][START_REF] Lasorne | Excited-state dynamics[END_REF][START_REF] Robb | Chapter three -in this molecule there must be a conical intersection[END_REF][START_REF] Gatti | Applications of Quantum Dynamics in Chemistry[END_REF][START_REF] Guan | Neural network based quasi-diabatic Hamiltonians with symmetry adaptation and a correct description of conical intersections[END_REF] . Alternatively, diabatic representations of the electronic states proa) Electronic mail: etienne.mangaud@sorbonne-universite.fr b) Electronic mail: benjamin.lasorne@umontpellier.fr c) Electronic mail: osman.atabek@u-psud.fr d) Electronic mail: michele.desouter-lecomte@u-psud.fr vide a smoother picture of the CI in the branching space [START_REF] Desouter-Lecomte | Nonadiabatic interactions in unimolecular decay. V. Conical and Jahn-Teller intersections[END_REF][START_REF] Van Voorhis | The diabatic picture of electron transfer, reaction barriers, and molecular dynamics[END_REF] . Even if a rigorous diabatisation raises issues when restricting the basis to a model made of two electronic states, as abundantly discussed in the literature [START_REF] Baer | Beyond Born-Oppenheimer: Electronic Nonadiabatic Couling Terms and Conical Intersections[END_REF][START_REF] Mead | Conditions for the definition of a strictly diabatic electronic basis for molecular systems[END_REF] , a fair approximation may be derived in practice so that the diagonalization of the diabatic model yields the two adiabatic energies of interest [START_REF] Guan | Neural network based quasi-diabatic Hamiltonians with symmetry adaptation and a correct description of conical intersections[END_REF][START_REF] Thiel | Proposal and numerical test of a simple diabatization scheme[END_REF] . Within a harmonic approximation, the diabatic potential energies are described by two paraboloids crossing along a line (3N -7 dimensional hyperline). The diabatic electronic coupling is often considered as varying linearly in the vicinity of conical intersections, where it is zero by construction. In the branching space, the degeneracy is lifted to first order, except at the CI point when diagonalizing the diabatic matrix provides adiabatic potential energies. The two main degrees of freedom are the tuning mode (TM), Q t , along which the diabatic energy gap varies to first order, and the coupling mode (CM), Q c , that induces the linear electronic coupling [START_REF] Davidson | Global topology of triatomic potential surfaces[END_REF][START_REF] Desouter-Lecomte | Nonadiabatic unimolecular reactions of polyatomic molecules[END_REF][START_REF] Atchity | Potential energy surfaces near intersections[END_REF] .

Simulating excited-state non-adiabatic dynamics in extended systems remains a challenge and an active topic. This requires a compromise between the computational scaling and the need to explore accurately, and in real time, many environmental degrees of freedom at any temperature since dissipation is responsible for electronic relaxation. Numerical methods may be classified according to the explicit or implicit description of the environmental degrees of freedom. The first strategy involves an explicit description of the elec- tronic system and of a large sample of vibrators. It has been addressed by various strategies ranging from fully quantum methods to semi-classical approaches, for instance multiconfiguration time-dependent Hartree (MCTDH) [START_REF] Worth | Using the MCTDH wavepacket propagation method to describe multimode non-adiabatic dynamics[END_REF][START_REF] Gonon | A generalized vibronic-coupling Hamiltonian for molecules without symmetry: Application to the photoisomerization of benzopyran[END_REF] , multilayer-MCTDH [START_REF] Wang | On regularizing the ML-MCTDH equations of motion[END_REF] , variational approaches with Gaussian basis sets (ab initio multiple spawning AIMS [START_REF] Ben-Nun | Ab initio quantum molecular dynamics[END_REF][START_REF] Levine | Isomerization through conical intersections[END_REF] , direct-dynamics variational multiconfiguration Gaussian DD-vmCG) [START_REF] Richings | Quantum dynamics simulations using Gaussian wavepackets: The vMCG method[END_REF] , tensorial matrix product approaches [START_REF] Schroder | Tensor network simulation of multi-environmental open quantum dynamics via machine learning and entanglement renormalisation[END_REF] or mixed quantum-classical trajectories [START_REF] Schaupp | A classical ride through a conical intersection[END_REF] , semiclassical surface hopping [START_REF] Tully | Mixed quantum-classical dynamics[END_REF][START_REF] Xie | Assessing the performance of trajectory surface hopping methods: Ultrafast internal conversion in pyrazine[END_REF] . Full quantum methods based on wave functions suffer from the scaling and from the difficulty of simulating high temperature. On the other hand, semi-classical strategies may lead to a questionable description of the quantum noise or of tunneling effect. An alternative description is dissipative dynamics [START_REF] Vega | Dynamics of non-Markovian open quantum systems[END_REF] , which treats only implicitly the environment at any temperature with statistical Schrödinger equations [START_REF] Ke | Hierarchy of stochastic Schrödinger equation towards the calculation of absorption and circular dichroism spectra[END_REF][START_REF] Wang | The hierarchical and perturbative forms of stochastic Schrödinger equations and their applications to carrier dynamics in organic materials[END_REF][START_REF] Nakamura | Hierarchical Schrödinger equations of motion for open quantum dynamics[END_REF] , multi-state quantum Fokker-Planck equation (MSQFPE) [START_REF] Ikeda | Probing photoisomerization processes by means of multi-dimensional electronic spectroscopy: The multi-state quantum hierarchical fokker-planck equation approach[END_REF][START_REF] Ikeda | Phase-space wavepacket dynamics of internal conversion via conical intersection: Multi-state quantum Fokker-Planck equation approach[END_REF][START_REF] Ikeda | Low-temperature quantum fokker-planck and smoluchowski equations and their extension to multistate systems[END_REF][START_REF] Ikeda | Modeling and analyzing a photo-driven molecular motor system: Ratchet dynamics and non-linear optical spectra[END_REF] , or master equations in Liouville space. Observables related to the selected active subspace can thus be obtained exactly in the framework of the harmonic baths with linear couplings by non-perturbative methods such as the "hierarchical equations of motion" (HEOM) formalism [START_REF] Tanimura | Time Evolution of a Quantum System in Contact with a Nearly Gaussian-Markoffian Noise Bath[END_REF][START_REF] Tanimura | Stochastic Liouville, Langevin, Fokker-Planck, and master equation approaches to quantum dissipative systems[END_REF] .

The partitioning between an active subsystem and a thermal bath determines not only the coupling strength and therefore the computational scaling and the tractable Markovian or non-Markovian master equations but also the accuracy of the information about some environmental modes. In a two-electronic state case, the two main strategies are the spin-boson or "star model" and the effective reaction coordinate model. In the first case (denoted two-level system TLS), the electronic states are coupled to the ensemble of environmental modes and in the second, one or more primary effective coordinates are retained in the active system and are coupled to a residual bath. In the CI situation, the TM and CM modes are included in the system denoted two-level system two-coordinate model (TLS-TC). The relative canonical transformation of the vibrational modes has been discussed for a long time [START_REF] Levine | Hamiltonian theory for vibrational dephasing rates of small molecules in liquids[END_REF][START_REF] Cederbaum | Short-time dynamics through conical intersections in macrosystems[END_REF][START_REF] Gindensperger | Short-time dynamics through conical intersections in macrosystems. I. Theory: Effectivemode formulation[END_REF][START_REF] Hughes | Effective-mode representation of non-Markovian dynamics: A hierarchical approximation of the spectral density. II. Application to environment-induced nonadiabatic dynamics[END_REF][START_REF] Pereverzev | Energy and chargetransfer dynamics using projected modes[END_REF][START_REF] Martinazzo | Unraveling a Brownian particle's memory with effective mode chains[END_REF][START_REF] Roden | Accounting for intra-molecular vibrational modes in open quantum system description of molecular systems[END_REF] and in particular the relation between the corresponding spectral densities has been given in a seminal paper about electron transfer [START_REF] Garg | Effect of friction on electron transfer in biomolecules[END_REF] . This has given rise to an abundant literature using one or the other partitioning with different dynamical treatments [START_REF] Tanimura | Optical stark spectroscopy of a brownian oscillator in intense fields[END_REF][START_REF] Chernyak | Collective coordinates for nuclear spectral densities in energy transfer and femtosecond spectroscopy of molecular aggregates[END_REF][START_REF] Cao | A phase-space study of Bloch-Redfield theory[END_REF][START_REF] Cao | A unified framework for quantum activated rate processes. II. The nonadiabatic limit[END_REF][START_REF] Thoss | Self-consistent hybrid approach for complex systems: Application to the spin-boson model with Debye spectral density[END_REF][START_REF] Tanaka | Quantum dissipative dynamics of electron transfer reaction system: Nonperturbative hierarchy equations approach[END_REF][START_REF] Tanaka | Multistate electron transfer dynamics in the condensed phase: Exact calculations from the reduced hierarchy equations of motion approach[END_REF][START_REF] Chenel | Exciton dissociation at donor-acceptor heterojunctions: Dynamics using the collective effective mode representation of the spin-boson model[END_REF][START_REF] Iles-Smith | Environmental dynamics, correlations, and the emergence of noncanonical equilibrium states in open quantum systems[END_REF][START_REF] Iles-Smith | Energy transfer in structured and unstructured environments: Master equations beyond the Born-Markov approximations[END_REF][START_REF] Wertnik | Optimizing co-operative multi-environment dynamics in a dark-state-enhanced photosynthetic heat engine[END_REF] .

The particular CI case has been frequently dealt for instance with a perturbative Markovian Redfield approach [START_REF] Kühl | Multilevel Redfield description of the dissipative dynamics at conical intersections[END_REF] , time nonlocal second-order master equation [START_REF] Qi | Track-ing an electronic wave packet in the vicinity of a conical intersection[END_REF] , Wigner distribution [START_REF] Ikeda | Phase-space wavepacket dynamics of internal conversion via conical intersection: Multi-state quantum Fokker-Planck equation approach[END_REF] , hybrid quantum master equation [START_REF] Schile | Simulating conical intersection dynamics in the condensed phase with hybrid quantum master equations[END_REF] , and recently beyond perturbative regime by HEOM [START_REF] Chen | Dissipative dynamics at conical intersections: Simulations with the hierarchy equations of motion method[END_REF][START_REF] Chen | Mapping of wave packet dynamics at conical intersections by time-and frequency-resolved fluorescence spectroscopy: A computational study[END_REF] . One of the advantages is the accurate monitoring of the dissipative wave packet in the branching space that can be mapped by time and frequency resolved fluorescence spectroscopy [START_REF] Chen | Mapping of wave packet dynamics at conical intersections by time-and frequency-resolved fluorescence spectroscopy: A computational study[END_REF] . In particular, it enables to get information about the vibrational coherence that influences the dynamics through the CI. However, the dimension of the basis set for a two-dimension two-or many-electronicstate system may become prohibitive for HEOM simulations since all the so-called auxiliary matrices of the corresponding coupled equations have the same dimension as that of the system. An interesting alternative has been suggested recently [START_REF] Duan | Quantum mechanical wave packet dynamics at a conical intersection with strong vibrational dissipation[END_REF][START_REF] Dijkstra | Simulation of photo-excited adenine in water with a hierarchy of equations of motion approach[END_REF] . In this approach, the system describes only two (or more) electronic states (two-level system TLS) while the effective TM and CM are included in two highly structured baths. In most of the cases, this requires non-perturbative, non-Markovian treatments and, therefore, HEOM with a high level of hierarchy. The number of auxiliary matrices may be very large but they have a small dimension (only two by two in the TLS). This scaling is thus favorable mainly when several electronic states are considered [START_REF] Dijkstra | Simulation of photo-excited adenine in water with a hierarchy of equations of motion approach[END_REF] .

The main point addressed in the present work is to show that HEOM techniques enable us to access information about environmental processes even in the minimal TLS partition. We define tools able to extract the probability distribution of the effective modes of the structured spectral densities in real time from the auxiliary matrices. In this way, we can visualize the interplay between the electronic transition and the vibrational motion leading to a strong restructuring of the statistical distribution, which does not remain a passive unperturbed Gaussian distribution during the non-adiabatic evolution. We extend for two baths the interesting strategy proposed by Shi et al. [START_REF] Zhu | Explicit system-bath correlation calculated using the hierarchical equations of motion method[END_REF][START_REF] Liu | Reduced quantum dynamics with arbitrary bath spectral densities: Hierarchical equations of motion based on several different bath decomposition schemes[END_REF] so as to extract this statistical distribution of the collective TM and CM from the auxiliary matrices. The time evolution of the distribution, in particular the average position or the distortion of the initial Gaussian profile before reaching the asymptotic regime provide valuable information about the damped vibrational motion at any temperature.

We illustrate the general trends induced by typical changes made on the system or on the environmental spectral densities even if it is not possible to entirely explore the large parameter space. Some parameters related to the system could be controlled by interaction with laser fields such as the energy gap and the position of the CI by a Stark shift [START_REF] Choi | Selective photodissociation in diatomic molecules by dynamical Stark-shift control[END_REF][START_REF] Scheit | Control scheme of nonadiabatic transitions with the dynamical of potential curve crossing[END_REF][START_REF] Saab | Full-dimensional control of the radiationless decay in pyrazine using the dynamic Stark effect[END_REF][START_REF] Richter | Sub-laser-cycle control of coupled electron-nuclear dynamics at a conical intersection[END_REF] or the preparation of the initial wave packet by modifying the initial superposition of the electronic states during the ultra fast excitation. Some parameters are also molecule-dependent, and exploring the main behaviors is a timely issue in molecular engineering to design new generations of quantum-assisted molecular systems for energy transfer [START_REF] Levine | Conical intersections at the nanoscale: Molecular ideas for materials[END_REF] . Indeed it is possible to vary the parameters entering the model Hamiltonian (potential energies and couplings) according to the chemical nature of the chromophore within the molecular system. They will also be sensitive to the presence of covalent bridges and substituents (intramolecular interactions) and to the surrounding medium (intermolecular interactions with a cage, a protein, or a solvent). For instance, a crossing between a locally excited state and a charge-transfer state can be displaced both in energy and position upon playing with the stabilisation or destabilisation of the charge separation with different types of electroactive substituents or with solvents of various polarities. The relative positions, energies, and frequencies of the two minima will also be affected by the extent of the reorganization of the bonding pattern induced by the electronic excitation: a large Stokes shift between absorption and emission spectra will be the signature of a strong rearrangement. One can choose a system in which one electronic state has similar parameters with respect to the ground state, while the other is strongly shifted. This is likely to be the case when a state is very bright and the other is dark. Finally, modulating the intrastate coupling can be achieved upon changing the nature of the bridge between the donor and the acceptor in a system experiencing charge or excitation energy transfer between two sites. Through-bond and through-space processes should correspond to strong and weak couplings, respectively.

The paper is organized as follows. In Sec.IIA we define the main parameters the TLS and TLS-TC models and their relations. The spectral densities of both models and the correlation functions are given in Sec.IIB. The methodology to extract the statistical distribution of the baths is summarized in Sec.IIC while the HEOM equations for two uncorrelated baths are given in the Appendix. Sec.IID presents one possible non-Markovianity signature used in this work. Results of numerical simulations are presented in Sec.III where we illustrate the interplay between TM and CM collective modes and the dissipative electronic wave packet around a CI in a TLS model.

II. METHODS

A. System-environment partitioning

We consider three electronic states with a ground state |0 decoupled from a doublet of coupled excited diabatic states |1 and |2 . In our model, the ground state |0 and the first excited state |1 have the same equilibrium positions but the second one |2 is spatially displaced. By restricting the electronic basis set to two excited states, the full vibronic Hamiltonian in the diabatic representation reads

H = ∑ 2 n=1 |n H n n| + (|1 H 12 2| + h.c.) (1)
where h.c denotes Hermitian conjugation H n is the vibrational Hamiltonian associated with each electronic state |n and treated within a harmonic level of approximation. By assuming that the frequencies are the same in both states but the equilibrium positions are different, and by separating the N t modes (label t standing for tuning) that tune the energy gap from the N c modes (c) that induce a variation of the electronic coupling, H n reads in mass-weighted coordinates

H n = ε n + 1/2 ∑ N t k=1 p 2 t,k + ω 2 t,k q t,k -q (n) t0,k 2 +1/2 ∑ N c k=1 p 2 c,k + ω 2 c,k q c,k -q (n) c0,k 2 . (2)
The reference point that defines the equilibrium geometry of the initial vibrational bath for the partition is the equilibrium position of the ground state. It is chosen as q t,k = 0 and q c,k = 0. The equilibrium positions of the tuning modes are also equal to zero in the first electronic state, q

t0,k = 0, and those of the second state will be denoted q (2) t0,k = q 0k . The equilibrium positions of the coupling modes are the same in both states and equal to zero q (1)

c0,k = q (2)
c0,k = 0. The electronic coupling is varying linearly from the CI position,

H 12 = ∑ N c k=1 c k q c,k -q CI,k . (3)
We consider a symmetrical case for which the CI position corresponds to the equilibrium positions of the coupling modes, q CI,k = q c0,k = 0. We now summarize the operational relations of the two TLS or TLS-TC partitions for the CI model. They are illustrated in Fig. 2.

TLS partition

When partitioning so as to isolate the electronic subsystem, the full Hamiltonian adapted to the CI case is recast as H = H S + H Bt + H Bc + S t B t + S c B c . According to the chosen reference point the system Hamiltonian in the basis set |1 , |2 is :

H S = ∆ 2 σ z + 1 2 (I -σ z ) λ t +W 12 σ x (4)
where σ are the Pauli matrices, I the unity matrix, ∆ = ε 2ε 1 is the energy gap between the two minima of the diabatic potentials, W 12 is the electronic coupling at the reference point. When it is chosen at the equilibrium position of the coupling modes, q c0,k , one has W 12 = 0. λ t = ∑ k ω 2 t,k q 2 0k /2 is the renormalization energy (also called reorganization energy) of the tuning bath. With the chosen reference point, it is added only to the energy of state |2 . The two baths are

H Bt = 1 2 ∑ N t i=1 p 2 t,i + ω 2 t,i q 2 t,i and H Bc = 1 2 ∑ N c i=1 p 2 c,i + ω 2 c,i q 2 c,i . The system-bath interaction involves collective coordinates B t = ∑ N t k=1 d k q t,k with d k = q t0,k ω 2 t,k and B c = ∑ N c
k=1 c k q c,k (c k being defined through Eq.( 3)). B t is diagonally coupled to the electronic system through the operator S t = 1 2 (σ z -I) (since only state |2 is coupled to the bath due to the choice made for the reference point) while B c is off-diagonally coupled by the operator S c = σ x .

TLS-TC partition

In order to incorporate the largest part of the vibronic coupling into the enlarged system, a unitary transformation of the bath oscillators is made to separate in each case the effective coordinate that captures the whole coupling [START_REF] Cederbaum | Short-time dynamics through conical intersections in macrosystems[END_REF][START_REF] Gindensperger | Short-time dynamics through conical intersections in macrosystems. I. Theory: Effectivemode formulation[END_REF][START_REF] Hughes | Effective-mode representation of non-Markovian dynamics: A hierarchical approximation of the spectral density. II. Application to environment-induced nonadiabatic dynamics[END_REF][START_REF] Pereverzev | Energy and chargetransfer dynamics using projected modes[END_REF][START_REF] Martinazzo | Unraveling a Brownian particle's memory with effective mode chains[END_REF][START_REF] Roden | Accounting for intra-molecular vibrational modes in open quantum system description of molecular systems[END_REF][START_REF] Cao | A phase-space study of Bloch-Redfield theory[END_REF][START_REF] Cao | A unified framework for quantum activated rate processes. II. The nonadiabatic limit[END_REF][START_REF] Thoss | Self-consistent hybrid approach for complex systems: Application to the spin-boson model with Debye spectral density[END_REF][START_REF] Pollak | Theory of activated rate processes: A new derivation of Kramers' expression[END_REF] ,

q t,i ⇒ Q t , Q j j = 2, N t and q c,i ⇒ Q c , Q j j = 2, N c where Q t = 1 D t ∑ N t k=1 d k q t,k with D 2 t = ∑ N t j=1 d 2 j and Q c = 1 D c ∑ N c k=1 c k q c,k with D 2 c = N c ∑ j=1 c 2 j .
The unitary transformation is built so that the residual modes are coupled to the effective one but are not coupled with each other. The transformed Hamiltonian (Eq.1) is by accounting for the reference point

H = ∑ 2 n=1 |n Hn n| + |1 W12 2| + h.c. (5)
with

H1 = ε 1 + 1/2(P 2 t + Ω 2 t Q 2 t ) + 1/2(P 2 c + Ω 2 c Q 2 c ) + HB (6) H2 = ε 2 +λ 2 +1/2(P 2 t +Ω 2 t Q 2 t )-D t Q t +1/2(P 2 c +Ω 2 c Q 2 c )+ HB (7) HB = 1/2 ∑ N t k=2 P 2 k + Ω2 k (Q k + dk Q t / Ω2 k ) 2 (8) + 1/2 ∑ N c k=2 P 2 k + Ω 2 k (Q k + ck Q c / Ω 2 k ) 2 (9) H12 = D c Q c . (10)
The TLS-TC partitioning retains in the system Hamiltonian HS all the terms of H related to the two-dimensional subspace and the counter-terms coming from Q 2 t and Q 2 c in HB [START_REF] Iles-Smith | Environmental dynamics, correlations, and the emergence of noncanonical equilibrium states in open quantum systems[END_REF][START_REF] Duan | Quantum mechanical wave packet dynamics at a conical intersection with strong vibrational dissipation[END_REF][START_REF] Tanimura | Reduced hierarchy equations of motion approach with Drude plus Brownian spectral distribution: Probing electron transfer processes by means of two-dimensional correlation spectroscopy[END_REF] .

For the TM, the coupling constant D t determines the displacement D t /Ω 2 t between the equilibrium position of the two effective diabatic parabolas while the coupling D c fixes the slope of the linear electronic coupling in the diabatic representation. The residual tuning and coupling baths gather all the new oscillators Q j or Q j and the system-bath couplings now involve cross terms d j Q t Q j or c j Q c Q j , respectively. The effective frequencies resulting from the transformed Hessian matrix 54 are given by

Ω 2 t = 1 D 2 t ∑ N t i=1 d i 2 ω 2 t,i
and

Ω 2 c = 1 D 2 c ∑ N c i=1 c i 2 ω 2 c,i .

B. Spectral densities and correlation functions

In open quantum systems, the system-bath interaction is characterized by the spectral density, which has a generic expression:

J(ω) = π 2 ∑ k f 2 k ω -1 k δ (ω -ω k ) (11)
where f 2 k is the strength of the linear coupling at each frequency. In the TLS, the spectral densities of the tuning and coupling baths are J t/c (ω) with f k = d k and f k = c k , respectively, while the spectral densities of the residual baths in the TLS-TC are J res t/c (ω) with f k = dk and f k = ck . The spectral density has been introduced here in a discrete representation as it should be the case when the displacements come from an ab initio normal mode analysis. When the frequency density is high, one may switch to the continuous representation by broadening the delta distribution by a Lorentzian function δ (ωω k ) → 1 π Γ (ω-ω k ) 2 +Γ 2 as suggested for instance in refs [START_REF] Martinazzo | Effective spectral densities for system-environment dynamics at conical intersections: S2-S1 conical intersection in pyrazine[END_REF][START_REF] Tamura | Quantum dynamics of ultrafast charge transfer at an oligothiophene-fullerene heterojunction[END_REF] . Note that a continous spectral density is directly provided by experiments or by molecular dynamics simulations giving the correlation of the energy gap [START_REF] Valleau | On the alternatives for bath correlators and spectral densities from mixed quantum-classical simulations[END_REF][START_REF] Aghtar | Different types of vibrations interacting with electronic excitations in phycoerythrin 545 and Fenna-Matthews-Olson antenna systems[END_REF] . A discretization taking into account the frequency density is then used for simulations of MCTDH type [START_REF] Wang | Systematic convergence in the dynamical hybrid approach for complex systems: A numerically exact methodology[END_REF] .

The temperature-dependent spectral density

C(ω) = J (ω) e β ω -1 -1 with β = 1/k B T is the Fourier transform of the correlation function of the collective mode C(t) = T r B B(t)B(0)ρ eq B where ρ eq B = exp (-β H B ) /T r B [exp (-β H B )]
is the Boltzmann equilibrium density matrix of the bath and

B(t) = exp (iH B t) B exp (-iH B t) (with h = 1 a.u.) C(t) = 1/π +∞ -∞ dω C(ω)e iωt . ( 12
)
The relation between the spectral density J res (ω) of the one-effective mode two-level model with J(ω) of the twolevel model (spin-boson model) has been derived in ref. [START_REF] Garg | Effect of friction on electron transfer in biomolecules[END_REF] . When J res (ω) is Ohmic with a friction coefficient η given by the slope for ω = 0, the spin-boson spectral density is a Lorentzian function whose width is proportional to this friction. A small friction in one-effective mode model leads to a highly peaked J(ω) in TLS model.

Here, we assume Lorentzian functions for the spectral densities of the tuning and coupling baths given by the Tannor-Meier expression 90

J t/c (ω) = p t/c ω ω + ω t/c 2 + Γ 2 t/c ω -ω t/c 2 + Γ 2 t/c (13)
where ω t/c gives the peak position, Γ t/c is the width and p t/c is a parameter fixing the coupling strength. The corresponding residual spectral density of the TLS-TC model is computed by the following relation, 62,65-68,91

J res t/c (ω) = D 2 t/c J t/c (ω) J 2 t/c (ω) +W 2 0,t/c (ω) (14)
where

W 0,t/c (ω) = 1 π PV +∞ -∞ dω J t,c (ω ) ω -ω (15)
with PV the principal value distribution of the Hilbert transform of the TLS spectral density. By using Eq.( 14), the residual spectral density is an Ohmic function with a friction coefficient η t/c = 2Γ t/c . With this parametrization (Eq.( 13)), the correlation function (Eq.( 12)) takes the form of a sum of decreasing complex exponentials well adapted for applying the HEOM algorithm, 77

C t/c (t) = ∑ n cor,t/c k=1 α t/c,k e iγ t/c,k t . (16)
As shown in ref. [START_REF] Pomyalov | The importance of initial correlations in rate dynamics: A consistent non-Markovian master equation approach[END_REF] , the complex conjugate may be expressed with the same γ t/c,k but modified α t/c,k coefficients

C * t/c (t) = ∑ n cor,t/c k=1 αt/c,k e iγ t/c,k t . ( 17
)
n cor,t/c contains two terms for each two-pole Lorentzian function plus a number of Matsubara terms corresponding to the poles of the Bose function. Analytical expressions of α, α, and γ are given in refs. [START_REF] Chenel | Exciton dissociation at donor-acceptor heterojunctions: Dynamics using the collective effective mode representation of the spin-boson model[END_REF][START_REF] Pomyalov | The importance of initial correlations in rate dynamics: A consistent non-Markovian master equation approach[END_REF] .

In the continuous limit, the renormalization energy λ t/c , the effective mode frequencies Ω t/c , and the coupling constant

D t/c read λ t/c = 1/π ∞ 0 dωJ t/c (ω)ω -1 (18)
D 2 t/c = 2 π +∞ 0 dωJ t/c (ω)ω (19
)

Ω 2 t/c = 2 πD t/c 2 +∞ 0 dωJ t/c (ω)ω 3 . ( 20
)
The linear couplings of the TLS-TC model may be related to the renormalization energies by the relation [START_REF] Wertnik | Optimizing co-operative multi-environment dynamics in a dark-state-enhanced photosynthetic heat engine[END_REF][START_REF] Duan | Quantum mechanical wave packet dynamics at a conical intersection with strong vibrational dissipation[END_REF] given here in mass weighted coordinates

D t/c = 2λ t/c Ω t/c . ( 21
) C. Statistical distribution of B t/c
The hierarchical equations of motion are one of the reference dynamical methods for open quantum systems modelled with a bath of harmonic oscillators inducing stochastic modulation of the diabatic energy gap and of the electronic coupling. The equations giving the evolution of the reduced density matrix, ρ(t) = T r B [ρ tot (t)], were originally derived for a Drude-Lorentz spectral density in the high-temperature limit from the Kubo stochastic Liouville equation [START_REF] Kubo | Stochastic Liouville Equations[END_REF] and the Feynman-Vernon influence functional formalism [START_REF] Tanimura | Time Evolution of a Quantum System in Contact with a Nearly Gaussian-Markoffian Noise Bath[END_REF][START_REF] Tanimura | Stochastic Liouville, Langevin, Fokker-Planck, and master equation approaches to quantum dissipative systems[END_REF] . They have been extended for lower temperature [START_REF] Ishizaki | Quantum Dynamics of System Strongly Coupled to Low-Temperature Colored Noise Bath: Reduced Hierarchy Equations Approach[END_REF][START_REF] Shi | Efficient hierarchical Liouville space propagator to quantum dissipative dynamics[END_REF] and for fermion and boson grand-canonical bath ensembles [START_REF] Jin | Dynamics of quantum dissipation systems interacting with fermion and boson grand canonical bath ensembles: Hierarchical equations of motion approach[END_REF] . The efficiency of the algorithm relies on the expansion of the correlation function C(t) as a sum of n cor complex exponential functions associated to pseudo decay modes of decay. Different expansions have been proposed [START_REF] Liu | Reduced quantum dynamics with arbitrary bath spectral densities: Hierarchical equations of motion based on several different bath decomposition schemes[END_REF] among which the Tannor-Meier one leading to Eqs.16 and 17 [START_REF] Meier | Non-Markovian evolution of the density operator in the presence of strong laser fields[END_REF][START_REF] Pomyalov | The importance of initial correlations in rate dynamics: A consistent non-Markovian master equation approach[END_REF] , the Padé approximation [START_REF] Hu | Padé spectrum decompositions of quantum distribution functions and optimal hierarchical equations of motion construction for quantum open systems[END_REF] , and the Fano spectrum decomposition scheme [START_REF] Cui | Highly efficient and accurate sum-over-poles expansion of Fermi and Bose functions at near zero temperatures: Fano spectrum decomposition scheme[END_REF] . A direct decomposition of the bath correlation function based on Chebyshev polynomials and Bessel functions has also been proposed leading to the so called C-HEOM scheme [START_REF] Rahman | Chebyshev hierarchical equations of motion for systems with arbitrary spectral densities and temperatures[END_REF] .

The method has been applied to several processes, for instance to simulate excitation energy transfer in photosynthetic light harvesting complexes 100-104 , electron transfer [START_REF] Tanaka | Multistate electron transfer dynamics in the condensed phase: Exact calculations from the reduced hierarchy equations of motion approach[END_REF][START_REF] Tanimura | Reduced hierarchy equations of motion approach with Drude plus Brownian spectral distribution: Probing electron transfer processes by means of two-dimensional correlation spectroscopy[END_REF]105,106 , quantum transport 107 , heat transport 108 , charge separation dynamics 109 . We assume here an initial factorization between the system and the bath, ρ tot (0) = ρ(0)ρ eq B valid for an ultrafast electronic transition (vertical Franck-Condon transition) 100 . The equilibrium initial bath is then the Boltzmann distribution of vibrational states in the ground electronic state projected in the bright excited state, and this initial condition is valid when the vibrational relaxation is not faster than the electronic dynamics. Initial correlation could be taken into account with different strategies 105,110,111 . The operational equations adapted for two uncorrelated baths and to the chosen expansion of C(t) are summarized in the Appendix. Here, we focus on the way the bath information is extracted from the auxiliary matrices in the two-bath case (i.e. the CI TLS model) by generalizing the methodology originally derived by Shi et al. [START_REF] Zhu | Explicit system-bath correlation calculated using the hierarchical equations of motion method[END_REF][START_REF] Liu | Reduced quantum dynamics with arbitrary bath spectral densities: Hierarchical equations of motion based on several different bath decomposition schemes[END_REF] and recently illustrated in a single bath and three-state case [START_REF] Chin | Visualising the role of non-perturbative environment dynamics in the dissipative generation of coherent electronic motion[END_REF] .

The auxiliary operators are matrices of the same dimension as the system density matrix. They are denoted by a collective index n = {n t , n c } giving the excitation number in the pseudo modes of each correlation function:

n j = n 1, j , • • • , n n cor , j
where j = t, c and n k is the quantum number of the k th pseudo mode. The system density matrix ρ(t) is given by the top row, i.e.

n = {0, • • • , 0} hence ρ(t) = ρ I (t) = ρ {0,••• ,0} (t).
The level of the hierarchy is equal to the sum of all the quantum numbers of the two baths. For each bath, one can also define the partial level to which belongs a given auxiliary matrix by summing the quantum numbers of only n t or n c . Extending the Shi et al methodology [START_REF] Zhu | Explicit system-bath correlation calculated using the hierarchical equations of motion method[END_REF] , one may obtain the high-order moments of the collective bath coordinate B t/c in each electronic state from an operator in the system subspace defined for the n th order by

X (n) t/c (t) = T r B t/c B n t/c ρ tot (t) ] (22)
In particular, each diagonal element X t/c,αα of X

(1)

t/c gives the expectation value of the collective coordinate B t/c in the corresponding electronic state α. This first moment can be obtained from the matrices contributing to the partial first level for each bath,

X (1) t/c (t) = -∑ n ρ n (t) (23)
where the sum runs over all the collective index n for which ∑ k n k,i = 1 and ∑ k n k, j = 0 with i = t/c and j = c/t respectively. Each diagonal element X t/c,αα may be renormalized to have the dimension of the effective coordinate Q t/c . These diagonal elements X j,αα /D j with j = t, c will be simply denoted X t/c /D t/c without notation of the electronic state when the latter is clearly specified. This quantity may be compared to the average position Q t/c (t) of the effective mode computed in the present work in an undamped two-dimensional effective space (i.e., with the system Hamiltonian HS . The first moment already provides interesting information about the TM vibrational relaxation and the evolution towards equilibrium. However, in the particular case of a CI with a symmetric preparation in the bright state, the average of the CM vanishes and the first moment is not sufficient to analyze the relaxation.

Deeper insight is provided by the full statistical distribution of the collective modes B t/c or of the collective coordinate x t/c = B t/c /D t/c in each electronic state. We first discuss the projection on the collective coordinates P α (x t/c ,t) where α denotes the electronic state. P α (x t/c ,t) may be compared to the projection of the two-dimensional undamped wave packet in each electronic state

P α Q t/c ,t = dQ c/t |Ψ α (Q t , Q c ,t)| 2 (24)
P α (x t/c ,t) is expanded in a Hermite polynomial basis set,

P α (x j ,t) = ∑ n a α, j,n (t)Φ n (x j ) (25)
with

Φ n (x j ) = D j 2πC j (0) √ 2 n n! H n D j x j / 2C j (0) e - D j x 2 2C j (0) (26
) where j = t/c, H n is the Hermite polynomial of order n and C j (0) is the initial value of the correlation function (Eq.( 16)), i.e. the integral of the temperature-dependent spectral density Cj (ω). The expression is adapted from refs. [START_REF] Zhu | Explicit system-bath correlation calculated using the hierarchical equations of motion method[END_REF][START_REF] Liu | Reduced quantum dynamics with arbitrary bath spectral densities: Hierarchical equations of motion based on several different bath decomposition schemes[END_REF] to ensure a correct normalization of the distribution and to map from B j to coordinate x j = B j /D j . The procedure may be understood from an early interpretation of the auxiliary matrices as time dependent coefficients of an expansion of the stochastic Liouville equation [START_REF] Kubo | Stochastic Liouville Equations[END_REF] in a basis set of the eigenvectors of the corresponding relaxation operator [START_REF] Cui | Highly efficient and accurate sum-over-poles expansion of Fermi and Bose functions at near zero temperatures: Fano spectrum decomposition scheme[END_REF]113 . The ratio D 2 j /C j (0)

in P(x j ) does not depend on the coupling strength. For each electronic state α the coefficient a α, j,n (t) is the diagonal element of a matrix built from the auxiliary matrices of partial level n for the bath j [START_REF] Thiel | Proposal and numerical test of a simple diabatization scheme[END_REF] where j = t, c. The sum runs over the index n for which the partial level of bath j is equal to n and that of the other bath is equal to zero. It is worth noting that a similar treatment has been proposed to extract the distribution of the moment conjugated to the collective mode in order to calculate thermal fluxes 108 .

A j,n (t) = (-1) n n!C n j (0) ∑ n n! ∏ k n k, j ! ρ n (t) ∑ k n k, j = n (
We now generalize the procedure to get the more informative two-dimensional distribution P α (x t , x c ,t) that may be compared to the square modulus of the undamped two-dimensional wave packet |Ψ α (Q t , Q c ,t)| 2 . The coherence χ α,β (x t , x c ,t) corresponds to the cross product of the wave packets in each electronic state

ξ αβ (Q t , Q c ,t) = Ψ * α (Q t , Q c ,t)Ψ β (Q t , Q c ,t).
We build an array of matrices of the same dimension as ρ(t)

M n,m (t) = (-1) n D t n!C n t (0) (-1) m D c m!C m c (0) ∑ n t n! ∏ k n k,t ! ∑ n c m! ∏ k n k,c ! ρ n t ,n c (t) ∑ k n k,t = n ∑ k n k,c = m. (28)
The two-dimensional distribution P α (x t , x c ,t) in the electronic state α is then given by

P α (x t , x c ,t) = ∑ n,m M αα n,m (t)Φ n (x t )Φ m (x c ) (29)
where M αα n,m is the corresponding diagonal element. Moreover, one can extract information about the electronic coherence as a function of the two coordinates from the off diagonal elements:

χ αβ (x t , x c ,t) = ∑ nm M αβ n,m (t)Φ n (x t )Φ m (x c ). (30)
The quantities extracted from the auxiliary matrices are here related to the diabatic description of the dynamics. The adiabatic electronic coherence during relaxation through a CI has been recently computed by the MSQFPE method [START_REF] Ikeda | Phase-space wavepacket dynamics of internal conversion via conical intersection: Multi-state quantum Fokker-Planck equation approach[END_REF] .

The distribution P α (x t , x c ,t) and the coherence χ α,β (x t , x c ,t) are computed here in the diabatic representation. Transforming to the adiabatic representation is possible by using the parameters of the corresponding two-dimensional TLS-TC model. At each point (x t , x c ) one may compute the real orthogonal diabatic-to-adiabatic transformation matrix by diagonalizing the local two-by-two diabatic electronic matrix. The adiabatic electronic states are then given by :

ψ ad k (x t , x c ) = ∑ k =1,2 S(x t , x c )ψ dia k (x t , x c ). (31)
As the distribution P α (x t , x c ,t) maps the product

ψ * dia α (x t , x c )ψ dia α (x t , x c ),
D. Non-Markovianity signature

Definition and measure of non-Markovianity has given rise to an abundant literature in the last decade 114,115 . The first crude estimation is based on the comparison between a characteristic timescale of the system dynamics and the bath correlation time. Assuming a delta-function for C(t) and therefore a flat spectral density (white noise) leads to memoryless dynamics which may be treated by Markovian Lindblad or Redfield approaches (at least in perturbative regime) 116 . In the TLS model of a CI, the spectral density is expected to be strongly structured, thus leading to long correlation times so that CI electronic dynamics described by this partition should be strongly non-Markovian. It is commonly said that non-Markovianity means back-flow from the bath to the system [START_REF] Vega | Dynamics of non-Markovian open quantum systems[END_REF] . In the electronic dynamics, this involves a boost of electronic average energy from the vibrational modes, in other words a non-monotonic electronic relaxation. We use here a non-Markovianity witness given by the volume of the accessible states in the Bloch sphere 117 . This geometric measure is not a very sensitive method to detect the non-Markovianity 118 mainly at long times but it is easy to compute and is independent of the initial conditions of the system.

For an N-level system, HEOM defines a non-Markovian dynamical map, ρ(t) = φ t [ρ(0)], that may be expressed in a basis set of N 2 operators in the Liouville space including the identity G 0 = I/ √ N and the N 2 -1 generators of SU(N), G i (i = 1, ..., N 2 -1) 119,120 . These generators are the Pauli spin matrices when N = 2. In matrix form the map reads F m,n (t) = Tr (G m φ t [G n ]) and the volume of accessible states is obtained from the determinant of the matrix of the map,

V (t) = det(F(t)). (33)
A non-monotonic decrease of this volume is a signature of non-Markovianity. The V (t) decay of a TLS results from two processes, the dephasing and the electronic population relaxation. In a pure dephasing case for which the electronic coupling vanishes, the populations remain constant and the volume is closely related to the decoherence function given by 116,121,122 ν(t) = exp(-Λ(t)) with

Λ(t) = ∞ 0 dωJ(ω) coth(β ω/2) (1 -cos(ωt)) /ω 2 . (34)
The bumps in the non-monotonic decay of ν(t) or V (t) occur with the same period as the correlation function C(t) and depends only on the spectral density. The volume decays without bumps when C(t) does not exhibit any oscillation, i.e. when the spectral density is very flat and broad. The population exchange when the electronic coupling is active, accelerates the decay of the volume and reduces or possibly flattens the bumps.

The volume evolution depends on the return of the collective mode back to its initial position leading to the bumps and on the non-adiabatic transition rate. Switching to the two-level-system one-coordinate representation, the volume behaves in a similar way as the autocorrelation function of the initial wave packet, which exhibits oscillations if the coordinate is not overdamped and has a decreasing norm due to the electronic transition. The decay of the initial wave packet norm depends on several factors, the time for reaching the crossing region and its usual Landau-Zener characteristics. The non-adiabatic dynamics is thus expected to be intrinsically non-Markovian when the collective mode is not overdamped (very thin spectral density) and the electronic decay is not too fast to avoid the blurring of the bumps.

In the CI case with a reference point where the electronic coupling vanishes, removing the off-diagonal coupling induced by the CM should lead to a pure dephasing case for the TM. A non-monotonic profile should be obtained as long as it is not overdamped. The population decay that flattens the bumps is due to the coupling bath. In the TLS-TC representation, it depends on the non-adiabatic evolution of a wave packet centered at the equilibrium position of the CM (since we consider here a symmetrical case).

III. RESULTS

We consider two diabatic excited states |1 and |2 and we assume that the system is initially prepared in the upper state |2 by an ultra fast excitation from a ground state |0 not included in the basis set. Modifying the parameters of the TLS model or of its spectral densities generates different shapes of the corresponding TLS-TC diabatic and adiabatic potential energy surfaces, in particular the CI position and the slopes at this point. This strongly influences the behavior of the electronic wave packet. The highly structured environment of the TLS is characterized by two peaked spectral densities J t/c (ω) (Eq.( 13)) with three main parameters (position, height, and width). Typical spectral densities used in the simulations are displayed in Fig. 3a and the corresponding real part of the correlation functions (Eq.( 12)) for T = 298K are given in Fig. 3b.

The positions of the spectral density peaks fixed by the ω k parameters (Eq.( 13) nearly correspond to the effective frequencies of the TM and CM in the TLS-TC model when the functions are sharply peaked. We choose ω t = 1000 cm -1 and ω c = 500 cm -1 (different frequencies are considered in the supporting information 123 ). The widths Γ t/c are proportional to the friction coefficient of the Ohmic spectral density of the residual bath in the TLS-TC model. In other words, this determines the damping of the TM or CM in this TLS-TC model. Two delta-function spectral densities ( with Γ t/c = 0) should lead to undamped TLS-TC dynamics, i.e., simple twodimensional dynamics. In most simulations, we take Γ t = Γ c = 60 cm -1 but we will also illustrate the effect of a variation of 12)) for T = 298K. The parameters of the spectral densities (Eq.( 13)) are p t = 1.95 × 10 -11 a.u., p c = p t /2.5 ; ω t = 4.5639 × 10 -3 a.u. (1000 cm -1 ), ω c = ω t /2; Γ t = Γ c = 2.7188 × 10 -4 a.u. (60 cm -1 ). λ t = 188 cm -1 and λ c = 300 cm -1 . The p t/c are taken as parameters to modify the renormalization energy λ t/c and therefore the linear couplings.

Γ c . The weight p t (Eq.( 13)) with constant Γ t determines the renormalization energy λ t , which is added to the energy difference ∆ to built H S of the TLS (Eq.( 4)). In the TLS-TC model, modifying p t alters the distance D t /Ω 2 t (Eq.( 21)) between the minima of the diabatic parabolas. Varying p c modifies λ c , which does not directly appear in the TLS system Hamiltonian but affects the linear coupling D c through the same relation (Eq.( 21)) and therefore the shape of the adiabatic potential energy curves along the CM. In summary, frequencies, spatial displacement and slope of the electronic coupling depend on the spectral densities. The relative position in energy is fixed by the difference ∆ between the minima corrected by λ t . Finally, the reference point defining the initial equilibrium position of the bath determines the Franck-Condon (FC) region. We assume here that the equilibrium geometry of the ground electronic state also corresponds to that of the lower dark excited state |1 while the bright state |2 is displaced. When the reference point is located at the equilibrium position of the CM, the corresponding electronic coupling W 12 of the TLS model vanishes. By varying its position one could also explore situations where the electronic coupling W 12 is not zero in the center of the FC region.

We will present below different examples focusing on parameters which could be controlled by external laser fields (energy gap and initial state). We also illustrate the influence of the damping of the coupling mode, i.e., the impact of the residual friction coefficient η c (linked to the width via η c = 2Γ c ) of the CM in the corresponding TLS-TC situation. The supporting information analyzes the variation of other molecule-dependent parameters (linear couplings and effective frequencies) 123 .

In the HEOM simulations, the correlation function of each bath is expanded on three exponential functions. Two sets of coefficients in Eqs.( 12) and ( 17) come from the poles of the chosen spectral density (Eq.( 13)). A single Matsubara term is enough to ensure convergence at room temperature. A high hierarchy level is necessary for analyzing the statistical distribution even if the populations converge at a lower level. Convergence of the expansions (Eqs.( 25) and ( 29)) may require about 20 or 25 terms and therefore the same hierarchy level. A typical computation at levels 20/25 involves 230230/736281 matrices respectively. The integration is carried out by the Cash-Karp Runge-Kutta algorithm with adaptative time step 124 . The wave packet is computed on a two-dimensional grid by the non-adiabatic split-operator algorithm 125 using only the system Hamiltonian of the TLS-TC model HS . For each dimension, the grid contains 2 10 points. The time step is 10 a.u.

A. Variation of the energy gap Fig. 4 presents cuts in the TLS-TC potential energy surfaces along the TM (left panel) and CM (right panel) at the reference position (Q t = 0 and Q c = 0) for two typical situations in which the adiabatic lower curve along the CM has either a single minimum or a barrier at the equilibrium position.

Environmental statistical distribution

Fig. 5 gives iso-value contours of the projected distribution P α (x t/c ,t) (Eq.( 25)) on the collective bath coordinate in each electronic state for the two contrasting cases illustrated in Fig. 4. In each case, the initial distribution is a Gaussian one in the upper excited state. When the electronic gap is ∆ = 1000 cm -1 (upper panels of Fig. 4), the environmental dynamics of the TM presented in the left panels of Fig. 5 are similar to the behavior already obtained in a previous work in a one-dimensional model [START_REF] Chin | Visualising the role of non-perturbative environment dynamics in the dissipative generation of coherent electronic motion[END_REF] . One observes the initial displacement of the average position towards the equilibrium position of the upper state after the preparation in the FC region and during about 50fs the average position is similar to the trajectory of a Gaussian wave packet with a decreasing norm due to the non-adiabatic transition induced by the CM. This can be related to the so-called "Zeno time" as analyzed in ref. [START_REF] Chen | Mapping of wave packet dynamics at conical intersections by time-and frequency-resolved fluorescence spectroscopy: A computational study[END_REF] . The component appearing in the lower state also exhibits damped oscillations. The equilibrium distribution in the final state is reached after 250 in agreement with the correlation time scale (Fig. 3b). For the coupling mode (right panels), the Gaussian distribution smoothly melts at the initial equilibrium position. The density appears in the lower state with two symmetric lobes, fulfilling for a short time (about 10fs) the expected node at the equilibrium position since the electronic coupling is zero at this point. The damping then evolves towards a uni-modal distribution in agreement with the shape of the potential energy curves with a single minimum.

The distribution P α (x t ,t) (Eq.( 25)) of the tuning bath in the ∆ = 0 case (lower panels of Fig. 4) evolves during the first fifty femtoseconds as in the previous example but the final norm is not zero since the final population is the same in both states. The asymptotic distribution is Gaussian in both states and centered at the respective equilibrium positions. When the energy gap decreases, the profile of the diabatic and adiabatic lower potential energy curves along the coupling mode are very different, with either a minimum or a barrier at the equilibrium position. Interestingly, the asymptotic distribution of the coupling bath becomes bi-modal. Even if for the reduced electronic system Hamiltonian H S the diabatic and adiabatic basis sets coincide for the chosen reference point, the full entangled system-bath dynamics provide a vibrational asymptotic distribution that correctly "feels" the adiabatic potential of the coordinate representation. This reveals the efficiency of the information extracted from the auxiliary matrices.

Fig. 6 compares the full statistical distribution P α (x t , x c ,t) (Eq.29) in the two diabatic electronic states at two times with the square modulus of the two-dimensional undamped wave packet |Ψ α (Q t , Q c ,t)| 2 in the diabatic representation. We illustrate only the case with ∆ = 1000 cm -1 (first two rows in Fig. 5 and corresponding TLS-TC model in the first row of Fig. 4). The first chosen time t = 25 fs is shorter than the Zeno time so that the behaviors are expected to be very close. The only noticeable difference concerns the nodal line in the wave packet map that is already blurred at this short time in the damped case. On the contrary, for t = 120fs, the damped distribution has nearly reached its asymptotic shape around the equilibrium position of the lower electronic state while the undamped wave packet remains fully delocalized.

For the same example with ∆ = 1000 cm -1 , Fig. 7 gives the imaginary part of the coherence (Eq.( 30)) in the diabatic representation for the damped and the undamped wave packet. At short time t = 25 fs, this confirms the efficiency of the extraction of information from the off diagonal elements M αβ n,m (t) (Eq.( 28)). At longer time t = 120 fs, the coherence becomes two orders of magnitude smaller in the damped case while it remains at the same order of magnitude in the undamped case, as expected.

We compare in Fig. 8 the two-dimensional distribution in the lower excited state in the diabatic (left panels) and in the adiabatic (right panels) representations, for the two gaps ∆ = 1000 cm -1 and ∆ = 0 illustrated in Fig. 4 and for two times, a short ( t = 25 fs) and a long one (t =240 fs). The adiabatic representation is obtained through Eqs.(32) using the parameters of the TLS-TC model to build the diabatic-to-adiabatic transformation matrix at every point of the grid.

Electronic dynamics

Fig. 9 gives the population of the TLS model in the upper excited state when the energy gap ∆ decreases from 1500 cm -1 to zero in the situation where the spectral densities are centered at Ω c = 500 cm -1 for J c (ω) and Ω t = 1000 cm -1 for J t (ω) with renormalization energies λ c = 150 cm -1 and λ t = 188 cm -1 . The population evolution for the undamped TLS-TC is shown in Fig. 9b. One observes a perfect agreement during the first thirty femtoseconds before the damping operates. This corresponds to a period of the tuning mode as may be seen in Fig. 9c where is drawn the first moment of the collective mode X t /D t of the tuning bath in the upper electronic state. The damped TLS decay curves show markedly different profiles from a smooth overdamped decay (long dashes in Fig. 9) for the large gap ∆ = 1500 cm -1 to damped oscillatory decay with variable periods of partial recurrence and asymptotic values. The different recurrence periods and initial decay yields during the first oscillation are already observed in the undamped TLS model (Fig. 9b).

These general trends may be explained by the way the FC wave packet will explore the CI region. Varying the diabatic gap modifies the location of the CI and the slopes at the crossing. In the present examples, the CI is of the sloped type (see upper panels in Fig. 4), except for the case ∆ = 0 leading to a peaked case (see lower panels in Fig. 4) [START_REF] Atchity | Potential energy surfaces near intersections[END_REF] . For the largest gap ∆ = 1500 cm -1 , the CI is above the FC region and the smooth decay occurs via tunneling. In the case ∆ = 500 cm -1 on the contrary, the CI is close to the turning point. In the ∆ = 0 case, the CI is very close to the FC region, thus leading to early Stueckelberg oscillations 126 .

Another way to explain the different oscillation timescales in the population evolution is the way the initial wave packet is expanded in the eigen vibronic basis set. This superposition depends on the gap ∆ and involves different energy gaps between the main eigenstates and therefore different oscillation periods. The variation of ∆ also leads to different asymptotic populations (Fig. 9a) and different average positions of the collective mode (Fig. 9c). The final equilibrium Botzmann mixture concerns the eigen vibronic states and not only the electronic states [START_REF] Iles-Smith | Environmental dynamics, correlations, and the emergence of noncanonical equilibrium states in open quantum systems[END_REF] .

Finally Fig. 9d gives the volume of the accessible states in the Bloch sphere, which is a signature of non-Markovianity when the decay is not monotonic. The non-Markovianity signature is given rigorously by the integrals of the region where the volume increases but the shape of the bumps provides a qualitative information. In the context of non-adiabatic dynamics, we have already discussed [START_REF] Chin | Visualising the role of non-perturbative environment dynamics in the dissipative generation of coherent electronic motion[END_REF] that this volume is closely linked to the correlation function of the wave packet in the corresponding effective mode model. It reveals the tendency of both the electronic and vibrational wave packets to return towards the initial position. The black solid curve in Fig. 9d gives the volume for the pure dephasing case of the only TM. The bumps merely correspond to the recurrences of the bath correlation function of the TM alone C(t) (see Fig. 3b) and should match with the oscillations of the wave packet in the effective coordinate model. When the CM is taken into account, the electronic transition (see populations in Fig. 9b) increases the decay rate and blurs the bumps. Electronic dynamics of a CI described by a TLS is non Markovian in the present example even if the population decay seems smooth (however non exponential) as it is the case for ∆ = 1500 cm -1 (long dashes). x x 10 -3

x x 10 -3 FIG. 7. Electronic coherence χ(x t , x c ,t) (Eq.( 30)) computed by HEOM in the TLS model and wave packet coherence ξ (Q t , Q c ,t) for two times, t = 25 fs and t = 120 fs. The corresponding TLS-TC model is represented in the upper panels of Fig. 4. The gap between the two diabatic minima is ∆ = 1000 -1 , λ t = 188 cm -1 and λ c = 300 cm -1 . T = 298K.

B. Initial electronic state

We now examine the difference of behaviors related to the initial preparation of the electronic system. The energy gap is ∆ = 1000 cm -1 and the renormalization energies are λ t = 188 cm -1 and λ c = 150 cm -1 . Control strategies with ultrafast designed laser fields could modify the initial superposed electronic state. For the examples in the previous section, the initial state is prepared as |2 . As the center of the FC region is located at Q c = 0, the electronic coupling vanishes and the diabatic and adiabatic electronic states coincide. The initial state is thus a pure state, which does not involve electronic coherence. The relaxation of the coupling vibrational bath and the evolution towards equilibrium in the final state maintain the symmetry as it has been illustrated in the previous examples.

When the initial state is a superposed state, electronic coherence is created and persists during the timescale of the correlation function (about 200 fs). Panels (a) and (b) of Fig. 10 give two examples for the initial superposed states

|χ 1 = |1 + √ 9 |2 / √ 10, and |χ 2 = (|1 + |2 ) / √ 2
, respectively. In addition to the initial electronic coherence, the new behavior is the breaking of symmetry in the vibrational relaxation and evolution towards equilibrium of the coupling bath. This is illustrated for P(x c ,t) in state |1 . This unexpected vibrational effect is already observed in the evolution of the undamped two-dimensional wave packet as can be seen in panel (d) of Fig. 10 where we compare the average position of CM computed by the undamped two-dimensional wave packet and by HEOM. The curves are labeled UNDAMPED1 and DAMPED1 for |χ 1 , and UNDAMPED2, DAMPED2 for |χ 2 . The dissymmetry comes from the interferences of the left or right parts of the upper wave packet when it makes a transition due to the antisymmetrical electronic coupling. The initial wave packet in the lower state has everywhere the same sign while the decaying parts get the opposite sign due to the coupling. The resulting wave packet in the lower state becomes asymmetrical. This effect is illustrated in panel (c) of Fig. 10 where we give the square modulus of the undamped wave packet in the lower state at a short time t = 12 fs.

C. Variation of Γ c Fig. 11 compares iso-contours of the projection P (Q t ,t) (Eq.( 24)) in the lower electronic state |1 of the localization density probability for the undamped two-dimensional wave packet dynamics with the statistical distribution P(x c ,t) (Eq.25) for increasing width of J c (ω), i.e., increasing the friction coefficient of the residual bath in the TLS-TC model. The energy gap is ∆ is 1000 cm -1 and the two renormalization energies are λ t = 188 cm -1 and λ c = 150 cm -1 . The adiabatic lower potential energy curve has a single minimum along the CM. When the width Γ c of the spectral density J c (ω) is modified, the weight p c is adapted to impose the same λ c and therefore the same linear coupling D c . The damping begins after 30fs and one sees the blurring of the characteristic nodal line of the distribution and the evolution towards the final distribution is faster with increasing Γ c .

IV. CONCLUSIONS

The main result of this work is to illustrate the ability of the powerful tool provided by HEOM to get an interesting information about the statistical distribution of the position of the collective mode associated with different baths from the auxiliary matrices. We have extended for two uncorrelated baths the methodology suggested by Shi et al. [START_REF] Zhu | Explicit system-bath correlation calculated using the hierarchical equations of motion method[END_REF][START_REF] Liu | Reduced quantum dynamics with arbitrary bath spectral densities: Hierarchical equations of motion based on several different bath decomposition schemes[END_REF] by deciphering the partial level related to each bath for each auxiliary matrix. We have derived the projected and the twodimensional distributions of the tuning and coupling modes of a conical intersection from the minimal TLS electronic partitioning. This strategy already proposed in the literature [START_REF] Duan | Quantum mechanical wave packet dynamics at a conical intersection with strong vibrational dissipation[END_REF][START_REF] Dijkstra | Simulation of photo-excited adenine in water with a hierarchy of equations of motion approach[END_REF] could provide a better computational scaling than the more familiar two-state two-coordinate representation [START_REF] Chen | Dissipative dynamics at conical intersections: Simulations with the hierarchy equations of motion method[END_REF] . The present work shows how to complement this approach by visualizing the evolution of the environmental distributions in each electronic state that correspond to the square modulus of the dissipative wave packets. When the electronic dynamics induces a strong restructuring of the environment, one may visualize the departure from the initial Gaussian equilibrium distribution (within the initial factorization assumption). By using the parameters of the associated TLS-TC model, we have also illustrated how to derive the adiabatic distributions from those computed in the diabatic representation in the TLS system. In the present work, we have considered only a very symmetrical model and an initial reference point located in a region where the electronic coupling vanishes so that diabatic and adiabatic states coincide at this point. As a consequence, the minimal electronic system will not give any information on the pro- 29)) and P adia 1 (x t , x c ) (Eq.( 32)) computed by HEOM in the TLS model and using the parameters of the TLS-TC model to obtain the adiabatic map. The gaps are ∆ = 1000 cm -1 and ∆ = 0 (see Fig. 4). The distributions are given for two times, t = 25 fs and t =240 fs. The color code gives P(x t , x c ) × 10 3 . λ t = 188 cm -1 and λ c = 300 cm -1 . T = 298K. jected electronic coherence if the initial state is a single state. However, maps of the real and imaginary parts of the coherence corresponding to the cross product of the wave packet components in the two electronic states may be nicely derived from the auxiliary matrices.

Working in the TLS model with HEOM is the opposite strategy to the projection of electronic dynamics from a full dimension wave packet and it is particularly rewarding to see the efficiency of this HEOM technique. However, high levels of hierarchy are required when the distribution becomes bi-or multi-modal since the distribution is expanded on a Hermite polynomial basis set and convergence may require higher levels than the simple population computation. We have used a home made code parallelized with OpenMP. However other implementations on new architectures such as GPU could be useful 127 .

The inclusion of the TM and CM in their respective bath are expected to generate two highly structured spectral densities and therefore long correlation times and strongly non-Markovian dynamics as confirmed by the non monotonic decay of the volume of accessible states in the Bloch sphere in our example. As already discussed in a previous work [START_REF] Chin | Visualising the role of non-perturbative environment dynamics in the dissipative generation of coherent electronic motion[END_REF] , the volume evolution is closely linked to the autocorrelation function of the nonadiabatic wave packet of the corresponding TLS-TC model. The back flow of the electronic energy is induced by the return of the collective mode towards its initial position. Another measure of non-Markovianity during exciton-phonon dynamics has also been extracted from the auxiliary matrices 128 .

It is worthwhile noting that even if we have treated the specific example of a conical intersection, the method by itself has a large panel of applications to analyze the electronvibration dynamics of complex systems such as excitonphonon dynamics. HEOM offers the advantage of easily accounting for high temperature and this paves the way towards interesting applications for instance by considering baths at different temperatures [START_REF] Wertnik | Optimizing co-operative multi-environment dynamics in a dark-state-enhanced photosynthetic heat engine[END_REF]129,130 . Moreover, the possibility to include radiative couplings and exploit laser induced conical intersection as recently suggested 131,132 in order to control the electronic dynamics is also a possible future application. Optimal control with laser fields has recently been implemented with HEOM 133 and the new tool could provide an interesting information about the controlled environmental dynamics.

SUPPLEMENTARY MATERIAL

The supplemental material presents some complementary examples of electronic dynamics by varying the linear couplings D t/c of the TLS-TC model extracted from the TLS model (Eq.( 19)) or the peak positions of the spectral densities fixed by the parameters ω t/c (Eq.( 13)).

FIG. 1 .

 1 FIG. 1. Left panel: Schematic view of a conical intersection in the adiabatic representation along the branching space of the tuning Q t and coupling Q c modes. Right panel: the diabatic representation in which the electronic coupling is varying linearly along the coupling mode.

FIG. 3 .

 3 FIG. 3. (Color online) Panel (a): Spectral densities for the tuning (solid blue line) and coupling (dashed red line) baths of the TLS model. The renormalization energies (Eq.18 ) are λ t = 188 cm -1 and λ c = 300 cm -1 , respectively. Panel (b): real parts of the corresponding normalized bath correlation functions C t/c (t) (Eq.(12)) for T = 298K. The parameters of the spectral densities (Eq.(13)) are p t = 1.95 × 10 -11 a.u., p c = p t /2.5 ; ω t = 4.5639 × 10 -3 a.u. (1000 cm -1 ), ω c = ω t /2; Γ t = Γ c = 2.7188 × 10 -4 a.u. (60 cm -1 ). λ t = 188 cm -1 and λ c = 300 cm -1 . The p t/c are taken as parameters to modify the renormalization energy λ t/c and therefore the linear couplings.

FIG. 4 .

 4 FIG. 4. (Color online)Potential energy curves of the TLS-TC model corresponding to the TLS with the spectral densities given in Fig.3at the reference point Q t = 0 and Q c = 0. Top panels: diabatic energy gap ∆ = 1000 cm -1 ; lower panels: ∆ = 0. Left panels: tuning mode Q t (the diabatic and adiabatic curves coincide). The curves correspond to an inverted Marcus-type crossing when ∆ = 1000 cm -1 and to a normal Marcus-type crossing when ∆ = 0. Right panels: coupling mode Q c , dashed lines diabatic curves, solid lines adiabatic curves. The renormalization energy λ t = 188 cm -1 is added to ε 2 . λ c = 300 cm -1 .

FIG. 5 .FIG. 6 .

 56 FIG.5. Projected statistical distributions P α (x t/c ,t) (Eq.25) on the collective coordinates computed in the TLS model for the corresponding TLS-TC model represented in Fig.4. The color code gives P(x t/c ,t) × 10. The left/right panels show the tuning/coupling baths, respectively. The diabatic energy gaps are ∆ = 1000 cm -1 and ∆ = 0. λ t = 188 cm -1 and λ c = 300 cm -1 . The spectral densities are given in Fig.3a. T = 298K.

FIG. 8 .

 8 FIG.8. Statistical distributions in the lower excited state P dia 1 (x t , x c ) (Eq.(29)) and P adia

  FIG. 9. (Color online) Variation of dynamical quantities with the electronic energy gap between the diabatic minima ∆. Upper panels : Population in the initial upper electronic state; panel (a)TLS damped HEOM dynamics, panel (b) undamped TLS-TC wave packet dynamics (b). Panel (c) : first moment of the collective TM in the upper electronic state extracted from the auxiliary matrices. Panel (d) Volume of the accessible states in the Bloch sphere. The parameters of the spectral densities J t (ω) and J c (ω) are Ω t = 1000 cm -1 , λ t = 188 cm -1 and Ω c = 500 cm -1 , λ c = 150 cm -1 , T = 298K.

FIG. 10

 10 FIG. 10. (Color online) Statistical distribution P(x c ,t) projected on the coupling mode in the lower excited state |1 for two initial superposed states. Panels (a) and (b): Preparation in states |χ 1 = |1 + √ 9 |2 / √ 10 and |χ 2 = (|1 + |2 ) / √ 2, respectively (The color code gives P(x c ) × 10); panel (c) distribution of the undamped wave packet in the lower excited state at time t = 12 fs for state |χ 2 (the color code gives |Ψ| 2 × 10 3 ; panel (d): Average position of the effective mode Q c (t) computed with the undamped twodimensional wave packet and of the collective coordinate x c (t) obtained from the damped HEOM. λ t = 188 cm -1 and λ c = 150 cm -1 , ∆ = 1000 cm -1 , T = 298K.

FIG. 11

 11 FIG. 11. (Color online) Comparison of the undamped twodimensional wave packet projected dynamics P(Q c ,t) (panel (a)) and the projected statistical distribution for the CM P(x c ,t) in the lower electronic state for increasing width Γ c and the same linear coupling D c (panels (c) and (d)). The color code gives P × 10. Panel (b) gives the decay of the population for the undamped wave packet and for the HEOM dynamics with increasing width. T = 298K, ∆ = 1000 cm -1 , λ t = 188 cm -1 and p c is adjusted to give the same λ c = 150 cm -1 when Γ is modified.

  2S 11 S 21 ℜ e (χ 12 ) P ad 2 = S 2 12 P dia 1 + S 2 22 P dia 2 + 2S 12 S 22 ℜ e (χ 12 )

	P ad 1 = S 2 11 P dia 1 + S 2 21 P dia 2 +

one easily gets at each point (x t , x c ) :

ACKNOWLEDGMENTS

We warmly acknowledge Dr A.X. Chin for interesting discussions. One of us (EM) acknowledges support from the ANR-DFG COQS, under Grant No. ANR-15-CE30-0023-01. This work has been performed within the French GDR 3575 THEMS.

Appendix A: HEOM for two uncorrelated baths

Solving the hierarchical equations of motion for a conical intersection requires set of equations for two In this derivation, we mainly assume that they are uncorrelated. Such expressions have been derived in several works, for particular spectral densities [START_REF] Duan | Quantum mechanical wave packet dynamics at a conical intersection with strong vibrational dissipation[END_REF][START_REF] Dijkstra | Simulation of photo-excited adenine in water with a hierarchy of equations of motion approach[END_REF] . Here, we essentially aim at recalling the derivation of the coupled equations and the assumptions that we used.

Auxiliary matrices in the two-uncorrelated-bath case

The reduced density matrix ρ S,I (t) for the system in interaction representation ρ S,I (t) = e iH 0 t ρ S (t)e -iH 0 t with H 0 = H S + H B , , is given by the partial trace over the bath B of the time evolution of the total density matrix,

where T (+) is a time-ordering operator and a factorization is assumed at the initial time t = 0. L(t) is the Liouvillian in interaction representation, which reads

(A2) The partial trace over B is a successive trace over bath B t and bath B v :

] and ρ eq B = ρ eq B t ρ eq B c . In the case of harmonic oscillators with linear coupling, it has been shown that the second-order cumulant expansion of Eq. (A1) is exact [START_REF] Kubo | Stochastic Liouville Equations[END_REF][START_REF] Ishizaki | Quantum Dynamics of System Strongly Coupled to Low-Temperature Colored Noise Bath: Reduced Hierarchy Equations Approach[END_REF]101,134 ,

Owing to this feature, A3) is a second-order memory kernel similar to the one found in second-order perturbation theory. It will be denoted K (τ) = τ 0 dt K [START_REF] Michl | Electronic Aspects Of Organic Photochemistry[END_REF] (τ,t ). First, one can easily derive from Eq.(A2) an expression for the two-bath case,

When the two baths are not correlated, cross-terms in Eq.(A4) are considered as negligible and discarded for the rest of the present derivation. For Eq.(A4), one can apply strategies similar to the one used for second-order perturbation theory. By using usual open quantum-system methods 135 , one can eliminate the bath and cast its whole influence in two-time correlation functions

] for the tuning and coupling mode respectively:

In this respect, using correlation function parametrizations C t/c (τt ) given in Eq.( 16) and its complex conjugate C * t/c (τ -t ) (Eq.( 17)) and following strategies proposed for instance in [START_REF] Pomyalov | The importance of initial correlations in rate dynamics: A consistent non-Markovian master equation approach[END_REF] and [START_REF] Ishizaki | Quantum Dynamics of System Strongly Coupled to Low-Temperature Colored Noise Bath: Reduced Hierarchy Equations Approach[END_REF]101,134 , equation (A5) is separable into a sum of operators K k={ j,l j } ,

where K k={ j,l j } (τ) is assigned to each exponential element l j of the correlation function of bath j = {t, c}. In addition, we decompose it as

where we have defined

and

with Θ k (t ) = α j,l j S j t -α j,l j S j t . (A10) By Fourier-Laplace transforming Eq.(A1) [START_REF] Tanimura | Time Evolution of a Quantum System in Contact with a Nearly Gaussian-Markoffian Noise Bath[END_REF]136 , one can recast Eq. (A3) as

ρ S,I (0).

(A11) Finally, we define auxiliary matrices 101,136 as

It is interesting to note that ρ S,I = ρ (0) 0,••• ,0 is the system density matrix we are looking for. To better clarify the notations, let us emphasize that • k = {l j , j} is a global index related to an exponential element l j of the j's bath correlation function and the j = {t, c} bath itself.

• n k is an integer that is a power of W k .

• K = n cor,1 + n cor,2 is the total number of exponential functions in the correlation function parameterization (Eq.16),

} is an index vector bound to a hierarchy matrix.

• (ν) is the hierarchy level of the matrix defined such that

Derivation of the hierarchical equations of motions

For any k = {l j , j}, one has the following relation,

ρ S,I (0).

(A13)

By performing an integration by parts, one obtains

δ is the Kronecker symbol (δ ν0 = 1 if ν = 0, 0 otherwise).

For auxiliary matrices, we assume that for all ν and vector {n

With an inverse Fourier-Laplace transform ( for a differentiable function f , +∞ 0 dte -iωt d dt f (t) = iω f [ω]f (0)), this leads to the set of hierarchical equations of motion:

(A15)


Amel Fraisse

Zheng Zhang

Shelley Fisher Fishkin

Ronald Jenn

Preserving Endangered European Cultural Heritage and Languages Through Translated Literary Texts

Keywords: under-resourced languages, parallel corpus, translated literary

We present the interdisciplinary ROSETTA project which consists in collecting all the translations worldwide of one fictional text in order to build multilingual parallel corpora for a large number of under-resourced languages. Building such corpora is vital to help preserve and expand language and traditional knowledge diversity. These corpora will be useful to handle under-resourced languages in a number of interconnected research fields such as computational linguistics, translation studies and corpus linguistics. Our project taps into a wealth of translated versions of a single fictional text spanning a period of over a century. It consists in collecting, digitizing, transcribing and aligning translations of this text. Our data collection process is based on volunteer work from the scientific and scholarly communities, the power of the crowd and national libraries and archives. Our first experiment was conducted on the world-famous and well-traveled American novel "Adventures of Huckleberry Finn" by the American author Mark Twain. This paper reports on the parallel corpus that are now sentence aligned pairing English with Basque.

1 The example of Mark Twain's text for under-resourced languages -Mark Twain's books are some of the most well-travelled texts on the planet.

-The advantage of using a work of fiction such as "Adventures of Huckleberry Finn", is that it uses a very broad vocabulary linked to every day life, which makes it a valuable asset for those languages that are currently lacking such computational resources. 1. Identifying existing translations : calling on the international community of Mark Twain scholars as well as Translation Studies scholars in order to identify existing translations in different languages.

Data Collection

2. Mining online digital libraries and national archives : using the title in the target languages, we crawled the web and mined online digital libraries and national archives in order to find the full texts.

3. Text transcription : We used the CrowdFlower platform to transcribe digital versions that came as images.

Experiments and Results

English -Basque Chapter 38 English Basque

Making them pens was a distressid-tough job, and so was the saw ; and Jim allowed the inscription was going to be the toughest of all. That's the one which the prisoner has to scrabble on the wall. But we had to have it ; Tom said we'd got to ; there warn't no case of a state prisoner not scrabbling his inscription to leave behind, and his coat of arms.

Idazteko lumak egitea oso lan nekagarria zen, eta zerra egitea ere bai, eta Jimek esan zuen inskripzioa egitea lanik gogorrena izango zela. Hori presoak horman egin behar zuen zirriborroa zen. Baina nahitanahiez izan behar genuen inskripzioa. Tomek esan zuen egin behar genuela ; ez zela estatuko preso bakar bat inskripzioa idatzita utzi ez zuenik, bere armarria eta guzti gainera. "Why, Mars Tom, I hain't got no coat o' arms ; I hain't got nuffn but dishyer ole shirt, en you knows I got to keep de journal on dat." -Baina, Tom jauna -esan zuen Jimek-, nik ez dut armarririk ; ez daukat alkandora zahar hau besterik, eta badakizu bertan egunkaria idatzi beharra dudala. "Oh, you don't understand, Jim : a coat of arms is very different." -Ez duzu ulertzen, Jim ; armarria bestelako gauza da.

FIGURE 1 -

 1 FIGURE 1 -The UNESCO Index Translationum

FIGURE 3 -

 3 FIGURE 3 -The global knowledge diagram representing a subset of existing translations of Adventures of Huckleberry Finnin in Basque and Catalan.

FIGURE 4 -

 4 FIGURE 4 -Pie chart of the chapter distribution of Basque translation compared with the original English version.

TABLE 1 -

 1 List

of languages "Adventures of Huckleberry Finn" was translated into.

TABLE 2 -

 2 Original text of the chapter 38 of "Adventures of Huckleberry Finn" and its translations in Basque