Prof Vuilleumier Rodolphe

Prof Giovanni Ciccotti

Coupled-trajectory mixed quantum-classical GD Gauge-dependent GI Gauge-invariant HPC High-performance computing IC Internal conversion ISC

BO

Born-Oppenheimer CT-MQC

Overview of my research activity

My research activity in the domain of molecular dynamics started during my PhD with the study of the transfer process of an excess proton in water. The particularly high mobility of an excess positive charge in water has always attracted much attention. A lone proton H + does not exist in water but it is weakly bonded, via a covalent bond, to a water molecule forming the hydronium cation H 3 O + . Hydronium's high mobility in water [1] can be justified qualitatively by invoking Grotthuss mechanism [2]: charge migration is the consequence of successive proton jumps between oxygen atoms, taking place, preferentially, along the hydrogen bonds. To describe Grotthuss mechanism in a molecular dynamics simulation one has to be able to account for the reorganization of chemical bonds during the proton-jump event, either using ab initio techniques or "reactive" force fields. Due to the system size limitations imposed by explicit calculations of the ground-state electronic structure in ab initio methods, I have focused on the second approach. The empirical valence bond (EVB) method allows to account for changes in the chemical-bond structure around the excess charge, by adapting the force fields used to describe different chemical species, i.e., water molecules and hydronium cation. By calculating infrared spectra of different protonated water clusters in the gas phase and of bulk protonated water, it was possible to identify [3,4] the vibrational signatures of the most stable structure formed in protonated water, i.e., the Eigen cation H + (H 2 O) 4 , and a transient structure, i.e., the Zundel cation H + (H 2 O) 2 , charecterizing the proton-transfer configuration. During my first postdoc I worked on the development of an analysis tool to characterize the local microscopic structure of liquid systems. The method is based on the statistical interpretation of the positions of particles, e.g., of atoms, in a liquid system, with the aim of identifying regions of space where those particles are more likely localized. Such regions are the maximum probability domains, defined as the regions where the probability of finding one, and only one, particle is maximum, for a given choice of the reference frame. I have developed and implemented the algorithm that searches for the maximum probability domains by computing the shape derivative of the domains and by maximizing the probability based on a geometric optimization method called the level set method. Maximum probability domains have been determined [5] in the case of water at ambient and high pressure, and for Na + cations diluted in water.

During my second postdoc I started working on excited-state molecular dynamics, which is still today the main focus of my research activity. An ultrashort visible-UV laser pulse can create coherent superpositions of electronic and vibrational states. As a consequence, excited electrons and vibrational modes exchange energy, HDR Excited-state molecular dynamics and the molecule can relax via radiationless channels at an ultrashort time scale. This kind of processes are called nonadiabatic because they involve changes of the electronic states of molecules -nonadiabatic transitions -induced by the strong coupling with nuclear motion. My goal is to be able to describe nonadiabatic processes accounting for various effects: (i) the changes in the electronic character occurring during the evolution of the system; (ii) the role of nuclear motion on electronic dynamics, accounting for the quantum nature of the nuclei and for anharmonic effects; (iii) electron-nuclear correlation beyond the weak-coupling regime necessary to capture effects such as quantum decoherence and coherence revivals. Simulation methods that fit all the requirements necessary to describe nonadiabatic phenomena triggered by light in complex molecular systems are trajectory-based quantum-classical schemes. Combining a classical description of nuclear motion with a quantum-mechanical description of electronic dynamics, trajectory-based approaches are well-suited for treating large systems at the atomistic level, since they allow for the on-the-fly calculation of electronic structure properties. The exact factorization of the electronnuclear wavefunction [6][START_REF] Agostini | TDDFT and Quantum-Classical Dynamics: A Universal Tool Describing the Dynamics of Matter[END_REF][START_REF] Agostini | Exact factorization of the electron-nuclear wavefunction: Theory and Applications[END_REF][START_REF] Agostini | [END_REF] offers a new perspective to analyze and simulate nonadiabatic phenomena. It will be extensively presented in the following sections. Initially, my work focused on the analysis of the properties of the evolution equations [10][11][12][13][14][15] derived from the exact-factorization Ansatz, with the aim to link the new picture provided by the exact factorization to a more standard picture used to interpret nonadiabatic processes [START_REF] Agostini | [END_REF]. Afterwards, I have worked on proposing various strategies to solve the exact-factorization equations in approximate ways, with the aim to describe nonadiabatic processes in molecules. Based on my expertise in the field of quantum-classical dynamics developed during my master degree [16,17], I followed this route to propose a new trajectorybased approach to excited-state dynamics [18][19][20][21][22][23][24]. Very recently those theoretical and numerical developments have led to the application of the exact factorization in simulating the photo-induced ring-opening process in oxirane [25,26] (see Fig. 2). An alternative route for approximations, applicable for weakly nonadiabatic processes, is to treat nonadiabaticity as a perturbation. Such formulation of the problem can be derived by showing that the adiabatic regime can be obtained starting from the exact-factorization equations as the limit for the electron-nuclear mass ratio µ = m e /M n going to zero [27]. The numerical method thus derived has been applied to compute electronic current densities in a model for proton-coupled electron transfer [28], for vibrational circular dichroism in chiral molecules [29], for high-precision vibrational spectroscopic of hydrogen-based molecules [30]. Recently, I have also been interested on geometric-phase effects in molecular systems. I have conducted preliminary studies [31,32] on model systems to understand how the standard picture, i.e., that of conical intersections, is translated into the language of the exact factorization. These preliminary studies have shown that none of the typical features related to conical intersections appear in the framework of the exact factorization.

In the next sections, I will present an in-depth analysis of my work on the exact factorization and on excited-state molecular dynamics, with major focus on the development of trajectorybased quantum-classical algorithms and on the properties of the exact factorization in the presence of conical intersections. These topics are currently the focus of my research activity, and future projects aim at pursuing these research interests.

Excited-state molecular dynamics

Ab initio molecular dynamics is nowadays a cornerstone in the fields of Theoretical Chemistry and Chemical Physics. Numerical simulations allow us to achieve a profound understanding of structural and dynamical properties of matter at the microscopic level, for a large variety of systems, from isolated molecules, or liquids, to proteins, and a large variety of phenomena, from equilibrium processes to ultrafast photo-activated relaxation. Ab initio molecular dynamics relies on the approximate numerical solution of the time-dependent Schr ödinger equation for systems of electrons and nuclei. Approximations are vital for studying complex systems due to the exponential scaling of the computational cost with the number of degrees of freedom of the exact solution. Major effort in the theory community is, therefore, devoted to the development of such approximate treatments of the full quantum-mechanical problem.

In 1927, Born and Oppenheimer [33] posed the basis for what is currently the standard way of visualizing motion in molecular systems: the electrons generate the potential(s) felt by the nuclei during their evolution. This picture is retained either within the adiabatic, Born-Oppenheimer, approximation, where the electrons are "statically" found in a single eigenstate, or in nonadiabatic situations, when excited states are populated along the dynamics and exchange amplitude. In what follows, the expression "Born-Oppenheimer representation" or "Born-Oppenheimer framework", to be distinguished from the "Born-Oppenheimer approximation", will be used to indicate the interpretation of molecular dynamics in terms of static electrons providing the potentials that guide nuclear dynamics.

The concept of electronic (adiabatic) potential energy surfaces lends itself to introducing the idea of classical-like nuclei evolving in time along those potential energy surfaces. In fact, perhaps the most successful approximation used in the field of quantum molecular dynamics simulations is the use of -quantum, semiclassical, or classical -trajectories to sample the nuclear configuration space and to mimic the quantum-mechanical nuclear motion. In the adiabatic approximation, a single electronic eigenstate (usually the ground state) is populated during the dynamics, thus trajectories are driven by the force determined from the corresponding potential. However, in nonadiabatic conditions, excited states become "available" during the dynamics, and a unique choice for determining the classical force that drives the trajectories seems not possible. Consequently, a commonly-agreed and generally-applied procedure to guide nuclear trajectories does not yet exist, and various algorithms have been proposed over the years. Trajectory surface hopping [34], where trajectories can hop between states switching from one adiabatic potential energy surface to another, or the Ehrenfest approach [35], where trajectories are evolved on a mean-field potential created by the electrons, mimic the quantum nuclear evolution using purely classical trajectories. Other approaches, instead, use the support of trajectories to follow nuclear motion in time. This is the case of multiple spawning [36,37], where Gaussian basis functions evolving in adiabatic states are spawned in regions of strong coupling between states and are coupled together, or of DD-vMCG (direct-dynamics variational multiconfiguration Gaussian) [38][39][40], where the trajectories obey non-classical equations determined from the quantum-mechanical variational principle.

In such a beyond-Born-Oppenheimer representation of the coupled motion of electrons and nuclei, curious features of the adiabatic potential energy surfaces emerge, which are often invoked to interpret relaxation processes undergone by photoexcited molecules. Those features are conical intersections, typical examples of the breakdown of the Born-Oppenheimer approximation. Conical intersections represent efficient funnels [START_REF] Michl | Electronic Aspects of Organic Photochemistry[END_REF][START_REF] Michl | [END_REF][43] for population transfer between electronic states, mediated by nuclear motion. They are regions of configuration space where the adiabatic potential energy surfaces are degenerate -thus producing a diverging coupling between the corresponding electronic states -and exhibit, within the so-called branching space, a double-cone shape. In the literature, conical intersections have been studied not only for their role in nonadiabatic processes, but also for the effect that the related Berry phase has on adiabatic phenomena [44][45][46][START_REF] Domcke | Conical Intersections: Electronic Structure, Dynamics & Spectroscopy[END_REF][START_REF] Worth | [END_REF][START_REF] Baer | Beyond Born-Oppenheimer: Electronic Nonadiabatic Coupling Terms and Conical Intersections[END_REF][START_REF] Levine | [END_REF][51][52][START_REF] Domcke | Conical Intersections: Theory, Computation and Experiment[END_REF][START_REF] Malhado | [END_REF][55][56]. The concept of conical intersections seems to be grounded in the Born-Oppenheimer representation of the electron-nuclear problem, and despite the fundamental and numerical challenges encountered when dealing with conical intersections, such as geometric or topological phases and singularities, this representation is still commonly employed to perform nonadiabatic ab initio molecular dynamics simulations.

