Validating Fewest-Switches Surface Hopping in the Presence of Laser Fields

The capability of fewest-switches surface hopping (FSSH) to describe non-adiabatic dynamics of small and medium sized molecules under explicit excitation with external fields is evaluated. Different parameters in FSSH and combinations thereof are benchmarked against multi-configurational time dependent Hartree (MCTDH) reference calculations using SO$_2$ and 2-thiocytosine as model, yet realistic, molecular systems. Qualitatively, we find that FSSH is able to reproduce the trends in the MCTDH dynamics with (and without) an explicit external field; however, no set of FSSH parameters is ideal. An adequate treatment of the overcoherence in FSSH is identified as the driving factor to improve the description of the excitation process with respect to the MCTDH reference. Here two corrections were tested, the augmented-FSSH (AFSSH) and the energy-based decoherence correction. A dependence on the employed basis is detected for the AFSSH algorithm, performing better when spin-orbit and external laser field couplings are treated as off-diagonal elements instead of projecting them onto the diagonal of the Hamilton operator. In the presence of an electric field, the excited state dynamics was found to depend strongly on the vector used to rescale the kinetic energy along after a transition between surfaces. For SO$_2$, recurrence of the excited wave packet throughout the duration of the applied laser pulse is observed for long laser pulses (>100~fs), resulting in additional interferences not captured by FSSH and only visible in variational multi-configurational Gaussian when utilizing a large amount of gaussian basis functions. This feature essentially vanishes when going towards larger molecules, such as 2-thiocytosine, where this effect is barely visible in a laser pulse with a full width at half maximum of 200~fs.

The capability of fewest-switches surface hopping (FSSH) to describe non-adiabatic dynamics of small and medium sized molecules under explicit excitation with external fields is evaluated. Different parameters in FSSH and combinations thereof are benchmarked against multi-configurational time dependent Hartree (MCTDH) reference calculations using SO 2 and 2-thiocytosine as model, yet realistic, molecular systems. Qualitatively, we find that FSSH is able to reproduce the trends in the MCTDH dynamics with (and without) an explicit external field; however, no set of FSSH parameters is ideal. An adequate treatment of the overcoherence in FSSH is identified as the driving factor to improve the description of the excitation process with respect to the MCTDH reference. Here two corrections were tested, the augmented-FSSH (AFSSH) and the energy-based decoherence correction. A dependence on the employed basis is detected for the AFSSH algorithm, performing better when spin-orbit and external laser field couplings are treated as offdiagonal elements instead of projecting them onto the diagonal of the Hamilton operator.
In the presence of an electric field, the excited state dynamics was found to depend strongly on the vector used to rescale the kinetic energy along after a transition between surfaces.
For SO 2 , recurrence of the excited wave packet throughout the duration of the applied laser pulse is observed for long laser pulses (>100 fs), resulting in additional interferences not

I. INTRODUCTION
Femtosecond time-resolved spectroscopy has progressed drastically throughout the last decades, 1 challenging the computational excited state dynamics simulations to include explicitly laser pulses. [2][3][4][5][6][7][8][9][10][11][12] Following a laser excitation to some high-lying electronic state, a wave packet can evolve through different potential energy surfaces (PESs), ultimately deactivating to the electronic ground state by radiationless or radiative processes. The characterization of these dynamical processes requires to consider coupled nuclear-electronic motion, reaching beyond the Born-Oppenheimer approximation. If one also includes the interaction with the laser pulse, the ensuing excited state dynamics will be influenced according to the pulse amplitude and duration, thereby further challenging the calculations as compared to the case of dynamics in the absence of explicit external fields.
Even without explicit laser pulses, the exact quantum treatment of all degrees of freedom in non-adiabatic dynamics is a considerable burden for systems containing more than few atoms.
From the large large number of methods that proliferated, 13 the so-called mixed quantum-classical methods employ trajectories as basis functions, moving according to the classical laws of motion to describe the nuclei. In this category are Ehrenfest, 14 ab initio multiple spawning 15 (AIMS) or fewest-switches surface hopping 16,17 (FSSH) dynamics. These methods show a favorable scaling with size at the price of sacrificing some quantum effects and thus accuracy, in comparison to methods that more readily converge to the exact result like multi-configuratonal timedependent Hartree 18 (MCTDH), variational multi-configurational Gaussian 19,20 (vMCG) or full multiple spawning 21,22 (FMS). Contrary to Ehrenfest dynamics and AIMS, FSSH is not directly derived from first principles; instead, it is based on straightforward assumptions that allow a classical trajectory to describe movement and transitions within and between different electronic states.
Although recent investigations employing exact factorization [23][24][25] and the quantum-classical Liouville equation 26,27 have shed light on the nature of FSSH, a proper derivation of FSSH does not exist. Therefore, systematic improvements of the FSSH algorithm are scarce and known issues, such as the inherent overcoherence, 28,29 are treated by ad hoc corrections that while remedying the issue at hand in some test cases, fail in others. The resulting abundance of available options with different parameters to choose from both redeems and curses FSSH: On the one hand, an appropriate choice of parameters enables at least qualitative accuracy. On the other hand, most of these parameters are devised on one-or low-dimensional models, far away from reality, and when modified in real systems can dramatically affect dynamics. 30 The inclusion of explicit laser pulses in the Hamiltonian is straightforward in wave packet dynamical calculations, but due to the exponential scaling, such simulations are done only in reduced dimensionality. Extensions beyond are offered by modified versions of MCTDH, 10 vMCG, 9 ab initio multiple cloning, 7 or AIMS, 6 showcasing the capabilities of these methods to cope with this additional coupling taking into account many degrees of freedom. The implementation of explicit laser pulses in FSSH is also possible. 3,31,32 However, recent comparisons of FSSH to exact quantum results have revealed a strong dependence on the chosen representation of the interaction with the laser pulse. Wrong representations were found to fail even for H + 2 5, 33 where only a set of Floquet states was able to correctly follow the trends of the quantum results. Similar shortcomings were found for the dissociation dynamics of LiF 6 and for single and dual avoided crossing problems employing a set of Floquet states, 12 showing that no general methodology could be devised so far to treat the coupling correctly. Overall, it seems that including a laser field into FSSH simulations provides an additional uncertainty as to which parameters are more suitable to use. This paper is spurred on bridging recent investigations 33 that highlighted errors present in FSSH and more apparently in the presence of laser fields and the more pragmatic aim of striving to increase the comparability with experimental data by including laser fields in realistic molecular studies. With this in mind, we use two model systems, SO 2 and 2-thiocytiosine, to test the influence of different FSSH parameters dealing with decoherence, representations, and rescaling options of the kinetic energy after a change between PESs or a frustrated hop. These tests are conducted both using an explicit laser pulse to excite the system and the much more common approach of neglecting any external field and just placing the wave packet directly in the optical bright state at the beginning of the dynamics. As a reference, MCTDH (and vMCG) calculations are used to verify the validity of the FSSH simulations, and estimate errors for a given set of parameters, thus showing the grade of applicability of FSSH in realistic molecules.
The remainder of the paper is organized as follows: In chapter II, the framework of FSSH and the various options available therein are presented. Dynamics methods that go beyond FSSH are introduced in chapter III. Chapter IV describes the nomenclature employed and the methodological details. Finally, the numerical results are presented in chapter V followed by the conclusions in chapter VI.

