Nonadiabatic dynamics in multidimensional complex potential energy surfaces

Abstract

Despite the continuous development of theoretical methodologies for describing nonadiabatic dynamics of molecular systems, there is a lack of approaches for processes where the norm of the wave function is not conserved, i.e., when an imaginary potential accounts for some irreversible decaying mechanism. Current approaches rely on building potential energy surfaces of reduced dimensionality, which is not optimal for more involving and realistic multidimensional problems. Here, we present a novel methodology for describing the dynamics of complex-valued molecular Hamiltonians, which is a generalisation of the trajectory surface hopping method. As a first application, the complex surface fewest switches surface hopping (CS-FSSH) method was employed to survey the relaxation mechanisms of the shape resonant anions of iodoethene. We have provided the first detailed and dynamical picture of the π*/σ* mechanism of dissociative electron attachment in halogenated unsaturated compounds, which is believed to underlie electron-induced reactions of several molecules of interest. Electron capture into the π* orbital promotes C=C stretching and out-of-plane vibrations, followed by charge transfer from the double bond into the σ* orbital at the C–I bond, and, finally, release of the iodine ion, all within only 15 fs. On-the-fly dynamics simulations of a vast class of processes can be envisioned with the CS-FSSH methodology, including autoionisation from transient anions, core-ionised and superexcited states, Auger and interatomic coulombic decay, and time-dependent luminescence.

---

1 Introduction

Computational methods for nonadiabatic dynamics simulations comprise a fundamental theoretical tool for surveying the relaxation of electronically excited molecules. Basically, an approximate solution for the electronic problem provides the potential energy surfaces (PESs) that guide the nuclear dynamics. At the same time, the coupling between the motion of light and heavy particles is achieved in one of several ways, depending on the particular mixed quantum-classical approach. Applications of these methodologies are countless, as can be attested by three recent reviews of the topic.^1–3

Despite their success, the current formulations cannot be applied to open quantum systems, where the set of discrete states of the system interacts with a continuum of scattering states.^4,5 This class of problems requires the use of non-Hermitian quantum mechanics, which is a vast ever-growing topic on itself,^5,6 going beyond the molecular dynamics subject that concerns us here. But in short, the effect of the continuum of states can be modelled with imaginary potential terms, which mimic the decay from the discrete states and account for the decrease of the system wave function. Within the non-Hermitian formalism, the effective Hamiltonian eigenstates become resonances, with complex-valued energies E_r − iΓ/2, where E_r is the resonance energy and Γ the resonance width, whose inverse gives the lifetime of the metastable state.⁵ For example, in the context of molecular systems, this usually means localized discrete states, which are embedded into a continuum of autoionizing scattering states. Then, the PESs that drive the dynamics also become complex-valued, and Γ represents the rate of electron autoionization (or autodetachment).

Even though there are implemented methodologies based on quantum wave packet propagation on coupled and complex PESs,^7–10 in practical terms, they are restricted to problems with few degrees of freedom, hence limiting their potential applicability. In another direction, we have presented a computationally feasible formulation that accounts for the multidimensional and complex character of the PESs,¹¹ yet it missed the inclusion of nonadiabatic effects. It is also worth mentioning recent theoretical efforts in search of critical points in complex PESs.^12,13 Counting with quantum chemistry methods adapted with complex absorbing potentials,¹⁴ it is now possible to locate minimum energy crossings between anion and neutral states,¹² as well as exceptional points in complex PESs¹³ (which are the complex analogous of conical intersections of real PESs^15,16).

Despite these relevant developments, there is still not a computationally feasible method for full multidimensional molecular dynamics simulations that account for both non-Hermitian and nonadiabatic effects. Here, we propose a novel theoretical methodology that fills this important gap. It comprises a generalization of the trajectory surface hopping (TSH) method¹ to the case of complex-valued PESs, while employing the popular fewest switches surface hopping (FSSH) flavor¹⁷ for the hopping probabilities. As such, we name it complex surface fewest switches surface hopping, or CS-FSSH for short. Our main goal is to present the method and establish its potential for performing dynamics simulations of electronically metastable molecular states. We have first validated it by comparing to numerically exact quantum dynamics simulations on model analytical PESs. Then, we demonstrate its potential by surveying the dynamics of a resonant state in an actual molecular system. We mention an independent development¹⁸ where the FSSH approach was applied to complex-valued PESs. There, the authors have employed a one-dimensional model to describe cavity losses in the polaritonic dynamics of azobenzene.

