Marco Schirò

Florian G Eich

Federica Agostini
email: federica.agostini@universite-paris-saclay.fr

Federica Agostini Quantum

Marco Schir

Quantum-classical nonadiabatic dynamics of Floquet driven systems

published or not. The documents may come

I. INTRODUCTION

The development of theoretical approaches to treat light-matter interactions is nowadays a very active and productive field in various domains, from chemical physics to condensed-matter physics, from biophysics to material sciences. The most innovative studies are perhaps phenomena exhibiting strong coupling between the radiation field and matter, and requiring special frameworks even beyond the standard classical description of light [START_REF] Ebbesen | Hybrid light-matter states in a molecular and material science perspective[END_REF][START_REF] Schäfer | Ab initio nonrelativistic quantum electrodynamics: Bridging quantum chemistry and quantum optics from weak to strong coupling[END_REF][START_REF] Ruggenthaler | From a quantumelectrodynamical light-matter description to novel spectroscopies[END_REF][START_REF] Abedi | Shedding light on correlated electron-photon states using the exact factorization[END_REF][START_REF] Hoffmann | Light-matter interactions via the exact factorization approach[END_REF][6][START_REF] Lacombe | Exact potential energy surface for molecules in cavities[END_REF] . In other situations, for instance when the external field initiates and drives the highly non-equilibrium dynamics of molecular systems or solids, it is necessary to develop non-perturbative approaches to adequately describe light-matter interactions [START_REF] Sukharev | Optics of exciton-plasmon nanomaterials[END_REF][START_REF] Sato | Floquet states in dissipative open quantum systems[END_REF][START_REF] Giusti-Suzort | Dynamics of H + 2 in intense laser fields[END_REF][START_REF] Posthumus | The dynamics of small molecules in intense laser fields[END_REF][START_REF] Khosravi | Electronic non-adiabatic dynamics in enhanced ionization of isotopologues of hydrogen molecular ions from the exact factorization perspective[END_REF][START_REF] Khosravi | Exact potential driving the electron dynamics in enhanced ionization of H + 2[END_REF][START_REF] Suzuki | Electronic Schr ödinger equation with nonclassical nuclei[END_REF][START_REF] Lauvergnat | A simple and efficient evolution operator for time-dependent Hamiltonians: The Taylor expansion[END_REF] .

In photochemistry, the coupling of a molecular system with light is routinely used to steer or to drive chemical reactions [START_REF] Solá | Strong field laser control of photochemistry[END_REF] . Molecular dynamics simulations of photochemical reactions, however, not always consider the external field explicitly in the calculations, even though some examples have been reported [START_REF] Richter | SHARC: ab Initio molecular dynamics with surface hopping in the adiabatic representation including arbitrary couplings[END_REF][START_REF] Mignolet | Communication: XFAIMS -eXternal Field Ab Initio Multiple Spawning for electron-nuclear dynamics triggered by short laser pulses[END_REF][START_REF] Mignolet | Excited-state molecular dynamics triggered by light pulses -Ab initio multiple spawning vs trajectory surface hopping[END_REF][START_REF] Kelkensberg | A semi-classical model of attosec-ond electron localization in dissociative ionization of hydrogen[END_REF][START_REF] Marquetand | Nonadiabatic ab initio molecular dynamics including spin-orbit coupling and laser fields[END_REF] . In condensed matter physics an emerging frontier is the optical control of material properties [START_REF] Oka | Floquet engineering of quantum materials[END_REF] achieved by shining strong and frequency selective light fields on solid state materials [START_REF] Nicoletti | Nonlinear light matter interaction at terahertz frequencies[END_REF] and leading to a variety of light-induced collective phenomena, from transient superconductivity [START_REF] Fausti | Light-induced superconductivity in a stripe-ordered cuprate[END_REF] to non-trivial topological phases [START_REF] Mciver | Lightinduced anomalous hall effect in graphene[END_REF] .

In recent years, in order to simulate and interpret these phenomena, some interest has been devoted to the integration of Floquet theory for time-periodic Hamiltonians with techniques for excited-state, nonadiabatic dynamics. Floquet treatment of the periodicity induced on a quantum system by an external field effectively maps the time-dependent Schrodinger equation into an eigenvalue problem for Floquet states and their associated quasi-energies, analogously to the Bloch theorem for electrons in periodic potentials [START_REF] Floquet | Sur les équations différentielles linéaires à coefficients périodiques[END_REF][START_REF] Sambe | Steady states and quasienergies of a quantum-mechanical system in an oscillating field[END_REF][START_REF] Grifoni | Driven quantum tunneling[END_REF] . In the Floquet picture, time-dependent processes with periodic drive are analyzed in an extended vector space of "physical" electronic states and harmonics of the driving frequency. Therefore, electronic states can be interpreted as dressed states by the drive harmonics. The concept of static Floquet potential energy surfaces is, thus, introduced, which clearly generalizes the concept of static adiabatic, or Born-Oppenheimer, po-tential energy surfaces that are routinely used for interpreting field-free nonadiabatic processes [START_REF] Agostini | Different flavors of nonadiabatic molecular dynamics[END_REF] . Such a Floquet picture has been often employed in strong-field molecular physics, for interpretation purposes of bond-softening [START_REF] Bucksbaum | Softening of the H + 2 molecular bond in intense laser field[END_REF] and bond-hardening [START_REF] Giusti-Suzor | Vibrational trapping and suppression of dissociation in intense laser fields[END_REF] processes, photo-dissociation [START_REF] Giusti-Suzo | Above-threshold dissociatioionn of H + 2 in intense laser fields[END_REF][START_REF] Posthumus | Slow protons as a signature of zero-photon dissociation of H + 2 in intense laser fields[END_REF][START_REF] He | Laser-induced resonances in molecular dissociation in laser pulses[END_REF][35] , dynamical alignment [START_REF] Uhlmann | Molecular alignment of fragmenting H + 2 and H 2 in strong laser fields[END_REF] and anti-alignment [START_REF] Frasinski | Counterintuitive alignment of H + 2 in intense femtosecond laser fields[END_REF] . Recently, and due to the analogy between Born-Oppenheimer states/surfaces and Floquet dressed states/surfaces, the Floquet picture has been employed for actual simulations of driven time-dependent processes in molecules [START_REF] Fiedlschuster | Surface hopping in laser-driven molecular dynamics[END_REF][START_REF] Subramanian | Choice of the electronic basis for field-induced surface hopping[END_REF][START_REF] Zhou | Nonadiabatic dynamics in a laser field: Using Floquet fewest switches surface hopping to calculate electronic populations for slow nuclear velocities[END_REF][START_REF] Makhova | Floquet Hamiltonian for incorporating electronic excitation by a laser pulse into simulations of non-adiabatic dynamics[END_REF][START_REF] Fiedlschuster | Floquet surface hopping: Laser-driven dissociation and ionization dynamics of H + 2[END_REF][START_REF] Bajo | Mixed quantum-classical dynamics in the adiabatic representation to simulate molecules driven by strong laser pulses[END_REF] . In particular, trajectory-based schemes for excited-state dynamics, where the concept of electronic potential energy surfaces that guide nuclear dynamics is of utmost importance, lend themselves, naturally, to be combined with the Floquet picture.