An alternative perspective on the electron-nuclear time-dependent problem is described here: the exact factorization of the electron-nuclear wavefunction [6,[START_REF] Agostini | TDDFT and Quantum-Classical Dynamics: A Universal Tool Describing the Dynamics of Matter[END_REF][START_REF] Agostini | [END_REF]57]. The Born-Oppenheimer representation is abandoned, and the molecular wavefunction is expressed as a single product of a nuclear wavefunction and an electronic, conditional wavefunction with parametric dependence on the nuclear positions. These two components of the molecular wavefunction are time-dependent and evolve according to coupled equations, that are derived from the full Schr ödinger equation. Nuclear dynamics is generated by a "standard" time-dependent Schr ödinger equation, where the effect of electronic excited states is encoded in a time-dependent vector potential and a time-dependent potential energy surface [10][11][12][13][14]31]. Starting from the nuclear Schr ödinger equation of the exact factorization, a trajectory-based solution can be derived [18][19][20][21][23][24][25][26]58], where the time-dependent potentials are employed to uniquely define the force [19,20] used to generate classical-like nuclear trajectories. Due to their dependence on time, the vector and scalar potentials are able to encode electronic nonadiabatic effects and to drive nuclear dynamics beyond the Born-Oppenheimer approximation. The vector and scalar potentials can be interpreted as time-dependent generalizations of similar, static, quantities also appearing within the Born-Oppenheimer framework, i.e., the nonadiabatic coupling vectors and the adiabatic potential energy surface(s). A formal connection can be established [27,59], since in the limit of infinite nuclear mass M → +∞, or in the limit of the electron-nuclear mass ratio tending to zero µ = m e /M → 0, the adiabatic quantities are recovered. However, as the theory stands at the moment, for finite values of the nuclear masses, neither the time-dependent vector potential presents singularities at conical intersections in analogy to the nonadiabatic coupling vectors, nor the time-dependent potential energy surface shows peculiar conical shapes related, once again, to conical intersections.

The feature of the exact factorization that suggests promising applications in the field of ab initio molecular dynamics is, as just discussed above, the appearance of the time-dependent potentials, free of signatures strictly related to conical intersections, 1 and indicating a unique route to determine the classical force for a trajectory-based solution of the coupled electronnuclear dynamics. In what follows, after presenting the exact factorization, the time-dependent vector potential and the time-dependent potential energy surface will be discussed, supporting the theoretical analysis with numerical data. The algorithm proposed as a trajectory-based solution of the coupled electronic and nuclear equations of the exact factorization is presented as well, along with its application to the photo-induced dynamics of Oxirane.

The time-dependent molecular problem

A system of interacting electrons and nuclei is described, in the non-relativistic limit, by the molecular Hamiltonian

Ĥ(r, R) = Nn ν=1 -2 2M ν ∇ 2 ν + Te (r) + Vee (r) + Vnn (R) + Ven (r, R) (1)
= Nn ν=1 -2 2M ν ∇ 2 ν + ĤBO (r, R). (2)
The first term on the right-hand side is the nuclear kinetic energy operator, with ∇ ν the gradient with respect to the position of the νth nucleus, and M ν its mass; is the reduced Planck's constant. The symbols r and R are used to indicate the sets of 3N el and 3N n electronic and nuclear coordinates, respectively. The electronic kinetic energy operator is indicated as Te (r), whereas the electron-electron, nucleus-nucleus, and electron-nucleus interaction potentials are indicated as Vee (r), Vnn (R), and Ven (r, R), respectively. The sum of the electronic kinetic energy and of the interaction potentials is denoted ĤBO (r, R), i.e., the Born-Oppenheimer (BO)

Hamiltonian. In what follows, time-dependent Hamiltonians will not be considered, where, for instance, a classical time-dependent field acts on the electrons, Vext (r, t), and/or on the nuclei, Vext (R, t). However, such contribution can be easily included.

The time-dependent Schr ödinger equation (TDSE)

Ĥ(r, R)Ψ(r, R, t) = i ∂ t Ψ(r, R, t), (3)
determines the evolution of the molecular wavefunction, Ψ(r, R, t). As shown in Refs. [6,[START_REF] Agostini | Exact factorization of the electron-nuclear wavefunction: Theory and Applications[END_REF]57], the molecular wavefunction can be factored, without making any approximation, as

Ψ(r, R, t) = χ(R, t)Φ R (r, t), (4)
with χ(R, t) the nuclear wavefunction, and Φ R (r, t) a conditional electronic wavefunction that depends parametrically on the nuclear positions R. Existence and uniqueness, up to a gauge transformation, can be proven by imposing the partial nornalization condition [6,57]

dr |Φ R (r, t)| 2 = 1 ∀ R, t (5)
on the electronic wavefunction. It follows from Eq. ( 5) that the exact nuclear density is

dr |Ψ(r, R, t)| 2 = |χ(R, t)| 2 . (6)
Therefore, the nuclear wavefunction can be written as

χ(R, t) = exp i S (R, t) dr |Ψ(r, R, t)| 2 . ( 7
)
Here, S (R, t) is a real function of nuclear positions and time, which will be discussed below. The total molecular wavefunction is invariant under an (R, t)-dependent gauge transformation,

χ(R, t) → χ(R, t) = exp - i θ(R, t) χ(R, t) (8)
Φ R (r, t) → Φ R (r, t) = exp i θ(R, t) Φ R (r, t), (9)
therefore, uniqueness of Eq. ( 4) is guaranteed upon fixing the gauge freedom with a suitable choice of θ(R, t). The gauge is ultimately related to the phase S(R, t) of the nuclear wavefunction, since choosing S(R, t) is fully equivalent to choosing the gauge. This point will be clarified later on, where possible ways of fixing the gauge freedom are discussed.

As illustrated in great detail in Ref. [23], inserting the factored form of the molecular wavefunction, Eq. ( 4), in the TDSE (3), and using the partial normalization condition (5) [60,61], the evolution equations

ĤBO (r, R) + Ûen [Φ R , χ] -(R, t) Φ R (r, t) = i ∂ t Φ R (r, t) (10)
Nn ν=1

[-i ∇ ν + A ν (R, t)] 2 2M ν + (R, t) χ(R, t) = i ∂ t χ(R, t), (11)
can be derived. The nuclear equation has the form of a standard TDSE, where a time-dependent vector potential and a time-dependent scalar potential represent the effect of the electrons on nuclear dynamics. Since no approximation has been invoked so far in the derivation of Eqs. (10) and (11), those potentials can incorporate nonadiabatic, excited-state effects. The nuclear contribution to Eq. ( 4) evolves according to a TDSE, yields by definition the exact nuclear N-body density (Eq. ( 6)), and, using Eq. ( 11), it can be shown that it also yields the exact nuclear N-body current density. Therefore, χ(R, t) can be considered a genuine nuclear wavefunction. In Eqs. ( 10) and ( 11) the terms responsible for the dynamical, i.e., beyond the adiabatic BO approximation, coupling between electronic and nuclear motion are the electron-nuclear coupling operator Ûen [Φ R , χ] [6, 12, 13, 27-30, 57, 62]

Ûen [Φ R , χ] = Nn ν=1 1 M ν [-i ∇ ν -A ν (R, t)] 2 2 + -i ∇ ν χ χ + A ν (R, t) (-i ∇ ν -A ν (R, t)) , (12)
the time-dependent potential energy surface (TDPES) (R, t) [6,10,11,13,14,23,31,63,64],

(R, t) = Φ R (t)| ĤBO + Ûen -i ∂ t |Φ R (t) r , (13)
and the time-dependent vector potential A ν (R, t) [31,32],

A ν (R, t) = Φ R (t)| -i ∇ ν Φ R (t) r . (14)
The symbol • r indicates an integration over electronic coordinates only. Under the gauge transformations ( 8) and ( 9), the scalar potential and the vector potential transform as

˜ (R, t) = (R, t) + ∂ t θ(R, t) (15) Ãν (R, t) =A ν (R, t) + ∇ ν θ(R, t) . ( 16
)
The evolution equations ( 10) and ( 11) are form-invariant under such transformation.

A large amount of literature can be found on different topics related to the exact factorization of the electron-nuclear wavefunction. For instance, still in the time-dependent context, studies have focused on proposing trajectory-based solutions [18-22, 24-26, 58] of the coupled electronic (10) and nuclear (11) equations, on introducing a perturbation-theory framework to treat situations of weak nonadiabatic coupling between electronic and nuclear motion [27][28][29][30] and to cure the inconsistencies of the BO approximation, on deriving a density-functional theory of the coupled electron-nuclear problem [65,66], on exploiting the "inverse" exactfactorization formalism to derive the exact electronic TDSE with non-classical nuclei [67][68][69], on reformulating the dynamical problem for systems different from electrons and nuclei, as purely electronic systems [70] or electron-photon systems [71][72][73], on analyzing the nature of electron-nuclear entanglement [15]. The static formulation of the exact factorization, relying on the time-independent Schr ödinger equation, has been also subject of extensive studies that focus on the properties of the wavefunctions and of the potentials [74][75][76][77][78][79][80][81][82][83][84][85][86], or, interestingly, on tackling the problem of geometric phases in molecular problems [59,87,88].

In the following sections, the focus will be put on two aspects of the exact factorization: (i) the relation with the, more standard, BO representation of the electron-nuclear problem, based on the analysis of the time-dependent potentials of the theory, (ii) a procedure to solveapproximately -the coupled electronic and nuclear equations.

Born-Oppenheimer framework and exact factorization

Standard approaches to formulate and solve the electron-nuclear dynamical problem rely on the Born-Huang expansion of the time-dependent molecular wavefunction,

Ψ(r, R, t) = l χ (l) BO (R, t)ϕ (l) R (r), (17)
which previously has been referred to as BO framework. Here, {ϕ

R (r)} l=1,... are the eigenstates of the (electronic) BO Hamiltonian ĤBO , determined for each value of the parameter R, and χ (l) BO (R, t) are expansion coefficients. The corresponding eigenvalues are denoted (l) BO (R), the adiabatic, or BO, potential energy surfaces (PESs). We introduce also the matrix of nonadiabatic coupling vectors,

d kl,ν (R) = ϕ (k) R |∇ ν ϕ (l)
R r , which are three-dimensional vectors for each value of the nuclear index ν, to be used in the following. Since the adiabatic states are orthonormal, the nuclear density following from Eq. ( 17) is

|χ(R, t)| 2 = l χ (l) BO (R, t) 2 . ( 18
)
The quantity χ

