The capability of fewest-switches surface hopping (FSSH) to describe non-adiabatic dynamics under explicit excitation with external fields is evaluated. Different FSSH parameters are benchmarked against multi-configurational time dependent Hartree (MCTDH) reference calculations using SO_{2} and 2-thiocytosine as model, yet realistic, molecular systems. Qualitatively, FSSH is able to reproduce the trends in the MCTDH dynamics with (also without) an explicit external field; however, no set of FSSH parameters is ideal. The adequate treatment of the overcoherence in FSSH is revealed 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) correction and the energy-based decoherence correction. A dependence on the employed basis is detected in AFSSH, 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 laser pulses (>100 fs), resulting in additional interferences missed by FSSH and only visible in variational multi-configurational Gaussian when utilizing a large number of Gaussian basis functions. This feature vanishes when going toward larger molecules, such as 2-thiocytosine, where this effect is barely visible in a laser pulse 200 fs long.

## I. INTRODUCTION

Femtosecond time-resolved spectroscopy has progressed drastically throughout the past decades,^{1} challenging the computational excited state dynamics simulations to explicitly include laser pulses.^{2–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 the consideration of 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 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. Present in this category are the Ehrenfest,^{14} *ab initio* multiple spawning^{15} (AIMS), or fewest-switches surface hopping^{16,17} (FSSH) dynamics. These methods show a favorable scaling with the molecular size at the cost of some quantum effects and thus accuracy in comparison to methods that more readily converge to the exact result like multi-configurational time-dependent 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–25} and the quantum–classical Liouville equation^{26–28} have shed light on the nature of FSSH and make clear the use of the non-adiabatic coupling (NAC) vector as the proper direction for adapting the nuclear velocities at a hopping event, 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,^{29,30} are treated by *ad hoc* corrections that while remedying the issue at hand in some test cases fail in others. The resulting abundance of the available options with different parameters to choose from—some of which should be more accurate while others aim to extend the applicability of this method to larger and more complex systems—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 applied in real systems, can dramatically affect dynamics.^{31}

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.^{2,3,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 $H2+$^{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-thiocytosine, 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 (also 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 this paper is organized as follows: In Sec. II, the framework of FSSH and the various options available therein are presented. Dynamics methods that go beyond FSSH are introduced in Sec. III. Section IV describes the nomenclature employed and the methodological details. Finally, the numerical results are presented in Sec. V followed by the conclusions in Sec. VI.

## II. FEWEST-SWITCHES SURFACE HOPPING BACKGROUND

Since its original formulation^{16} as an improvement to the existing surface hopping methodologies,^{34,35} FSSH has seen a multitude of adjustments to overcome some of its inherent limitations.^{17,36–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.^{30} 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 background 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

where *c*_{α} are the coefficients for each electronic state. Their time dependence is given by

Here, **h**_{βα}(**R**) is the non-adiabatic coupling (NAC) vector $\u27e8\Phi \beta (r,R)|\u2202\u0124el\u2202R|\Phi \alpha (r,R)\u27e9$ that indicates the change in the electronic wave function with variation of the nuclear coordinates. $H\beta \alpha d$ is the matrix element of a complete Hamiltonian that contains any arbitrary coupling and can be written as

The dipole operator $\mu ^\beta \alpha (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 the 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 that 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 PESs when the active state changes. Multiple algorithms to determine 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 *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 such as 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 upside 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 representations—are 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. Solving the electronic Schrödinger equation with the MCH yields a set of diagonal and non-crossing states within each multiplicity. Additional coupling elements such as SOCs or external fields can then be included as off-diagonal elements in the MCH picture. This MCH set of states can be transformed to a new non-diagonal set of states by a diabatization.^{43} The PES of diabatic states can be chosen to be smoothly varying and can cross without showing avoided crossings. For a polyatomic molecule, no unique diabatic transformation exists^{44} and the transformation should be chosen as to minimize kinetic-energy coupling and retain a set of chemically relevant states obtained at a reference geometry or to other relevant observable. A third representation is the one considered in the surface hopping including arbitrary couplings (SHARC) approach:^{32,42} here, the complete Hamiltonian **H**^{d} with matrix elements $H\alpha \beta d$ as in Eq. (3) is diagonalized. 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 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 represents 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 $H2+$ 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, 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 with both 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 past 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 *i*th 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 ($pi\u2032$) 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 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 the 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 in 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. The use of the **h** vectors as the best possible option has subsequently been confirmed through work connecting the FSSH algorithm to the quantum–classical Liouville equation.^{27,28} In the case where NAC vectors are either unavailable within the employed electronic structure method or computationally limited, 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. It has been found that the differences between using **h** and **g** are negligible in the case of large nuclear momenta in some systems but become more apparent when these are small,^{53,54} highlighting the more rigorous derivation of the **h** vector to adjust the kinetic energy after a hop. An even more approximate method is to rescale the full nuclear momentum, which offers an intriguing simplicity where no additional properties need to be computed.

Hops to a lower PES result in an increase in the 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.^{55,56} 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 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 no 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.^{31,56} Additional care has to be taken when using the NAC vectors to adapt the nuclear momenta when a complex Hamiltonian is employed. In this case, additional measures have to be employed to guarantee a meaningful choice of the phase, which is not arbitrary anymore.^{57}

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.^{58} 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 labeled field-induced and, thus, no rescaling of the velocities is performed.

## III. REFERENCE METHODS

### A. Multi-configurational time-dependent Hartree

Introduced in 1990, the MCTDH^{18,59,60} method is one of the most versatile methods to simulate non-adiabatic processes due to the inherent possibility to converge toward 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, 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 therefore relies on a functional form of the PES. Subsequent, alternative formulations have been proposed, resulting in the development of the vMCG^{19,20} method, where frozen Gaussians are used as basis functions instead of SPFs.

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) $|\chi j(\alpha )\u3009$ and time-dependent expansion coefficients $Aj(s)$. As in MCTDH, the equations of motion can be derived variationally and result in the propagation of both the expansion coefficients and the parameters of the GBFs. The coupled motion of both the coefficients and the GBF parameters results in a “quantum” movement that goes beyond simple classical motion. Also similar to MCTDH, in the limit of infinite basis functions (SPFs for MCTDH and GBFs in vMCG), the complete space is covered in basis functions and the exact dynamics is obtained. The vMCG algorithm has the advantage that it can be employed in an on-the-fly fashion without relying on precomputed PESs.^{61}

## IV. COMPUTATIONAL DETAILS

### 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:

The possible options for each of the keywords are listed in Table I and explained as follows:

Representation: two representations for the propagation are used, the MCH or the completely diagonal picture (DIAG) previously introduced in Sec. II A.

Decoherence: to correct for overcoherence, we use AFSSH or EDC (see Sec. II B). Additionally, the NONE option indicates that no decoherence correction is included.

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, therefore violating the energy conservation within each trajectory (none).

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 canceled and termed frustrated. 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 labeled “+.”

Representation . | Decoherence . | Rescaling . | Frustrated hops . |
---|---|---|---|

MCH | AFSSH | v | −v |

DIAG | EDC | h | −h |

NONE | g | −g | |

none | + |

Representation . | Decoherence . | Rescaling . | Frustrated hops . |
---|---|---|---|

MCH | AFSSH | v | −v |

DIAG | EDC | h | −h |

NONE | g | −g | |

none | + |

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 $AFSSHh+DIAG$. Note that not all combinations of these parameters result in stable setups for FSSH dynamics: for example, in the MCH representation, the non-adiabatic coupling vector **h** between singlet and triplet states is zero and therefore cannot be used as a direction to rescale the velocities along in simulations including singlet and triplet states. Apart from these incompatibilities, all possible combinations between options (i)–(iii) have been used. Throughout this work, a consistent pairing of the used vectors in (iii) and (iv) has been used: This means 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, and *p*_{i,t} and *p*_{i,t,ref} are the populations of state *i* taken from the investigated and the reference dynamics at time *t*, respectively. If the populations of the reference and the investigated dynamics are identical for all time steps, *ε* equals zero. The maximum value of this error is achieved only if for every time step of the dynamics, all of the population in the investigated dynamics is found in states that are not populated at all in the reference dynamics. For example, if all of the population in the reference dynamics is in state *a*, while the investigated dynamics always run in state *b*, the finale error equates to $|pa\u2212pa,ref|+|pb\u2212pb,ref|+\u2211i\u2260a,bnstates\u22122|0\u22120|=2$. 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 ($\epsilon S0$) 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 $pi,t1\u2212pS0,t$ or $pi,t,ref1\u2212pS0,t,ref$ 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 (also 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. 62 has been employed.

### C. Dynamical propagations

The FSSH simulations have been carried out with the SHARC^{39,63} program suite using pySHARC^{64} to drastically reduce I/O overhead. A set of 1000 nuclear initial conditions obtained from a ground state Wigner sampling^{65} 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 and 0.002 fs in the field-free and field-including case, respectively.^{66} 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 linear vibronic coupling (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.^{67} In MCTDH, the Adams–Bashforth–Molton predictor–corrector integrator of sixth 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 and, 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 and time steps, where ⟨*n*⟩ and ⟨*dn*⟩ are the position expectation value and the corresponding standard deviation, respectively.

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 fourth 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 lowest-energy 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 supplementary material.

All FSSH, MCTDH, and vMCG simulations were run for 400 fs.

### D. Definition of the laser pulse

The coupling with the external field is described within the semi-classical dipole approximation and can be written as $\u2212\mu ^\beta \alpha (R)\epsilon (t)$, neglecting any further interaction terms. The pulse has a Gaussian shape, defined as

with field amplitude $\epsilon tp0$, 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 corresponding molecules: 4.49 eV for SO_{2} and 3.92 eV for 2-thiocytosine. Seven different values for *t*_{p} are used with the corresponding centers of the pulse (*t*_{0}) in parentheses: 2 fs (*t*_{0} = 10 fs), 10 fs (*t*_{0} = 30 fs), 17 fs (*t*_{0} = 40 fs), 30 fs (*t*_{0} = 70 fs), 50 fs (*t*_{0} = 90 fs), 100 fs (*t*_{0} = 140 fs), and 200 fs (*t*_{0} = 200 fs).

The field amplitude of the laser field for a given *t*_{p}, $\epsilon tp0$, was varied for different lengths of the laser pulse according to $\epsilon tp0=\epsilon 170\u22c517tp$. In 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 were set to $\epsilon 170=0.03$ a.u. (15.44 GV/m) for SO_{2} and $\epsilon 170=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.^{68} 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. 64, which also demonstrates that the LVC approximation is able to reproduce the main characteristics of the corresponding on-the-fly dynamics. Note that all spin orbit coupling elements entering the employed LVC model Hamiltonian have no imaginary parts.

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 past decade, illustrating a complex interplay between singlet and triplet states.^{69–71} In this work, we use a LVC Hamiltonian that contains four singlet and three 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.^{64} 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.

#### 1. 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. 1(a)], which will serve as a reference throughout, shows ultrafast population transfer from the initially populated ^{1}B_{1} state to the ^{1}A_{2} state. Within the first 50 fs, there is almost complete depletion of the ^{1}B_{1} population, which then oscillates on par with that of the ^{1}A_{2} during the rest of the propagation. This oscillatory behavior is due to the closeness of the respective minima in both energy and phase space. Therefore, only limited stabilization of one diabatic state with respect to the other can take place, resulting in the repeating pattern that gets more complex over time. During the first 400 fs, about 10% of the population crosses to the triplet manifold, where the ^{3}B_{2} state is populated almost exclusively—in line with the also MCTDH simulations performed on *ab initio* potentials by Lévêque *et al*.^{69}

The vMCG dynamics using 75 GBFs [Fig. 1(b)] show only small differences compared to the MCTDH reference. The excellent agreement of the vMCG singlet populations with the MCTDH ones illustrates the strength of the method. The oscillating behavior between the ^{1}B_{1} and the ^{1}A_{2} 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 toward the exact result (see Sec. S1 C of the supplementary material 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. 1(c) only the populations for a single parameter set: $EDCv+DIAG$. 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 with this simulation gives *ε* = 0.248, with the damped oscillations and the transfer to the ^{3}B_{2} state being the main sources of error. Unsurprisingly, FSSH has a larger error than vMCG (approximately three times larger).

In order to investigate the influence of different FSSH parameters, all sensible combinations of parameters according to Sec. IV A have been used to run FSSH dynamics. Figure 2 collects all the deviations against the MCTDH reference for each set. See Tables SI–SIII of the supplementary material for a list of all obtained *ε* values.

The *ε* ranges from 0.240 to 0.383, meaning that an unsuitable combination of surface hopping parameters can increase the error with respect to the best possible set by 50%. Each column of Fig. 2 gives the stacked error summed over all the states, singlet and triplet states, respectively, while each row depicts the influence of a particular option in the FSSH algorithm.

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 to the active state. Yet, both of these parameters 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. This is surprising, as adapting the kinetic energy along the NAC vectors was considered the option closest to real quantum dynamics, both from a theoretical point of view^{27,28} and in a set of practical applications.^{31,56} A possible explanation 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 the option of rescaling along the complete velocity vector and, of course, the option of not adapting 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 vs 3833 frustrated hops is found for $EDCg+MCH$, while 7484 hops against only 1041 frustrated hops are found for $EDCv+MCH$ on the complete set of 1000 trajectories. Thus, the higher available kinetic energy in the complete kinetic energy vector enables a lot of jumps that are otherwise forbidden when the non-adiabatic coupling or the gradient difference vector for rescaling is used. This increased number of frustrated hops might account for the discrepancies detected when rescaling the velocities along **h** or **g**, despite the fact that these options are both more physically sound and should have resulted in better adapted velocities after a successful hop than the full velocity vector or not rescaling the velocities at all. A possible remedy is the use of a modified FSSH algorithm that employs a time-uncertainty scheme, which can allow hops otherwise forbidden—as proposed by Truhlar and co-workers.^{72} We noted that switching to the diagonal basis increases the total number of observed hops drastically (39 908 for $EDCg+DIAG$ and 41 738 for $EDCv+MCH$) 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 ($EDCg+DIAG$) and 859 ($EDCv+DIAG$). 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 that “few frustrated hops equal 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.

The findings can be summarized as follows: (i) The MCH picture results in a stable description of the triplet population. (ii) No real favorite decoherence correction emerges. (iii) Rescaling along one of the physically more sound options, **h** or **g** results in larger deviations from the MCTDH dynamics, especially when paired with the parameter of not reflecting frustrated hops. Good results were achieved when using **v** for rescaling the kinetic energy combined with a continuation along the current velocity when encountering frustrated hops. (iv) A very bad combination of parameters for the description of the triplet states emerges in the form of using a decoherence correction in the diagonal representation and rescaling along **h** or **g** using non-reflected frustrated hops.

#### 2. 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 are 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, note that laser pulses of different lengths employ different field amplitudes $\epsilon tp0$ [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 a *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 ($\epsilon tp0$) 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 toward 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 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 and *t*_{p} = 200 fs lasers reveals that the diverging behavior is caused by 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 acts as a probe, detecting whether the returning wave packet causes interference terms. The resulting populations for different *τ* delay times are shown in Fig. 4(a). 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, and 97 fs. From this, it can be concluded that after ∼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 shown in Fig. 3(b), the S_{0} populations do not show this interference term anymore. Instead, the laser pulse induces a 55% S_{0} population inversion.

Figure 3(c) demonstrates that vMCG dynamics is in 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 and *t*_{p} = 200 fs lasers is 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. 4(b)]. 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 for agreement. However, in simulations employing long laser pulses with *t*_{p} = 100 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 created in the pump–probe setup. Exploratory calculations employing larger numbers of GBFs indicated that the reversal of the population flow occurring for the S_{0} population of the *t*_{p} = 200 fs dynamics needs more than 125 GBFs to be even visible. Unfortunately, calculations using this and larger number of GBFs could not be converged.

When comparing the MCTDH and vMCG dynamics to the FSSH dynamics using the $EDCv+DIAG$ set of parameters [see Fig. 3(d)], a general underestimation of about 5%–10% in the amount of excited population is observed for laser lengths up to *t*_{p} = 50 fs. The reduced amount of excitation is mainly due to the employed decoherence correction. Different from other decoherence corrections such as AFSSH, the EDC acts instantaneously, meaning that in every time step of the simulations, the electronic coefficients of non-active states will be damped and the coefficient of the active state increased. During the excitation process, the applied external field will periodically and continuously increase the electronic coefficient of excited states and decrease the coefficient in the ground state, thus enabling the chance for the trajectory to perform a surface hop. The EDC therefore directly counteracts the influence of the laser field, resulting in a lower amount of excited population. Interestingly, applying the same amount of energy over a longer time did not decrease the performance of the EDC correction further. Similar to vMCG, the S_{0} populations in FSSH do not show any significant irregularities when moving to the *t*_{p} = 100 fs and the *t*_{p} = 200 fs lasers and instead correspond closely to the populations observed in MCTDH when employing a CAP.

A screening for all combinations of surface hopping parameters considered in Sec. V A 1 has been conducted using the *t*_{p} = 30 fs laser pulse, where no additional interference due to the 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, $\epsilon S0$, 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 into the contributions of singlets and triplet states. The total calculated error is then just the sum over these three contributions: $\epsilon =\epsilon S0+\epsilon singr+\epsilon tripr$. All calculated *ε* values for all states, the S_{0}, the excited singlet state, and the excited triplet state, are presented in Fig. 5 and listed in 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 $\epsilon S0$, 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 toward 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 $\epsilon S0$ between 0.058 and 0.076 and therefore show the lowest $\epsilon S0$ values. However, when using AFSSH in combination with the diagonal representation, the errors essentially double to a range of 0.123–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 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 $\epsilon singr$ as compared to the dynamics in the absence of a laser field. Almost all parameters give identical averages, with these averages also being smaller than those 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. 1(a) 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 shown in Fig. S1 of the supplementary material. 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 a 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 is the driving factor 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 improve 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 the absence of a laser field: $NONEv\u2212vMCH$ (best average *ε* and *ε*_{trip}) and $EDCv\u2212vDIAG$ (best *ε*_{sing}). For dynamics initiated using a *t*_{p} = 30 fs laser, the combinations $EDCv\u2212vMCH$ (best *ε*), $EDCg\u2212gMCH$ ($\epsilon S0$), and $EDCv+MCH$ ($\epsilon singr$) were added, as the on-average best performing options resulted again in $NONEv\u2212vMCH$ for $\epsilon tripr$. These sets were then used to simulate dynamics using different laser lengths for excitation varying in size from *t*_{p} = 2 to *t*_{p} = 50 fs and $\epsilon S0$ 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 and *t*_{p} = 17 fs before decreasing for lasers *t*_{p} = 30 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 representation. 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 ($NONEvvMCH$). 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,^{73} 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 toward the triplet states after 500 fs, but to maintain all other essential features.^{64} With this in mind, in order to test the effect of different FSSH settings, here, we consider simulation times of up to 400 fs.

#### 1. 2-Thiocytosine dynamics in the absence of a laser field

Figure 7(a) shows FSSH dynamics including all normal modes in conjunction with the $EDCv+DIAG$ set of parameters in the absence of a laser field. The ground state Wigner distribution of geometries has been excited into the bright diabatic S_{2} state from where subsequent dynamics takes place. A strong coupling of this bright state with S_{0} results in a fast oscillating transfer between these two states, recovering population in the S_{0} state. At the same time, internal conversion to the diabatic S_{1} and intersystem crossing to the triplet manifold is observed. After 400 fs, both the S_{1} and the S_{2} have similar populations, and a higher population is found in the diabatic T_{1} state.

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 SHARC-gym.^{74} 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. 7(a) 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 Sec. 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.^{75,76} Following the SHARC-gym, we selected the ten most important modes that can reproduce the 33-dimensional dynamics of 2-thiocytosine. This reduced ten-dimensional set will be now used for the remainder of the simulations. As shown in Fig. 7(b), the FSSH dynamics using this ten-dimensional Hamiltonian captures the overall behavior of the full-dimension very reasonably.

With this reduced model at hand, MCTDH calculations have been converged yielding the populations presented in Fig. 7(c). Compared to the $EDCv+DIAG$ presented set of parameters, the MCTDH simulations predict larger S_{1} and smaller T_{2} population. The amount of S_{0} population is lower in MCTDH due to the inability of EDC to describe fast oscillations as it tries to end up in pure states. Overall, the agreement between both methods is nevertheless quite reasonable with an error *ε* = 0.214 for the $EDCv+DIAG$ set.

Using the ten-dimensional 2-thiocytosine model, we obtain the distribution of *ε* values presented in Fig. 8 (see also Tables SXII–XIV). The errors are spread wide, ranging from 0.163 to 0.488, with the largest accumulation of parameter sets giving very good results. Similar to SO_{2}, the MCH picture gives a better description of the triplet populations, while most simulations using the diagonal picture perform worse. The EDC decoherence treatment performed well both for the singlet and the triplet populations. Using no decoherence correction resulted in somewhat surprisingly good results when compared to the AFSSH, which in turn gave larger *ε* values on average than the other treatments. When including the effect of rescaling of the kinetic energies after a hop, the driving factor in the current dynamics becomes apparent: rescaling along the **h** and the **g** vectors is associated with the largest errors, while neglecting the rescaling completely or rescaling along the velocity vector gives very good agreement with the MCTDH populations. Interestingly, there are two exceptions to this, as the $EDCg+MCH$ and the $EDCggMCH$ 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 the population of the T_{1} state, thus giving a bad description. This effect is drastically enhanced when rescaling along **h** or **g**.

#### 2. 2-thiocytosine dynamics in the presence of a laser field

Following the framework presented in Sec. V A 2, a set of seven 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 9(a) 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 $EDCv+DIAG$ setup [see Fig. 9(b)] but in a much weaker extent. This is, in part, due to 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. 9(c) 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 for the explicit modes) are used to probe the interaction with the *t*_{p} = 200 fs laser pulse. With very few modes, almost no population remains in the excited state, as the recurrence of the wavepacket occurs very fast and further excitation is hampered and even reversed, similar to what occurred in SO_{2}. The more modes are added to the system, the more excited state pathways open up, increasing the time needed until the excited wave packet returns to the Franck–Condon region, reducing the transfer back to the S_{0} population as the laser duration is shorter than the recurrence time. The increase in recoherence time shows that this effect is only affecting simulations of very small systems and it is almost undetectable for a ten-mode system with a rather long laser pulse of *t*_{p} of 200 fs.

The *t*_{p} = 30 fs pulse duration was selected to estimate the performance of the different FSSH parameter combinations on the ten mode 2-thiocytosine model. As shown in Fig. 10, a very broad distribution of total errors is obtained for the overall error (see Tables SXV–XVIII in the supplementary material for the complete lists). This is mostly due to the high deviation observed in the S_{0} state population, yielding $\epsilon S0$ 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 dependent 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 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.^{58} 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 ground states. 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 one hand, 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 the 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, $\epsilon singr$ 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 $\epsilon S0$ 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 shown in Fig. S3 in the supplementary material. 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 behavior is the opposite of what was observed in the 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 $EDCg\u2212gDIAG$ and $EDCg+DIAG$ 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 non-adiabatic coupling vector is still found to decrease agreement with the MCTDH reference.

Summarizing this section, modeling the excitation of S_{0} population by an explicit laser pulse with *t*_{p} = 30 fs in the ten-dimensional 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 toward 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 durations.

Previous work^{5,6,12,33,58} has reported the inability of FSSH to correctly describe the excitation process and the subsequent dynamics in the presence of explicit light–matter interactions and 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 molecular systems 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 three and up to ten normal modes, are beyond the commonly one-dimensional test systems employed to investigate the interaction with an external field. The consideration of both singlet and triplet states poses an additional difficulty for FSSH as both very weak delocalized coupling between singlets an triplets and the time-dependent coupling via an external field have to be treated correctly.

A set of four parameters with different options available in each of those have been tried and tested. The various employed options span from exact procedures (such as adapting the nuclear velocities after a hop along the non-adiabatic coupling vector) to extremely crude approximations (e.g., not adapting the energies at all at a hopping event). This incomplete variation in parameters does not aim at providing a single set of “best” performing options, of which the conducted comparison of a single electronic property is not capable of. Instead, it should illustrate how diverse the field of surface hopping parameters and users thereof is, and call for caution when choosing and combining surface hopping parameters.

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 of reproducing these interference terms, and vMCG can only 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 ten-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 and 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 S_{0} population as compared to AFSSH in the fully diagonal representation. For the EDC, no such strong basis dependence was identified for the excitation process, but larger disagreements for the triplet states are observed. For the observed dynamics in the excited singlet states, most parameters were found with almost identical errors. When using a reasonable combination of FSSH parameters, similar deviations are obtained for the excitation process for different lengths of the laser pulse until the wave packet recurrence is observed. Finally, the importance of rescaling the kinetic energy after a hopping event (not influenced by the laser field) is evident in the excited state dynamics, with rescaling along both the non-adiabatic coupling vector and the gradient difference vector showing high deviations for the triplet state dynamics. The performance of rescaling along both the NAC and the gradient difference vectors to mimic the MCTDH dynamics suffers from a large amount of frustrated hops but would most probably be improved by including additional nuclear properties in the comparison. Overall, the best FSSH parameter sets are associated with ∼20% deviation in the amount of excited S_{0} population with respect to MCTDH, while the subsequent population evolution is described with a similar level of accuracy as without the laser pulse.

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 summary, 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 that paves the way for the 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 short laser pulses is encouraged due to the small amount of coupled electronic-nuclear motion during the pulse duration, thereby reducing the laser field to an almost purely electronic effect where FSSH yields good results. Additionally, quantum interferences are drastically reduced when employing short pulses, especially for larger 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.

## SUPPLEMENTARY MATERIAL

See the supplementary material for examples of population evolution in the presence of laser fields, highlighting the renormalization conducted within Eq. (10); sorted lists containing all calculated *ε* values used in Figs. 2, 5, 8, and 10; all *ε* values calculated for SO_{2} using vMCG with different numbers of basis functions (which are presented in the PDF file); input and operator files used to simulate all MCTDH and vMCG calculations in QUANTICS; and Molden files and LVC template files for SO_{2} and 2-thiocytosine (full-dimensional and one- up to ten-dimensional) use in SHARC (which are presented in the ZIP file).

## ACKNOWLEDGMENTS

The authors thank the University of Vienna for providing access to the Vienna Scientific Cluster where part of the presented simulations was computed. The authors also thank S. Mai and P. Marquetand for discussions on the various flavors of FSSH.

## DATA AVAILABILITY

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

## REFERENCES

^{+}with D

_{2}

_{2}

^{3}B

_{2}state in the dynamics of sulfur dioxide

_{2}. III. An ab initio quantum study of single- and multi-photon ionization

_{2}. II. The role of triplet states in the bound state dynamics studied by surface-hopping simulations