In this paper, we focus on the combination of the exact factorization of the electronic-nuclear wavefunction [START_REF] Agostini | Different flavors of nonadiabatic molecular dynamics[END_REF][START_REF] Agostini | Exact factorization of the electron-nuclear wavefunction: Theory and applications[END_REF][START_REF] Abedi | Exact factorization of the time-dependent electronnuclear wave function[END_REF][START_REF] Agostini | TDDFT and quantum-classical dynamics: A universal tool describing the dynamics of matter[END_REF] with the Floquet formalism [START_REF] Fiedlschuster | Surface hopping in laser-driven molecular dynamics[END_REF][START_REF] Fiedlschuster | Laser-driven molecular dynamics: an exact factorization perspective[END_REF] , devoting particular attention to the coupledtrajectory mixed quantum-classical (CT-MQC) algorithm [START_REF] Min | Coupled-trajectory quantum-classical approach to electronic decoherence in nonadiabatic processes[END_REF][START_REF] Agostini | Quantum-classical non-adiabatic dynamics: Coupled-vs. independent-trajectory methods[END_REF][START_REF] Min | Ab initio nonadiabatic dynamics with coupled trajectories: A rigorous approach to quantum (de)coherence[END_REF][START_REF] Gossel | Coupled-trajectory mixed quantum-classical algorithm: A deconstruction[END_REF][START_REF] Curchod | CT-MQC -A coupled-trajectory mixed quantum/classical method including nonadiabatic quantum coherence effects[END_REF][START_REF] Agostini | An exact-factorization perspective on quantum-classical approaches to excitedstate dynamics[END_REF] . CT-MQC is the numerical scheme allowing to solve the exact-factorization equations based on the quantum-classical approximation [START_REF] Abedi | Mixed quantum-classical dynamics from the exact decomposition of electron-nuclear motion[END_REF][START_REF] Agostini | Classical nuclear motion coupled to electronic nonadiabatic transitions[END_REF][START_REF] Agostini | Nuclear quantum effects in electronic (non)adiabatic dynamics[END_REF][START_REF] Talotta | Quantum trajectories for the dynamics in the exact factorization framework: a proof-of-principle test[END_REF] of the nonadiabatic electron-nuclear problem, by introducing a trajectory-based solution of nuclear dynamics as formulated within the exact factorization. In the exact factorization, the nuclear wavefunction evolves according to a standard time-dependent Schr ödinger equation where the dynamic, fully nonadiabatic, effect of the electrons is represented by a time-dependent vector potential and a time-dependent potential energy surface (TDPES). It has been shown [START_REF] Fiedlschuster | Surface hopping in laser-driven molecular dynamics[END_REF][START_REF] Fiedlschuster | Laser-driven molecular dynamics: an exact factorization perspective[END_REF] in various situations where nonadiabatic effects are induced by an external time-dependent field, either laser pulse or continuous wave laser, that such TDPES very much resembles Floquet surfaces, rather than Born-Oppenheimer surfaces or quasi-static surfaces [START_REF] Kelkensberg | A semi-classical model of attosec-ond electron localization in dissociative ionization of hydrogen[END_REF][START_REF] Suzuki | Laser-induced electron localization in H + 2 : Mixed quantum-classical dynamics based on the exact time-dependent potential energy surface[END_REF] . Therefore, it seems natural to employ the Floquet representation of electronic dynamics driven by an external time-periodic field in combination with CT-MQC. Note that, despite the fact that CT-MQC has been derived from the exact-factorization electronic and nuclear equations, electronic dynamics is solved in a basis. The choice of the "most suitable" electronic basis is, thus, crucial. It is worth mentioning here that other trajectory-based approaches to excited-state dynamics have been combined with the Floquet formalism in order to explicitly include the effect of the photo-exciting or driving field, as the quantum-classical Liouville equation [START_REF] Horenko | A theoretical model for molecules interacting with intense laser pulses: The Floquet-based quantum-classical Liouville equation[END_REF] , ab initio multiple spawning [START_REF] Kim | Ab initio multiple spawning on laser-dressed states: a study of 1,3-cyclohexadiene photoisomerization via light-induced conical intersections[END_REF] , ab initio multiple cloning [START_REF] Makhova | Floquet Hamiltonian for incorporating electronic excitation by a laser pulse into simulations of non-adiabatic dynamics[END_REF] , and trajectory surface hopping [START_REF] Zhou | Nonadiabatic dynamics in a laser field: Using Floquet fewest switches surface hopping to calculate electronic populations for slow nuclear velocities[END_REF][START_REF] Fiedlschuster | Floquet surface hopping: Laser-driven dissociation and ionization dynamics of H + 2[END_REF][START_REF] Bajo | Mixed quantum-classical dynamics in the adiabatic representation to simulate molecules driven by strong laser pulses[END_REF][START_REF] Chen | On the proper derivation of the Floquet-based quantum classical Liouville equation and surface hopping describing a molecule or material subject to an external field[END_REF] .

The present work develops and applies CT-MQC using the Floquet formalism (F-CT-MQC).

To this end, in Section II, we briefly recall the exact factorization, introduce the elements of the Floquet theory necessary for F-CT-MQC, and derive F-CT-MQC equations. In Section III, F-CT-MQC is applied to simulate the periodically-driven, nonadiabatic dynamics of a model system, for which the exact solution is available, allowing us to test the performance of the algorithm. Our conclusions are presented in Section IV.

II. EXACT FACTORIZATION FOR PERIODICALLY DRIVEN SYSTEMS: THE FLOQUET

PICTURE

We study a system of interacting electrons and nuclei subject to an external time-dependent classical field V (r, R, t). The Hamiltonian describing the system is

Ĥ(r, R, t) = Tn (R) + Ĥel (r, R) + V (r, R, t) (1)
with electronic positions labelled by r and nuclear positions labelled by R. The field-free molecular Hamiltonian is the sum of the nuclear kinetic energy Tn (R) and of the electronic Hamiltonian Ĥel (r, R) which is the sum of the electronic kinetic energy and the interaction potential.

The solution of the time-dependent Schr ödinger equation (tdSE), Ψ(r, R, t), with Hamiltonian (1) can be factored as the product of a nuclear wavefunction, χ(R, t), and an electronic conditional factor, Φ(r, t; R), that parametrically depends on R, namely [START_REF] Agostini | Exact factorization of the electron-nuclear wavefunction: Theory and applications[END_REF][START_REF] Abedi | Exact factorization of the time-dependent electronnuclear wave function[END_REF] Ψ(r, R, t) = χ(R, t)Φ(r, t; R)

The evolution equations for χ(R, t) and Φ(r, t; R) are derived from the full tdSE [START_REF] Agostini | Nuclear quantum effects in electronic (non)adiabatic dynamics[END_REF][START_REF] Abedi | Correlated electron-nuclear dynamics: Exact factorization of the molecular wave-function[END_REF][START_REF] Alonso | Comment on "Correlated electron-nuclear dynamics: Exact factorization of the molecular wave-function[END_REF][START_REF] Abedi | Reply to comment on "Correlated electron-nuclear dynamics: Exact factorization of the molecular wave-function[END_REF] and are

i ∂ t χ(R, t) = 1 2 M -1 [-i ∇ + A(R, t)] 2 + (R, t) + ext (R, t) χ(R, t) (3)
i ∂ t Φ(r, t; R) = Ĥel (r, R) + V (r, R, t) + Û [Φ, χ] -(R, t) -ext (R, t) Φ(r, t; R) (4)
where the symbol M stands for the diagonal (constant) mass tensor and ∇ for the spatial derivative with respect to nuclear positions. The time-dependent potentials of the exact factorization mediate the coupling between electrons and nuclei, and are the time-dependent vector potential 65-68

A(R, t) = Φ(r, t; R)| -i ∇Φ(r, t; R) (5)
and the time-dependent scalar potential, or time-dependent potential energy surface (TDPES) [START_REF] Lacombe | Exact potential energy surface for molecules in cavities[END_REF][START_REF] Abedi | Dynamical steps that bridge piecewise adiabatic shapes in the exact time-dependent potential energy surface[END_REF][START_REF] Agostini | Mixed quantum-classical dynamics on the exact time-dependent potential energy surfaces: A novel perspective on non-adiabatic processes[END_REF][START_REF] Min | Is the molecular berry phase an artefact of the born-oppenheimer approximation?[END_REF][START_REF] Agostini | The exact forces on classical nuclei in non-adiabatic charge transfer[END_REF][START_REF] Curchod | An exact factorization perspective on quantum interferences in nonadiabatic dynamics[END_REF] , which we decompose into two contributions, namely

(R, t) = Φ(r, t; R)| Ĥel (r, R) + Û [Φ, χ] -i ∂ t |Φ(r, t; R) (6)
and

ext (R, t) = Φ(r, t; R)| V (r, R, t) |Φ(r, t; R) (7)
The electron-nuclear coupling operator Û

[Φ, χ] = Û [Φ(r, t; R), χ(R, t)] is 54,55,74 Û [Φ, χ] =M -1 1 2 [-i ∇ + A(R, t)] 2 + -i ∇χ(R, t) χ(R, t) + A(R, t) • -i ∇ -A(R, t) (8)
and depends explicitly on the nuclear wavfunction χ(R, t), and implicitly on the electronic factor Φ(r, t; R), via its dependence on the vector potential. The integration operation, indicated as • in previous equations, will be discussed below.

While the derivation just presented is valid in general situations, we focus here on the case of an external drive V (r, R, t), with constant amplitude and periodic in time with frequency Ω = 2π/T (continuous wave (cw) laser). In such case, we will generalize the trajectory-based approach developed to solve Eqs. ( 3) and ( 4), and dubbed coupled-trajectory mixed quantum-classical (CT-MQC) algorithm, to the Floquet formalism.

A. Elements of Floquet theory for CT-MQC

The implementation of CT-MQC relies on the expansion of the electronic wavefunction Φ(r, t; R) on an electronic basis. Since the system is subject to a time-periodic Hamiltonian it seems natural to expand on a basis that takes already into account, to some extent, the driven nature of the problem. As we discuss below this can be done using Floquet theory and will naturally lead to two different classes of states that can be used as electronic basis, the Floquet adiabatic and diabatic (or dressed) states, that we briefly introduce here.

The electronic time-dependent periodic Hamiltonian is Ĥel (r, R) + V (r, R, t), at fixed nuclear positions. As such, according to the Floquet theorem [START_REF] Floquet | Sur les équations différentielles linéaires à coefficients périodiques[END_REF][START_REF] Sambe | Steady states and quasienergies of a quantum-mechanical system in an oscillating field[END_REF][START_REF] Grifoni | Driven quantum tunneling[END_REF] , a complete set of solutions of the electronic tdSE, if there were no nonadiabatic coupling to the nuclei, takes the form e iEα(R)t φ α (r, t; R), where φ α (r, t; R) are the so called Floquet adiabatic states which are constructed as the eigenstates of the Floquet Hamiltonian

ĤF l (r, R, t) = Ĥel (r, R) + V (r, R, t) -i ∂ t (9)
i.e., they satisfy the eigenvalue problem