Various kinds of molecular states can undergo irreversible decay to the continuum, and thus can be modelled within a non-Hermitian formalism. Transient anions (shape resonance) are formed by direct attachment of free electrons to a neutral species, which then decay either by electron autodetachment (such that the inverse of Γ provides the detachment lifetime) or by molecular dissociation. The latter relaxation channel, dissociative electron attachment (DEA), is ubiquitous in molecules of various types, even though its underlying working mechanisms are far from understood.¹⁹ Transient anions can also be populated by photoexcitation of an electronically bound anion,^20,21 possibly accessing different regions of the PESs and leading to distinct decaying channels. Related processes include dissociative recombination²² and associative detachment, the latter being the inverse process of DEA. Metastable electronic states that lie above other ionization thresholds can also decay by electron detachment, of which we cite multiply charged anions,²³ superexcited,^24,25 core-excited and core-ionized states.²⁶ The latter kind of state can also decay by two-particle processes such as Auger²⁷ and interatomic coulombic decay.^28,29 Furthermore, metastable states can be formed in bimolecular collisions, which include charge transfer in atom-molecule,³⁰ ion-molecule collisions,³¹ and collisional detachment.

From all the above-mentioned processes that can be formulated in terms of a complex-valued Hamiltonian, we have chosen to survey a DEA reaction as the first application of our proposed methodology. DEA is known to play important roles, to a greater or lesser degree, in DNA damage induced by ionizing radiation,^32,33 in the action of radiosensitizing drugs,^33,34 in the chemical evolution of interstellar media³⁵ and planetary atmospheres,^36,37 in plasma technologies,³⁸ among many other examples.^19,39 Of particular relevance herein, low energy electrons can very efficiently promote DEA reactions in halogen-containing unsaturated compounds. This process is believed to account for the radiosensitization properties of modified nucleobases,^34,40,41 and also for the degradation of persistent halogenated compounds.⁴² The general picture is that electron attachment forms a π* shape resonance (associated with the unsaturated region), which undergoes efficient vibronic coupling to σ* resonances (dissociative along the carbon–halogen bond).^43–46

While this π*/σ* mechanism is widely accepted to take place in halogenated species, being often evoked to interpret experimental data,^40,41,47 it has never been theoretically demonstrated, to the best of our knowledge. More importantly, though, the overall understanding about it is rather qualitative, as the underlying mechanisms that drive the reaction are not yet figured out. In particular, there is no information about the sequence of steps of the DEA pathway, its timescales, and the couplings among electronic states. On top of that, it is still not clear how DEA triggered by π* resonances can present cross sections extending over several orders of magnitude.

As a first application of the CS-FSSH methodology, we have investigated the low energy electron induced dynamics of iodoethene. The molecule is expected to display a prototypical π*/σ* indirect DEA mechanism,⁴⁶ which should also be encountered in biologically relevant halogen-containing molecules.^40,41 Our goals with this particular application are twofold. First, to provide theoretical evidence for the existence of the π*/σ* mechanism, and second, to unveil the underlying pathway that culminates in molecular dissociation. Some of our findings might start to shed some light on how DEA reactions take place in larger halogenated molecules.

2 Methodology

The construction of the CS-FSSH methodology follows the same steps of the usual FSSH approach.¹⁷ In short, nuclei are propagated classically on the real-valued PESs computed on-the-fly with electronic structure theory, while the electronic wave function is propagated with the time-dependent Schrödinger equation (TDSE) on the complex PESs. The FSSH prescription mediates the coupling between electronic and nuclear degrees of freedom.

Let us consider a complex-valued electronic Hamiltonian H = H^R + iH^I, with a real component H^R and an imaginary component H^I. The electronic wave function Φ(r,R,t) is governed by the TDSE

where r denotes the collective electron variables, and R those of the nuclei, with corresponding classical trajectory R(t) assumed to be known. Then, the time-dependent electronic wave function Φ(r,R,t) can be expanded in a set of orthonormal basis functions {ψ_j}, which depend explicitly on r and implicitly on R(t),

Introducing the above expansion into the TDSE, multiplying by and integrating in r leads to a set of coupled equations for the expansion coefficients c_j,

where the off-diagonal terms of the Hamiltonian are the diabatic couplings, 〈ψ_j|H|ψ_k〉 = H_jk = H^R_jk − iΓ_jk/2, while the diagonal terms, H_jj = E_j − iΓ_j/2, provide the energies E_j and widths Γ_j. We have introduced σ^NAC_jk ≡ 〈ψ_j|∂ψ_k/∂t〉 = F_jk·v, where F_jk = 〈ψ_j|∇_R|ψ_k〉 are the nonadiabatic coupling vectors and v = dR/dt is the classical nuclei velocity. The deduction of eqn (3) is exactly the same as in usual TSH,¹⁷ with the key difference that here we acknowledge that energies and couplings might be complex-valued. Decoherence effects should still be included, and we do it with an adapted version of the simplified decay of mixing⁴⁸ (equations are shown in the ESI, Sec. S1†).