(l)
BO (R, t) is sometimes indicated as the adiabatic contribution of the nuclear wavepacket in state l. Being a proper electronic basis, the adiabatic states can be used to expand the electronic wavefunction of the exact factorization (4) as well,

Φ R (r, t) = l C l (R, t)ϕ (l) R (r), (19)
where the coefficients satisfy the relation

l |C l (R, t)| 2 = 1 ∀ R, t (20)
by virtue of the partial normalization condition (5). It is interesting to note that, even though the nuclear wavepacket might have contributions, i.e., the χ (l) BO (R, t), in different electronic states, the single-product representation of the molecular wavefunction still holds. In fact,

Ψ(r, R, t) = e i S(R,t) l χ (l) BO (R, t) 2 l C l (R, t)ϕ (l) R (r) (21)
which follows from Eqs. ( 7) and (19), where the first term in parenthesis is the nuclear wavefunction χ(R, t) and the second term is the electronic wavefunction Φ R (r, t).

HDR

Excited-state molecular dynamics

The identity relation between Eqs. ( 4) and (17),

χ(R, t) l C l (R, t)ϕ (l) R (r) = l χ (l) BO (R, t)ϕ (l) R (r) (22)
yields a relation between the expansion coefficients

χ(R, t)C l (R, t) = χ (l) BO (R, t). (23)
Eqs. ( 18), ( 20), (21), and (23) show the relation between the BO representation of the electronnuclear wavefunction and the exact factorization. The key difference between the two representations is that in the latter, the electronic wavefunction acquires an explicit time dependence. By contrast, in the BO framework, the electrons appear in the problem as the eigenstates of the electronic BO Hamiltonian, providing the static potentials driving nuclear dynamics. Therefore, within a BO perspective of nonadiabatic dynamics, processes are visualized employing static adiabatic PESs that are coupled in regions of nuclear space where they are close or degenerate, i.e., avoided crossings or conical intersections. The nuclear wavepacket has adiabatic contributions evolving "on" the different surfaces, and transfers amplitude to other states in the regions of coupling. Such perspective is revisited within the exact factorization.

In the following sections, two situations will be presented and analyzed. First, a onedimensional model system, the Shin-Metiu model [89], will be employed to discuss the properties of the TDPES. It will be shown that a single, uniquely defined up to a gauge [6,10,11,13,57,63,64], PES that changes with time is able to guide the nuclear dynamics in different ways in different portions of nuclear configuration space. It is able to achieve this result simply by changing its shape over time: it has a diabatic shape in the region where the BO PESs are close in energy [10,11,13]; it develops steps connecting adiabatic shapes, forcing the wavepacket to spatially split onto BO PESs of different slopes [10,11,13]; it has a mean-field character when the adiabatic contributions of a nuclear wavepacket are localized in the same region [14]; it presents oscillations to account for interference effects [14]. Clearly, the TDPES encodes all features of (electronic) dynamics, and this information is transferred to the nuclear TDSE. Second, the observations on the TDPES will be generalized to a two-dimensional case. Additionally, the time-dependent vector potential will discussed. The vector potential couples to the momentum in the nuclear TDSE [13], and it can be related to the nuclear velocity field [31,32]. The feature highlighted in this section is that, in general, the time-dependent vector potential is not irrotational, meaning that it cannot be written as the gradient of a scalar function of R and t. The consequence is that the vector potential cannot be gauged away in general situations, a property that is important to keep in mind when approximations to the exact evolution equations, Eqs. ( 10) and (11), are to be developed.

Time-dependent potential energy surface

To proceed with the analysis of the TDPES and its relation to the static BO PESs, it is instructive to decompose the TDPES into gauge-invariant (GI) and gauge-dependent (GD) components, (R, t) = GI (R, t) + GD (R, t), where

GI (R, t) = Φ R (t)| ĤBO |Φ R (t) r + Nn ν=1 2 2M ν ∇ ν Φ R (t)|∇ ν Φ R (t) r - Nn ν=1 A 2 ν (R, t) 2M ν (24)
and

GD (R, t) = Φ R (t)| -i ∂ t |Φ R (t) r . (25)
The second and third terms on the right-hand side of Eq. ( 24) are obtained from the action of the electron-nuclear coupling operator of Eq. ( 12) on the electronic wavefunction. The GI part of the TDPES, GI , is invariant under the gauge transformation (9): ˜ GI (R, t) = GI (R, t); the GD part, on the other hand, transforms as

˜ GD (R, t) = GD (R, t) + ∂ t θ(R, t).
The analysis of the TDPES will be based, in this section, on a one-dimensional model for proton-coupled electron transfer [89], whose details are given in Refs. [10][11][12][13]. The model consists of two positively-charged ions fixed at a distance L, a positively-charged ion of mass M (in this case, the proton mass is used) moving in one dimension between the two fixed ions, thus interacting with the fixed ions via a bare Coulomb potential, an electron moving in one dimension and interacting with all ions via a soft-Coulomb potential. The schematic representation of the model is given in Fig. 3. Since we are dealing with a one-dimensional problem, we can choose a gauge such that the time-dependent vector potential is always zero. Inverting Eq. ( 4) to express Φ R (r, t) in terms of Ψ(r, R, t) and of χ(R, t), and inserting such expression in the definition of the vector potential (14), we obtain a relation between the vector potential itself, the nuclear velocity field, and the phase S(R, t) of the nuclear wavefunction

A ν (R, t) = Im Ψ(t)| ∇ ν Ψ(t) r |χ(R, t)| 2 -∇ ν S(R, t). ( 26
)
Adopting a one-dimensional representation,

A(R, t) = Im Ψ(t)| ∂ R Ψ(t) r / |χ(R, t)| 2 -∂ R S(R, t),
and imposing A(R, t) = 0, leads to an expression that defines the phase of the nuclear wavefunction,

S(R, t) = R dR Im Ψ(t)| ∂ R Ψ(t) r |χ(R , t)| 2 , ( 27
)
in this gauge. Since the time-dependent vector potential is identically zero all along the dynamics, only the TDPES appears in the nuclear TDSE (11). Therefore, the TDPES is the only potential affecting nuclear dynamics and carrying information about electronic dynamics. The dynamics is initiated by setting a Gaussian wavepacket, with zero average velocity, on the first excited adiabatic state S 1 , centered at R 0 = -4 bohr and with variance σ = 1/ √ 2.85 bohr. Fig. 3 shows the BO PESs corresponding to the ground S 0 and first-excited S 1 states. The time evolution of the molecular wavefunction is generated by solving the full TDSE by using the split-operator technique [90] with time step dt = 2.4 10 -3 fs (0.1 a.u.). Using Eq. ( 7) and fixing the gauge as in Eq. ( 27), the nuclear and electronic wavefunctions of Eq. ( 4) are determined. Once the electronic wavefunction is known, Eqs. ( 24) and ( 25) allow to determine the TDPES. A qualitative analysis of the TDPES and of its relation to the BO PESs is presented here, whereas a more quantitative discussion can be found in Refs. [10][11][12][13][14]. Fig. 3 shows different snapshots along the dynamics of the nuclear density (blue curves) and of the TDPES, decomposed in its GI (black lines) and GD (grey lines) components. The nuclear dynamics is initiated in the first excited state, corresponding to the green curve, in the negative R-region. Then, the nuclear wavepacket slides down towards the avoided crossing, where the energy gap between S 0 (corresponding to the red curve) and S 1 closes. While crossing the region of coupling, amplitude is transferred to the ground state, producing a two-component wavepacket evolving with different speeds depending on the shapes of the corresponding BO PES: in the positive R-region the red curve accelerates the wavepacket, whereas the green curve slows it down.

At early times (t = 4.84 fs), the dynamics is driven simply by the excited state, as the GI part of the TDPES completely lies on the green curve; later on (t = 14.52 fs), GI has a diabatic shape, connecting the adiabatic surfaces smoothly through the avoided crossing, thus allowing amplitude to "flow" from S 1 to S 0 ; finally (t = 24.20 fs and t = 31.46 fs), the GI component of the TDPES presents different adiabatic shapes in different portions of R-space connected by steps, clearly suggesting that the nuclear wavepacket is driven by different forces on the two sides of the steps. The spatial splitting of the nuclear density is, therefore, attributable to the appearance of the steps. However, one should keep in mind that also the GD component of the TDPES has an effect on the dynamics. As shown in Fig. 3, GD is either constant or piecewise constant, with steps appearing at the same positions of the steps in GI . Therefore, GD does not modify the slope of the GI component, it only reduces the size of the steps. These observations have been supported by an analytical justification in Refs. [10,13], which proves their general validity. According to Eq. ( 18), the nuclear density can be decomposed into adiabatic contributions, and this property is verified in Fig. 4. Initially at time t = 4.84 fs, the nuclear density has only the contribution corresponding to S 1 (green curve). When the wavepacket crosses the nonadiabatic coupling region at time t = 14.52 fs, the ground-state and excited-state contributions are localized in the same region. The spatial splitting appears at later times, t = 24.20 fs and t = 31.46 fs, and in Fig. 4 it is evident that the S 0 component (red curve) is localized in the same region where the GI part of the TDPES has the shape of the lower BO PES, while the S 1 component is localized on the other side of the step in the TDPES.

The coefficients C k (R, t) for k = 1, 2 (corresponding to S 0 and S 1 , respectively) of the expansion of the electronic wavefunction on the adiabatic basis, Eq. ( 19), have a similar behavior, that is analyzed in Fig. 5 by reporting their squared moduli. The quantity |C k (R, t)| 2 is determined from Eq. ( 23), thus it is only shown in regions where the nuclear density |χ| 2 is (numerically) larger than zero. In general, |C k (R, t)| 2 is equal to unity if the corresponding portion of the nuclear density is fully localized in one electronic state (see for instance the green curve at time t = 4.84 fs or the red curve at t = 24.20 fs for R > 2 bohr). The plots of the coefficients at the final times of the dynamics present a sigmoid shape, and cross when the nuclear wavepacket changes character, i.e., when it splits into different components evolving on different BO surfaces. This observation is fundamental for accounting for the

Φ R (t)| ĤBO |Φ R (t) r = k |C k (R, t)| 2 (k) BO (R).
This expression is simply a weighted average of the BO PESs, however if this average is computed at time t = 31.46 fs, it yields a form of the GI part of the TDPES that is far from a smooth potential, as the steps in |C k (R, t)| 2 translate to the steps in the time-dependent potential.