ĤF l (r, R, t)φ α (r, t; R) = E α (R)φ α (r, t; R) (10)
In the equation above the eigenvalues E α (R) are called Floquet quasi-energies and do not depend on time (but they depend on R as effect of the parametric dependence of the electronic Hamiltonian on the nuclear positions); φ α (r, t; R) = φ α (r, t+T ; R) has the same periodicity of the external drive. Given the periodicity of the Floquet states φ α (r, t; R), we can expand them in harmonics of the drive, as

φ α (r, t; R) = n=+∞ n=-∞ e iωnt φ α,n (r; R) (11)
with ω n = nΩ and n an integer. Inserting this expression into Floquet equation (10) and projecting onto the harmonic m we get

n Ĥel (r, R) + ω m Î δ mn + Vmn (r, R) φ α,n (r; R) = E α (R)φ α,m (r; R) (12)
where

Vmn (r, R) = 1 T T 0 dt e i(ωn-ωm)t V (r, R, t) . (13)
We recognize in Eq. ( 12) an eigenvalue problem in an extended space, including both the electronic degrees of freedoms as well as the harmonics of the drive. Note that, here, the scalar product in the space of the harmonics is defined via a time integral over a period of the drive and satisfies the orthonormality condition 1 T T 0 dt e i(ωn-ωm)t = δ nm

The Floquet adiabatic states obtained solving the eigenproblem (12) have a mixed electronic-field nature since Vmn (r, R) couples different harmonics, while the operator Ĥel (r, R) + ω m Î is diagonal in the space of harmonics. We can now expand the electronic wavefunction of the exact factorization in the basis of Floquet adiabatic states as

Φ(r, t; R) = α C α (R, t)φ α (r, t; R) = α n C α R, t e iωnt φ α,n (r; R) (14)
where the sum over α runs over the complete set of Floquet states, solution of Eq. (10), and the expansion coefficients C α (R, t) depend on the nuclear positions as well as on time. Note that, while the Floquet adiabatic states φ α (r, t; R) have periodicity T , this is in general not the case for the electronic wavefunction Φ(r, t; R), which is obtained from the full electron-nuclear wavefunction of the problem (for which one could also use Floquet theorem) through the exact factorization ansatz in Eq. ( 2). This observation justifies the use of time-dependent coefficients in the expansion in Eq. ( 14).

As we have seen so far, constructing the Floquet adiabatic states requires the solution of an eigenvalue problem in an extended space, which can be a burden for realistic applications beyond model systems. A different set of states, called Floquet diabatic states, can be obtained considering the electronic field-free Floquet Hamiltonian

Ĥ(0) F l (r, R, t) = Ĥel (r, R) -i ∂ t (15)
Following from Floquet theorem [START_REF] Sambe | Steady states and quasienergies of a quantum-mechanical system in an oscillating field[END_REF] , the eigenvalue equation

Ĥ(0) F l (r, R, t)ϕ α (r, t; R) = E (0) α (R)ϕ α (r, t; R) (16)
yields the Floquet eigenmodes ϕ α (r, t; R) and the Floquet quasi-energies

E (0) α (R).
As for the adiabatic modes, the diabatic eigenmodes are periodic in time ϕ α (r, t + T ; R) = ϕ α (r, t; R) with the same period T of the drive. However, as compared to the Floquet states introduced previously, we expect the diabatic states to have a much simpler structure in the space of harmonics, since by construction the time-dependent drive is not included in the Floquet Hamiltonian of Eq. (15). In fact, we can write explicitly the expression of these Floquet diabatic states, in terms of the standard Born-Oppenheimer basis {ψ k (r; R)} k=1,...,+∞ , defined as the eigenstates of the electronic Hamiltonian Ĥel (r, R) for each nuclear configuration R and satisfying the condition

Ĥel (r, R)ψ k (r; R) = E k (R)ψ k (r; R).
A direct calculation shows that the state

ϕ α (r, t; R) = e iωnt ψ k (r; R), α = (k, n) (17)
is eigenstate of the drive-free Floquet Hamiltonian (15) with quasi-energy

E (0) α (R) = E k (R)+ ω n .
The Floquet diabatic potential energy surfaces (PESs) are simply the standard Born-Oppenheimer PESs rigidly shifted by the energy of the corresponding harmonic.

Similarly to Eq. ( 14), we can choose the Floquet diabatic states as electronic representation for the electronic wavefunction Φ(r, t; R), and write

Φ(r, t; R) = k n C k,n R, t e iωnt ψ k r; R (18)
where the two sums run over the electronic physical states k and over the harmonics n. As observed previously, we allow an explicit time-dependence for the coefficients C k,n R, t in order to reproduce the dynamics of the electronic wavefunction which does not need to be periodic in time with the period of the drive.

While a priori the Floquet adiabatic and diabatic bases could be used in trajectory-based calculations, and one might expect the Floquet adiabatic basis to perform better since it contains already information on the external drive, we will argue that for the practical implementation of F-CT-MQC, the Floquet diabatic basis is preferred. For this reason in the next section we will derive F-CT-MQC equations only for this basis, although the derivation for the Floquet adiabatic basis can be obtained via an obvious generalization.

B. F-CT-MQC with Floquet diabatic states

We introduce the trajectory-based representation of nuclear dynamics [START_REF] Talotta | Internal conversion and intersystem crossing with the exact factorization[END_REF] that follows from the exact-factorization equations ( 3) and ( 4). The nuclear wavefunction is written in polar represen-

tation as χ(R, t) = |χ(R, t)| exp [(i/ )S(R, t)],
such that the coupled evolution equations

-∂ t S(R, t) = M -1 [∇S(R, t) + A(R, t)] 2 2 + (R, t) + ext (R, t) (19)
∂ t |χ(R, t)| 2 = -∇ • M -1 (∇S(R, t) + A(R, t)) |χ(R, t)| 2 (20)
are derived from the nuclear time-dependent Schr ödinger equation (3). The equation for the phase S(R, t) is given in the classical limit by neglecting in Eq. ( 19) the quantum potential term.

The first equation can be solved with the method of characteristics, introducing a set of ordinary differential equations -the characteristic equations -yielding the values of the field S(R, t) ∀ R, t.

The characteristic equations are Hamilton-like equations that describe the evolution in time of the "variables" R(t) and P(t) ≡ ∇S(R(t), t) + A(R(t), t) appearing in Eq. ( 19), i.e.,

Ṙ(t) = M -1 P(t) (21)
Ṗ(t) = -∇ R(t), t + ext R(t), t + Ṙ(t) • A R(t), t + Ȧ R(t), t (22)
The procedure to derive Eqs. ( 21) and ( 22) has been illustrated in detail in Refs. [START_REF] Agostini | Nuclear quantum effects in electronic (non)adiabatic dynamics[END_REF][START_REF] Talotta | Quantum trajectories for the dynamics in the exact factorization framework: a proof-of-principle test[END_REF] , therefore, we just mention here that the first term on the right-hand side of Eq. ( 22) is the gradient of the (pseudo-classical) Hamiltonian defined on the right-hand side of Eq. ( 19). The time-dependent potentials R(t), t , ext R(t), t and A R(t), t are evaluated along each characteristic R(t) by solving the electronic equation (4) along the same characteristic. Solving the partial differential equation (4) along the flow of trajectories -the characteristics -requires to switch from the Eulerian frame to the Lagrangian frame, but before discussing the electronic equation, let us discuss the continuity equation (20).

The evolution of the nuclear density, Eq. ( 20), is described by a "standard" continuity equation, that can be solved coupled to Eqs. ( 21) and ( 22). This would allow us to reconstruct the nuclear wavefunction, the only approximation being that we dropped the quantum potential in Eq. ( 19).

Neglecting the quantum potential has the effect of decoupling the evolution of the phase from the evolution of the density, while the continuity equation still depends on the phase. Therefore, as done in previous work [START_REF] Agostini | Quantum-classical non-adiabatic dynamics: Coupled-vs. independent-trajectory methods[END_REF][START_REF] Agostini | Nuclear quantum effects in electronic (non)adiabatic dynamics[END_REF] , we will not solve the continuity equation, and we will only reconstruct a classical-like nuclear density from the distribution of the trajectories, working in the hypothesis that for short enough times, an ensemble of trajectories evolving with Eqs. ( 21) and ( 22) will sample portions of nuclear configuration space with high probability density.

Equations ( 21) and ( 22) are the basis of the CT-MQC algorithm, that has been derived and tested in Refs. [START_REF] Min | Coupled-trajectory quantum-classical approach to electronic decoherence in nonadiabatic processes[END_REF][START_REF] Agostini | Quantum-classical non-adiabatic dynamics: Coupled-vs. independent-trajectory methods[END_REF] , applied to study the relaxation dynamics through conical intersections in

Refs. [START_REF] Marsili | Quantum and quantum-classical studies of the photoisomerization of a retinal chromophore model[END_REF][START_REF] Pieroni | Relaxation dynamics through a conical intersection: Quantum and quantum-classical studies[END_REF] , and in Refs. [START_REF] Min | Ab initio nonadiabatic dynamics with coupled trajectories: A rigorous approach to quantum (de)coherence[END_REF][START_REF] Curchod | CT-MQC -A coupled-trajectory mixed quantum/classical method including nonadiabatic quantum coherence effects[END_REF] in combination with time-dependent density functional theory to study the photo-induced dynamics in oxirane, and thoroughly analyzed in Refs. [START_REF] Gossel | Coupled-trajectory mixed quantum-classical algorithm: A deconstruction[END_REF][START_REF] Agostini | An exact-factorization perspective on quantum-classical approaches to excitedstate dynamics[END_REF] .