The first term in the right-hand side of eqn (3) controls the phase of the amplitude c_j, according to the energy E_j. The second term accounts for the decrease of the complex amplitudes, whose rate is given by Γ_j/ℏ. The populations ρ_jj are defined as the diagonal elements of the density matrix , and the total population p as the sum over the contribution of each PES, . In this way, the total population of each classical trajectory will not be conserved, but will lie between 0 and 1 instead. Finally, the third term in eqn (3) is responsible for the population transfer between states and has three contributions. Diabatic H^R_jk and nonadiabatic σ^NAC_jk coupling terms promote the direct interaction between states j and k. Additionally, the states may interact indirectly, via Γ_jk, which represents a coupling through the continuum.

While the electronic wave function is propagated quantum mechanically, nuclei are assumed to be classical particles, which thus follow Newton's equations,

where a single PES j provides the forces, and M_α is the mass of nucleus α.

The decision of which is the current PES j that guides the nuclei dynamics is based on Tully's FSSH scheme.¹⁷ Within the adiabatic representation for the set {ψ_j}, the formula for the hopping probabilities becomes

The above hopping formula differs from the usual FSSH one¹⁷ by the indirect coupling term Γ_jk/2ℏ.

Whether a hopping event takes place or not is decided stochastically, by comparing the hopping probabilities to a random number between 0 and 1, at each time step of the TDSE. When the hopping algorithm indicates a hop to an upper PES, we must first ensure that total energy is conserved after the velocity vector is rescaled along the nonadiabatic coupling vector.⁴⁹ If that is not possible, the hopping is said to be frustrated, and the dynamics resumes on the same PES.

Once the dynamics simulation is finished, any time-dependent observable of interest ŝ(t) can be computed by performing a weighted average over the contribution from each trajectory, as

where the sum runs over a total of N_t trajectories, s_i(t) represents the observable and p_i(t) the total population, for trajectory i and at time t. While in usual FSSH, statistical averages are performed with the same expression, there is a key conceptual difference to the CS-FSSH methodology. In the latter, the total populations are different from 1, such that the contribution from each trajectory is weighted by the factor .

A limitation of the methodology is that the coupling between discrete and continuum states is unidirectional, i.e., only hops from the former to the latter can take place, which is a reasonable approximation for most situations of interest. In the particular case of transient anions, it means that the molecule cannot recapture the detached electron. These so-called nonlocal effects should only be relevant at the lowest collision energies and for very broad resonances,⁵⁰ which do not apply to our present case. The CS-FSSH methodology should also present the same limitations as in usual FSSH, concerning the lack of zero-point energy, interference effects, and quantum tunnelling. At the same time, the quality of the results is intimately related to the employed level of theory for describing the real and imaginary components of the PESs.

3 Validation against quantum dynamics

Before employing the CS-FSSH methodology to a real molecular system, we have first benchmarked it against numerically exact quantum propagation results on model PESs, shown in Fig. 1.

Fig. 1 Results for the 2RHE PESs, for stronger (left) and weaker (right) diabatic couplings. On the top panels, energies (full lines) and widths (dashed lines) for the lower (blue) and upper (red) adiabatic resonant states, and for the state that defines the decaying threshold (black). Adiabatic populations are shown on the other panels, when the dynamics starts from the upper state (middle panels) or lower state (bottom panels). Populations for the lower (p₁, blue) and upper (p₂, red) states were obtained with the quantum dynamics (QD, darker dashed lines) and with the CS-FSSH methodology (lighter full lines).

Our diabatic potentials comprise two coupled resonances, one exponential and one harmonic (2REH model). An additional state is uncoupled and defines the decaying threshold where the resonance becomes a bound state, and was modelled as

The exponential (dissociative) and the harmonic (non-dissociative) states are coupled, allowing nonadiabatic transitions between them in the adiabatic representation. Their real component diabatic potentials are given as

V₂(R) = [V₂(0) − D]e^−αR + D, (9)

and they interact through the diabatic coupling

V₁₂(R) = V₁₂(R_c)e^{−β(R − R_c)2}, (10)

where R_c is the crossing point between V₁ and V₂. The coupling via the continuum was neglected, Γ₁₂ = 0. Resonance widths were modelled as being proportional to the resonance energy,

The vibrational ground state of V₀ was employed as the initial wavepacket and as the sampling distribution for the ensemble of initial conditions.

The parameters of the 2RHE model, as well as the computational details concerning the CS-FSSH and quantum propagation, are presented in Sec. S2 of the ESI.† While different profiles for the PESs have been considered, the most important aspect to be assessed in our validation would be the strength of the coupling. Thus, we decided to compare results for two diabatic coupling strengths for the same set of diabatic PESs.