Time-dependent potential energy surface and time-dependent vector potential

The purpose of this section is to generalize to higher dimensions the observations on the TDPES and to analyze the time-dependent vector potential in a situation where it cannot be set to zero by any choice of gauge. To this end, we will be dealing with a two-dimensional model and a more general choice of gauge than in the previous section will be made. We choose the gauge such that the nuclear wavefunction χ(R, t) is real and non-negative, i.e. χ = |χ| ∀ R, t, achieved by imposing S(R, t) = 0. The time-dependent vector potential is a two-dimensional vector field depending on the nuclear coordinates R = X, Y and time t, and with this choice of gauge its expression is

A α (R, t) = Im Ψ(t)|∂ α |Ψ(t) r |χ(R, t)| 2 . ( 28
)
Here, α = X, Y , and Im Ψ(t)|∂ α |Ψ(t) r (divided by the nuclear mass) is the Cartesian αcomponent of the nuclear velocity field. The additional term of Eq. ( 26) depending on the gradient of the nuclear phase S(R, t) is identically zero in the chosen gauge. The details of the model Hamiltonian used in this section can be found in Ref. [32].

In the diabatic representation, this model describes the crossing of two similar parabolas, one being slightly displaced both in the X direction and in energy. The adiabatic PESs present a conical intersection at (3.0, 0.0) bohr. The initial nuclear wavepacket is taken as Gaussian, with widths σ X = 0.15 bohr and σ Y = 0.197 bohr, and is centered at R 0 = (X 0 , Y 0 ) = (2.0, 0.0) bohr. The dynamics is initiated in the second adiabatic state (S 1 ).

Fig. 6 shows snapshots of the nuclear wavepacket in the X, Y -plane. It is initialized on S 1 (red contour lines), slides down towards the conical intersection (red dot in the figure), transferring amplitudes to S 0 (grey contour lines), and finally separates into two wavepackets, one accelerating on the lower surface (grey contour lines) and one moving slower on the more confining upper surface (red contour lines). The dynamics is similar to the one-dimensional case studied above, but here it takes place in a two-dimensional space and it involves a funnelling process through the conical intersection. The analysis reported here will focus on two main points: first, we are interested in understanding if the TDPES shows peculiar features related to the conical intersection; second, and more importantly, we look for the appearance of any singular behavior of the time-dependent vector potential. Fig. 7 shows two views of the TDPES at times t = 24 fs (left panel) and t = 48 fs (right panel). At t = 24 fs only a cut of the TD-PES is reported along the Y -axis for X = 3 bohr, the position of the conical intersection. At this time, the nuclear wavepacket is localized around the conical intersection, and in fact we observe two contributions, denoted |χ BO |2 , corresponding to the lower BO PES S 0 and to the upper BO PES S 1 , respectively. In this analysis, the GI part of the TDPES has been decomposed into three contributions GI = GI1 + GI2 + GI3 , that are the three terms in the expression (24). The cut of GI1 lies between the adiabatic surfaces, and in this case does not present steps, in fact the two portions of the wavepacket are localized in the same region: the TD-PES resembles an average potential. Additionally, we observe a peculiar feature of |χ BO | 2 , namely peaks at the position of the conical intersection. At this position, there is an infinite coupling between the states S 0 and S 1 (the nonadiabatic coupling diverges here). However, the peaks exactly cancel out, yielding a very smooth [91] (total) nuclear density

|χ| 2 = |χ (1) BO | 2 + |χ (2)
BO | 2 . In Fig. 7, it appears that GI2 + GI3 exactly sum up to zero: this is not a general feature of GI2 and GI3 , however, we have observed in various situations that their sum is a smooth function of R. Such observation should be kept in mind, since those terms will be completely neglected [20,21,25,58] in the approximations introduced to construct the trajectory-based solution of the exact equations ( 10) and (11). Finally, with the current choice of gauge, the GD part of the TDPES is mostly constant and would lead to a rigid shift of the GI parts of the TDPES.

In Fig. 7 (right panel) a color-map plot of the TDPES at t = 48 fs is shown, 2 along with the contour lines representing the nuclear density (black) and its two adiabatic contributions. The steps observed in the one-dimensional case discussed previously, leading to a splitting of the nuclear wavepacket and explained -in a BO framework -by the S 1 nuclear component separating from the S 0 nuclear component, are clearly visible from the TDPES. The TDPES is composed of two main regions (3.7 < X < 4.6 bohr and 4.6 < X < 5.5 bohr), with a pronounced change of behavior at their interface. A sharp repulsive potential in the x direction composes the first region, supplemented by a central repulsive component at Y = 0 bohr. Conversely, the second region only shows a slowly increasing potential towards larger X value. This strong variation of the TDPES is responsible for the splitting of the nuclear wavepacket into two components.

In the chosen gauge, the time-dependent vector potential equals the nuclear velocity field. Therefore, it encodes information about how the nuclear wavepacket moves. The vector potential is shown in Fig. 8 at times t = 24 fs (left panel) and t = 48 fs (right panel). The color-code indicates the modulus of the vector field, whereas the arrows indicate the direction. For reference, the position of the nuclear wavepacket is represented as the black contour-lines.

The time-dependent vector potential exhibits a simple behavior when the wavepacket reaches the nonadiabatic region: it mainly points towards larger X, with a strength increasing along X, since the nuclear wavepacket itself, after being initiated on S 1 , moves rapidly towards the conical-intersection region without spreading significantly in the Y -direction. After the passage through the conical intersection (t = 48 fs), the time-dependent vector potential still mostly points towards larger X, but now with some variations in the region 3.5 < X < 4.5 bohr. In the region 4.0 < X < 4.5 bohr, the magnitude of A(R, t) is lower than at larger values of X, reflecting the change in behavior of the nuclear wavepacket: the nuclear wavepacket component on S 0 overtakes the one on S 1 at t = 48 fs, which feels the repulsive potential of S 1 . Therefore, the nuclear component evolving on S 1 , localized in the region 4.0 < X < 4.5 bohr, moves slower than the one on S 0 , that spreads along larger values of X. In the region 4.5 < X < 6.0 bohr, the magnitude of A(R, t) starts decreasing, reflecting how the wavepacket slows down when it enters a more repulsive part of the PES.

Computing the circulation of the vector potential along a closed path C,

γ ν (C) = C A ν (R, t) • dR ν = S curl A ν (R, t) • ds (29)
one finds that, in general, the circulation integral (here denoted γ ν (C)) is non-zero and depends on the path. The line integral along C can be replaced by a surface integral according to Stokes' theorem, yielding the flux of the vector field curl A ν (R, t) through S, with S the surface enclosed by C. If γ ν = 0, it means that curl A ν (R, t) = 0: the time-dependent vector potential A ν (R, t) has non-zero curl, therefore it cannot be written as the gradient of a scalar function (since the curl of a gradient is identically zero). Recalling the transformation expression of

A ν (R, t) under a change of gauge, Ãν (R, t) = A(R, t) + ∇ ν θ(R, t), it is clear that one can set Ãν (R, t) = 0 if A ν (R, t) = -∇ ν θ(R, t
). However, if in a general situation A ν cannot be written as the gradient of a scalar function θ because it has a non-zero curl, then in general the vector potential cannot be gauged away.

Trajectory-based quantum-classical dynamics

The coupled electronic (10) and nuclear (11) evolution equations derived from the exact factorization (4) of the molecular wavefunction have been solved in Refs. [20,21] based on the representation of the nuclear dynamics in terms of trajectories. Note that, as discussed in Ref. [23], the nuclear time-dependent Schr ödinger equation ( 11) can be exactly reformulated as two coupled equations, i.e., a Hamilton-Jacobi equation for the phase of the nuclear wavefunction χ and a continuity equation for the nuclear density |χ| 2 . The former, a partial differential equation that determines the phase of χ, can be solved by characteristics. When the coupling between the Hamilton-Jacobi and the continuity equations is neglected, an approximate solution of Eq. ( 11) is achieved. This procedure will be denoted "quantum-classical", as the electronic equation, solved at the quantum-mechanical level, is employed to determine the time-dependent potentials, that, in turn, are used to compute the force that generates classical, thus approximate, trajectories. The algorithm has been dubbed coupled-trajectory mixed quantum-classical (CT-MQC). CT-MQC has been employed to simulate nonadiabatic dynamics in model systems [20][21][22][23][24], and implemented in a developer version of the electronic-structure package CPMD [92], based on (time-dependent) density functional theory (TDDFT), to simulate the photo-induced dynamics in Oxirane [START_REF] Agostini | TDDFT and Quantum-Classical Dynamics: A Universal Tool Describing the Dynamics of Matter[END_REF]25,26]. The algorithm is also available as a module of the E-CAM Software library [93].