The exact-factorization product form of the solution of the tdSE is invariant under a (R, t)dependent phase transformation of the nuclear and electronic wavefunctions. Therefore, the gauge freedom has to be fixed by choosing the gauge. As previously done, in CT-MQC the gauge

is chosen such that R(t), t + Ṙ(t) • A R(t), t = 0.
With this choice of gauge, the characteristic equation ( 22) for the nuclear momentum yields

Ṗ(t) = -∇ ext R(t), t + Ȧ R(t), t (23)
which gives the expression of the classical force Ṗ(t) = F(t) that drives F-CT-MQC trajectories.

The explicit form of each term will be given below, after determining the evolution equation for the electronic wavefunction.

The electronic equation ( 4) is expressed along the characteristics R(t) as

i Φ(r, t; R(t)) = Ĥel r, R(t) + V r, R(t), t + Û [Φ, χ] -R(t), t -ext R(t), t +i Ṙ(t) • ∇ Φ(r, t; R(t)) (24)
Switching from the Eulerian to the Lagrangian frame, only total time derivatives can be evaluated along the flow, that is why the symbol Φ(r, t; R(t)) has been introduced. Furthermore, we

used the relation Φ(r, t; R(t)) = ∂ t Φ(r, t; R(t)) + Ṙ(t) • ∇Φ(r, t; R(t))
to replace the partial time derivative of Eq. ( 4) with the total time derivative, and thus obtain the last term in square brackets on the right-hand side of Eq. ( 24). The expansion of the electronic wavefunction in the Floquet diabatic basis is expressed as well along a trajectory R(t)

Φ(r, t; R(t)) = k n C k,n R(t), t e iωnt ψ k r; R(t) (25)
and is inserted into Eq. ( 24). When using it in the left-hand side of Eq. ( 24), we have

Φ(r, t; R(t)) = k,n Ċk,n R(t), t ψ k (r; R(t)) + C k,n R(t), t ψk (r; R(t)) +iω n C k,n R(t), t ψ k (r; R(t)) e iωnt (26)
Applying the total time derivative to the electronic state in the second term on the right-hand side, we find that the partial time derivative is zero because ψ k (r; R(t)) depends on time only via its dependence on the trajectory, then

Φ(r, t; R(t)) = k,n Ċk,n R(t), t + C k,n R(t), t Ṙ(t) • ∇ + iω n C k,n R(t), t ψ k (r; R(t))e iωnt (27)
With the aim to derive the evolution equations for C k,n R(t), t , the expansion in Floquet diabatic states is used as well in the right-hand side of Eq. ( 24). The whole procedure is detailed in Appendix A, and we discuss here only the main differences if compared to the standard derivation of CT-MQC equations.

In order to isolate in Eq. ( 24) a term Ċl,m R(t), t , it is necessary to project onto the state e iωmt ψ l r; R(t) . Working in the Floquet basis, the projection operation involves an integral over time, as well as an integral over the electronic configuration space. To do this, we introduce the change of variable t → s in e iωnt (and e iωmt ). This is done to separate the time dependence in the electronic equation into a periodic contribution induced by the drive, and indicated with s, and a non-periodic contribution, indicated with t. The time s is only used when the integral over the period is computed, in order to distinguish it from the non-periodic time dependence t. After the projection operation, the dependence on s disappears because it has been integrated out over a period, but the dependence on t remains. The projection operation onto the state e iωms ψ l r; R(t) is performed in the extended harmonic-electronic space, with the time integration over a period and the spatial integration over the whole electronic configuration space. Therefore, the instantaneous expectation value on the electronic wavefunction Φ(r, t; R(t)) of a general electron-nuclear operator Ô(r, R(t)) that does not depend explicitly on time is computed as

Φ(r, t; R(t)) Ô r, R(t) Φ(r, t; R(t)) = 1 T T 0 ds dr k,n l,m Ck,n R(t), t e -iωns ψk (r; R(t)) Ô r, R(t) C l,m R(t)
, t e iωms ψ l (r; R(t)) (28) showing that the periodic time dependence carried by the harmonics and the variable s is integrated out, whereas the time dependence in t remains. This operation is applied in the definition of the TDPES and of the time-dependent vector potential. Note that, this separation of times is crucial to keep working in a time-dependent picture while using the Floquet formalism. As discussed in the Introduction, the Floquet picture maps a periodic time-dependent problem into an effective time-independent problem in an enlarged harmonic-electronic vector space, thus reducing the solution of the time evolution to a diagonalization in this extended space. When separating the -full -molecular time-periodic process into a coupled electron-nuclear process to apply trajectory-based quantum-classical schemes, the time periodicity is only partially accounted for in the definition of the Floquet diabatic states, because, in general, the electronic and the nuclear dynamics are not time periodic (only the full molecular wavefunction is). It should be mentioned that the benefit of this procedure has been identified in previous work on the combination of the Floquet formalism with the quantum-classical Liouville equation [START_REF] Horenko | A theoretical model for molecules interacting with intense laser pulses: The Floquet-based quantum-classical Liouville equation[END_REF] . The identification of two time scales allows one to account for the oscillatory time dependence induced by the drive, while the remaining part of the time evolution is explicitly accounted for via the evolution of the expansion coefficients (as in Eq. ( 25)). Essentially, the operation described in Eq. ( 28) allows us to obtain time-dependent quantities after their dependence on the driving frequency of the laser has been averaged out, when considering their effect on the nuclei; the effect of the driving frequency on the electrons is treated in Fourier space, as it is treated in terms of Floquet diabatic states and transitions among them. A similar separation of time scales has been used in Floquet surface hopping [START_REF] Fiedlschuster | Floquet surface hopping: Laser-driven dissociation and ionization dynamics of H + 2[END_REF] with the aim to extend the Floquet picture to processes that are not strictly periodic, as for instance nonadiabatic processes initiated by a laser pulse.

Henceforth, the dependence on the trajectories R(t) will be replaced by a superscript ν, which stands for a trajectory index -this is done mainly to simplify the notation used in the equations below. The solution of Eq. ( 19) with the method of characteristics requires to solve Eqs. (21) and (22) for a large number of initial conditions, ν = 1, . . . , N tr . The index ν is, thus, used to label the trajectories, or characteristics. Similarly to the original form of the electronic evolution equation of CT-MQC [START_REF] Agostini | Quantum-classical non-adiabatic dynamics: Coupled-vs. independent-trajectory methods[END_REF][START_REF] Gossel | Coupled-trajectory mixed quantum-classical algorithm: A deconstruction[END_REF][START_REF] Agostini | An exact-factorization perspective on quantum-classical approaches to excitedstate dynamics[END_REF][START_REF] Talotta | Internal conversion and intersystem crossing with the exact factorization[END_REF] , the final expression of the time derivative of the Floquet diabatic coefficients, derived in Appendix A, is

Ċν l,m (t) = Ċν l,m (t) EH + Ċν l,m (t) QM + Ċν l,m (t) EXT (29)
which has to be solved along a trajectory ν. The three contributions identified in Eq. ( 29) are: an Ehrenfest-like term (EH)

Ċν l,m (t) EH = - i (E ν l + ω m ) C ν l,m (t) -Ṙν (t) • k,n d ν lk C ν k,n (t)δ mn (30)
a term depending on the quantum momentum (QM)

Ċν l,m (t) QM = 1 M -1 P ν (t) • f ν l,m -A ν (t) C ν l,m (t) (31)
and a term with explicit dependence on the external field

Ċν l,m (t) EXT = - i ν ext (t)C ν l,m (t) - k,n i V ν mn,lk C ν k,n (t) (32)
Note that, the only difference between the expression given in ( 29) and the field-free CT-MQC electronic equation resides in the additional terms depending on the external field (EXT), via ν ext (t) and V ν mn,lk given in Eq. (32). The contribution ν ext (t) to the TDPES in the Floquet diabatic basis reads

ν ext (t) = k,l n,m Cν k,n (t)C ν l,m (t)V ν mn,lk (33)
where the matrix elements of the driving field are

V ν mn,lk = 1 T T 0 ds dr e i(ωn-ωm)s ψν l (r) V ν (r, s)ψ ν k (r) (34)
As done above, all quantities depending on R(t) have acquired a superscript ν. The "periodic" time dependence has been indicated with the symbol s, which is integrated over a period in the projection operation. Combining Eqs. ( 33) and (34) shows that the expression of ν ext (t) easily follows from its definition in Eq. ( 7). Since only the periodic time dependence is integrated out, ν ext (t) still depends on time t. The quantity indicated as V ν mn,lk stands for the Fourier transform of the matrix elements of the external drive, from Eq. (34). The term V ν mn,lk induces population transfer along the trajectory ν between different electronic states and different harmonics as effect of the action of the external field. In particular, if only one driving frequency is considered, as in the present case, the non-zero matrix elements of V ν mn,lk are those satisfying the condition m -n = ±1.