The adiabatic PESs and the evolution of the populations are shown in Fig. 1, as computed with both theoretical approaches. Along with the conventional population transfer between states, here the total population (sum of lower and upper populations, not shown) is reduced because of the coupling to the continuum. The main finding is that both the CS-FSSH methodology and the reference quantum propagation deliver quite consistent results for both weaker and stronger coupling strengths. Our methodology correctly reproduces the decay of the population and its transfer among the two resonant states.

We have surveyed four cases, for dynamics starting at the lower and upper states, and for weaker and stronger couplings. In all of them, the initial population drops by around two orders of magnitude, due to the positive resonance widths (notice the scale of the y-axis in Fig. 1). Then, population transfer takes place, which is more pronounced in the case of a weaker diabatic coupling, as would be expected. In the meantime, the population keeps decaying before the resonance becomes a bound state (blue-black crossing). Therefore, when dynamics starts in the upper state, the final total population is larger for weaker diabatic couplings. In contrast, a stronger coupling minimises the decay of population when the lower state is initially populated. Finally, the trajectories (or the wavepacket, for the quantum propagation) eventually dissociate along the lower PES, and the total population stabilises.

When starting at the lower state (bottom panels of Fig. 1), the final total population is somewhat superestimated in the CS-FSSH dynamics (with respect to the quantum results) for both coupling strengths. This effect is probably related to the presence of frustrated hoppings: the trajectory should hop to the upper state and lose an extra population, but instead, it remains on the lower one, where population decay stops shortly after. This issue could be minimised with alternative prescriptions when a frustrated hopping is encountered⁵¹ or with methods that impose energy conservation only at the ensemble level.^52,53

4 Application to a realistic case

4.1 Computational details

For the dynamics propagation, both anion and neutral electronic states of iodoethene were described at the multireference configuration interaction (MRCI) level of theory. Two roots have been computed for the anion, which in the neutral equilibrium geometry corresponded to a lower-lying state (σ* for short) and a higher-lying π* state. At other geometries, they mix their diabatic characters, and we refer to the lower and upper states of the anion. Inspection of the orbitals and energies along the trajectories confirmed that our computational protocol for the electronic problem avoided undesired pseudo-continuum states.

The resonance widths were modelled based on the ab initio electron scattering calculations, performed with the Schwinger multichannel method.⁵⁴ Additionally, these calculations have corroborated the MRCI results, as the same valence anion states were obtained in both types of calculations.

The dynamics simulations were performed with a development version of the Newton-X package,^55,56 interfaced with the Columbus package.⁵⁷ Further details about the bound state, scattering, and dynamics calculations can be found in Secs. S3.1–S3.3 of the ESI.†

It is worth emphasizing that the CS-FSSH methodology is not restricted to this particular choice of methods for solving the electronic problem. It could be combined with any methodology that provides the required ingredients for the dynamics (energies, nuclear gradients, nonadiabatic couplings and resonance widths).

4.2 Modeling the complex PES

While the real-valued PESs were computed on-the-fly with the MRCI method, for the imaginary component, we have employed a model that combined scattering and bound state results.

First, some test dynamics simulations on real-valued-surfaces revealed which were the most active coordinates upon electron attachment. Based on that information, ab initio electron scattering calculations were performed for a couple of selected geometries, which provided resonance energies and widths. A cubic interpolation for this set of results provided our model for Γ_j as a function of the corresponding energy E_j. Their profiles and the fitting parameters are shown in Fig. S1 and Table S2 of the ESI.† At each time step of the simulation, resonance energies were obtained from the bound state calculations (MRCI), which fed the model that provided the instantaneous widths.

In this first application, we neglected the nondiagonal terms Γ_kj, which couple the two anionic states through the continuum. This contribution should be much smaller than the direct coupling viaF_kj. Qualitatively, the former kind of coupling means the extra electron undergoes two hops (discrete to continuum, and then to discrete). This should be much less likely than the latter kind, which involves a single hop (discrete to discrete).

Furthermore, we have introduced a systematic correction to the real-valued PESs, based on the energies obtained with the scattering calculations. Additional details about this correction and the construction of the width model are given in Sec. S3.5 of the ESI.†

4.3 Electron-induced dynamics

In the previous subsections, we presented the crucial computational aspects of our calculations and a practical way of computing Γ_k on-the-fly, which is used with the usual real quantities to integrate eqn (3). In the present subsection, we discuss our main findings for the electron-induced dynamics of iodoethene, whose dissociation mechanisms are schematically depicted in Fig. 2.

Fig. 2 Schematic PESs for the neutral ground state (grey), and diabatic σ* (blue) and π* (red) anion states of iodoethene, along dissociating (x) and coupling (y) coordinates. Electron attachment into the σ* state promptly induces cleavage into the iodine ion, whereas a π*/σ* coupling offers a route for dissociation when the π* state is initially formed.