In what follows, the focus will be mainly put on two phenomena occurring in nonadiabatic events, the amplitude transfer between electronic states and quantum decoherence. The former is mainly induced by the nonadiabatic coupling vectors, thus occurring in regions of configuration space where the adiabatic states are close in energy. The latter [58, is usually defined as the decay of the off-diagonal elements of the reduced density matrix of the electronic subsystem, consequence of the coupling with the nuclear environment. Within the BO framework, the electronic subsystem is expected to collapse onto a selected adiabatic state after the dynamics has led the full system through a region of nonadiabatic coupling. It will be argued that decoherence is related to the spatial splitting, or delocalization, of the nuclear wavefunction.

CT-MQC: The approximations

Trajectories will be employed to mimic the dynamics of the nuclear wavefunction χ(R, t), that evolves according to the TDSE (11). The trajectories can be interpreted as a moving grid, thus, switching to the Lagrangian reference frame, we calculate time derivatives "along the flow": partial time derivatives are replaced by total derivatives, using the chain rule d/dt = ∂ t + ν V ν • ∇ ν . The velocity of the moving grid point, i.e., of the trajectory, is

V ν = P ν /M ν .
Writing the nuclear wavefunction in polar form, χ(R, t) = |χ(R, t)|e (i/ )S(R,t) , the real part of Eq. ( 11) yields

∂ t S(R, t) = - Nn ν=1 [∇ ν S(R, t) + A ν (R, t)] 2 2M ν -(R, t) - Nn ν=1 -2 2M ν ∇ 2 ν |χ(R, t)| |χ(R, t)| , (30)
a Hamilton-Jacobi equation in the presence of the time-dependent vector potential A ν (R, t), the TDPES (R, t), and of a potential term (last term in Eq. ( 30)) known in the framework of Bohmian dynamics as quantum potential [23]. The imaginary part of Eq. ( 11) yields a continuity equation for the nuclear density. Neglecting the quantum potential, Eq. ( 30) decouples from the continuity equation, and becomes a (standard) classical Hamilton-Jacobi equation. Then, the linear momentum is defined as M ν V ν = ∇ ν S(R, t) + A ν (R, t) = P ν , the full time derivative of S(R, t) is introduced, and the spatial derivative ∇ ν is applied on both sides. Eq. ( 30) reduces to Newton's equation,

Ṗν (t) R (I) (t) = -∇ ν (R, t) + Nn ν =1 A ν (R, t) • P ν (R, t) M ν + Ȧν (R, t) R (I) (t) , (31)
i.e., to one of the characteristic (ordinary differential) equations associated to the Hamilton-Jacobi equation (30) without the quantum potential (as shown in Ref. [23]). At time t all quantities are evaluated at R (I) (t). The gauge is imposed so as to set to zero the term in parenthesis, (R, t) + ν A ν (R, t) • P ν (R, t)/M ν = 0. Eq. ( 31) thus becomes Ṗν (t) = Ȧν (t).

The expression of the TDPES is simplified within CT-MQC by neglecting the terms GI2 and GI3 (as observed previously, their sum has a smoother R-dependence than GI1 ). It is worth noting that both terms have to be neglected in order to maintain gauge invariance. In order to allow for the interface of CT-MQC with standard quantum-chemistry codes, the BO representation of the electronic wavefunction is introduced, such that energies, forces and nonadiabatic coupling vectors appear in Eq. ( 10) and can be computed on the fly. The partial time derivative of the electronic wavefunction in Eq. ( 10) is replaced by the total time derivative,

∂ t Φ R (r, t) R (I) (t) = ΦR (r, t) -Nn ν=1 Pν (t) Mν • ∇ ν Φ R (r, t) R (I) (t)
. These two operations allow us, starting from Eq. ( 10), to derive a set of coupled partial differential equations for the coefficients C l (R, t) → C (I) l (t) -that become functions of the trajectory R (I) (t), and are thus indicated with a superscript (I). In the expression of the electron-nuclear coupling operator Û coup en of Eq. ( 12), a dependence on the nuclear wavefunction appears explicitly. Employing its polar representation,

-i ∇ ν χ (I) (t) χ (I) (t) + A (I) ν (t) = ∇ ν S (I) (t) + A (I) ν (t) + i -∇ ν χ (I) (t) χ (I) (t) = P (I) ν (t) + iQ (I) ν (t), (32)
we express such dependence in terms of P ν (t), the signature quantity of the exact factorization that we have dubbed "quantum momentum".

The only approximation introduced so far in the electronic equation (10) has been the neglect of GI2 and GI3 in the TDPES. Additionally, all terms that contain products of the quantum momentum and the nonadiabatic coupling vectors will neglected, and the spatial derivatives of the expansion coefficients will expressed as ∇ ν C (I)

l (t) = (i/ )f (I) l,ν (t)C (I) l (t) (with f (I) l,ν (t) = - t dt ∇ ν (l),(I)
BO ). The first approximation is justified if the nonadiabatic region is sufficiently localized; the second approximation is discussed in detail in Ref. [21] and it is of a more "phenomenological" nature, as it has been devised based on previous analysis [10][11][12][13]19] of the expansion coefficients. This second approximation is essential to transform the set of partial differential equations for the coefficients in a set of ordinary differential equations.

In conclusion, CT-MQC equations can be cast in a very simple form, that is

Ċk (t) = ĊEh. k (t) + Ċqm k (t), (33)
F ν = F Eh. ν + F qm ν , (34)
where the first terms on the right-hand sides are Ehrenfest-like terms (Eh.), whereas the second terms depend on the quantum momentum (qm). The Ehrenfest-like terms in Eqs. (33) and (34) are

Ċk (t) = - i (k) BO C k (t) - Nn ν=1 V ν • l d kl,ν C l (t) (35)
F Eh. ν = - k |C k (t)| 2 ∇ ν (k) BO - k,l C * l (t)C k (t) (k) BO - (l) BO d lk,ν . (36)
The additional terms, instead, are The adiabatic BO PESs are shown in all plots in black, and the positions of the trajectories moving the TDPES (the TDPES is approximated in CT-MQC as described in the main text) are represented as orange dots. The classical distributions are more localized than the quantum-mechanical, issue that can be corrected by including the quantum potential (last term in Eq. ( 30)) as discussed in Ref. [23].

Ċqm k (t) = Nn ν=1 Q ν M ν • f k,ν - l |C l (t)| 2 f l,ν C k (t) (37)
F qm ν = k |C k (t)| 2 Nn ν =1 2Q ν M ν • f k,ν f k,ν - l |C l (t)| 2 f l,ν . (38)
All R-dependent quantities have to be evaluated along a trajectory (I): Eqs. ( 36) and (38) generate the trajectory, that is coupled to n evolution equations (k = 1, . . . , n) representing the evolution of the electronic system according to Eqs. (35) and (37).

As clear from Eq. ( 32), the quantum momentum tracks the spatial variation of the nuclear density, as it contains its spatial derivative. At each time step of the simulated dynamics, the nuclear density has to be reconstructed, for instance by computing a histogram from the distribution of classical trajectories (plus smoothing, by using Gaussian functions localized at the positions of the trajectories [25]). Such calculation requires that at the end of each step of dynamics, the trajectories "communicate" -all at the same time -information about their positions, in order to compute the quantum momentum. Once Q (I) ν (t) is known, the trajectories can perform a new step of dynamics. On-the-fly calculation of the quantum momentum is possible only if the trajectories are propagated all at the same time, that is why the underlying algorithm has been dubbed "coupled-trajectory"-MQC. It has been shown [20-22, 24-26, 58] that inclusion of the quantum momentum is essential to reproduce quantum decoherence effects.

As an example of the importance of the spatial delocalization of a nuclear wavepacket, the nuclear density reconstructed from the distribution of CT-MQC trajectories for a two-state onedimensional model of extended nonadiabatic coupling with reflection [20-22, 24, 34] is shown here. The adiabatic PESs are plotted in black in Fig. 9. A Gaussian-shaped wavepacket centered at -15 bohr is prepared "on" the lower surface and launched with a positive initial momentum of 20 bohr -1 towards the region of nonadiabatic coupling between the two adiabatic states. When the coupling region is overcome, the lower BO PES has a well, whereas the upper BO PES has a barrier. For the chosen initial momentum, the contribution propagating on the lower state is transmitted, while the contribution propagating on the upper state is reflected. In Fig. 9 the first crossing of the nonadiabatic region of the incoming wavepacket (left panel), along with the splitting of the transmitted and reflected portions of the wavepacket (central and right panels) are represented at different time-steps. Here, the histograms constructed from the distribution of classical trajectories (blue lines) are compared with the solution of the full time-dependent Schr ödinger equation (red lines). The orange dots show the actual positions of the trajectories that move on the (CT-MQC-approximate form of the) TDPES.

CT-MQC: Photo-induced ring opening in Oxirane

Simulations of the photo-dynamics of Oxirane start at the time in which, upon absorbing a photon, the molecule is excited from its ground electronic state S 0 to the lowest-lying bright state, S 2 . The photo-excitation induces rearrangements of the molecular structure that drive the system through a conical intersection between S 2 and S 1 within 7 to 15 fs [25,26]. This process is illustrated in Fig. 10, which shows (in the upper panel) the average population of the electronic states, i.e., ρ k (t) = N -1 traj I |C (I) k (t)| 2 , with N traj = 100 the number of trajectories (k = 1, 2, 3 or, equivalently, S 0 , S 1 and S 2 ). When the molecule reaches the region of the conical intersection, population is transferred from S 2 (fully populated at time t = 0) to the first excited state S 1 . After that, the dynamics carries on until a conical intersection between S 1 and S 0 is eventually reached (this second event is not discussed here).

Additional information about the S 2 -to-S 1 transfer process can be extracted from the analysis of the decoherence indicator

η 12 (t) = N -1 traj I |C (I) 1 * (t)C (I) 2 (t)| 2 ,
presented in the lower panel of Fig. 10. The quantity η 12 (t) (black line in Fig. 10) shows two pronounced peaks, indication that two groups of trajectories funnel through the S 1 /S 2 conical intersection at subsequent times. In order to understand if the time delay between different passages through the conical intersection has an effect on the overall dynamics, the indicator of decoherence has been decomposed in different contributions, associated to different final molecular configurations. We observe four final product structures, (i) a right-open ring structure (with probability 36%), (ii) a left-open ring structure (with probability 47%), (iii) a CC-extended bond structure (with probability 10%), and (iv) a closed-ring structure (with probability 7%). All reported values are computed as the ratio between the number of trajectories ending up in the target configuration and the total number of trajectories (N traj = 100). Structures (i) and (ii) are equivalent, and we expect that the observed difference in percentage can be reduced by improving the statistics, thus by increasing the number of simulated trajectories. In these two cases, the ring-opening of Oxirane is obtained through the breaking of one of the two CO bonds. In structure (iii), the Oxirane ring opens via the elongation of the CC bond. A few trajectories, identified as structure (iv), are not reactive and remain close to the original molecular configuration. The decoherence indicator is decomposed into contributions associated to the "reactive trajectories". Therefore, we observe that the first peak (cyan and magenta curves in Fig. 10) between 6 and 12 fs is produced by trajectories that lead to the breakage of one of the two CO bonds. However, these curves do not decay monotonically. Instead, the curves corresponding to the C 1 O and C 2 O groups both contribute to the second peak (between 12 fs and 17.5 fs), indicating that the first group of trajectories (associated to the first peak) is reached by a second group while funnelling through the conical intersection.

The main contribution to the second peak (blue line in Fig. 10) between 12 and 16 fs is produced by trajectories yielding an extended CC bond. These trajectories clearly encounter the nonadiabatic region with some delay when compared to the sets of trajectories (i) and (ii).

The different reaction channels are clearly a consequence of the topology of the TDPES in configuration space. Therefore, we will now analyze the TDPES for the reactive trajectories of type (i) and (iii). It is worth recalling that, even though the adiabatic basis has been used to expand the electronic wavefunction of the exact factorization, the nuclear dynamics is still governed by the TDPES and by the time-dependent vector potential of Eq. (11), in their approximate quantum-classical form. The electronic adiabatic basis has been used merely for convenience.