In Eqs. (30) and (31), the "standard" terms of the electronic evolution equation of F-CT-MQC contain: the time-dependent vector potential A ν (t), the Floquet diabatic force f ν l,m accumulated along the trajectory ν, the quantum momentum P ν (t), the classical velocity of the trajectory Ṙν (t), and the nonadiabatic coupling vectors d ν lk . In Eq. ( 31), the time-dependent vector potential is

A ν (t) = 1 T T 0 ds dr k,l n,m Cν k,n (t)e -iωns ψν k (r) (-i ∇) C ν l,m (t)e iωms ψ ν l (r) (35)
following its definition in Eq. ( 5). Also in this case, the average operation is performed over the periodic time dependence, such that A ν (t) still depends on time t. The quantity indicated with the symbol f ν l,m in Eq. ( 31) is the force from the Floquet diabatic PES E ν l + ω m accumulated over time up to time t along the trajectory ν. Since the Floquet diabatic PESs are parallel to the adiabatic PESs, such an accumulated force can be computed from the gradient of E ν l as in standard CT-MQC. Such an accumulated force is used in the F-CT-MQC-approximate expression of the time-dependent vector potential, namely

A ν (t) = l,m C ν l,m (t) 2 f ν l,m (36)
which is related to the gradient of the harmonic-electronic coefficients C ν l,m (t) of Eq. (35) (see for details Refs. [START_REF] Agostini | Quantum-classical non-adiabatic dynamics: Coupled-vs. independent-trajectory methods[END_REF][START_REF] Talotta | Internal conversion and intersystem crossing with the exact factorization[END_REF][START_REF] Requist | Exact factorization-based density functional theory of electrons and nuclei[END_REF] ). In Eqs. (30) and (31), the nuclear velocity Ṙν (t) and the quantum momentum P ν (t) appear, via the contribution of the electron-nuclear coupling operator (8) that depends on the nuclear wavefunction (as discussed in Appendix A). The nuclear velocity in Eq. ( 30) couples to the nonadiabatic coupling vectors d ν lk , also known as derivative couplings, driving the population transfer between the electronic states k and l that belong to the same harmonic (δ mn ), as effect of the nuclear displacement operator ∇. They can be evaluated numerically or analytically via the Hellmann-Feynman formula, exactly as it is done usually in the the field-free case.

To finalize F-CT-MQC equations using the Floquet diabatic states, we give the explicit expression of the classical force of Eq. ( 23) for the trajectory ν. In particular, we can compute the gradient of ν ext (t) from Eq. ( 33), and the time derivative of the vector potential. Therefore, as in Eq. ( 29), we identify three contributions to the classical force,

F ν (t) = F EH ν (t) + F QM ν (t) + F EXT ν (t) (37)
where 'EH' and 'QM' indicate the -standard -Ehrenfest-like [START_REF] Agostini | Electron-nuclear entanglement in the timedependent molecular wavefunction[END_REF] and quantum-momentum terms, whereas 'EXT' stands for the term depending on the external drive. The Ehrenfest-like contribu-tion is

F EH ν (t) = k,n (-∇E ν k ) C ν k,n (t) 2 + k,n l,m Cν l,m (t)C ν k,n (t) (E ν l -E ν k ) d ν lk δ mn (38)
the quantum-momentum contribution, that couples the trajectories, is

F QM ν (t) = 2 k,n C ν k,n (t) 2 M -1 P ν (t) • f ν k,n f ν k,n -A ν (t) (39)
and the part that depends on the external field is

F EXT ν (t) = k,n l,m -∇V ν mn,lk Cν l,m (t)C ν k,n (t) + 1 k,n l,m Im Cν l,m (t)C ν k,n (t) V ν mn,lk f ν k,n -f ν l,m (40)
Note that the last term in the definition of F EXT ν (t) is formally similar to the contribution due to spin-orbit coupling recently introduced in G-CT-MQC [START_REF] Talotta | Internal conversion and intersystem crossing with the exact factorization[END_REF][START_REF] Talotta | Spin-orbit interactions in ultrafast molecular processes[END_REF] .

Before concluding this section, it is worth mentioning the partial normalization condition on the electronic wavefunction that is used to derive Eqs. ( 3) and ( 4). The partial normalization condition still holds, and using the expansion in harmonic-electronic eigenstates of the field-free Floquet Hamiltonian it reads

1 = 1 T T 0 ds dr k,n l,m
Ck,n (R, t)C l,m (R, t)e i(ωm-ωn)s ψk (r; R)ψ l (r; R)

(41a) = k,n |C k,n (R, t)| 2 ∀ R, t. (41b)
Equations ( 21), ( 29) and ( 37) define the F-CT-MQC algorithm in the Floquet representation, that will be used in the following section to study the dynamics of a two-electronic-state onedimensional model system subject to a cw laser of different intensities.

III. NUMERICAL STUDIES

The one-dimensional model employed for the numerical studies using F-CT-MQC is defined by the field-free electronic Hamiltonian

Ĥel =   H 11 (R) H 12 (R) H 12 (R) H 22 (R)   =   1 2 K(R -R 1 ) 2 γe -α(R-R 3 ) 2 γe -α(R-R 3 ) 2 1 2 K(R -R 2 ) 2 + ∆   (42)
given in the diabatic basis and depending on one nuclear coordinate R. The parameters are chosen as: K = 0.02 au, ∆ = 0.01 au, γ = 0.01 au, α = 3.0 au, R 1 = 6.0 au, R 2 = 2.0 au, R 3 = 3.875 au.

The diabatic PESs are parabolas, one with minimum in R 1 and the other with minimum in R 2 < R 1 ; the parabolas are, thus, displaced in R, and in energy as well by the amount ∆; the diabatic coupling has a Gaussian shape, with maximum in R 3 , and strength regulated by the parameter γ. The energy difference between the minimum of the H 22 (R) potential, at R = R 2 au, where the initial wavepacket will be prepared for the dynamics described in the following, and

H 11 (R 2 ) is ∆E = 0.15 au.
The interaction with the laser field is described in the dipole approximation as Since the Hamiltonian explicitly depends on time, the stability of the numerical integration has been confirmed by monitoring the norm of the wavefunction. The diabatic-to-adiabatic change of basis has been performed to compute the population of the adiabatic states S 0 and S 1 , for comparison with F-CT-MQC calculations. While quantum dynamics is performed to determine the full wavefunction as function of time, its exact-factorization form can be determined as well, in order to compute the TDPES (that will be shown below). To this end, a choice of gauge has to be made. In the quantum dynamics case, the gauge can be set by imposing that the time-dependent vector potential is zero.

V (t) = -μE(t) = -μE 0 cos(Ωt),
Trajectory-based F-CT-MQC calculations are performed in the Floquet diabatic basis. N tr = 100 trajectories are evolved according to Eqs. ( 21) and ( 37) with the velocity-Verlet algorithm starting with Wigner-sampled initial conditions (the Wigner distribution is determined as the Wigner transform of the initial nuclear probability density of quantum calculations). The Fourier series used to represent the Floquet eigenmodes is truncated to include 2N max + 1 components.

In the calculations presented below, convergence is tested by increasing the size of the basis with N max = 1, 2, 3, 4, 5 (calculations with N max = 6, 7 have been performed as well but they are not shown because convergence is reached with N max = 5). The electronic initial condition is chosen by setting equal to one the harmonic-electronic coefficient corresponding to k = S 0 and n = 0 for all trajectories. We recall that the index k indicates the electronic state, in this case k can be S 0 or S 1 , and n indicates the harmonic, in this case n ∈ [-N max , N max ]. With this initial condition, the electronic equation ( 29) is integrated with the forth-order Runge-Kutta algorithm. The time step for the nuclear and the electronic integration is dt = 0.1 au.

F-CT-MQC calculations provide direct access to the populations of the Floquet diabatic states along a trajectory ν. In order to compare trajectory-based results to the reference results, the populations ρ k (t) of the adiabatic states k = S 0 , S 1 as functions of time are determined as

ρ k (t) = 1 N tr Ntr ν=1 n=Nmax n=-Nmax C ν k,n (t)e iωnt 2 (43)
Similarly, it is possible to interpret the dynamics in terms of the number of photons exchanged between the system and the external field. To this end, the photon population can be estimated as

ρ n (t) = 1 N tr Ntr ν=1 k=S 0 ,S 1 C ν k,n (t) 2 (44)
where n ∈ [-N max , N max ].

A. Case of a weak cw laser

In the weak-field case, some of the density that is initially prepared in the vibrational ground state of the left well, in S 0 , is transferred to the electronic excited state S 1 . While the field keeps driving population back and forth between the two states, with the main portion of the nuclear density still centered around R = 2 au in S 0 , the excited wavepacket evolves towards the avoided crossing located at R ∼ 4 au, and transfers back population to the ground state. The dynamics just described is detailed in Fig. 1 where the nuclear density at three different times is shown, along with its decomposition in terms of Born-Oppenheimer contributions and the gauge-invariant part of the TDPES given as the sum of Eqs. ( 6) and ( 7) (left panels).