Key observables from the simulations are presented in Fig. 3. Formation of the σ* anion state (left panels) is followed by a significant loss of population due to electron autodetachment, but also by rapid stretching of the C–I bond. This stretching, in turn, provides stability against detachment, until the resonance becomes an electronically bound state (at 8.3 fs in the ensemble average, see purple-black curve crossing). From there, cleavage of the bond resumes and the iodine ion is formed. The process is barrierless, as the ensemble average anion energy (which virtually always coincides with the lower state energy), decreases monotonically. At longer times, the total population stabilises at 0.20.

Fig. 3 Results for dynamics starting at the lower (left) and upper (right) anion states. On the top panel, lower (blue), upper (red) and total (purple) adiabatic populations. On the middle, lower (blue), upper (red) and average (purple) anion energies, as well as average neutral energies (black), where 0 eV corresponds to the neutral at its equilibrium geometry. And on the bottom, ensemble average of the C–I (purple) and C=C (pink) stretching from the neutral equilibrium geometry.

At electron impact energies around 1 eV, the π* resonance is formed (Fig. 3, right panels), and its relaxation is considerably more involved, yet extremely fast. Within 5 fs, the C=C bond is stretched (which works as the tuning coordinate), and the hydrogens move out of the molecular plane (coupling coordinate), while electron autodetachment reduces the total population to around 0.25. These geometrical changes bring the system close to a conical intersection, where several trajectories hop to the lower PES, mostly between 5 and 15 fs. This feature can be seen in the profile of the average anion energy, which initially matches with the upper state energy, and later approaches the lower one. Meanwhile, the anion changes its character from π*, to a mixed π*/σ* and finally to σ*. Once the latter dominates, the dynamics unfolds as described before, culminating into dissociation of the iodine ion.

In addition, strong vibrational excitation of the C=C mode takes place. Between 5 fs and 15 fs, the anion becomes electronically stable for most of the trajectories, such that the average ensemble population decreases much more slowly than in the beginning. The anion takes an average of 8.3 fs to become more stable than the neutral, which coincidentally is the same as when the σ* initiates the dynamics. Trajectories that missed the first opportunity to hop to the lower PES eventually do so, but with increasingly smaller populations as time passes by, since they still remain on regions of the PES with positive resonance widths. The total population finally converges to around 0.17.

An important finding concerns the observation of an energy barrier when the π* resonance triggers the dynamics (see middle right panel of Fig. 3). Since it represents an ensemble quantity, it does not mean that every trajectory encounters a barrier, but rather that it manifests itself in an average way. In the present case, the top of the barrier lies below the initial energy, and therefore it does not hinder the dissociation pathway. This result is consistent with what has been discussed for the similar chloroethene molecule,¹³ which is believed to display an analogous π*/σ* DEA pathway.⁴⁶ In this previous study, however, the discussion was supported by one-dimensional PESs interpolated between critical points, including the minimum energy crossing for the relevant anion states and between anion and neutral. In contrast, our conclusion about the energy profile has been obtained dynamically. Importantly, we have also found that between 5 fs and 10 fs, the ensemble neutral energy rises even faster than the anion energy, such that the resonance energy actually decreases and drops to zero before the top of the barrier is reached. If that were not the case, the presence of the barrier would leave extra time for electron autodetachment, thereby disfavouring decay by DEA.

In Fig. 4, the computed DEA cross sections (see Sec. S3.4 in the ESI†) are compared to the single available set of experimental data.⁴⁶ The latter was not presented in absolute units, and as such the data has been normalised to our most intense peak. The main finding from this comparison is that our calculations are able to reproduce the relative contribution of each state to the DEA cross section. The difference in magnitude is not related to the final populations, which are quite similar (0.20 and 0.17). Instead, it is simply due to the 1/E kinetic factor that appears in the DEA cross sections expression (eqn (S6) in the ESI†).

Fig. 4 Calculated DEA cross section for iodoethene (purple), when starting from the lower (blue) or the upper (red) anion states, compared to available experimental data (grey dots).⁴⁶

Part of the observed discrepancy of around 0.25 eV for the position of the first peak can be traced back to inaccuracies in the scattering calculations (details in the ESI, Sec. S3.5†). Additionally, nonlocal and zero-point energy effects (not encompassed in the methodology) could play a role, given the proximity of this resonance to the detachment threshold.

5 Conclusions and outlook