II. FEWEST-SWITCHES SURFACE HOPPING BACKGROUND
Since its original formulation 16 as an improvement to existing surface hopping methodologies, 34,35 FSSH has seen a multitude of adjustments to overcome some of its inherent limitations. 17,[36][37][38] The most severe ones concern energy conservation and the fact that a single independent classical trajectory is unable to describe branching and interference correctly. 29 Instead, the whole wave packet is piggybacking on a single trajectory, resulting in a phenomenon termed overcoherence.
In the following, a short overview of the FSSH methodology with its deficiencies and remedies is given, to provide a backgroud for the parameters that will be used in the subsequent dynamical simulations.
At the core of each FSSH scheme is a classical propagation of the nuclei coupled to a quantum propagation of the electronic wave function bound to the classical trajectory. In both cases, the propagation is based on a set of electronic states |Φ(r, R) α that are commonly obtained from solving the electronic Schrödinger equation. Surface hopping is not restricted to the adiabatic basis and different sets of states can serve as a basis 39 (see section below). For now, a general non-diagonal basis will be assumed although it has been argued that the adiabatic representation is most fitting for the FSSH methodology. 40 The electronic wave function (|Ψ(r, R,t) ) along each trajectory can then be written as |Ψ(r, R,t) = ∑ α c α (t)|Φ α (r, R) (1) where c α are the coefficients for each electronic state. Their time-dependence is given by Here, h β α (R) is the non-adiabtic coupling (NAC) vector Φ β (r, R)| ∂Ĥ el ∂ R |Φ α (r, R) , that indicates the change of the electronic wave function with variation of the nuclear coordinates. H d β α is the matrix element of a complete Hamiltonian that contains any arbitrary coupling and can be written as H d β α (R,t) = H β α (R) −μ β α (R)ε(t) + H SOC β α (R).
The dipole operatorμ β α (R) mediates coupling between states α and β with an external laser field ε(t) and H SOC couples states of different multiplicity (here singlet and triplet states) via relativistic spin-orbit coupling (SOC). H β α (R) is the matrix element of the electronic Hamiltonian in absence of laser and SO couplings. In the case of an adiabatic basis, H β α (R) equals zero for α = β .
Contrary to the coefficients of the electronic wave function that can be distributed over multiple states at once and will fluctuate across a simulation, the nuclei are restricted to move on only one of the PES in each time step, termed the active state. Non-adiabatic effects are included in FSSH simulations via instantaneous switches between PES when the active state changes. Multiple algorithms of when the active state should be switched have been proposed -most of them adhering to the concept that the number of switches that occur during a simulation run of a single trajectory should be minimized, thus coining the term of fewest-switches. 16,41 Since a single trajectory is only able to follow one distinct nuclear rearrangement at a time, swarms of trajectories are employed to mimic a nuclear wave packet and obtain meaningful branching ratios or excited state deactivation times.
One of the flaws of FSSH can be readily seen in the use of a set of independent classical trajectories, which prevents the simulation of nuclear quantum phenomena like tunneling or interference. The advantages of surface hopping, however, rely on its on-the-fly application as a time step in every trajectory only needs to evaluate properties that can be obtained from a single quantum chemistry calculation. The up-side is that the algorithm itself scales well with the size of the system, only depending on the cost of the corresponding quantum chemistry calculation.
As running FSSH simulations necessitates the calculation of multiple trajectories at once, but the independent trajectory approximation inherent to FSSH allows for all trajectories to be computed independently, it is trivial to parallelize.

A. The Choice of Representation
When performing FSSH simulations, different sets of electronic basis states -termed representationsare available. While exact wave packet quantum dynamics in a complete basis is invariant to the choice of representation, FSSH is not. The nuclei of each trajectory are propagated on the PES of the active electronic state and thus changing the definition of the electronic states will change the PES the nuclei evolve on and in turn the observed dynamics.
The most accessible representation is the so-called molecular Coulomb Hamiltonian (MCH), 32,39,42 in which no coupling between states of different multiplicity and no external field is considered. In this completely diagonal picture (DIAG), the amount of small coupling regions is reduced in favor of more strongly coupled avoided crossings. In this work, this DIAG representation will include both, coupling with an external field and SOC, although it is conceivable to include only one of those couplings in the diagonalization while keeping the other as an off-diagonal elements.
A fourth representation that has been employed when it comes to include external fields is the Floquet representation. 4,12,45 In the Floquet representation a set of time-independent states is obtained for a continuous wave by diagonalizing the Floquet Hamiltonian. 46,47 For each state in the original basis a set of infinite new states is created in the Floquet picture that represent the original PES shifted by [-n,-n+1,...,n] times hω where ω corresponds to the frequency of the applied field. These surfaces are not simply shifted by this amount but show additional avoided crossings at the intersection of states with different photon numbers. Performing surface hopping simulations employing a reduced set of Floquet states has been shown to result in surfaces that can coincide with those obtained from exact factorization 33 and thus give the best description for H + 2 in the presence of laser fields -a system where all other representations were found to fail. 5 Floquet theory, however, is not guaranteed to give the best FSSH results 12 and is only exact in the regime of continuous external fields, and breaking down in the regime of ultra short few-cycle pulses. Thus, the Floquet representation will not be applied in the current work.

B. Electronic Decoherence
A wave packet traversing a conical intersection branches into a part continuing on the upper state and another propagating on the lower state. The two parts can reach different regions of phase space or could interact at a later time.
When it comes to mixed quantum-classical simulation methods, description of a passage through a conical intersection is one of the most decisive steps. The FSSH formalism mimics the splitting of the wave packet into two parts while traversing a single conical intersection as a swarm of trajectories that split in two sets following one or the other state. 48 Interestingly, the challenge of FSSH to represent the quantum behavior does not lie on the hopping itself, but on the subsequent evolution. To picture this, we take a look at the evolution of a single trajectory.
When passing through the strong coupling region, the classical trajectory will end up in one of two states, randomly selected based on the coupling strength. While the classical part of this trajectory is subject to a binary choice, namely which gradient the nuclei will follow, the electronic part tells another story. The propagation of the electronic coefficients leads to a ratio of state occupations that resembles the branching ratio of the complete wave packet at the conical intersection.
This distribution of electronic occupations is necessary to ensure a binary hopping choice that resembles the real wave packet so that the swarm of trajectories undergoes a reasonable splitting although every trajectory has no information on any other trajectory. Following the conical intersection, every trajectory is subject to a discrepancy between the classical population (100% in the currently active state) and the electronic populations, which are a distribution between the two states involved in the strong coupling. Doing this, the electronic population of a single trajectory mimics both parts of the wave packet at once although the nuclear movement is dictated by the active state and therefore only reminiscent of this single branch of the wave packet. Dragging this "wrong" part of the wave packet along is termed overcoherence. When another coupling region is encountered, there is interaction both with the "right" and the "wrong" parts of the wave packet and wrong hopping probabilities will be predicted as the presence of this second part of the wave packet is nonphysical.
In the last decades, a plethora of modifications (decoherence methods) to the plain surface hopping have been presented to remedy overcoherence, i.e. trying to adapt the electronic populations in a more or less physical way to resemble the quantum dynamical results. In this work, we will work with two decoherence methods: the energy-based decoherence correction (EDC) 37 and the augmented fewest switches surface-hopping (AFSSH). 49,50 EDC is based on the simple assumption that the branched part of the wave packet in another state will dephase and move into another spatial part where the interactions between these two wave packet parts vanish. For a single trajectory this means that electronic population in non-active states should slowly decay because any branched part of the wave packet in another state will dephase and move into another spatial part where the interactions between these two wave packet parts vanish. This decay is realized by modifying the electronic populations 51 (p i ) of every non-active state in every time step via where E i and E α are the energy of the ith non-active and the currently active state, respectively.
E kin is the kinetic energy and C is a parameter commonly set 52 to 0.1 E h . Any population that is reduced from non-active states is added to the population of the active state to keep the overall population constant. The modified populations (p ′ i ) are then used subsequently. The AFSSH mechanism 50 is more intricate as it tries to track where any part of the wave packet in a non-active state is moving to. In every time step, gradients in the non-active states are collected and auxiliary trajectories propagated in those states. If the trajectories deviate too far or too fast from the active trajectory, the population in this non-active state is considered dephased and set to 0.
Both EDC and AFSSH decoherence corrections have been shown to be improve dynamics over the case of using no decoherence correction at all in a set of test cases without including laser fields. 37,49,50 The EDC has been applied to the excited state dynamics of LiF in the presence of different laser pulses, where no significant improvement over the basic FSSH algorithm without any treatment of the overcoherence was observed when compared to exact quantum simulations. 6 Yet, one should keep in mind that the presented decoherence corrections only explicitly tackle the problem of overcoherence but do not improve FSSH results when it comes to a wave packet recombination event. These recombination events can be important in very specific systems but have minor influence in most cases 17 .