The gauge-invariant part of the TDPES does not contain the term with the partial time derivative of Eq. ( 6). It has been shown in previous work [START_REF] Agostini | Mixed quantum-classical dynamics on the exact time-dependent potential energy surfaces: A novel perspective on non-adiabatic processes[END_REF][START_REF] Agostini | The exact forces on classical nuclei in non-adiabatic charge transfer[END_REF][START_REF] Curchod | An exact factorization perspective on quantum interferences in nonadiabatic dynamics[END_REF] , and it is the case here as well, that this gauge-dependent contribution to the TDPES nearly mirrors the shape of the gauge-invariant part, and it is a piecewise constant function of R; its main effect is to reduce the height of the steps, as the step observed at t = 856 au in Fig. 1 at about R = 3 au, but leaving unchanged the slope of the gauge-invariant TDPES outside the steps. It should be noted that, in the present case, where an external time-dependent field is applied, the features of the gauge-invariant and gauge-dependent parts of the TDPES are not as neat as in the field-free case. Figure 2 shows the populations of the Floquet diabatic states for the simulation where N max = 5. These quantities are estimated as the average over the trajectories of the squared moduli of the coefficients C ν k,n (t). The analysis of the populations reported in Fig. 2 gives information about the number of photons exchanged between the system and the external field during the dynamical process. Figure 2 shows that only few Floquet diabatic states are substantially populated during the dynamics. The state labelled k = S 0 , n = 0 starts with full occupation, but rapidly transfers some population to k = S 1 , n = -1, suggesting that a one-photon process is taking place to excite a portion of the ground-state wavepacket. Since the nonadiabatic couplings mediate population transfer between different electronic states within the same harmonic, and since in the region where the nuclear density is prepared these couplings are mostly zero, population transfer at the initial times is only driven by the external field. A small amount of population is also transferred to k = S 1 , n = +1 at the initial time, once again suggesting a one-photon process, but the population remains close to zero (< 0.02 ∀ t), as well as all other populations but k = S 0 , n = 0 and

k = S 1 , n = -1.
Using Eq. ( 43), we can estimate the populations of the physical electronic states as functions of time from the Floquet-based calculations. The population of the ground state S 0 is shown in Fig. 3. 45), from exact quantum-dynamics calculations, and from F-CT-MQC calculations with N max = 1, 2, 3, 4, 5. The colore code is the same as in the upper panel. 44) and averaged over a period of the driving field with an expression equivalent to Eq. ( 45). Only the populations with values > 0.05 ∀ t are shown.

In the upper panel, the estimate from Eq. ( 43) with k = S 0 is used, with different values of N max for F-CT-MQC calculations, to test convergence with the number of harmonics by comparing with exact results. In the lower panel, the oscillations in the populations are averaged out by computing a moving average over a period, namely

ρ k (t) T = 1 T t+T t dτ 1 N tr Ntr ν=1 n=Nmax n=-Nmax C ν k,n (τ )e iωnτ 2 ( 45
)
and similarly for quantum calculations, which allows us to focus on the overall behavior of the population as function of time. As predicted from the analysis of Fig. 1 of the TDPES, N max = 4

is sufficient to converge F-CT-MQC, even though not to the exact value of the S 0 population.

However, for the simulated dynamics, the relative error on the population remains within 10% of the reference (estimated from the lower panel of Fig. 3). Finally, note that, in the upper panel of Fig. 3, the amplitude of the oscillations of F-CT-MQC results is smaller if compared to exact results, probably as effect of the average operation over the trajectories.

Using Eq. ( 44), we can estimate the photon population to interpret the dynamics in terms of single-or multi-photon absorption and emission processes. In order to average out the fast oscillations in photon populations that are similar to the oscillations observed in the upper panel of Fig. 3, we show in Fig. 4 ρ n (t) T , which is a quantity analogous to the one derived in Eq. ( 45), but applied to the photon population. Figure 4 shows that the dynamics is dominated by zero- photon processes, with a small contribution from one-photon absorption (negative values of n) mainly appearing after 2000 au.

In Fig. 5, we compare the nuclear density as function of time with the distribution of F-CT-MQC trajectories for N max = 2 (gray circles) and N max = 5 (black circles). As discussed above, the nuclear density mainly remains localized around R = 2 au at all times, even though already at about 500 au, a small portion of the density -evolving on the excited state -moves towards large values of R. With N max = 2, classical trajectories follow closely the main portion of the nuclear density, but the diverging branch is missed. With N max = 5, instead, some trajectories diverge from the main bundle localized around R = 2 au, even though the time scale is not as in the reference results, and the trajectories do not go as far as the quantum wavepacket within the simulated dynamics. FIG. 6. Strong-field case. The panels are similar to those in Fig. 1 and the same color-code is used.

The nuclear density and the TDPES are shown at times t = 867, 1443, 1874 au.

B. Case of a strong cw laser

Even though the amplitude of the external field is only doubled in going from the weak-field case to the strong-field case, the dynamics in this second example is very different with respect to the previous one, which can be confirmed by comparing nuclear dynamics from Fig. 1 encounters the excited-state wavepacket, producing complex interference patterns at long times.

Clearly, this dynamics is very difficult to capture based on classical trajectories without accounting for interference effects.

Analysis of the Floquet diabatic populations shows that multi-photon processes take place during the dynamics. As in the previous case, the state labelled k = S 0 , n = 0 is fully populated initially, but it rapidly transfers population mainly to k = S 1 , n = -3, to k = S 1 , n = -1 and to k = S 1 , n = +1 . At later times, after 1000 au, other states become populated, even if their population remains below 0.1.

Applying Eq. ( 43), the populations of S 0 and S 1 states are determined from the populations of the Floquet diabatic states, and the S 0 population is shown in Fig. 8. As mentioned above, convergence in achieved for N max = 4, but the results slightly differ from the reference calculations.

In particular, the low frequency oscillations of the S 0 population are not quite captured, as it is evident in the plot of the populations averaged over a period of the drive (lower panel). Indeed, inclusion of more harmonics, for instance from N max = 2 to N max = 4, allows to better reproduce the decay between 500 au and 1500 au, even though quantitative agreement is missing all along the simulated dynamics. Also in this case, and probably stronger than in the weak-field case, the amplitude of the high-frequency oscillations of the population (upper panel) is suppressed in F-CT-MQC calculations, once again as effect of the average over the trajectories.

Similarly to the weak-field case, the populations of the Floquet states allows for the calculation of photon populations. Figure 9 shows that at short times the dynamics is dominated by a zero-photon process, since the population ρ n=0 (t) is much larger that other contributions for t < 1000 au. By contrast with the weak-field case, however, a single-photon absorption process has non-zero contribution since the beginning of the simulated dynamics. After 1000 au, however, multi-photon absorption processes are observed, as well as single-and two-photon emissions (ρ n=+1 (t) and ρ n=+2 (t) are small but non-zero after 2000 au).

Finally, we compare in Fig. 10 the nuclear probability density with the distribution of F-CT-MQC trajectories for N max = 2 (gray circles) and N max = 5 (black circles). In general, we observe that the density spreads over time and moves towards large values of R. For N max = 2 the trajectories only follow one portion of the nuclear density, whereas for N max = 5 the whole nuclear configuration space spanned by the quantum wavepacket is sampled as well by the trajectories.

In this case, and in contrast with the weak-field case, the maximum number of harmonics 44) and averaged over a period of the driving field with an expression equivalent to Eq. (45). Only the populations with values > 0.05 ∀ t are shown. included in F-CT-MQC calculations provides satisfactory results for nuclear dynamics, while it misses some features in the electronic populations. Conversely, in the weak-field case, nuclear dynamics was not correctly captured by the trajectories, but a better agreement was achieved for the electronic populations.

C. General discussion on F-CT-MQC

The exact factorization, together with its trajectory-based treatment, can be formulated employing the Floquet formalism in different, equivalent, ways. In Section II, we have already discussed two possibilities, which we referred to as Floquet adiabatic and Floquet diabatic picture.

Nonetheless, the explicit derivation of F-CT-MQC equations and the numerical results have been presented only in the Floquet diabatic representation. The choice has been made based on practical reasons: (i) the preparation of the initial state for the dynamics, and (ii) the calculation of energy gradients and derivative couplings. The initial state of a simulation is usually chosen to be one of the "physical" states of the system, e.g., the ground state, which is easily identified in the Floquet diabatic picture by selecting the harmonic n = 0. However, in the Floquet adiabatic picture, this identification is not straightforward. In addition, quantum-chemistry codes usually provide gradients of the Born-Oppenheimer PESs, which are identical to the gradients of Floquet It is important to mention that, F-CT-MQC is readily suitable for on-the-fly calculations in combination with ab initio PESs and derivative couplings. The only additional elements used in Section III in comparison to field-free simulations are the electronic transition dipole moment of the molecule and its spatial derivatives, which have to be computed at each nuclear configuration visited by the trajectories, and are both accessible based on quantum-chemistry calculations [START_REF] Sisto | Atomistic non-adiabatic dynamics of the LH2 complex with a GPU-accelerated ab initio exciton model[END_REF] .