This contribution provides a novel and general theoretical framework for describing multidimensional nonadiabatic dynamics of molecular systems where irreversible decaying mechanisms take place. The CS-FSSH methodology is a natural generalisation of the FSSH method for complex-valued potential energy surfaces. As such, it assimilates the underlying concepts of an ensemble of classical trajectories and hopping between PESs, but here the total population may decay with time, and the hopping probability is modified by imaginary couplings terms. The comparison with respect to quantum dynamics results for our one-dimensional 2RHE model is outstanding. Both theoretical approaches provided quite similar behaviours for the populations, in both strong and weak couplings regimes. The method has been implemented into a developing version of the Newton-X package, which is going to be freely available to the community in a future release of the code.

As a first realistic application of the CS-FSSH methodology, we have surveyed the dissociation of iodoethene induced by low-energy electron attachment. Our simulations provided a detailed picture of the underlying relaxation mechanisms of the transient anion states. Formation of the σ* resonance promptly induces cleavage into the iodine ion. In contrast, DEA from the π* resonance was found to be intrinsically multidimensional, as it first relaxes by elongation of the C=C bond and out-of-plane vibrations, which promote the coupling to the σ* state and C–I bond breaking. This result is the first theoretical confirmation of the long proposed π*/σ* coupling mechanism in DEA reactions of halogenated molecules. The change of character between anion states takes place mostly between 5 fs and 15 fs, while electronic stability is reached within 10 fs on average. Furthermore, the dissociation pathway encountered a low energy barrier, already past the autodetachment region, which therefore did not quench the DEA channel. However, the presence of such barriers might play a role in explaining the large variations in the efficiency of electron-induced reactions. As a final validation of the methodology and calculations, the shape of the computed DEA cross sections agreed quite well with the available measured data.

The exploitation of the methodology in the context of transient anions is expected to push forward our understanding of the very rich processes induced by low-energy electrons. This, in turn, would provide important insights to various fields, ranging from astrochemistry, to radiobiology and material sciences.

In addition, the methodology could also be employed to simulate the dynamics of molecular systems that present arbitrary decaying mechanisms to a continuum, which would be modelled by imaginary potentials. As stated in the Introduction, this includes Auger and interatomic coulombic decay, and electron detachment from neutral super-excited and core-excited states and transient anion states, the latter being formed by electron impact, photoexcitation or bimolecular collisions. Dynamics simulations of these processes might become feasible by employing the CS-FSSH methodology in conjunction with proper modelling of the decaying widths. This would be accomplished with approaches that combine different types of calculations, as we have done here, with scattering and bound state methods. Alternatively, one could use adapted quantum chemistry methods that directly provide information about the widths.^14,58,59

As one final potential application of the methodology, we put forward the idea that it could be employed to incorporate the possibility of radiative decay into nonadiabatic dynamics of excited states. Although there are available computational codes for mixed quantum-classical dynamics that place internal conversion and intersystem crossing on an equal footing,⁶⁰ we are not aware of implementations that also account for decay by fluorescence. Since its typical time scales (from hundreds of ps to μs) are a lot longer than the ultrafast processes (typically below 1 ps) that are usually probed in nonadiabatic dynamics studies, it is safe to ignore radiative decay for most purposes. Nevertheless, as simulations advance toward longer time scales,⁶¹ decay by fluorescence will have to be acknowledged, and the present methodology enables exactly that. Since radiative decay is an irreversible process, the emission rate would be incorporated into the TDSE as an imaginary potential, suitable for the CS-FSSH method. The same reasoning applies to the case of phosphorescence.

Overall, our proposed methodology should largely extend the range of phenomena that can be surveyed with nonadiabatic dynamics simulations. Finally, we hope the present work is going to foment further theoretical developments on the topic.

Conflicts of interest

There are no conflicts of interest to declare.

Acknowledgements

FK and MB thank the support of the Excellence Initiative of Aix-Marseille University (A*MIDEX) and the project Equip@Meso (ANR-10-EQPX-29-01), both funded by the French Government “Investissements d’Avenir” program. They also acknowledge funding from the WSPLIT project (ANR-17-CE05-0005-01).

References