C. Energy Conservation
Propagating a swarm of non-interacting trajectories raises an important issue when it comes to energy conservation throughout the dynamics. The total energy contained within the system should not change in the course of the dynamics. This energy can be distributed into kinetic and potential energy in all nuclear and electronic degrees of freedom. When a hop between PES occurs, the potential energy of the trajectory undergoes an instantaneous change because the active state is switched. This behavior would cause discontinuities in the total energy of the single trajectory and possibly even that of the total ensemble. To ensure energy conservation of the total ensemble one typically enforces total energy conservation along each individual trajectory. For this, every change in the potential energy of the trajectory induced by a surface hop is compensated by adapting the kinetic energy correspondingly. In Tully's original prescription, 16 only the nuclear momenta parallel to the NAC vector h αβ are rescaled. If NAC vectors are not available in a simulation, the gradient difference vector g can be used as a substitute.
with F α and F β being the gradient in the active state and the non-active state, respectively. Even simpler is to rescale the full nuclear momentum, which offers an intriguing simplicity as no NACs vectors need then to be computed. Hops to a lower PES result in an increase in kinetic energy while transitions to higher-lying PESs are accompanied by a reduction in the kinetic energy of the system. When trying to enforce any of the mentioned energy conservation procedures, a special case can be encountered if the FSSH algorithm wants to switch the active state to a higher energy PES but the kinetic energy in the system is insufficient to bridge the energy gap to this PES. Such an attempt at hopping is rejected by most algorithms and fittingly coined as a frustrated hop. The presence of frustrated hops has been noted early on and it was found that they are necessary to retain detailed balance within surface hopping simulations. 53,54 However, the number of observed frustrated hops can vary significantly, depending on which of the rescaling schemes is employed. This is due to the different amounts of available kinetic energy in e.g. rescaling along h or rescaling along the full velocity vector. More energy is available in the full velocity vector as the velocity component along the NAC vector is included therein. Using the full kinetic energy rescaling scheme, i.e. enlarging a system, e.g. by including additional non-interacting molecules at large distances, increases the total kinetic energy thus enabling the possibility to jump to higher-lying states than when non-interacting molecules are absent. If rescaling along h is used, inclusion of additional non-or weakly interacting molecules has none or limited influence on the NAC vector and thus the same amount of energy as without the additional molecules is available for hopping to surfaces of higher energy. Rescaling along the NAC vectors has been found to give results that are in best agreement with quantum dynamics and is then treated as the preferred option when NAC vectors are available. 30,54 Unfortunately, the concept of strict energy conservation throughout the dynamics does not hold in the presence of an external field, which stimulates absorption or emission of photons, adding or subtracting energy to the system. When an external field is applied, two options are available to deal with this flow of energy: One is to suspend the conservation of the total energy for the duration of the pulse. Then, the nuclear momenta will not be adapted at a surface hopping event and the trajectory will be propagated on a different PES using the same nuclear momenta. Therefore, any changes in the potential energy are also found in the total energy. Another option is to find criteria to verify whether a hop is laser-induced. This can be realized by tracking the change in the electronic coefficients due to the laser coupling as opposed to changes in the coefficients due to NACs. 55 Alternatively, one can set an energy interval around the center frequency of the external field 32 and if the energy gap at a hopping event falls within this interval, the hop is labelled fieldinduced and thus no rescaling of the velocities is performed.

A. Multi-configurational time-dependent Hartree
Introduced in 1990, the MCTDH 18,56,57 method is one of the most versatile methods to simulate non-adiabatic processes, due to the inherent possibility to converge towards the exact solution.
Briefly, MCTDH is based on an expansion of the wave function using a time-dependent basis of single-particle functions (SPFs), |ϕ : where the coordinates of Ψ are f explicit degrees of freedom (q f ), which can be combined to form a set of actual coordinates Q i . Each Q i then represents one or more degrees of freedom at once, in a process called "mode combination". Φ J is a single Hartree-product preceded by its corresponding coefficient A J . Through the Dirac-Frenkel variational principle, the best suited expansion coefficients and SPFs can be determined directly for each time-step. Therefore, MCTDH can be understood as a method that allows describing the evolving wave packet in a reduced number of used basis functions for each time-step, and thereby minimizing the computational effort while maintaining a maximum amount of accuracy.

B. Variational multiconfigurational Gaussian
A drawback of MCTDH is that the SPFs themselves still use a static grid and the whole method In short, the wave function ansatz for Ψ then reads: where the electronic wave function for a state α, |ϕ (α) (r; R) is multiplied by a set of Gaussian basis functions (GBFs) |χ

A. Nomenclature
In order to distinguish between the different combinations of FSSH parameters employed in this work, the following short-hand notation will be used: REPRESENTATION DECOHERENCE frustrate hops rescaling (8) The possible options for each of the keywords are listed in Table I and explained as follows.
(i) Representation. Two representations for the propagation are used, the MCH or the completely diagonal picture (DIAG) previously introduced in Section II A (ii) Decoherence. To correct for overcoherence we use AFSSH or EDC, see Section II B. Additionally, the NONE option indicates that no decoherence correction is included.
(iii) Rescaling after a surface hop. For the treatment of velocity rescaling after a transition between states, we rescale the velocity vector after a hopping event along the full velocity vector (v), the non-adiabatic coupling vector of the involved states (h) or the gradient difference vector (g), or we do not rescale the momenta at all (none).
(iv) Frustrated hops. If an insufficient amount of kinetic energy is available to compensate for a hop to a higher-lying PES during the rescaling process, the hop is cancelled and termed frustrated. REPRESENTATION  At such frustrated hops, the velocities of the trajectory can be reflected. The same set of vectors as in (iii) is available to reflect the momenta, resulting in the options -v, -h and -g, where the minus indicates the reflection event. Alternatively, the velocities can be kept at a frustrated hop without any reflection, which is labelled "+".
Using this notation, a FSSH setup in the diagonal representation, using the AFFSH decoherence correction and rescaling along the NAC vectors at hopping events and no reflection at frustrated hops would be denoted as DIAG AFSSH + h . Note that not all combinations of these parameters result in stable setups for FSSH dynamics: e.g. in the MCH representation, the non-adiabatic coupling vector h between singlet and triplet states is zero and therefore cannot be used to rescale in simulations including singlet and triplet states. Apart from these incompatibilities, all possible combinations between options (i), (ii) and (iii) have been used. Throughout this work, a consistent pairing of the used vectors in (iii) and (iv) has been used: This means for example, that when using the gradient difference vector in (iii), reflection of frustrated hops in (iv) cannot be conducted along -v or -h but has to occur along -g. Using the option to not rescale the kinetic energy at all in (iii) leads to zero frustrated hops and is therefore only paired with "+" for the treatment of frustrated hops. The option "+" in (iv), however, can be paired with all rescaling vectors (iii) to shed light on the difference between dynamics where reflection takes place at frustrated hops compared to dynamics, where the same trajectory is not reflected at this point.