In the literature [START_REF] Zhou | Nonadiabatic dynamics in a laser field: Using Floquet fewest switches surface hopping to calculate electronic populations for slow nuclear velocities[END_REF][START_REF] Fiedlschuster | Floquet surface hopping: Laser-driven dissociation and ionization dynamics of H + 2[END_REF][START_REF] Bajo | Mixed quantum-classical dynamics in the adiabatic representation to simulate molecules driven by strong laser pulses[END_REF][START_REF] Chen | On the proper derivation of the Floquet-based quantum classical Liouville equation and surface hopping describing a molecule or material subject to an external field[END_REF] , the Floquet diabatic picture has been preferred over the Floquet adi-abatic one, especially in combination with the trajectory surface hopping algorithm. In fact, as shown in Figs. 1 and6, the Floquet adiabatic PESs can present many avoided crossings, as well as trivial crossings. These features potentially pose challenges for the fewest switches surface hopping algorithm, because the trajectories hop often from one state to another; in addition, a large number of trajectories is required for convergence in order to sample satisfactorily the "hopping space". Perhaps, when working with dense manifolds of coupled Floquet adiabatic states, an Ehrenfest-like approach would be preferable. This suggests that F-CT-MQC in the Floquet adiabatic representation might perform well, similarly to the case of Floquet diabatic representation, due to its Ehrenfest-like form -provided that the issues (i) and (ii) above are efficiently circumvented.

An alternative strategy for deriving F-CT-MQC can be envisaged: rather than invoking the Fourier representation of Floquet eigenmodes, and introducing the harmonics to account for the periodicity of the external drive in the electronic basis, one could directly work with Floquet eigenmodes, as, for instance, those determined by solving the eigenvalue problem (10). This avenue has not been explored in this work, but it would be interesting to investigate how the two formulations compare in the trajectory-based scheme. Clearly, if the exact equations were solved, the two formulations would yield the same result, however, this is not always true when approximations are invoked.

An additional approach to include the effect of the external drive has been tested, but the results have not been reported in the present work. Such an approach simply includes the external drive in CT-MQC equations, where the used electronic representation is the "standard" Born-Oppenheimer representation. CT-MQC yields, in this case, directly the populations of the physical electronic states, and interpretation in terms of photon absorption or emission processes is not possible. For the model studied here, trajectory-based results are very close to exact and Floquet-based results in both weak-and strong-field cases. Therefore, even though such strategy to include the external time-dependent field in CT-MQC has not been yet systematically investigated, it is a promising route for future theoretical and numerical developments.

IV. CONCLUSIONS

Trajectory-based excited-state simulations of systems subject to an external drive are nowadays very challenging. Standard trajectory-based schemes might yield different results depending on the used electronic basis, or electronic representation, since the approximations they are based on are derived in a particular electronic basis. The choice of electronic representation becomes, thus, of paramount importance. In order to easily generalize standard schemes, one might be tempted to use the so-called quasi-static representation, where the electronic states are defined at each time of the propagation as the eigenstates of the electronic problem for a given nuclear configuration (as in standard field-free situations) including the effect of the external time-dependent field. However, unless that external field varies slowly in time, the electronic quasi-static states are not descriptive of the state of the system, and numerical simulations yield unphysical results 83 . The Floquet formalism is, instead, adequate to treat periodically driven systems because the periodic time dependence is somehow treated in a static way -in Fourier space -by extending the electronic space to the space of harmonics. The drawbacks are, clearly, the fact the electronic Hamiltonian becomes unbounded, and that convergence studies on the number of harmonics to be included have to be carried out (and depending on the intensity of the drive, a large number of harmonics need to be considered).

As it has been shown in previous work [START_REF] Abedi | Shedding light on correlated electron-photon states using the exact factorization[END_REF][START_REF] Hoffmann | Light-matter interactions via the exact factorization approach[END_REF][START_REF] Lacombe | Exact potential energy surface for molecules in cavities[END_REF][START_REF] Khosravi | Electronic non-adiabatic dynamics in enhanced ionization of isotopologues of hydrogen molecular ions from the exact factorization perspective[END_REF][START_REF] Khosravi | Exact potential driving the electron dynamics in enhanced ionization of H + 2[END_REF][START_REF] Suzuki | Electronic Schr ödinger equation with nonclassical nuclei[END_REF][START_REF] Fiedlschuster | Surface hopping in laser-driven molecular dynamics[END_REF][START_REF] Abedi | Exact factorization of the time-dependent electronnuclear wave function[END_REF][START_REF] Suzuki | Laser-induced electron localization in H + 2 : Mixed quantum-classical dynamics based on the exact time-dependent potential energy surface[END_REF] , the exact-factorization formalism naturally lends itself to the treatment of dynamics in the presence of an external time-dependent field. In addition, the Floquet formalism can be used together with the exact factorization to interpret and to justify common approaches based on trajectories [START_REF] Fiedlschuster | Surface hopping in laser-driven molecular dynamics[END_REF] We presented the first application of F-CT-MQC to the treatment of excited-state processes with explicit inclusion of a time-dependent external field. The algorithm has been adapted to treat periodically-driven systems with the support of the Floquet formalism, and tested on a model systems subject of an external field with different intensities. The results are promising, especially aiming at the combination with quantum-chemistry approaches to compute electronic-structure properties, but there is clearly room for the development of refined approximation strategies to F-CT-MQC that will be the focus of future studies.

In order to derive an explicit expression for Ċk,n R(t), t only, from Eqs. (A1), we project onto a Floquet diabatic state e iωmt ψ l r; R(t) . To this end, as discussed in the main text, we introduce the change a variable t → s in the harmonics, and we perform a time integration over a period of the driving field only in the time variable s. This operation allows us to integrate out the periodic time dependence due to the external drive and carried by the Floquet diabatic states, while keeping a time-dependent perspective in the "non-periodic" variable t. The projection in the electronic physical space r is performed via the integral over the whole configuration space.

In the left-hand side of Eqs. (A1), the projection in time and space yields

The first term in parenthesis of Eq. (A8a) is the velocity of the trajectory, Ṙ(t), from the characteristic equation ( 21), whereas the second term contains the quantum momentum [START_REF] Min | Coupled-trajectory quantum-classical approach to electronic decoherence in nonadiabatic processes[END_REF] , P(R(t), t), that has been defined in previous work on CT-MQC (but it has a longer history in the context of the quantum-trajectory formalism [START_REF] Garashchuk | Energy conserving approximations to the quantum potential: Dynamics with linearized quantum force[END_REF][START_REF] Garashchuk | Quantum dynamics with Bohmian trajectories: energy con-serving approximation to the quantum potential[END_REF] ). The quantum momentum induces quantum decoherence effects by tracking the spatial delocalization over time of the nuclear density (or, equivalently, of its modulus) [START_REF] Gossel | Coupled-trajectory mixed quantum-classical algorithm: A deconstruction[END_REF][START_REF] Agostini | An exact-factorization perspective on quantum-classical approaches to excitedstate dynamics[END_REF] . Note that the definition (8) of Û [Φ, χ] contains two terms, but the first term has been shown [START_REF] Eich | The adiabatic limit of the exact factorization of the electron-nuclear wave function[END_REF][START_REF] Scherrer | On the mass of atoms in molecules: Beyond the born-oppenheimer approximation[END_REF][START_REF] Schild | Electronic flux density beyond the born-oppenheimer approximation[END_REF][START_REF] Scherrer | Nuclear velocity perturbation theory for vibrational circular dichroism: An approach based on the exact factorization of the electron-nuclear wave function[END_REF] to be smaller if compared to the second term, and we will thus neglect it henceforth.

In addition, such term contains second-order spatial derivatives of Floquet diabatic states, which are usually not available in quantum-chemistry codes to compute electronic-structure properties.

The corresponding term in the TDPES (6) is neglected as well to maintain gauge invariance of the exact-factorization equations in the trajectory-based formulation [START_REF] Min | Coupled-trajectory quantum-classical approach to electronic decoherence in nonadiabatic processes[END_REF][START_REF] Agostini | Quantum-classical non-adiabatic dynamics: Coupled-vs. independent-trajectory methods[END_REF] .

Evaluating the term -i ∇ -A R(t), t of Eq. (A8b) for the remaining part of Eq. (A1b) that contains the quantum momentum, we get

-i ∇C l,m R(t), t -A R(t), t C l,m R(t), t -i k,n C k,n R(t), t d lk R(t) δ mn = (A9) k,n 1 T
T 0 ds e i(ωn-ωm)s dr ψl r; R(t) -i ∇ -A R(t), t C k,n R(t), t ψ k r; R(t)

Here, the nonadiabatic coupling vectors appear, namely d lk R(t) δ mn = 1 T T 0 ds e i(ωn-ωm)s dr ψl r; R(t) ∇ψ k r; R(t) (A10) that only contribute to population transfer between states k and l that belong to the same harmonic, due to the presence of δ mn in the first line of Eq. (A9).

Putting together the results from Eq. (A2) to Eq. (A9), we get the time evolution of the expansion coefficients of the electronic wavefunction in the Floquet diabatic basis along the trajectory ν

Few additional steps are necessary to arrive at the final expression given in Eqs. (29), namely:

(i) we set the gauge freedom such that ν (t) + Ṙν (t) • A ν (t) = 0 (see discussion in Section II B);

(ii) we neglect the term containing the product quantum momentum P ν (t) and of the nonadiabatic couplings d ν lk (see Refs. [START_REF] Agostini | Quantum-classical non-adiabatic dynamics: Coupled-vs. independent-trajectory methods[END_REF][START_REF] Min | Ab initio nonadiabatic dynamics with coupled trajectories: A rigorous approach to quantum (de)coherence[END_REF][START_REF] Talotta | Internal conversion and intersystem crossing with the exact factorization[END_REF] );