1. R. Crespo-Otero and M. Barbatti, Chem. Rev., 2018, 118, 7026–7068.
2. B. F. E. Curchod and T. J. Martínez, Chem. Rev., 2018, 118, 3305–3336.
3. T. R. Nelson, A. J. White, J. A. Bjorgaard, A. E. Sifain, Y. Zhang, B. Nebgen, S. Fernandez-Alberti, D. Mozyrsky, A. E. Roitberg and S. Tretiak, Chem. Rev., 2020, 120, 2215–2287.
4. I. Rotter, J. Phys. A, 2009, 42, 153001.
5. N. Moiseyev, Non-Hermitian Quantum Mechanics, Cambridge University Press, 2011.
6. R. El-Ganainy, K. G. Makris, M. Khajavikhan, Z. H. Musslimani, S. Rotter and D. N. Christodoulides, Nat. Phys., 2018, 14, 11–19.
7. D. J. Haxton, C. W. McCurdy and T. N. Rescigno, Phys. Rev. A: At., Mol., Opt. Phys., 2007, 75, 012710.
8. D. J. Haxton, T. N. Rescigno and C. W. McCurdy, Phys. Rev. A: At., Mol., Opt. Phys., 2007, 75, 012711.
9. A. Moradmand, D. S. Slaughter, D. J. Haxton, T. N. Rescigno, C. W. McCurdy, T. Weber, S. Matsika, A. L. Landers, A. Belkacem and M. Fogle, Phys. Rev. A: At., Mol., Opt. Phys., 2013, 88, 032703.
10. S. Felicetti, J. Fregoni, T. Schnappinger, S. Reiter, R. de Vivie-Riedle and J. Feist, Photoprotecting uracil by coupling with lossy nanocavities, 2020.
11. F. Kossoski, M. T. d. N. Varella and M. Barbatti, J. Chem. Phys., 2019, 151, 224104.
12. Z. Benda and T.-C. Jagau, J. Chem. Theory Comput., 2018, 14, 4216–4223.
13. Z. Benda and T.-C. Jagau, J. Phys. Chem. Lett., 2018, 9, 6978–6984.
14. T.-C. Jagau, K. B. Bravaya and A. I. Krylov, Annu. Rev. Phys. Chem., 2017, 68, 525–553.
15. H. Estrada, L. S. Cederbaum and W. Domcke, J. Chem. Phys., 1986, 84, 152.
16. S. Feuerbacher, T. Sommerfeld and L. S. Cederbaum, J. Chem. Phys., 2004, 120, 3201.
17. J. C. Tully, J. Chem. Phys., 1990, 93, 1061.
18. P. Antoniou, F. Suchanek, J. F. Varner and J. Foley, 2020, https://chemrxiv.org/articles/preprint/The_Role_of_Cavity_Losses_in_Polaritonic_Chemistry/12673259.
19. I. I. Fabrikant, S. Eden, N. J. Mason and J. Fedor, Advances in Atomic, Molecular and Optical Physics, Elsevier Inc., 1st edn, 2017, vol. 66, pp. 545–657.
20. C. S. Anstöter, J. N. Bull and J. R. Verlet, Int. Rev. Phys. Chem., 2016, 35, 509–538.
21. C. S. Anstöter, G. Mensa-Bonsu, P. Nag, M. Ranković, R. Kumar T. P., A. N. Boichenko, A. V. Bochenkova, J. Fedor and J. R. R. Verlet, Phys. Rev. Lett., 2020, 124, 203401.
22. A. Florescu-Mitchell and J. Mitchell, Phys. Rep., 2006, 430, 277–374.
23. A. Dreuw and L. S. Cederbaum, Chem. Rev., 2002, 102, 181–200.
24. Y. Hatano, Phys. Rep., 1999, 313, 109–169.
25. W. Li, R. R. Lucchese, A. Doyuran, Z. Wu, H. Loos, G. E. Hall and A. G. Suits, Phys. Rev. Lett., 2004, 92, 083002.
26. P. Norman and A. Dreuw, Chem. Rev., 2018, 118, 7208–7248.
27. P. Auger, J. Phys. Radium, 1925, 6, 205–208.
28. L. S. Cederbaum, J. Zobeley and F. Tarantelli, Phys. Rev. Lett., 1997, 79, 4778–4781.
29. T. Jahnke, J. Phys. B: At., Mol. Opt. Phys., 2015, 48, 082001.
30. D. Almeida, F. Ferreira da Silva, G. García and P. Limão-Vieira, Phys. Rev. Lett., 2013, 110, 023201.
31. M. S. Huq, L. D. Doverspike and R. L. Champion, Phys. Rev. A: At., Mol., Opt. Phys., 1983, 27, 785–794.
32. B. Boudaïffa, P. Cloutier, D. Hunting, M. A. Huels and L. Sanche, Science, 2000, 287, 1658–1660.
33. E. Alizadeh, T. M. Orlando and L. Sanche, Annu. Rev. Phys. Chem., 2015, 66, 379–398.
34. H. Abdoul-Carime, M. A. Huels, E. Illenberger and L. Sanche, J. Am. Chem. Soc., 2001, 123, 5354–5355.
35. C. R. Arumainayagam, R. T. Garrod, M. C. Boyer, A. K. Hay, S. T. Bao, J. S. Campbell, J. Wang, C. M. Nowak, M. R. Arumainayagam and P. J. Hodge, Chem. Soc. Rev., 2019, 48, 2293–2314.
36. X. D. Wang, X. F. Gao, C. J. Xuan and S. X. Tian, Nat. Chem., 2016, 8, 258–263.
37. V. Vuitton, P. Lavvas, R. Yelle, M. Galand, A. Wellbrock, G. Lewis, A. Coates and J.-E. Wahlund, Planet. Space Sci., 2009, 57, 1558–1572.
38. S. Samukawa, M. Hori, S. Rauf, K. Tachibana, P. Bruggeman, G. Kroesen, J. C. Whitehead, A. B. Murphy, A. F. Gutsol, S. Starikovskaia, U. Kortshagen, J.-P. Boeuf, T. J. Sommerer, M. J. Kushner, U. Czarnetzki and N. Mason, J. Phys. D, 2012, 45, 253001.
39. C. R. Arumainayagam, H.-L. Lee, R. B. Nelson, D. R. Haines and R. P. Gunawardane, Surf. Sci. Rep., 2010, 65, 1–44.
40. F. Kossoski, M. H. F. Bettega and M. T. d. N. Varella, J. Chem. Phys., 2014, 140, 024317.
41. F. Kossoski and M. T. D. N. Varella, Phys. Chem. Chem. Phys., 2015, 17, 17271–17278.
42. J. Luo, J. Hu, X. Wei, L. Fu and L. Li, Chemosphere, 2015, 131, 17–33.
43. J. P. Johnson, L. G. Christophorou and J. G. Carter, J. Chem. Phys., 1977, 67, 2196–2215.
44. P. Burrow, A. Modelli, N. Chiu and K. Jordan, Chem. Phys. Lett., 1981, 82, 270–276.
45. R. Kaufel, E. Illenberger and H. Baumgärtel, Chem. Phys. Lett., 1984, 106, 342–346.
46. J. K. Olthoff, J. A. Tossell and J. H. Moore, J. Chem. Phys., 1985, 83, 5627.
47. A. Scheer, K. Aflatooni, G. Gallup and P. Burrow, Phys. Rev. Lett., 2004, 92, 068102.
48. G. Granucci and M. Persico, J. Chem. Phys., 2007, 126, 134114.
49. P. Pechukas, Phys. Rev., 1969, 181, 174–185.
50. W. Domcke, Phys. Rep., 1991, 208, 97–188.
51. A. W. Jasper, C. Zhu, S. Nangia and D. G. Truhlar, Faraday Discuss., 2004, 127, 1–22.
52. A. V. Akimov and O. V. Prezhdo, Phys. Rev. Lett., 2014, 113, 153003.
53. C. C. Martens, J. Phys. Chem. Lett., 2016, 7, 2610–2615.
54. R. F. da Costa, M. T. d. N. Varella, M. H. F. Bettega and M. A. P. Lima, Eur. Phys. J. D, 2015, 69, 159.
55. M. Barbatti, M. Ruckenbauer, F. Plasser, J. Pittner, G. Granucci, M. Persico and H. Lischka, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2014, 4, 26–33.
56. M. Barbatti, G. Granucci, M. Ruckenbauer, F. Plasser, R. Crespo-Otero, J. Pittner, M. Persico and H. Lischka, NEWTON-X: A package for Newtonian dynamics close to the crossing seam. Version 2, 2016, http://www.newtonx.org, accessed July 1, 2020.
57. H. Lischka, R. Shepard, T. Müller, P. G. Szalay, R. M. Pitzer, A. J. A. Aquino, M. M. Araújo do Nascimento, M. Barbatti, L. T. Belcher, J.-P. Blaudeau, I. Borges, S. R. Brozell, E. A. Carter, A. Das, G. Gidofalvi, L. González, W. L. Hase, G. Kedziora, M. Kertesz, F. Kossoski, F. B. C. Machado, S. Matsika, S. A. do Monte, D. Nachtigallová, R. Nieman, M. Oppel, C. A. Parish, F. Plasser, R. F. K. Spada, E. A. Stahlberg, E. Ventura, D. R. Yarkony and Z. Zhang, J. Chem. Phys., 2020, 152, 134110.
58. V. Averbukh and L. S. Cederbaum, J. Chem. Phys., 2005, 123, 204107.
59. N. Sisourat, S. Engin, J. D. Gorfinkiel, S. Kazandjian, P. Kolorenč and T. Miteva, J. Chem. Phys., 2017, 146, 244109.
60. S. Mai, P. Marquetand and L. González, Int. J. Quantum Chem., 2015, 115, 1215–1231.
61. J. Westermayr, M. Gastegger, M. F. S. J. Menger, S. Mai, L. González and P. Marquetand, Chem. Sci., 2019, 10, 8100–8107.

---

Footnote

† Electronic supplementary information (ESI) available: Details about the methodology, computational details about the propagation on the model potential energy surfaces, and computational details about the electronic structure, scattering and dynamics calculations. See DOI: 10.1039/d0sc04197a

---