B. Error quantification
To quantify the difference between a reference MCTDH calculation and a FSSH calculation with a particular set of parameters, we define an error ε calculated as where max t is the final simulated time, ∆t is the employed time step, p i ,t and p i,t,re f are the The error is further split into contributions from the singlet and triplet states, by summing the errors over all singlet (ε sing ) and all triplet states (ε trip ), respectively.
The error measure in Eq. 9 is not well suited to describe the error in the dynamics if a laser pulse is present: As mentioned above, the maximum error is 2, achieved if all the population is found in differing states for both the reference and the investigated dynamics. However, when we imagine an extreme case, where a weak laser excites only 1% of the ground state population to higher-lying states in both the reference and the investigated dynamics, one would obtain a maximum error of 0.02 due to the agreement in the 99% of population staying correctly in the S 0 ground state. Therefore, we use two differently calculated errors for simulations that include laser fields, in order to provide errors that can be compared to errors from dynamics without a laser field: First, the deviation in the ground state population (ε S 0 ) is calculated according to Eq. 9 only using the ground state population; this will indicate the capability of a given set of FSSH parameters to describe the initial excitation process. To factor out the differences in the excitation process itself from the subsequent excited state dynamics, the error in the excited state populations (ε r ) is calculated using a renormalized excited state population. This is done by renormalizing the population in the excited states in every time step to 1 for both the reference and the FSSH dynamics yielding ε r (the superscript thus indicates the renormalization within this error measure: To avoid numerical instabilities where the denominator in p i ,t 1−p S 0 ,t or p i,t,re f 1−p S 0 ,t,re f is close to zero, ε r is only calculated for time steps where the S 0 populations in both dynamics are below 0.98. To compare the diabatic MCTDH (and vMCG populations) to FSSH populations and to obtain diabatic populations for FSSH dynamics that contain both information about the distribution of active surfaces and the corresponding electronic populations at the same time, the methodology applied in Ref 59 has been employed.

C. Dynamical Propagations
The FSSH simulations have been carried out with the SHARC 39,60 program suite using pySHARC 61 to drastically reduce I/O overhead. A set of 1000 nuclear initial conditions obtained from a ground state Wigner sampling 62 is used for both SO 2 and 2-thiocytosine. The simulations employed a time step of 0.5 fs for the nuclear propagation in the field-free case and a nuclear time step of 0.05 fs in the presence of an external field. The electronic wavefunction is propagated in a locally diabatic basis with nuclear time steps of 0.02 fs and 0.002 fs in the field-free and field-including case, respectively. 63 The reduction of step size when including an electromagnetic field is necessary to capture the rapidly oscillating field.
As explained above, various FSSH options and modifications will be tested and compared to reference quantum dynamical results. Note however, that the representation in which FSSH is performed is never the purely diabatic representation of the LVC model (see below). Two different sets of initial electronic coefficients have been employed: For the simulations without an explicit laser field, the initially active electronic state at the start of the FSSH dynamics is set to the MCH state that has the largest overlap with the diabatic bright state of SO 2 and 2-thiocytosine for each initial condition. After setting the initial electronic state, the electronic coefficients are adapted correspondingly so that the initial electronic population at the start of the dynamics amounts to 1 for the diabatic bright state. This methodology represents an instantaneous excitation of the complete ground state wave packet to a single diabatic bright state, as it would follow an ideal instantaneous δ -pulse. For simulations in the presence of a laser field, the initial electronic state was set to the lowest energy state and no modification of the initial electronic coefficients was done.
The QUANTICS package has been employed to run MCTDH and vMCG dynamics. 64 In MCTDH, the Adams-Bashforth-Molton predictor-corrector integrator of 6th order and the multi-set formalism have been used. Convergence of MCTDH dynamics was deemed to be reached if both of the following criteria have been met: First, the weight of the last SPF assigned to a degree of freedom did not exceed a value of 0.001. Second, the number of grid points in a mode was taken to be sufficient if it was larger than n + 3 · dn for all states.
vMCG calculations were run in a single-set formalism which was found to be computationally more efficient than the multi-set formalism and converged faster to the MCTDH results as well, which is in line with previous observations. 9 Integrals between gaussian wave packets have been calculated up to 4th order. The Runge-Kutta integrator of fifth order has been used. The number of considered gaussian basis functions in vMCG ranged from 10 to 100. For both, vMCG and MCTDH, initial conditions were obtained by populating the vibrational ground state of the lowestenergy electronic state of all normal modes considered. In simulations without an external field, this initial wave packet was set to start in the diabatic bright state while for simulations in the presence of an external field, no additional steps were taken. The input and operator files used can be found in the supporting information.

D. Definition of the laser pulse
The coupling with the external field is described within the the semi-classical dipole approximation and can be written as −μ β α (R)ε(t) neglecting any further interaction terms. The pulse has a Gaussian shape, defined as with field amplitude ε 0 t p , carrier frequency ω i , carrier envelop phase η and pulse duration t p equivalent to the full width at half-maximum (FWHM). We assume a linearly polarized pulse along the x-direction (e i = x) of the transition dipole moment and a phase η=0. The frequency is set to be in resonance with the brightest state of the the corresponding molecules The field amplitude of the laser field for a given t p , ε 0 t p , was varied for different lengths of the laser pulse according to ε 0 t p = ε 0 17 · 17 t p . This way the pulse energy (area of intensity throughout the pulse duration) for different laser pulse lengths is kept constant. The pulses with a t p of 17 fs serve as reference pulses for the determination of the field amplitudes for other pulse lengths and where set to ε 0 17 = 0.03 a.u. (15.44 GV/m) for SO 2 and ε 0 17 = 0.01 a.u. (5.15 GV/m) for 2-thiocytosine. Using these amplitudes, about half of the S 0 population in the FSSH simulations was excited in both molecules, which was deemed a good tradeoff between minimizing non-linear behavior and Rabi-oscillations while still exciting enough population in FSSH to yield reasonable statistics.

E. LVC model
The PES employed for the dynamical simulations are parameterized using a linear vibronic coupling (LVC) model. 65 Vibronic coupling models are diabatic representations of the electronic states close to a reference point and capable of describing dynamics and conical intersections close to the reference point. A suitable reference point capable of describing interactions with all excited states directly after initial photoexcitation is the ground state equilibrium geometry. The spatial dependence of the excited states is then cast into mass-frequency scaled normal mode coordinates using the normal modes of the ground state as a basis. When truncating this Taylor expansion at first-order, only linear terms that depend on a single normal mode displacement Q i are obtained, resulting in where H LVC is the diabatic LVC Hamiltonian and H (0) contains the zero-order harmonic potential approximations to the PES Here, V 0 is the harmonic ground state potential along every normal mode i. The first order terms in W (1) are state-specific and consist of electronic energy shifts ε, intrastate gradients κ and couplings between two states λ : Truncating the Hamiltonian after the first-order terms results in a rather crude PES able to accurately capture the form of potential only near the reference point. Parameters for the systems investigated here were taken from Ref 61, that also demonstrates that the LVC approximation is able to reproduce the main characteristics of the corresponding on-the-fly dynamics.
For use in the remainder of the work, the diabatic states and properties were transformed to either the MCH or completely diagonal representation during the propagation.

V. RESULTS AND DISCUSSION
A. SO 2 The excited state dynamics of SO 2 after irradiation has been extensively investigated in the last decade, illustrating a complex interplay between singlet and triplet states. [66][67][68] In this work, we use a LVC Hamiltonian that contains 4 singlet and 3 triplet states obtained with multi-reference configuration interaction including single excitations, which is able to reproduce the main features of the full-dimensional excited state dynamics. 61 In order to find an optimal set of surface hopping parameters able to describe the excited states dynamics of SO 2 in the presence of an explicit laser pulse, we first validate different parameter sets for the ensuing excited state dynamics in the absence of the laser field.

SO 2 dynamics in the absence of a laser field
To simulate the non-adiabatic dynamics of SO 2 without an explicit laser pulse excitation, the ground state nuclear wave packet is vertically placed in the bright 1 B 1 state. Figure 1 shows the results obtained with MCTDH, vMCG and FSSH dynamics. The MCTDH dynamics (Fig. 1a), which will serve as a reference throughout, shows ultrafast population transfer from the initially states is very well reproduced, with minor deviations at later times. The differences in the triplet states, however, are more pronounced. There is a continuous transfer to the triplet states, mostly to the 3 B 2 state, which is not observed in the MCTDH dynamics, in which the population in this state stagnates after 250 fs. To quantify the differences between vMCG and MCTDH dynamics, we calculate the error ε following Eq. 9. This results in the average deviation between the reference populations and the dynamics in question in each time step. The so-derived error is ε=0.069 and it is attributed almost equally to the 1 B 1 , the 1 A 2 , and the 3 B 2 states. Increasing the number of used GBFs is expected to decrease the error until convergence towards the exact result (see Section S1 C in the supporting information for the ε values for simulations using different numbers of GBFs).
The FSSH simulations yield different results depending on the choice of the parameters discussed above (representation, decoherence correction, rescaling of the kinetic energy after a hopping event and reflection of the kinetic energy after a frustrated hop). For the sake of briefness, we show in Fig. 1c only the populations for a single parameter set: DIAG EDC + V . Compared to the MCTDH reference, the FSSH dynamics start with a similar transfer between the 1 B 1 and the 1 A 2 state, matching the time needed for this transfer. However, after the initial 50 fs the oscillatory transfer between both states is far less pronounced than it was in MCTDH. The damped oscillations in FSSH follow the general trend of the MCTDH populations, with a slow decline in 1 A 2 population at later times. The triplet state population is reminiscent of the vMCG dynamics, with a continuous transfer to the 3 B 2 state. The error associated to this simulation gives ε=0.248, being the the main sources of error the damped oscillations and the transfer to the 3 B 2 state. Unsurprisingly, FSSH has a larger error than vMCG (approximately 3 times larger).
In order to investigate the influence of different FSSH parameters all sensible combinations of parameters according to Section IV A have been used to run FSSH dynamics. When comparing the influence of the representation on the errors, it is found that using the MCH representation gives smaller deviations than its completely diagonal counterpart. While this does not hold for ε sing , where the average errors obtained using either representation are almost identical, big differences can be observed in ε trip : all MCH simulations yield an almost identical error while only half of the simulations employing the diagonal representation yield comparable errors. The difference observed in the average error due to the representation is therefore due to the better performance of the MCH representation to describe the triplet state populations. The decoherence correction inflicts very small differences among the three different parameters investigated. However, the FSSH dynamics on the investigated SO 2 model system is very sensitive to the choice of the vector along which rescaling of the kinetic energy should be conducted. Here, two clear favorites emerge of which rescaling along the total velocity vector outperforms the crude option of not changing the velocities at all after switching the active state. Yet, both of these pa-All states rameters are found with significantly smaller average error than rescaling along the non-adiabatic coupling vector (h) or the gradient difference vector (g), which gives the highest average error of all investigated options for all parameters. The bad performance of the h and g parameters is mainly due to the triplet states, where they are the cause for almost all large deviations from the MCTDH reference. Interesting to note here, is that choosing to rescale along h or g significantly reduces the amount of available kinetic energy that can be used to reach a higher-lying surface, as compared to option of rescale along the complete velocity vector and of course the option of not to adapt the energies at all. When comparing the number of successful and frustrated hops for two sets of parameters that differ only in the rescaling vector, e.g. v and g, a total of 3629 executed-versus 3833 frustrated hops is found for MCH EDC + g while 7484 hops against only 1041 frustrated hops are found for MCH EDC + v on the complete set of 1000 trajectories. The higher available kinetic energy in the complete kinetic energy vector therefore enables a lot of jumps that are forbidden when just using the gradient difference vector for rescaling. We noted that switching to the diagonal basis increases the total number of observed hops drastically (39908 for DIAG EDC + g and 41738 for MCH EDC + v ) due to the presence of more trivial and avoided crossings, as the triplet states are split into their respective components. Although the total number of hops increases in the diagonal picture, the amount of frustrated hops stays on a comparable level, with 4664 frustrated hops ( DIAG EDC + g ) and 859 ( DIAG EDC + v ). This might shed some light on why the largest errors are associated with using either of the two decoherence corrections in the diagonal representation along with rescaling the kinetic energy after a hop along h or g, as an increased number of hops slows down both types of decoherence corrections. In general, for the treatment of frustrated hops, reflection of the full velocity vector results in the lowest average error. As this option is always tied to using the full velocity vector for estimating if sufficient kinetic energy is present, which is the largest amount of available energy apart from not rescaling the energy at all, the number of frustrated hops is very low in these dynamics. The notion of "few frustrated hops equals better agreement between populations" is reinforced when realizing that most best-performing sets did not reflect frustrated hops at all. However, combinations of rescaling along h or g with a complete neglect of frustrated hops gave the highest disagreement for singlet and triplet populations, showcasing that a correct treatment of frustrated hops is needed in these cases.

SO 2 dynamics in the presence of a laser field
In the presence of a laser field, both the excitation and the subsequent excited state evolution of the population is influenced by the laser pulse. We shall investigate both processes separately, starting with the excitation process itself. The effect of applying longer pulses to excite population from the 1 A 1 ground state is shown in Fig. 3. Before discussing the details of the observed trends, the reader is reminded that laser pulses of different length employ different field amplitudes ε 0 t p (see Eq. 11) to carry the same total energy.
Within MCTDH, three different dynamical regimes are observed: i) the shortest pulse (t p =2 fs) excites around 30% of the population to higher-lying states, ii) laser pulses with a t p between 10 and 50 fs excite about 55% of the population, and iii) very long pulses beyond t p of 100 fs induce a diverging behavior, where the ground state is repopulated at later times where the laser is still active. The reduced excitation in the dynamics including the very short t p =2 fs laser pulse is due to the fact that the pulse carries only few cycles. This results in a higher uncertainty for the energy and interference effects are increased as the amplitude of the laser field (ε 0 t p ) is adapted to pack the same amount of energy as the longer pulses in this short laser pulse. When comparing the excitation process of the t p = 10 fs laser pulse to the t p = 17 fs, t p = 30 fs or t p = 50 fs dynamics, most of the effects hindering the excitation towards excited states present in the t p = 2 fs case are gone and an almost identical level of excited population is achieved. This excitation of an identical amount of population is reminiscent of an ideal non-interacting case where dynamics in the excited states initiated at the start of the laser pulse does not enhance or hamper further excitation during the duration of the laser pulse. When going to the t p = 100 fs and t p = 200 fs laser pulses, this ideal picture does not hold true anymore, and a more complex S 0 population behavior is observed.
Inspection of the nuclear wave function for the dynamics using t p = 100 fs and t p = 200 fs lasers reveals that the diverging behavior is caused by the return parts of the excited state wave packet that re-enter the Franck-Condon region. Once part of the wave packet is in the Franck-Condon region, additional interference terms arise as the returning excited wave packet and the remaining ground state wave function can interfere if the laser still couples the bright excited state with the ground state. To verify whether this is actually the cause for the observed behavior, new sets of MCTDH simulations have been carried out where two identical t p = 2 fs lasers separated by a time interval τ are employed. This way, one laser acts as a pump pulse and the other as probe, detecting whether the returning wave packet causes interference terms. The resulting populations for different τ delay times are shown in Fig. 4a. It can be seen that the probe laser just excites more S 0 population for small τ values. Upon reaching a time delay τ > 85 fs, a strong dependence on τ is found, which first enhances the excitation induced by the probe pulse but results in a reduced excitation for delays of 89, 93, 97 fs. From this it can be concluded that after ca. 90 fs, a recurrence of the excited wave packet takes place and further coupling with the laser field yields a new interference term, drastically altering the observed overall excitation process.
As the lasers up to t p = 50 fs are too short in time to be still active when the wave packet returns, the excitation follows the observed pattern. Only when going beyond this, a disturbance in the excitation pattern due to this effect is observed. To further support this argument, a complex absorbing potential (CAP) has been employed in the MCTDH calculations that destroys the excited state wave packet once it leaves the Franck-Condon region after initial excitation, thus hindering a possible return. Indeed, as seen in Fig. 3b, the S 0 populations do not show this interference term anymore. Instead, the laser pulse induces a 55% S 0 population inversion. Figure 3c demonstrates that vMCG dynamics shows an extremely good agreement with MCTDH for laser pulses up to t p = 50 fs. Therefore, for short to medium laser pulses, vMCG is able to capture the most important parts of the excitation process induced by explicit laser pulses.
However, for longer laser pulses, vMCG deviates strongly from the MCTDH reference as the behavior due to the interference with the returning wave packet observed for the S 0 population of the t p = 100 fs and t p = 200 fs lasers are not reproduced.
Conducting the same pump-probe laser sequence dynamics as for MCTDH using two t p = 2 fs pulses resulted in a perfect agreement between the vMCG and the MCTDH dynamics (see Fig. 4b).
This shows that vMCG is well suited to describe the additional interference in the pump-probe setup, especially considering that only 50 GBFs were necessary to obtain the agreement. However, in simulations employing long laser pulses with t p = 100 fs and t p = 200 fs, vMCG does not capture the interference that was observed in MCTDH. In these cases, the used GBFs are spread too thin, as the vMCG algorithm tries to describe the complete distributed excited state wave function at once -a much more difficult task than tracking the single excited state wave packet A screening for all combinations of surface hopping parameters considered in Section V A 1 has been conducted using the t p = 30 fs laser pulse, where no additional interference due to recurrence of the wave packet is observed. Note that the total calculated error with respect to the reference MCTDH dynamics is now split into an error in the S 0 population, ε S 0 , calculated according to Eq. 9 while the error in the excited states, ε r , is calculated using renormalization following Eq. 10. ε r is again split up the contributions of singlets and triplet states. The total calculated error is then just the sum over these three contributions: ε = ε S 0 + ε r sing + ε r trip . All calculated ε values for all states, the S 0 , the excited singlet and the excited triplet states are collected in Fig. 5 and Tables SIV-SVII. A very flat distribution of error values is observed for the total error, ranging from 0.265 to 0.462 showing an increase in overall error by 75% for the worst set of parameters with respect to the best set.
First, the capability of the different parameters to describe the initial excitation is discussed.
When investigating the influence of the representation on ε S 0 , the MCH picture is found to be better on average than the diagonal counterpart.
Regarding decoherence, the EDC performs best on average, giving very similar deviations from the MCTDH populations, irrespective of the chosen representation. This consistently good performance is rather surprising following the observations above as the EDC in general leads to an underestimation of the excited population due to its instantaneous nature. However, the ability to converge fast towards a trajectory that no longer contains considerable electronic S 0 population after the initial excitation is found to be beneficial here. This lowered population in the S 0 of an already excited trajectory reduces the chance of stimulated emission and therefore increases the overall agreement with the MCTDH reference. The use of AFSSH is found with a higher average error than the EDC. All sets of parameters that use AFSSH paired with the MCH representation result in a ε S 0 between 0.058 to 0.076 and therefore show the lowest ε S 0 values. However, when using AFSSH in combination with the diagonal representation, the errors essentially double to a range of 0.123 to 0.136, showing a strong sensitivity with the representation. Such differences are not observed in the more robust EDC, indicating that the diagonal representation is interfering with the AFSSH algorithm in SO 2 .
Using no decoherence correction at all is associated with a very bad agreement throughout all considered sets of parameters. This indicates that the correct decay of the ground state electronic coefficient plays a crucial role once the trajectory is excited. The performance in the absence of a treatment for the inherent overcoherence is slightly increased in the MCH representation and worsened in the diagonal one.
The different parameters for rescaling the kinetic energy after a hop and treating frustrated hops yielded almost identical averages. Indeed, this behavior is expected as these options manipulate the excited state dynamics after the initial excitation has occurred and are therefore of limited importance to the excitation process itself.
We now analyze the dynamics in the excited states and see that a much more compact error distribution in error values is observed for the renormalized ε r sing as compared to the dynamics in absence of a laser field. Almost all parameters give identical averages, with these averages also being smaller than for the dynamics without the laser field. The lower average errors are mostly due to loss of the fine structure of the population curves in the 1 B 1 and the 1 A 2 states: This fast and strongly oscillating behavior of the excited state populations depicted in Fig. 1a was one of the major sources of disagreement between MCTDH and FSSH, with FSSH predicting more damped oscillations. Applying a long laser pulse completely changes the observed picture, as the starting point of the dynamics in the excited states is now distributed across the duration of the complete laser pulse as can be seen in Fig. S1 of the SI. The delayed excitation results in a smearing of the transfer between the 1 B 1 and the 1 A 2 states that now show a simple oscillating behavior without fine structure. This behavior is easier to reproduce using FSSH, which also profits from the very forgiving nature of using populations to calculate the deviation from exact quantum dynamics.
For the triplet states, all parameters that were found to increase description of triplet states in the simulations without explicit laser fields are also found to be the dominating factors in the presence of laser fields. Hence, the MCH representation gives significantly improved results over the diagonal representation for the triplet state populations too. The EDC is associated with higher errors in the triplet populations than using no decoherence correction or AFSSH. For rescaling the energies after a hopping event and the treatment of frustrated hops, the velocity vector is the most prominent option followed by the complete neglect of adapting the kinetic energies.
Altogether, we find that an appropriate choice of the decoherence correction together with a compatible representation are the driving factors to capture the excitation process correctly. The EDC emerged here as a rather robust variant while AFSSH was found to perform better in the MCH basis. For the remaining dynamics, the MCH basis was found to increase agreement for the triplet populations with the MCTDH reference. Additionally, using the h and g vectors to rescale or reflect the velocity vector along after a real or frustrated hop was found to give larger deviations with respect to the reference.
There are multiple approaches to find the best set of surface hopping parameters in the presented case: either the best performing set of parameters for a given ε is taken or the parameter of each option that obtained the lowest average error is considered. In this work, the parameters yielding the lowest average errors have been taken for each option, as this represents a more robust error measure instead of cherry-picking a single combination that could well be the result of error compensation. Therefore, the following sets have been determined in absence of a laser field: MCH NONE −v v (best average ε and ε trip ), DIAG EDC −v v (best ε sing ). For dynamics initiated using a t p = 30 fs laser, the combinations MCH EDC −v v (best ε), MCH EDC −g g (ε S 0 ) and MCH EDC + v (ε r sing ) were added, as the on average best performing options resulted again in MCH NONE −v v for ε r trip . These sets where then used to simulate dynamics using different laser lengths for excitation varying in size from t p = 2 fs to t p = 50 fs and ε S 0 calculated with respect to the MCTDH dynamics. The resulting deviations from the reference dynamics are depicted in Fig. 6. Three different classes of surface hopping parameters can be clearly distinguished: (i) surface hopping sets using the MCH representation and the EDC clearly show the same trend for different t p of the laser, starting with a very small error that increases for t p = 10 fs and t p = 17 fs before decreasing for lasers t p = 30 fs and t p = 50 fs long. Essentially, no differences between the different rescaling or reflection parameters can be discerned for these three combinations. (ii) Employing the diagonal representation initially results in a larger deviation from the MCTDH S 0 population for short laser lengths, which can be attributed to the different form the laser coupling takes in the diagonal repre- sentation. For longer laser lengths this is evened out and no differences are observed between EDC in the diagonal and the MCH representation when looking at the S 0 populations. (iii) A completely different picture can be observed when no decoherence correction is used ( MCH NONE v v ). Here, the excitation due to very short laser pulses is well described but as soon as the excitation events and subsequent dynamics start occurring at the same time scale, the need to treat overcoherence is apparent. Longer laser pulses lead to a larger deviation from the MCTDH S 0 population resulting in very larger errors for this set of parameters.