Representative trajectories have been selected for the groups (i) (Fig. 11, left) and (iii) (Fig. 11, right). The populations of the electronic states and the adiabatic potential energy for each configuration visited along the trajectories are reported in Fig. 11 (upper and lower panels, respectively). The upper panels of Fig. 11 confirm that the region of strong coupling between states S 2 and S 1 is encountered by trajectories of type (i) at earlier times compared to trajectories of type (iii). In fact, the populations of the electronic states corresponding to group (i) sharply switch at the conical intersection at around 10 fs. This behavior is the consequence of the different shapes of the TDPES, represented as dotted lines in Fig. 11. The trajectories of group (iii) are driven by a TDPES that is initially flat for about 10-15 fs, following the shape of S 2 adiabatic state, until it smoothly approaches and then switches to S 1 . Later, these trajectories continue on S 1 without showing a clear tendency to approach the S 1 /S 0 conical intersection. By contrast, the TDPES sampled by the trajectories of group (i) follows a steeper path that brings to a fast closure of the S 2 /S 1 gap (within about 15 fs) and subsequently of the S 1 /S 0 gap (after about 25 fs), suggesting the presence of a funneling process that guides the trajectories to the ground state.

These calculations have been performed with the plane-waves based electronic structure package CPMD [92], employing the PBE [115] functional for ground-state and excited-state calculations. Linear-response TDDFT calculations [116][117][START_REF] Casida | Time-dependent density-functional response theory for molecules[END_REF] are based on the Tamm-Dancoff approximation [START_REF] Tamm | [END_REF]120]. The Kleinman-Bylander [121] pseudo-potential has been used for all atom species together with a plane-wave cutoff of 70 Ry. Initial conditions, i.e., positions and momenta, have been sampled from an ab initio ground-state trajectory of 2 ps at 300 K. The N traj = 100 trajectories are propagated with a time step of 0.12 fs (5 a.u.).

Perspectives Spin-orbit interactions in ultrafast molecular processes

Spin-orbit (SO) interactions are responsible for fascinating phenomena that often have direct implications for technological advances. SO coupling makes it possible to tune magnetic properties of materials with light [122,123], whereas the interaction between electronic charge and spin is exploited in the field of spintronics [124,125]; in organic light-emitting diodes SO effects allow to harness both singlet and triplet excitations for achieving high efficiency [126,127]. Assessing a priori the importance of SO coupling is difficult. Being a relativistic effect, it is essential in describing intersystem crossings (ISCs) in systems with heavy transition metals [128][129][130][131][132][133]. ISC denotes a spin-forbidden non-radiative electronic transition between states of different spin multiplicity, to be distinguished from the ubiquitous internal conversion (IC) that takes place between states of the same spin multiplicity. Nonetheless, ISCs have been observed in processes involving (light) organic species [134][135][136][137][138][139][140][141][142], thus proving that the strength of SO interactions depends on the shape of the electronic potential energy landscape and, consequently, on the molecular geometry. Therefore, a theory that is able of treating IC and ISC on equal footing, without an a priori knowledge of the relative importance of the two effects, is highly desirable, especially for the study of (photo-induced) ultrafast dynamics in complex molecular systems.

In this ongoing project, the reaction of ground-state oxygen O( 3 P) with ethylene C 2 H 4 , a prototype reaction [138,143] in atmospheric chemistry, will be studied based on the extension of CT-MQC to the treatment of ISC. In fact, from a fundamental perspective, O( 3 P)+C 2 H 4 is a paradigm reaction for excited-state molecular dynamics, as it involves coupling among electronic states of different spin multiplicity. Theoretically, the potential-energy landscape is complex, and in previous work [138] it has been reconstructed based on static high-level ab initio calculations via parametrized model potentials [143,144] in both singlet and triplet spin configurations. However, avoided crossings or conical intersections between states of the same spin multiplicity (IC) can be encountered along the dynamics, as well as ISCs. Furthermore, reaction pathways might become available only under the effect of the interplay between nuclear and electronic motion, thus the dynamical aspect of the problem is as important as the electronic-structure problem. The idea is, therefore, to perform on-the-fly ab initio calculations based on the extension of CT-MQC to treat SO interactions.

In general, the most common strategy to simulate ISC is to adapt molecular-dynamics schemes designed for IC: quantum wavepackets propagation techniques and trajectory-based approaches become, thus, readily available. Wavepacket propagation [134,139,145,146] requires to pre-compute electronic potentials as functions of nuclear configurations, whereas trajectories evolve under the effect of forces [147][148][149][150][151][152] determined based on on-the-fly electronicstructure calculations -and for this reason are better suited to access complex systems. Nowadays, a widely-used technique for the treatment of IC and ISC is trajectory surface hopping [153,154]. For IC, classical trajectories evolve "on" electronic potential energy surfaces, and hop from one surface to another to mimic an electronic transition preserving spin multiplicity. The extension of surface hopping to ISC has been proposed by different authors about a decade ago [148][149][150][151]. However, its use in this context remains still highly debated, concerning the choice of the most appropriate electronic representation to be used in the simulations [155], physically problematic, exhibiting issues in preserving rotational invariance of the coupling with states in the same multiplet [155,156], and subject to inconsistencies, in the use of different hopping rules for IC and ISC [148]. An intriguing and timely question thus arises, as to whether IC and ISC can be described within the same trajectory-based approach without incurring in the above-mentioned issues.

The aim of this project is to employ the exact factorization [6,[START_REF] Agostini | [END_REF] to treat in a consistent way, within the same theoretical construction IC and ISC. The exact factorization will be formulated accounting for relativistic SO interactions in the TDSE, and a numerical procedure for trajectory-based on-the-fly calculations will be derived to simulate ultrafast IC/ISC processes [157].

Ongoing collaborations: Sabine Morisset (Chargée de Recherche at the Institut de Sciences Moléculaires d'Orsay (ISMO) of the University Paris-Sud) and Nathalie Rougeau (Maître de Conférences at the ISMO) participate in this project, as well as David Lauvergnat (Directeur de Recherche at the LCP).

Quantum decoherence in light-induced charge transfer phenomena

Since 2007, after the first experimental observation of Fleming and coworkers [94], several studies have been reported on quantum coherent exciton transport in photosynthetic systems, even at room temperature and up to a few hundreds of femtoseconds. How is it possible that neither temperature nor the fluctuations of the environment are able to hinder the quantum mechanical nature of the process? Answering this question is key to understanding and controlling quantum coherence, with the aim to engineer devices that, similarly to photosynthetic systems, transform light energy into other forms of energy with high efficiency. For instance, what are the properties of materials able to harness light energy and convert it into electric current, but minimizing the energy loss onto the environment? In this context, π-conjugated donoracceptor complexes [103] are of central interest for the development of organic photovoltaic devices.

The purpose of this research project is to study quantum decoherence in light-induced charge transfer processes in π-conjugated molecular systems including the effect of temperature and of the environment. To address this problem, the idea is to develop theoretical and computational tools in the framework of the exact factorization to propose new, flexible algorithms for ab initio simulations. The starting point is CT-QMC, the algorithm derived from the exact-factorization equations. However, CT-MQC has been applied so far only to the study of small molecules in the gas phase, exploiting its implementation in CPMD, an electronicstructure package based on TDDFT. The limitations of TDDFT are several, as for instance it might incorrectly capture the excitation energy, thus yielding an incorrect order of the excited states [158], and it cannot properly capture conical intersections between ground and firstexcited state [26]. In addition, CPMD is a plane-wave-based code, property that limits its numerical performance (even though its excellent parallelization makes it suitable for HPC infrastructures). Therefore, it is clearly necessary to, first, circumvent these limitations, and, second, to adapt CT-MQC to the new challenge. These two goals will be achieved based on theoretical developments of the algorithm, on code developments, and on applications to systems of growing complexity.

Planned collaborations: Carine Clavaguéra (Chargée de Recherche at the LCP) and David Lauvergnat (Directeur de Recherche at the LCP) will contribute to this research project. An international collaboration with Ivano Tavernelli (IBM Research Laboratory in Zurich, Switzerland) is also envisaged, related to the new algorithms to be implemented in the code CPMD, maintained and released by IBM.

Dynamics at conical intersections

Since very recently we are witnessing a rebirth of the interest of the Physical Chemistry and Chemical Physics communities in conical intersections. Numerical simulations at different levels of theory [56], and experiments of different nature [START_REF] Yuan | Teaching activity Teaching activities • First-year Chemistry 2016-2019 -Chemistry introductory course: Exercises and Laboratory (∼ 110 hours)[END_REF] are being conducted to discover the effect of conical intersections on chemical reactions. The reason for this interest lies in the peculiar behavior of both electrons and nuclei at conical intersections, that can be observed in excited-state (nonadiabatic) processes and in ground-state (adiabatic) processes.

It has been pointed out [56] that geometric phases can have observable effects on nuclear quantum dynamics, for instance arising from destructive interferences of two wavepackets propagating around a conical intersection in opposite directions. The photodissociation of phenol constitutes a striking molecular example of such a process. In phenol, the H atom release upon light absorption involves a nonadiabatic tunneling process for which geometric phase effects seem critical to describe tunneling lifetimes adequately [56]. Also very recently the collision reaction H+HD→H 2 +D has been studied [START_REF] Yuan | Teaching activity Teaching activities • First-year Chemistry 2016-2019 -Chemistry introductory course: Exercises and Laboratory (∼ 110 hours)[END_REF]. A conical intersection between ground and first-excited state appears for all nuclear configurations yielding an equilateral H 2 D geometry. In the experiment the angular distribution of H 2 products has been measured, and it has been reproduced based on two sets of simulations. Firstly, simulations have been performed by forcing the electronic system to be in the ground state. However, the presence of the conical intersection, and thus of an infinite singular coupling with the first-excited state, clearly invalidates this approach, and a correction is needed to reproduce experimental results. The correction is the geometric phase. Secondly, ground and excited states have been explicitly considered. Experimental results have also been reproduced with this second approach, without invoking any effect related to the conical intersection.

In the context of the exact factorization it has been shown [31,32] for model systems that typical features related to conical intersections, e.g., the singularity of the potential energy surfaces and geometric phase effects, cannot be observed. A similar observation can be reported if one works in the (quasi)diabatic basis, defined as the representation where the nonadiabatic couplings are zero. Therefore, how can an experimental proof of observable effects of geometric phases be presented, as in the hydrogen-exchange reaction, if such effects appear to be rooted in the particular theoretical representation chosen to interpret experimental observations? Also, are there situations where typical features related to conical intersections are visible within the exact-factorization framework?

To address these questions, numerical studies based on quantum wavepackets propagation techniques will be carried out on various systems, and analyzed based on the exact factorization. Initially, two-and three-dimensional models will be investigated, followed by studies on molecular systems, such as pyrazine and the protonated Schiff base. The properties of the exact-factorization equations will be analyzed, to highlight possible singular behaviors related to geometric phases and conical intersections. Eventually, the possibility of solving the exactfactorization equations in a fully quantum-mechanical way will be investigated. The idea is that the absence of singularities in the exact factorization might simplify the design of a numerical scheme to address the problem of quantum dynamics at conical intersections.

Planned collaborations:

A collaboration is envisaged with Basile F. E. Curchod (Durham University, UK) as a continuation of the preliminary work on this topic initiated with Refs. [31,32].

Figure 1 :

 1 Figure 1: Colored areas define the maximum probability domains that have been computed around Na + . The separation between the two colored areas can be identified approximately as the first minimum of the radial distribution function of water surrounding the sodium ion.

Figure 2 :

 2 Figure 2: When oxirane is photo-excited from the ground to the second-excited electronic state, an ultrafast reorganization of the molecule is triggered, ending with the relaxation to the ground state within 35 fs.

Figure 3 :

 3 Figure 3: Left: Schematic representation of the model for proton-coupled electron transfer. Right: GI part (black) and GD part (grey) of the TDPES at times t = 4.84, 14.52, 24.20, 31.46 fs. The two lowest BO surfaces are shown for reference. The TDPES is shown in the regions where the nuclear density (blue lines in the panels) is larger than 10 -8 . The GD part of the TDPES is shifted in all panels to superimpose it to the plot of the GI part.

Figure 4 : 1 ) 2 )

 412 Figure 4: Nuclear density |χ(R, t)| 2 (blue lines), along with its BO contributions |χ (1) BO (R, t)| 2 (red lines) and |χ(2) BO (R, t)| 2 (green lines) at times t = 4.84, 14.52, 24.20, 31.46 fs.