(iii) we approximate ∇C ν l,m (t) such that we neglect the spatial dependence of the modulus of C ν l,m (t) and we only keep the spatial derivative of the phase, which we call f ν l,m , i.e., -i ∇C ν l,m C ν l,m f ν l,m (see Refs. [START_REF] Abedi | Dynamical steps that bridge piecewise adiabatic shapes in the exact time-dependent potential energy surface[END_REF][START_REF] Agostini | The exact forces on classical nuclei in non-adiabatic charge transfer[END_REF][START_REF] Talotta | Internal conversion and intersystem crossing with the exact factorization[END_REF] ).

These additional approximations allow us to derive Eq. ( 29).

  with μ the dipole moment operator. Neglecting the nuclear contribution to the total dipole moment, and setting to zero the diagonal elements of μ in the diabatic basis, the off-diagonal elements of the system-laser interaction Hamiltonian are -µ(R)E 0 cos(Ωt) = -βRE 0 cos(Ωt). The transition dipole moment µ(R) = βR is chosen to be a linear function of R, and β = 0.05 au. The frequency of the driving field is set as Ω = 0.05 au, and two cases will be studied to test F-CT-MQC, a weak field case with E 0 = 0.25 au and a strong field case with E 0 = 0.5 au.Quantum dynamics is initiated in the ground vibrational state of the left diabatic potential well (corresponding to S 0 ), centered at R = 2 au with average zero initial momentum. Nuclear mass has been set to the value M = 20000 au. Reference results are obtained by the integrating the full tdSE in the diabatic basis by using the split-operator technique[START_REF] Feit | Solution of the Schr ödinger equation by a spectral method[END_REF] with a time step of dt = 0.1 au.

FIG. 1 . 1 (

 11 FIG. 1. Weak-field case. Left panels: Floquet diabatic PESs (S 0 red lines and S 1 blue lines) with N max = 4. The thick red and blue lines identify the S 0 and S 1 PESs for the harmonic n = 0. Nuclear density (black thin lines), along with its decomposition in S 0 (dashed red lines) and S 1 (dashed blue lines) contributions, and the gauge-invariant part of the TDPES (black line-dots) are shown at times t = 856, 1585, 1962 au. Right panels: Floquet adiabatic PESs computed with N max = 4. Nuclear density (black thin lines) and the TDPES (black line-dots) are shown at the same times as in the left panels.

FIG. 2 .FIG. 3 . 1 (

 231 FIG. 2. Weak-field case. Populations of the Floquet diabatic states as functions of time. Only the non-zero populations are shown.

FIG. 4 .

 4 FIG.4. Weak-field case. Photon population as function of time, estimated from Eq. (44) and

FIG. 5 .

 5 FIG. 5. Weak-field case. Colored areas: nuclear density as function of R and t. Circle: distributions of classical trajectories for N max = 2 (gray circle) and N max = 5 (black circles) plotted every 100 au.

and fromFig. 6 . 4 .FIG. 7 .

 647 Fig. 6. However, similar general conclusions on the performance of F-CT-MQC can be drawn. Let us first analyze the dynamics shown in Fig. 6, based on the comparison between the TDPES and the Floquet diabatic PESs (left panels) or the Floquet adiabatic PESs (right panels). The TDPES is almost at all times comprised within the Floquet diabatic PESs shown in the figure, for N max = 4. In this case, the external field excites the nuclear wavepacket from S 0 to S 1 . While population transfer between the two states is continuously driven, the excited portion of the wavepacket moves towards the avoided crossing and transfer population back to the ground state. At the same time, the ground-state wavepacket is driven towards the avoided crossing itself, where it

FIG. 8 .

 8 FIG.8. Strong-field case. The panels are similar to those in Fig.3and the same color-code is used.

FIG. 9 .

 9 FIG. 9. Strong-field case. Photon population as function of time, estimated from Eq. (44) and

FIG. 10 .

 10 FIG.10. Strong-field case. Similar to the results reported in Fig.5with the same color-code.

  . In the present work, we took a step forward, and we employed the trajectory-based solution of the exact-factorization equations, i.e., the CT-MQC scheme, together with the Floquet formalism, to propose a new algorithm designed to treat periodically driven electron-nuclear systems in the presence of nonadiabatic effects. The CT-MQC algorithm has the advantage of being easily adapted to treat various physical situations, like standard nonadiabatic effects[START_REF] Min | Coupled-trajectory quantum-classical approach to electronic decoherence in nonadiabatic processes[END_REF][START_REF] Agostini | Quantum-classical non-adiabatic dynamics: Coupled-vs. independent-trajectory methods[END_REF][START_REF] Min | Ab initio nonadiabatic dynamics with coupled trajectories: A rigorous approach to quantum (de)coherence[END_REF] , spin-orbit coupling[START_REF] Talotta | Internal conversion and intersystem crossing with the exact factorization[END_REF][START_REF] Talotta | Spin-orbit interactions in ultrafast molecular processes[END_REF] (G-CT-MQC), or external time-dependent fields (F-CT-MQC). The fundamental equations do not have to be completely modified or adapted to include different effects, since such additional effects can be easily included starting from the -quantum-mechanical -formulation of the time-dependent Schr ödinger equation in its exact-factorization form.

  Ċk,n R(t), t ψ k (r; R(t))e iωns (A2)where we have used the normalization condition of the Floquet diabatic states1 T T 0 ds e i(ωn-ωm)s dr ψl r; R(t) ψ k (r; R(t)) = δ mn δ lk (A3)An analogous procedure allows us to deriveṘ(t) • ∇C l,m R(t), t = i(ωn-ωm)s dr ψl r; R(t) Ṙ(t) • ∇C k,n R(t), t ψ k r; R(t) (A4)from the first term on the right-hand side of Eq. (A1a), E l R(t) + ω m ] C l,m R(t), t = (A5)

Ċν

  ds e i(ωn-ωm)s dr ψl r; R(t) R(t), t+ ext R(t), t C k,n R(t), t ψ k r; R(t)from Eq. (A1c). In Eq. (A1b), the term containing the external drive is projected onto the Floquet diabatic state e iωmt ψ l r; R(t) , by performing the time-variable transformation t → s in V (r, R(t), t), before the integration over a period, such thatThe expression of the electron-nuclear coupling operator Û [Φ, χ] used in Eq. (A1b) is obtained by representing the nuclear wavefunction in terms of its modulus and phase, such thatÛ [Φ, χ] M -1 ∇S R(t), t + A R(t), t + i -∇ χ R(t), t= Ṙ(t) + iM -1 P R(t), t • -i ∇ -A R(t), t .

				2
				χ R(t), t	2	• -i ∇ -A R(t), t (A8a)
				1
			k,n	T
		1		T
	k,n	T	0

T 0 ds e i(ωn-ωm)s dr ψl r; R(t) Ĥel r, R(t) + ω n C k,n R(t), t ψ k r; R(t) from the first two terms in square brackets in Eq. (A1b), and

R(t), t + ext R(t), t C l,m R(t), t = (A6) k,n V mn,lk R(t) C k,n R(t), t = k,n

1 T

T 0 ds e i(ωn-ωm)s dr ψl r; R(t) V (r, R(t), s)ψ k r; R(t) C k,n R(t), t (A7)

  l,m (t) = -i E ν l + ω m -ν ext (t) -ν (t) + Ṙν (t) • A ν (t) C ν l,m (t) (A11a) + M -1 P ν (t) • -i ∇C l,m (t) ν -A ν (t)C Ṙν (t) • d ν lk δ mn C k,n (t) -M -1 P ν (t) • d ν lk δ mn C ν k,n (t) -

	ν l,m (t)	(A11b)
	-	i V ν mn,lk C ν k,n (t)
	k,n

J. Flick, H. Appel, M. Ruggenthaler, and A. Rubio, "Cavity Born-Oppenheimer approximation for correlated electron-nuclear-photon systems," J. Chem. Theory Comput. 13, 1616-1625 (2017).

DATA AVAILABILITY STATEMENT

The data that support the findings of this study are available from the corresponding author upon reasonable request.

Appendix A: Electronic equation of F-CT-MQC

The evolution equation for the coefficients C l,m R(t), t is derived from Eq. ( 24), using, in the left-hand side, the result of Eq. ( 27), and introducing, in the right-hand side, the expansion of Eq. (25). In this Appendix, we discuss in detail this procedure, ultimately leading to the derivation of Eq. ( 29).

The evolution of the electronic wavefunction ( 24) can be rewritten in the Floquet diabatic basis as

Here, we rearranged the terms of the original equation, keeping only the term Ċk,n R(t), t in the left-hand side. Moreover, the term containing the spatial derivative ∇ acting on the adiabatic states ψ k (r; R(t)) in the left-hand side of Eq. ( 24) exactly cancels out the corresponding term on the right-hand side, where the operator ∇ acts on C k,n R(t), t e iωnt ψ k r; R(t) . For this reason, only the first term on the right-hand side of Eq. (A1a) is left, which contains the spatial derivative of the coefficients (and not of the states).