B. 2-Thiocytosine
The substituted nucleobase 2-thiocytosine represents a much more challenging system for dynamics due to the increased number of normal modes with respect to SO 2 . Ab initio on-the-fly FSSH simulations and experimental results are available for 2-thiocytosine, 69 showing a significant triplet yield. Substituting the on-the-fly PES by the more rigid LVC Hamiltonian was found previously to result in a reduced transfer towards the triplet states after 500 fs but to maintain all other essential features. 61 With this in mind, in order to test the effect of different FSSH settings, here we consider simulation times of up to 400 fs.  In order to obtain a reference dynamics to assess different FSSH parameters using an affordable MCTDH computation, the amount of normal modes was reduced using the so-called SHARCgym. 70 In the SHARC-gym systematic calculations to identify the most important normal modes are carried out as follows. For every vibrational normal mode, an independent FSSH simulation is conducted with all gradients and coupling elements associated with that mode set to zero. The so-obtained populations are then compared to the reference populations in Fig 7a by calculating ε values according to Eq. 9. If the omission of a mode results in a large ε, the mode carries essential coupling terms for the overall dynamics and should be included in the reduced set (see the ε values associated with the neglect of specific normal modes in Section S2 A). Using this approach, it is possible to find the modes relevant for the deactivation dynamics without any selection bias that could derive if the selection would be conducted based on the strength of specific coupling elements -a common approach in low-dimensional systems carrying limited amounts of states. 71,72 Following the SHARC-gym we selected the ten most important modes that can reproduce the 33dimensional dynamics of 2-thiocytosine. This reduced 10-D set will be now used for the remainder of the simulations. As it can be seen in Fig. 7b   are two exceptions to this, as the MCH EDC + g and the MCH EDC g g sets of parameters are the best performing sets investigated. All but a single set of parameters using v gave better ε values than for any set using rescaling along h. Reflection of frustrated hops is found to be important; however, due to the pairing criterion of using the same vector for reflecting at a frustrated hop and to rescale along in the case of a hop, the errors are very similar to the errors observed for the corresponding rescaling option. This does not hold true for not reflecting hops at all which was combined with all possible rescaling options. The associated average error for not reflecting frustrated hops at all is found at almost the same value as the average over all calculated ε values , indicating that this parameter has no influence here. The analysis of the sets resulting in larger ε values reveals that AFSSH dynamics in the diagonal representation results in too fast transfer from the S 2 to the 1 and increases population of the T 1 state, thus giving a bad description. This effect is drastically enhanced when rescaling along h or g.