Figure 5 :

 5 Figure 5: Nuclear density |χ(R, t)| 2 (blue lines), along with the squared moduli of the coefficients |C 1 (R, t)| 2 (red lines) and |C 2 (R, t)| 2 (green lines) of the expansion (19) at times t = 4.84, 14.52, 24.20, 31.46 fs.

Figure 6 :

 6 Figure 6: Schematic representation of the nuclear wavepacket at the beginning of the dynamics (t = 0), during the passage through the conical intersection (t = 24 fs), and after the nonadiabatic region (t = 48 fs). The colormap and the black contour lines show |χ(R, t)|, whereas the red contour lines represent |χ

( 2 )

 2 BO (R, t)| and the grey contour lines |χ

( 1 )

 1 BO (R, t)|. A red filled circle indicates the position of the conical intersection. Figure adapted from Ref. [32].

( 1 )

 1 BO | 2 and |χ

Figure 7 :

 7 Figure 7: Left panel: Representation of all the components forming the TDPES ( GI1 , GI2 , GI3 , and GD ) along a cut in the Y direction at X = 3 bohr at t = 24 fs. The total nuclear probability density (|χ| 2 ), adiabatic nuclear probability densities (|χ (1) BO | 2 and |χ (2) BO | 2 ), and adiabatic PESs (

  superimposed for comparison. Right panel: Sum of all the GI contributions to the TDPES at t = 48 fs (colormap). The black contour lines show |χ(R, t)|, whereas the red contour lines represents |χ

( 2 )

 2 BO (R, t)| and the grey contour lines |χ

( 1 )

 1 BO (R, t)|. Figures adapted from Ref.[32].

Figure 8 :

 8 Figure 8: Time-dependent vector potential at t = 24 fs (left panel) and at t = 48 fs (right panel). A(R, t) is represented by arrows, and the color-map represents its magnitude. The black contour-lines show |χ(R, t)| and the red filled circle the position of the conical intersection. The plot of A(R, t) is restricted to the region of space where |χ(R, t)| 2 > 10 -10 . Figures adapted from Ref. [32].

ν

  (t), the classical nuclear momentum, and Q

Figure 9 :

 9 Figure9: Snapshots at t = 30, 48, 57 fs of the nuclear density constructed from the distribution of CT-MQC-evolved trajectories (blue lines) for a model of extended nonadiabatic coupling with reflection. CT-MQC results are compared with exact dynamics (red lines). The adiabatic BO PESs are shown in all plots in black, and the positions of the trajectories moving the TDPES (the TDPES is approximated in CT-MQC as described in the main text) are represented as orange dots. The classical distributions are more localized than the quantum-mechanical, issue that can be corrected by including the quantum potential (last term in Eq. (30)) as discussed in Ref.[23].

2 Figure 10 :

 210 Figure 10: Upper panel: electronic populations of S 0 (orange), S 1 (green) and S 2 (red) as functions of time. Lower panel: (normalized) S 1 /S 2 indicator of decoherence (black line), and its decomposition in contributions from the three sets of reactive trajectories. The trajectory sets labeled with C 1 O (cyan line) and C 2 O (magenta line) lead to a final configuration where the Oxirane ring opens via the breaking of one of the two equivalent CO bonds; the set of trajectories labeled C 1 C 2 (blue line) yields final configurations where the ring opens through the elongation of the CC bond. The error-bars are computed as the standard deviations of the data. Figure adapted from Ref. [26].

Figure 11 :

 11 Figure11: Upper panels: populations of the electronic states S 0 , S 1 , and S 2 as functions of time for two selected trajectories of type (i) (left) and of type (iii) (right). The color code is the same used in Fig.10. Lower panels: energy profiles (in eV) along the selected trajectories. The zero is set to be the value of the energy of S 0 at time t = 0. In the upper panels, Oxirane at the final time is shown. Figures adapted from Ref.[26].

  2 (red lines) and |C 2 (R, t)| 2 (green lines) of the expansion (19) at times t = 4.84, 14.52, 24.20, 31.46 fs.

	spatial nuclear splitting,
	in fact, if expansion (19)
	is inserted in the first
	term on the right-hand
	side of Eq. (24) on gets

It is worth stressing here that the properties of the time-dependent vector potential and time-dependent potential energy surface in the presence of conical intersections are currently under investigation. Previous work has not identified peculiar features of the potentials that are related to conical intersections. However, the possibility of encountering such features in future analysis cannot be excluded.

We do not include GD (R, t) in this sum as it only contributes a nearly constant negative contribution to TDPES.

Research Grants conical intersection" -COST Action MP1306 EUSpec (Project leader: F. Agostini, host of B. F. E. Curchod).

HDR

Curriculum Vitae

• 10/2015 COST Short-Time Scientific Mission (STSM): Project "Charge transfer in photovoltaic processes beyond Ehrenfest dynamics" -2.5k euros (COST Action MP1306 EU-Spec).

• 2009 -2010 Support for a research project abroad: Project "Maximum probability domains for the analysis of the microscopic structure of liquids"-21.6k euros over 1 year (University Sapienza, Rome, Italy).

• 2007 -2009 Support for double PhD programs (co-tutelle de thèse between Italy and France) from the Franco-Italian University: Project "Study of the transfer process of an excess proton in water via molecular dynamics simulations" -5k euros over 2 years. evaluation of proposals for the organization of workshops and schools. DFG, Deutsche Forschung Gemeinschaft (German Science Foundation

Teaching activity

Supervision of students

• two PhD students: one in co-tutelle between the University Paris-Saclay, Orsay, France and the University of Pisa (Italy) -2019-2022; one at the University Martin-Luther of Halle-Wittenberg, Halle, Germany (supervisors: Prof. D. Sebastiani and Prof. R. Vuilleumier) -2014-2016.

• one postdoc at the University Paris-Sud, Orsay, France -11/2018-10/2020.

• two master students at the University Paris-Sud, Orsay, France -04-08/2018 (M1) and 02-07/2019 (M2).

• four bachelor students at the University Paris-Sud, Orsay, France -03-05/2018 (L2) and 11-12/2019 (L2).

Administrative responsibilities

• 2019 -2022: part of the scientific board "Ultrafast dynamics" of the LabEx PALM (Physique, Atomes, Lumière, Matière)

• 2019 -present: elected member of the LCP Council.

Outreach activities


E Lamballais
email: eric.lamballais@univ-poitiers.fr

T Dairay

S Laizet

J C Vassilicos

Implicit/explicit spectral viscosity and large-scale SGS effects

Introduction

Subgrid-scale (SGS) modelling based on regularization has become a popular approach for Large Eddy Simulation (LES). When the regularization is driven by the numerical error or by an extra discrete operator like a filter, it is usual to refer to implicit LES in the sense that the discretization provides an artificial dissipation interpreted as a substitute of SGS modelling. Typically, it is expected that this artificial dissipation is inactive at very large scales thanks to the numerical convergence of the associated discretization. This assumption of large-scale dynamics virtually free from any artificial dissipation can even be intentionally extended on a wide range of scales through an optimal design of the associated discrete schemes. This idea can also be recovered in explicit SGS models with for instance the concept of Spectral Vanishing Viscosity (SVV) [START_REF] Maday | Analysis of the spectral vanishing viscosity method for periodic conservation laws[END_REF] and the Variational Multiscale (VMS) methods [START_REF] Hughes | Large eddy simulation and the variational multiscale method[END_REF].

The goal of this study is to assess this inviscid assumption at very large scales for a flow at high Reynolds number while using DNS results to estimate the exact energy transfers from large to SGS. These transfers are investigated in the challenging situation of a flow subjected to a complete transition up to a fully developed turbulent state. The corresponding benchmark is the Taylor-Green vortex problem at Re = 20,000. In a previous work [START_REF] Dairay | Numerical dissipation vs. subgridscale modelling for large eddy simulation[END_REF], using an implicit SVV associated to the differentiation errors of the viscous term, we have shown that very accurate results can be obtained by LES at Re = 10,000 with a reduction of the number of degrees of freedom (DOF) of 8 3 by reference to DNS. In this study, we want to investigate the ability of this type of SGS modelling (without any direct influence on large scales) at higher Reynolds number and with a stronger reduction of DOF. Note that both Reynolds numbers Re = (10,000; 20,000) correspond to fully turbulent conditions as suggested by the value of their counterparts based on the Taylor microscale Re λ ≈ (200; 300) obtained after the complete turbulence breakdown at t ≈ 13.

A priori analysis from DNS results

The DNS of reference as well as the LES are performed using the sixth-order flow solver "Incompact3d" which is kinetic energy conserving in the discrete and inviscid sense (up to the time advancement error). For the present high Reynolds number case Re = 20,000, 3456 3 mesh nodes are required for a computational domain of (2π) 3 but using some symmetries of the problem, the number of DOF is actually divided by 8. Here, the goal is to carry out counterpart LES where the number of DOF and computational cost are reduced by 16 3 and 16 4 respectively leading to a cutoff wavenumber of k c = 108 for the LES mesh against 1728 for the DNS one.

Following our conclusions in a previous work, the targeted LES solution is defined using a progressive spatial filter as illustrated in figure 1-left where raw and filtered DNS energy spectra are compared. The filter is obtained by solving the Lin equation using a simplified spectral Pao-like closure while taking the implicit SGS dissipation into account [START_REF] Dairay | Numerical dissipation vs. subgridscale modelling for large eddy simulation[END_REF]. An important remark is that this filter is applied once in each spatial direction and not in all directions as it would be for an isotropic filter. This 1D definition of the filter is believed to be more significant by reference to the actual anisotropy of the LES mesh. Using this specific filter applied on the DNS data, the time evolution of the supergrid scale kinetic energy Ēk can be computed a priori with its associated total dissipation ε = -d Ēk /dt. Then, it is easy to estimate the viscous large-scale dissipation ε LS and its complementary SGS part ε SGS such as ε = ε LS + ε SGS . These dissipations as well as the full DNS dissipation are presented in figure 1-right where it can be seen that this benchmark is very challenging with SGS dissipation ε SGS up to 90% of the total dissipation ε, this unequal distribution giving a major role to the SGS model.

To have a more detailed view of the kinetic energy transfer, a scale by scale analysis of the SGS dissipation can be done starting from the large-scale Lin equation decomposed as

∂ ∂t + 2νk 2 Ē(k,t) = T (k,t) + T SGS (k,t) (1)
where Ē(k,t) is the kinetic energy of the filtered solution, T (k,t) the transfer term involving only the filtered solution (i.e. explicitly computed in LES) and T SGS (k,t) the remaining term that describes transfers between the supergrid and subgrid scales.

In this formalism expressed in the Fourier space, T SGS (k,t) is simply the spectral density of ε SGS that leads to the introduction of the spectral eddy viscosity Using the DNS data and following the procedure of [START_REF] Domaradzki | An analysis of subgridscale interactions in numerically simulated isotropic turbulence[END_REF] adapted in the present context, ν t (k,t) is estimated from 200 snapshots distributed throughout the calculation.

ν t (k,t) = - T SGS (k,t) 2k 2 Ē(k,t) (2)
Figure 2-left presents 5 samples of ν t (k,t) with a normalization based on the molecular viscosity ν. This figure clearly exhibits the dominant transfers close to k c meaning that the "hyperviscous feature" is observed for this benchmark, especially in the early transition. Figure 2-left also reveals that more distant triad interactions result in high values of ν t (k,t) at small k, not only during the transition (for instance at t = 10) but also until the end of the simulation where a fully developed nonequilibrium turbulence is observed. Then, at least qualitatively, the "plateau-cusp" profile of the spectral eddy viscosity, as predicted by two-point closure theories at high Reynolds number [START_REF] Kraichnan | Eddy viscosity in two and three dimensions[END_REF][START_REF] Chollet | Parameterization of small scales of the three-dimensional isotropic turbulence utilizing spectral closures[END_REF][START_REF] Lesieur | Turbulence in fluids[END_REF], is well recovered in the present a priori analysis.

In figure 2-right, the average value ν t <kc/4 of ν t (k,t) for 2 < k < k c /4 (as an estimation of the "plateau" value) and its cutoff value ν t =kc at k = k c are plotted throughout the simulation. The values of ν t =kc are found to be very high in the early transition but exhibit a global decrease as the turbulence develops. More importantly, the significance of direct effect of SGS on the large scale dynamics is confirmed during the turbulence breakdown (with ν t <kc/4 that can be more than 4 times larger than ν) but also when the turbulence is fully developed where ν t <kc/4 is still of the same order as ν. This observation is against the lack of any direct dissipative effect at large scales in the SGS modelling as it is assumed in implicit LES, SVV or VMS.

It could be thought that despite the high values of ν t (k,t) for k < k c /4, the corresponding fraction of SGS dissipation

ε SGS <kc/4 = 2 k c /4 0 ν t (k,t)k 2 Ē(k,t) dk (3)
is negligible. In figure 1-right, it can be seen that this quantity can actually be about 20% of the full SGS dissipation ε SGS during the turbulence breakdown while remaining about 10% until the end of the calculation. In the next section, it will be examined whether this significant contribution can be ignored in practical LES. /ν ν t =k c /ν "optimal" a posteriori ν t <k c /4 /ν Fig. 2 Left: spectral eddy viscosity ν t (k,t) at t = 5, 10, 15 and 20. Right: time evolution of the average ν t <kc/4 and cutoff ν t =kc values of ν t (k,t).

A posteriori analysis of LES results

The same flow configuration is investigated by LES (i) without any SGS modelling;

(ii) with the standard/dynamic Smagorinsky model and (iii) with our implicit SVV [START_REF] Dairay | Numerical dissipation vs. subgridscale modelling for large eddy simulation[END_REF]. Figure 3 presents the time evolution of the total dissipation ε obtained for the different LES. The very unrealistic behaviour observed for the no-model case confirms the major role of the SGS modelling for the present high Reynolds number case where the LES are based on a coarse mesh by comparison to DNS. The strong overestimation of ε without SGS model corresponds to an almost complete thermalization of the flow due to the development of small-scale spurious oscillations during the transition, as it can be clearly observed by visualization and spectral analysis (not shown for conciseness). The standard and dynamic Smagorinsky models are found to lead to a partial thermalization resulting in an overdissipative behaviour in the early transition. The resulting damping of small-scale spurious oscillations has, in a second stage, a feedback effect with an underestimation of ε. The use of implicit SVV prevents any thermalization with a good prediction of ε in the early transition. However, the excellent agreement obtained in [START_REF] Dairay | Numerical dissipation vs. subgridscale modelling for large eddy simulation[END_REF] at Re = 10,000 with a less coarse mesh is not recovered for the present more demanding benchmark. In particular, the peak of dissipation cannot be captured (underestimation of ε) and as a subsequent feedback effect, a spurious secondary peak (overestimation of ε) can be clearly observed. A similar spurious secondary peak can be observed for the enstrophy ζ (see figure 4-right).

For the implicit SVV, the main discrepancy can be attributed to the poor reproduction of the main peak of ε that has the potential to spoil the flow any time thereafter. Even a very strong increase of the implicit SVV near k c is unable to capture this peak (not shown for conciseness). A spectral analysis shows that the kinetic energy is overestimated at large scales during this particular moment, especially in the range 10 < k < k c /2, as illustrated in figure 4-left. This overestimation is interpreted as the consequence of the quasi-inviscid cascade at large scales. For this type of LES, free from distant triad interaction modelling, the overestimation of Ē(k,t) in the range 10 < k < k c /2 can be interpreted as a bottleneck effect. This interpretation is consistent with the a priori analysis presented in previous section. This view can even be supported quantitatively by observing that the main peak of dissipation is underestimated by about 22%. Because the implicit SVV is essentially inactive for k < k c /4, this behaviour may be related to the a priori ratio ε SGS <kc/4 /ε SGS ≈ 20% reported in the previous section.

To restore the peak, an idea could be to combine an explicit model with our implicit SVV with the hope that the former can boost the total dissipation during the transition while the latter can avoid the unrealistic partial thermalization. This kind of mixed model approach has been tried with the standard Smagorinsky (see figure 3-right) but without any improvement for the prediction of the main peak of ε. However, it is worth noting that this mixed model can remove the spurious secondary peak.

Another attempt was to modify our implicit model to allow non-vanishing spectral viscosity while adjusting its "plateau" value in order to follow precisely the time evolution of ε. Such LES can be considered as "optimal" in terms of ability to predict the filtered kinetic energy Ēk . The time evolution of the resulting plateau value is presented in figure 2-right. Even if significantly higher levels of eddy viscosity are required by comparison with the a priori estimation of ν t <kc/4 , it confirms the ability of a direct influence of the SGS model on very large scales to ensure the correct energy dissipation. However, this "optimal" approach (that is only a test case in the sense that it requires to know the expected time evolution of Ēk ) is found to underestimate Ē(k,t) at small scales (see figure 4-left) with a resulting strong underestimation of enstrophy (see figure 4-right).

Conclusion

In order to investigate the scale-selective influence of SGS on the large scale dynamics, DNS and LES are performed for the Taylor-Green vortex problem. An a priori analysis confirms the interest of the hyperviscous feature at small scale as used in implicit LES, SVV and VMS. However, the assumption of zero SGS dissipation at very large scales is found unrealistic for the high Reynolds number and coarse LES mesh considered. A posteriori analysis shows that SGS modelling based on the assumption of an inviscid cascade leads to a bottleneck effect on the kinetic energy spectrum with a significant underprediction of the total SGS dissipation. The simple addition of a constant eddy viscosity, even targeted to be optimal in terms of SGS dissipation, is unable to give realistic results. To allow accurate predictions by LES, a specific closure that incorporates both the hyperviscous feature (i.e. regularisation) and the expected SGS dissipation at large scales has to be developed.

Fig. 1

 1 Fig. 1 Left: raw and filtered DNS energy spectra. Right: time evolution of the viscous large-scale, SGS, SGS for k < k c /4, filtered DNS and full DNS dissipations (ε LS , ε SGS , ε SGS <kc/4 , ε, ε).

Fig. 3

 3 Fig. 3 Time evolution of the total dissipation ε predicted by LES.

Fig. 4

 4 Fig. 4 Left: energy spectra Ē(k,t) at t = 10. Right: time evolution of the enstrophy ζ .

Acknowledgements This work was granted access to the HPC resources of TGCC under the allocation 2014/2015/2016-2a0912 made by GENCI.