2-thiocytosine dynamics in the presence of a laser field
Following the framework presented in Section V A 2, a set of 7 lasers that differ in their t p and their maximum amplitude have been applied to excite population from the ground state and initiate dynamics. First, the influence of different lengths of laser excitation has been investigated. Figure 9a presents the results obtained for MCTDH. Although the pulses are tuned to carry the same amount of overall energy, the amount of excited population increases with the pulse length. Accordingly, the longest pulse achieves almost ground state population inversion. The same trend is achieved by the FSSH dynamics using the DIAG EDC + v setup, see Fig 9b, but in a much weaker extent. This is in part due to the the EDC acting against the excitation process. The hampering effect of the EDC is more pronounced than in SO 2 because the strength of this decoherence correction is based on the kinetic energy of the system, which on average is much larger for 2-thiocytosine than SO 2 due to size.
Intriguingly, the MCTDH S 0 populations do not show any surprises for longer pulses -contrary to what was observed for SO 2 where the recurrence of the wave packet resulted in additional interferences. This is investigated in more detail in Fig. 9c where ten different model systems containing one to ten normal modes of 2-thiocytosine (sorted according to their importance for the overall dynamics, see Table SX  The t p = 30 fs pulse duration was selected to estimate the performance of the different FSSH parameter combinations on the 10 mode 2-thiocytosine model. As can be seen in Fig. 10, a very broad distribution of total errors is obtained for the overall error (see Tables SXV-XVIII in the supporting information for the complete lists). This is mostly due to the high deviation observed in the S 0 state population, yielding ε S 0 values larger than 0.2 for all sets.
The description of the S 0 population was found to be rather independent of the representation but strongly depending on the presence of a decoherence correction, with better results when this is considered. Rescaling the momenta after a surface hop impacts the S 0 population, with not rescaling the energies resulting in the bigger disagreement with MCTDH. The overall increase in deviation from the MCTDH S 0 population when compared to the deviations in SO 2 can partially be attributed to the use of the rather strong laser pulse that inverts 90% of the population in the MCTDH reference. Indeed, strong laser fields have been previously found to increase the deviation from quantum dynamics reference results. 55 In that paper it was reported that the excitation process in general results in a net vibrational cooling of the swarm of excited trajectories with respect to the average kinetic energy of the ground state population, resulting in slower movement of the FSSH trajectories away from the region where the laser couples the bright and the ground state.
In the presence of a strong laser field, dwelling in the region of strong coupling increases the chance of inducing radiative emission to the ground state, effectively decreasing the amount of excited population. This poses an important dilemma for carrying out FSSH dynamics with laser pulses when generating initial conditions: On the one hand side, the stronger the pulse, the more trajectories will be excited, thereby decreasing the amount of unexcited trajectories that are now obsolete for describing excited state dynamics and thus represent computation time essentially wasted. On the other hand, stronger lasers will lead to larger deviations from the exact dynamics due to the cooling effect and the increased interplay between different parts of the wave packet that are more strongly coupled.
For the excited state dynamics following excitation, ε r sing shows a rather narrow distribution. Surprisingly, a complete change in the observed order of importance for the various parameters can be detected with respect to the dynamics without any external field. Two effects are responsible for the presented errors: First, the influence due to the fast oscillating coupling element with the S 0 on the ε values is reduced as the S 0 component of this error is included in ε S 0 and therefore only taken into account via the S 2 state. Second, the fine structure is smoothed out, giving rise to simpler decay patterns as can be seen in Fig. S 3 in the supporting information. Still, it comes as a big surprise to find that the performance of the EDC is reduced and that rescaling of the kinetic energy after a hop along g in the MCH representation gives by far the best renormalized singlet populations. The latter behaviour is the opposite of what was observed in absence of a laser field, i.e. this rescaling option was found with the largest average ε sing . When including an explicit laser pulse, rescaling along g gives good agreement with MCTDH independent of the applied decoherence correction if the MCH representation is used. However, a strong dependence on the applied decoherence correction in the diagonal representation can be seen where DIAG EDC −g g and DIAG EDC + g are among the largest deviations from the MCTDH results for all sets observed.
The deviation in the populations of the triplet states follows the same trends found in the excited state dynamics without an external field, with a decrease in the importance of the chosen representation and an overall better performance of the AFSSH. Rescaling along the nonadiabatic coupling vector is still found to decrease the agreement with the MCTDH reference.
Summarizing this section, modelling the excitation of S 0 population by an explicit laser pulse with t p = 30 fs in the 10-D 2-thiocytosine model system using FSSH results in deviations of about 30%. Differences due to the representation are negligible while using a decoherence correction was found to be an important factor, similar to observations in SO 2 . For the excited state dynamics following excitation, rescaling of the energies after a hop was found to be important, this time favoring rescaling along the gradient difference vector for the dynamics in the singlet manifold.
Transfer of population towards the triplet states is found to be only slightly influenced by the decoherence treatment and the chosen representation.

VI. CONCLUSIONS
In this work, we benchmark the performance of different parameters that are needed in FSSH against a MCTDH reference for two molecular systems: SO 2 and 2-thiocytosine using parametrized LVC potentials. We investigate the effect of such parameters on simulations that do not explicitly include a laser field and compare with dynamics initiated by laser pulses of different time duration.
Previous work 5,6,12,33,55 has reported the inability of FSSH to correctly describe the excitation process and the subsequent dynamics in the presence of explicit light-matter interactions, as well as pointed out the dependence of the results on the chosen FSSH parameters. In this work, we quantitatively measure the deviation of FSSH from real wave packet dynamics, and by studying the system with and without the influence of external fields, we discern between effects coming from the dynamics itself or from the laser interaction. Furthermore, the two chosen models SO 2 and 2-thiocytosine, with 3 and up 10 normal modes, are beyond the commonly 1-dimensional test systems employed to investigate the interaction with an external field. The consideration of both singlet and triplet states pose an additional difficulty for FSSH as both very weak, delocalized coupling between singlets an triplets as well as the time-dependent coupling via an external field have to be treated correctly.
We find that FSSH has difficulties in low-dimensional models if long laser pulses (FWHM>100 fs) are present, where the excited state wave packet returns to the Franck-Condon region within the time scale of the applied laser pulse. MCTDH simulations in SO 2 indicate that the returning wave packet interacts with the active laser pulse, resulting in an additional interference, for laser pulses of at least t p = 100 fs. FSSH is not capable to reproduce these interference terms and vMCG can ony capture this effect when going beyond 100 GBFs. Moving to larger model systems decreases the recurrence time until it is not observed with the longest (t p =200 fs) laser pulse in the 10-mode model of 2-thiocytosine, indicating that these effects will be absent in systems with more degrees of freedom.
In SO 2 , no perfect set of parameters was found suitable to describe quantitatively all different transition processes triggered by an ultrafast (t p =30 fs) laser pulse or when the simulation starts directly from the bright state. The decoherence correction emerges as the most essential parameter to model the excitation process and therefore also the overall laser-induced dynamics in the presence of laser pulses longer than a few cycles. The MCH representation where both the coupling via a laser field as well as the spin-orbit coupling are treated as off-diagonal elements, is found beneficial to describe adequately both the S 0 and the triplet populations. This is mainly due to a better performance of AFSSH in the MCH representation to describe the change in the For 2-thiocytosine, larger differences in the amount of excited population are observed due to the strong laser applied, increasing the difference in excited population by 35% and more, almost irrespective of the chosen parameters. The use of a decoherence correction increases the overall agreement throughout the excitation process, similar to the observations in SO 2 . In the subsequent dynamics, rescaling of the kinetic energy after a hop reveals itself as the most important factor to obtain lower errors with respect to the MCTDH populations. Interestingly, the gradient difference vector is found to mostly result in large deviations from the reference if no laser is present but becomes the best performing parameter in the dynamics including the t p = 30 fs laser pulse.
In conclusion, this work shows the difficulties in choosing a universal set of parameters that guarantees quantitative agreement against quantum reference results, particularly in the presence of an electric field. However, and despite many difficulties, FSSH can qualitatively reproduce the dynamics ensuing after excitation using an exemplary laser pulse of FWHM=30 fs on two different test systems, with three and ten dimensions -an encouraging result to pave the way to use of explicit laser pulses in FSSH simulations. Based on the observations made in this paper, a few caveats have been identified. The use of very long laser pulses is discouraged due to the increased chance of quantum interferences occurring, especially for small systems. The same holds true for the use of very strong laser pulses, which will enhance problems inherent to FSSH without any laser pulse. Additionally, it is advised to use a decoherence correction (from the many available) because it is demonstrated to be in most cases the dominating factor to increase the agreement with the MCTDH reference, especially when considering explicit excitation. Caution is also advised when it comes to choosing a vector to adjust the momenta after a non-laser induced hop, as this parameter strongly influences the presented dynamics.          . An alternative solution would be the use of multi-layer MCTDH which has been employed successfully in the investigation of even larger molecules ? ? .

C. vMCG deviation
In the SHARC-gym approach, FSSH dynamics are used to determine vibrational normal modes that are essential to describe the excited state dynamics, which can then be used to form a lowerdimensional Hamiltonian still capable of describing the overall excited state dynamics. For this, a reference FSSH simulation using the full-D LVC model employing the DIAG EDC + v set of parameters is created, where the initially excited state was chosen to be the second excited state at the reference geometry. All initial electronic populations and initial states of each of the 1000 trajectories was adapted to start in the corresponding state for each initial geometry. Subsequently, 33 new simulations, each employing 1000 new trajectories starting in the same electronic state but each neglecting a different single vibrational mode (all couplings and gradients associated with this mode in the LVC template were set to zero) have been conducted. For each of these dynamics, ε values according to Equation 9 in the main manuscript with respect to the full-D reference simulation have been calculated as listed in TABLE X and can be seen in FIG 2. As can be seen there, most modes are associated with ε values below 0.05 and therefore result in almost identical dynamics. After sorting the modes based on the associated error that is introduced when this mode is omitted, an arbitrary cut-off of 10 modes has been chosen that will be included in the reduced

B. Population dynamics with laser fields
The influence of employing longer laser pulses to excite 2-thiocytosine are shown in FIG 3 where the top two rows show the population evolution for MCTDH (first row) and the DIAG EDC + v set of parameters (second row). The remaining columns show the influence of renormalization using the S 0 population in every time step (see Section S 1 A for a more detailed explanation).  where each excited state population is adapted according to p r i = p i 1−p 1 where p 1 denotes the ground state population and p i the population in the ith excited state in a-c). j-l) Renormalized excited state populations for the FSSH dynamics shown in d-f).