A trajectory surface hopping study of the vibration-induced autodetachment dynamics of the 1-nitropropane anion

Abstract

In this study, we investigate the autodetachment dynamics of the 1-nitropropane anion after vibrational excitation of the energetically lowest C–H stretching mode using our recently developed extended quantum classical surface hopping approach including the detachment continuum. Therein the detachment from an electronic bound anion state is treated as a nonadiabatic transition into discretized detachment continuum states for an ensemble of classical nuclear trajectories propagated on quantum mechanical potential energy surfaces. The initial ensemble is obtained by sampling a phase space distribution accounting for the vibrational excitation of the C–H stretching mode of the molecule to match the experimental conditions. The simulated kinetic energy distribution of the ejected electrons reproduces characteristic features of the available experimental data. Analysis of the nuclear dynamics points out that the approach to neutral-like geometries with decreased pyramidalization angle of the NO$_2$ group and reduced the N–O bond lengths are the crucial factors enhancing the ultrafast autodetachment process in vibrationally excited 1-nitropropane. This is facilitated when the dipole-bound first excited state of the anion is populated, which is structurally similar to the neutral system. Although only a small transient population of this state is observed, it acts as an efficient doorway to the detachment continuum and is responsible for a significant amount of the ejected electrons.

1 Introduction

Molecular anions have been of great interest for many years to both theoretical and experimental researchers because of their widespread occurrence in a wide variety of contexts such as surface chemistry [1], lower-atmosphere chemical reactions [2, 3], nucleophilic substitution reactions [3, 4], molecular biology [5,6,7] or acid–base chemistry [8], and simultaneously for the arising difficulties in observing and describing these systems accurately [9, 10]. An early spectroscopic method targeting the detection of negative ions for high-resolution data was autodetachment spectroscopy, in which anions were excited into a metastable resonance state high above the detachment threshold, enabling drastically increased signal resolution by measuring the ejected electrons rather than observing the anionic molecules directly [11]. The necessity of this detection method arises from the fact that bound molecular anions, in contrast to neutral molecules, seldom have more than a handful of bound electronic states, if any other than the ground state at all. Therefore, the use of detachment transitions represents a valuable spectroscopic alternative for these systems [12].

Nevertheless, molecular anions with bound excited states are known since the discovery of the excited state in C$_2^-$ [13, 14], and the study of such systems has since developed into an active area of modern research. Notable examples include their involvement in the formation of DNA strand breaks [5, 6], in a competing deactivation process to E$\rightarrow$Z isomerization in biochromophores [7], or the binding of electrons in the field of molecular multipoles [15,16,17,18,19], especially in the case of dipole-bound anion states [16, 20,21,22,23], which are formed when an electron is bound by the strong permanent dipole moment of a neutral molecular core with very small binding energies usually in the range of only a few up to at most some hundred meV [24]. Possible deactivation paths of such systems include the transition to vibrationally excited states of the electronic ground state, subsequently resulting in autodetachment on timescales as short as hundreds of femtoseconds [25,26,27].

Some molecular anions feature detachment energies that are small enough that even vibrational excitation of the electronic ground state is sufficient to eject electrons after a finite amount of time, which is a special case of vibration-induced autodetachment [28, 29]. There are several examples where the excitation of a vibrational mode induces electron loss [30], as, for instance, recently observed in the isomerization process of the vinylidene anion (whose neutral system is a high-energy isomer of acetylene) [31, 32].

The comprehensive theoretical treatment of this phenomenon is challenging due to a number of reasons: The low binding energies in weakly bound anions result in very diffuse electron distributions of the extra electron, requiring the inclusion of sufficiently large and diffuse basis sets in the quantum chemical treatment. Furthermore, the wave function of the ionized system requires the description of an unbound electron scattered at the molecular core, which poses a computationally very demanding problem that can only be tackled be introducing a number of approximations. Finally, the simulation of time-dependent processes requires an accurate treatment of the detachment continuum to define appropriate occupation probabilities for the continuum states. In an effort to gain more insight into the intricate dynamical processes of weakly bound anions upon vibrational excitation, we recently presented a novel surface hopping-based methodology for the quantum classical description of ultrafast vibration-induced autodetachment dynamics in such systems [33] and illustrated it on the example of the vinylidene anion, allowing us to contribute to the understanding of its autodetachment dynamics. We were able to establish a reasonable connection between the dynamic geometrical deformations, ultimately leading to acetylene formation upon excitation of specific vibrational modes of anionic vinylidene, and the efficiency of electron ejection, obtaining good agreement with the experimental data.

Our method has been implemented in a software package called HORTENSIA (Hopping Real-time Trajectories for Electron-ejection by Nonadiabatic Self-Ionization in Anions) that has been recently published on the GitHub platform [34, 35].

Opposite to vinylidene, whose anion does not have electronically excited states below the detachment threshold, the group of nitroalkane anions [36,37,38] exhibits a dipole-bound excited state involved in the autodetachment dynamics after vibrational excitation [39, 40]. Extensive research, experimental as well as theoretical, is available for the two smallest nitroalkane anions: nitromethane [40,41,42,43,44,45,46,47,48,49] and nitroethane [41, 50,51,52,53,54]. For the next larger homologue 1-nitropropane, theoretical data are much scarcer, but there are several experimental studies available, namely by Tsuda et al. [41] on the formation of the anion and by Weber et al., who recorded photoelectron spectra [37] and investigated the characteristics of the autodetachment process induced by infrared excitation of C–H stretching vibrations in comparison with other nitroalkanes [38]. Remarkably, in stark contrast to nitromethane, no vibrational structure was visible in the ejected electron distributions for the higher homologues, which was attributed to fast energy randomization following the population of the metastable vibrational state [38].

In the present contribution, we aim to deepen the molecular-level understanding of this process by investigating the autodetachment dynamics of 1-nitropropane following vibrational excitation of the energetically lowest lying C–H stretching mode, employing our above-mentioned quantum classical surface hopping approach. The remainder of the paper is structured as follows: In Sect. 2, the theoretical methodology is presented. Computational details are provided in Sect. 3, which is followed by the presentation and discussion of results in Sect. 4. Finally, conclusions are given in Sect. 5.

On the present occasion, we take the liberty to point out that our work in the field of quantum classical dynamics has over the years been profoundly inspired by seminal works of Maurizio Persico, such as Refs. [55,56,57,58,59].

2 Theoretical approach

The theoretical methodology employed here for the simulation of the autodetachment dynamics of anionic nitropropane has been presented in detail in Ref. [33] and will here be only sketched briefly. Our approach is based on the trajectory surface hopping method [60], which describes the nuclear motion classically while maintaining a quantum mechanical treatment of the electronic population dynamics. To account for detachment processes, we augment the manifold of bound electronic states by approximate continuum states and include the coupling between both manifolds. Specifically, we expand the general time-dependent molecular wave function, which is given along the classical trajectory ${\textbf{R}}(t)$, into a basis of discrete and continuous states as

$$\begin{aligned} \Psi \big ({\textbf {r}},{\textbf {R}}[t],t\big ) =&\sum _{m'} c_{m'}(t) \Phi _{m'}\big ({\textbf {r}},{\textbf {R}}[t]\big )\ \nonumber \\&+\sum _{m''} \int \! d^3{\textbf {k}}\ {\tilde{c}}_{m''}({\textbf {k}},t) {\tilde{\Phi }}_{m''}({\textbf {k}},{\textbf {r}},{\textbf {R}}[t]). \end{aligned}$$

(1)

The continuum part of the wave function is subsequently discretized according to

$$\begin{aligned} \int \! d^3{\textbf {k}}\ {\tilde{c}}_{m''}({\textbf {k}},t)&{\tilde{\Phi }}_{m''}({\textbf {k}},{\textbf {r}},{\textbf {R}}[t])\nonumber \\&\approx \sum _i c_{m''}({\textbf {k}}_i,t) \Phi _{m''}({\textbf {k}}_i,{\textbf {r}},{\textbf {R}}[t]), \end{aligned}$$

(2)

and the continuum state wave functions ${\tilde{\Phi }}_{m''}$ are obtained approximately as an antisymmetrized product of the electronic wave function of the neutral molecular core (with $N-1$ electrons) and a scattering state of the free electron:

$$\begin{aligned} {\tilde{\Phi }}_{n''}({\textbf {k}}_i) = \mathcal{A} \left( \Phi ^{\text {(n)}}_{n''} \cdot \psi ({\textbf {k}}_i) \right) , \end{aligned}$$

(3)

where $\mathcal A$ denotes an antisymmetrization operator and the free electron states can be approximated by a plane wave. The discrete bound states are obtained from quantum chemical calculations.

The wave function basis constructed in this particular way is not strictly adiabatic: The bound states are only approximations to the adiabatic eigenstates of the anionic system, and this approximation deteriorates with decreasing electron binding energy of the system because of the limited ability of the Gaussian basis sets commonly employed in quantum chemical calculations to mimic the increasingly diffuse electron distribution of a more and more weakly bound system. Moreover, also the continuum states are not adiabatic since they are constructed from plane waves rather than from true scattering states. Therefore, when inserting the discretized expansion (2) into the time-dependent Schrödinger equation, in the arising coupled set of equations for the expansion coefficients $c_m$,

$$\begin{aligned} {\dot{c}}_n(t) = \sum _m \left[ -\frac{i}{\hbar } H_{nm}({\textbf {R}}[t]) - D_{nm} ({\textbf {R}}[t]) \right] c_m(t), \end{aligned}$$

(4)

besides the nonadiabatic couplings $D_{nm}=\langle \Phi _n |d/dt|\Phi _m\rangle =\dot{{\textbf{R}}}\cdot \langle \Phi _n |\nabla _{\textbf{R}}|\Phi _m\rangle$ also nondiagonal matrix elements $H_{nm}$ of the electronic Hamiltonian remain to be present in general. For the sake of compact terminology, we refer to the latter as diabatic couplings, although we stress that they are not the result of an adiabatic-to-diabatic transformation of the electronic basis states. We note, however, that a similar structure of the coupled equations is obtained when constructing a quasi-diabatic basis from adiabatic states, where besides the diabatic couplings also residual nonadiabatic couplings persist, the inclusion of which is necessary to obtain accurate results [61].

Using the above-described basis implies that the electronic states considered as "bound" contain a certain amount of the true continuum states and the continuum states have partial bound character. The extent of this depends on the actual deviation of the states from the strictly adiabatic states and should, at least for the bound states, be the smaller the more strongly the system is bound, i.e., the larger the vertical electron detachment energy. Therefore, we assume that it is justified to take the initial bound state wave function obtained from a quantum chemical calculation as a good approximation to the actual adiabatic state, which subsequently, when the energy gap to the detachment continuum closes, gets more mixed with continuum states and thus less adiabatic in character. This effect will be less pronounced the more diffuse the basis sets employed in the bound state calculation are, and consequently, the importance of the diabatic couplings can be smaller or bigger. To substantiate these general considerations, we include in the Appendix numerical results obtained for an exactly solvable model system, where the exact adiabatic case can be compared with an approximate treatment and conditions can be found under which either the nonadiabatic or the diabatic couplings are dominant. It should be noted that when the approximate states are very similar to the adiabatic ones, the diabatic couplings are found to be small.

Although beyond the scope of the present paper, it might be valuable to further examine the relation between the strictly adiabatic and the approximate treatment also for the molecular case and to investigate whether conditions can be identified under which the computationally expensive diabatic couplings could be safely neglected.

While the nonadiabatic couplings are calculated employing a finite difference approximation for the time derivative following Ref. [62], the diabatic ones can be shown to depend entirely on two-electron integrals involving molecular orbitals of the neutral and anionic species as well as the unbound electron wave function. Specifically, if we assume detachment from an anionic wave function given by the configuration interaction expansion

$$\begin{aligned} \Phi _m^\mathrm {(a)} = \sum _I C_I^{(m)} \Theta _I, \end{aligned}$$

(5)

where $\Theta _I$ represents an individual Slater determinant formed from the anion’s molecular orbitals (MOs), to an electron-detached continuum state in which the neutral core is represented by the ground state Slater determinant $\Phi _0^\mathrm {(n)}$, the diabatic couplings can be written as

$$\begin{aligned} H_{0m}({\textbf{k}}_i)&=\Delta {\mathcal {V}}_k\sum _I C_I^{(m)} \left\langle {\tilde{\Phi }}_0({\textbf{k}}_i)|{\hat{H}}|\Theta _I\right\rangle \end{aligned}$$

(6)

$$\begin{aligned}&\equiv \Delta {\mathcal {V}}_k\sum _I C_I^{(m)} V^\textrm{dia}_{iI}({\textbf{k}}_i) \end{aligned}$$

(7)

In this expression, $\Delta {\mathcal {V}}_k$ is a discretized volume element in k-space, and $V^\textrm{dia}_{iI}({\textbf{k}}_i)$ denotes the part of the diabatic coupling due to the interaction between the Slater determinants of the electron-detached state, composed of the free electron and the neutral ground state, and the Ith Slater determinant contributing to the anion state m.

These terms are computed approximately, employing orthogonalized plane waves as the one-electron scattering states and approximating the molecular integrals involving plane waves based on a Taylor series expansion, as outlined in detail in Ref. [33], which results in

$$\begin{aligned} V_{iI}^{\textrm{dia}}({\textbf {k}}_i) =&\sum _n^{\textrm{occ}} \sum _u^{\textrm{virt}} \langle \chi _n|\phi _u\rangle \langle {\textbf {k}}_i \phi _u|{\hat{v}}| \nonumber \\&\times \sum _{q,p<q}^{\textrm{occ}} \left( \phi _p\phi _q-\phi _q\phi _p \right) \rangle \nonumber \\&\times (-1)^{n+p+q-1} \textrm{det}\,{\textbf{S}}_{0n,pq}, \end{aligned}$$

(8)

where $\hat{\nu}$ is the operator for electron-electron interaction, $|\phi _i\rangle$ and $|\chi _i\rangle$ denote the MOs of the anion and neutral molecule, respectively, while $|{\textbf {k}}_i\rangle$ is the orthogonalized plane wave describing the free electron and $\textrm{det}\,{\textbf{S}}_{0n,pq}$ denotes the minor determinant of the overlap matrix between bound and continuum state orbitals obtained by omitting the rows 0 and n (corresponding to the plane wave ${\textbf{k}}_i$ and the neutral MO $\chi _n$) as well as the columns p and q (corresponding to the anionic MOs $\phi _p$ and $\phi _q$). Expansion of the above equation into the atomic basis functions (denoted by Greek indices) subsequently leads to

$$\begin{aligned} V^{\textrm{dia}}_{iI}({\textbf {k}}_i)&= \sum _{\lambda \mu \nu } \Bigg [ A_{\lambda \mu \nu } \Big ( \langle {\textbf{k}}_i \lambda || \mu \nu \rangle -\sum _{\sigma } B_{\sigma } \langle \sigma \lambda || \mu \nu \rangle \Big ) \nonumber \\&\quad + {\bar{A}}_{\lambda \mu \nu } \Big ( \langle {\textbf{k}}_i \lambda | \mu \nu \rangle -\sum _{\sigma } B_{\sigma } \langle \sigma \lambda | \mu \nu \rangle \Big ) \Bigg ], \end{aligned}$$

(9)

where $\langle ab|cd\rangle$ denotes an electron–electron repulsion integral and $\langle ab||cd\rangle$ its antisymmetrized variant. The prefactors $A_{\lambda \mu \nu }$, ${\bar{A}}_{\lambda \mu \nu }$ and $B_{\sigma }$ read (cf. Ref. [33]):

$$\begin{aligned} A_{\lambda \mu \nu }&= \sum _n^{\textrm{occ},\alpha } \sum _{q,p<q}^{\textrm{occ},\alpha } (-1)^{n+p+q-1} \textrm{det}\ {\textbf{S}}_{0n,pq} \nonumber \\&\quad \times \left( c_\lambda ^{(n)} - \sum _u^{\textrm{occ},\alpha } c_\lambda ^{(u)} S_{nu} \right) c_\mu ^{(p)} c_\nu ^{(q)} \end{aligned}$$

(10)

$$\begin{aligned} {\bar{A}}_{\lambda \mu \nu }&= \sum _{{\bar{n}}}^{\textrm{occ},\beta } \sum _{p}^{\textrm{occ},\alpha } \sum _{{\bar{q}}}^{\textrm{occ},\beta } (-1)^{{\bar{n}}+p+{\bar{q}}-1} \textrm{det}\ {\textbf{S}}_{0{\bar{n}},p{\bar{q}}} \nonumber \\&\quad \times \left( c_\lambda ^{({\bar{n}})} - \sum _{{\bar{u}}}^{\textrm{occ},\beta } c_\lambda ^{({\bar{u}})} S_{{\bar{n}} {\bar{u}}} \right) c_\mu ^{(p)} c_\nu ^{({\bar{q}})} \end{aligned}$$

(11)

$$\begin{aligned} B_{\sigma }&= \sum _r^{\textrm{occ},\alpha } \sum _\rho c_\sigma ^{(r)} c_\rho ^{(r)} \langle {\textbf {k}}_i | \rho \rangle \end{aligned}$$

(12)

where $c_\lambda ^{(n)}$ denotes the expansion coefficient of atomic orbital $\lambda$ in MO n.

During the course of the simulation, vibration-induced autodetachment is described by nonadiabatic transitions (surface hops) into the continuum states in a stochastic manner after obtaining hopping probabilities from the changes of the electronic state coefficients as described in more detail in Refs. [33, 63, 64].

In addition, during the dynamics situations may arise in which the potential energy of the anionic state increases above that of the neutral state, i.e., the electronic system of the anion becomes unstable with respect to immediate electron loss. This would lead to an evolving free electron wave packet spreading rapidly in space. We model this adiabatic detachment process by including a gradual population loss determined from the 1s-type of a spherical electron distribution having the same electronic spatial extent $\langle {\textbf{r}}^2\rangle$ as the MO from which the electron is ejected. For the free-time propagation of an MO composed of Cartesian Gaussian basis functions of s, p and d type, the following analytic expression is obtained:

$$\begin{aligned} \varphi _\mu ({\textbf{r}},t) =&\,N_{l_xl_yl_z} \textrm{e}^{-\frac{\alpha }{1+i\beta t}\textrm{r}^2} \left[ -\Lambda \frac{i\beta t}{2\alpha } (1+i\beta t)^{-\frac{5}{2}} \right. \nonumber \\&\left. +(1+i\beta t)^{-\frac{3}{2} - \sum _j l_j} \prod _{j=x,y,z} (r_j - A_j)^{l_j}\right] . \end{aligned}$$

(13)

Here $l_i$ denotes the angular momentum quantum number for the ith spatial direction, ${\textbf {A}}$ is the center of the basis function and $\beta =\frac{2\hbar \alpha }{m_e}$. The constant $\Lambda$ is unity if one of the $l_i = 2$ and zero if all $l_i<2$. The population loss is then monitored within a sphere initially containing a population of 99% for the 1s-type electron distribution dispersing in time according to the evolution of the MO’s electronic spatial extent. This rather coarse procedure accounts for the finite lifetime of the electronically metastable anion state, leading to irreversible population decay and taking also into consideration that the electron ejection rate varies as a function of time within periods during which the anion state is electronically unbound.

More sophisticated approaches based on ab initio calculation of the decay lifetime have been presented recently, but are at the moment only available for the electronic ground state in the frame of the Hartree–Fock method [65].

3 Computational

Table 1 Comparison of adiabatic electron affinities (AEA), vertical detachment energies (VDE), anion vertical excitation energies ($\Delta$E) and anharmonic (harmonic) wave numbers of the energetically lowest lying CH stretching mode ($\nu _{27}$) of 1-nitropropane in its most stable conformer

The quantum chemical calculations on which the dynamics simulations are based are carried out using time-dependent density functional theory (TD-DFT) with the Gaussian09 program package [67] at the $\omega$B97XD level of theory, employing the 6-31+G** basis set augmented with 2 additional s- and p-functions each, located at the N atom. The exponents of these additional basis functions are generated from the most diffuse s- and p-functions within the underlying basis set with the geometric progression for a factor of 3.5 according to Ref. [66]. Table 1 provides a comparison of several quantum chemical methods for the electronic structure of the molecule carried out using the Gaussian16 program package [68]. Specifically the TD-DFT method employing the $\omega$B97XD [69], CAM-B3LYP [70] and LC-$\omega$PBE [71] functionals, as well as the coupled cluster approach using single and double excitations (CCSD) [72, 73] and additional perturbational treatment of triple excitations (CCSD(T)) [74], using the 6-31+G** [75], 6-311++G** [76] (with further addition of basis functions) as well as the aug-cc-pVDZ [77, 78], d-aug-cc-pVDZ [79] and d-aug-cc-pVTZ [79] basis sets are shown. It is evident that the CAM-B3LYP and LC-$\omega$PBE functionals either provide systematically too high electron affinities or are unable to produce a bound excited state and are therefore not suitable for the simulation. For the $\omega$B97XD functional, one can see that the system is described a lot more consistently with experimental data and the chosen combination of a small basis set augmented with diffuse basis functions is a good compromise between the accurate reproduction of experimental observables and the resource efficiency needed to carry out computationally expensive dynamics simulations. Furthermore, this type of basis set augmentation is a popular approach that has yielded good results in the description of weakly bound anions before [16, 20].

Fig. 1

Average absolute value of the dipole moment of neutral nitropropane along trajectory ensembles propagated in the anionic ground state. For normal mode $\nu _{27}$, excitation by a single quantum is taken into account by employing either the harmonic Wigner function for the first excited state (blue curve) or the product of squared wave functions in position and momentum space (orange curve, cf. Equation 14). A total number of 180 initial conditions of the lowest energy conformer (cf. Figure 2) is propagated at the $\omega$B97XD/6-31+G** + 2s2p level of theory, and a standard deviation is computed by calculating the average of four disjoint sub-ensembles of 45 trajectories for each set of initial conditions (areas shaded in light color)

The initial conditions for our dynamics simulations were obtained by sampling a phase space distribution function in terms of harmonic normal modes according to

$$\begin{aligned} P^{(i)}_\upsilon (Q_i,P_i) = |\chi ^{(i)}_\upsilon (Q_i)|^2 |{\tilde{\chi }}^{(i)}_\upsilon (P_i)|^2 \end{aligned}$$

(14)

with the harmonic oscillator wave functions of normal coordinate $\nu _i$ in position space, $\chi ^{(i)}_\upsilon (Q_i)$, and momentum space, ${\tilde{\chi }}^{(i)}_\upsilon (P_i)$. For $\upsilon =0$, this corresponds to the common Wigner function. In the example studied here, $\upsilon =1$ for the experimentally excited lowest energy C–H stretching vibration (mode 27 when enumerating all modes by increasing energy), while $\upsilon =0$ for all remaining normal modes. Notice that our choice for distribution function for $\upsilon =1$ differs from the Wigner distribution. However, when integrated over positions or momenta, it gives rise to the same marginal distributions and, at the same time, allows for higher computational efficiency since a single ensemble can be used to sample the distribution. For the excited state Wigner function, by contrast, due to the existence of positive and negative parts, two ensembles would need to be propagated in parallel, and for the calculation of observables, the difference of the averages obtained for these two ensembles has to be taken, which means that for a similar quality of the statistics, a higher number of trajectories is necessary. This is illustrated in Fig. 1 for the ensemble averages of the dipole moment of the neutral core along sets of anion trajectories propagated over 100 fs with vibrational excitation either incorporated by using the Wigner function for $v=1$ or the product of squared wave functions in position and momentum space (Eq. (14)). For both simulations, the averages are very similar (blue and orange curve). However, when calculating the standard deviation of the averages obtained for several smaller sub-ensembles, much larger values result for the Wigner ensemble (blue-shaded area) than for the ensemble generated by Eq. (14) (orange-shaded area).

After the generation of initial conditions, the nuclear trajectories were subsequently obtained by integration of Newton’s equations of motion using the velocity Verlet algorithm [80] with an integration time step of 0.2 fs over a total of 1 ps.

Concerning the electronic degrees of freedom, we included the bound anionic ground and first excited state as well as the detachment continuum discretized as outlined in Sect. 2 and initiated all trajectories in the anionic ground state. For the detachment continuum, we employed 96,000 discretized states, with the k-vectors chosen to account for 1000 evenly spaced kinetic energies between 0.0 and 1.5 eV and 96 different spatial orientations for each energy, generated using the Fibonacci sphere algorithm [81] to evenly cover a spherical surface in k-space. To include the impact of the molecular core on the wave function of the extra electron at least approximately, in each nuclear dynamics time step the plane waves were orthogonalized with respect to all occupied anionic molecular orbitals (MOs) using Gram–Schmidt orthogonalization, as has been done successfully in the past [82,83,84]. “Occupied” refers here to all MOs included in electronic configurations with significant contribution to the overall wave function of an electronic state. For the first excited state, which within TD-DFT is described by a CIS-like expansion into excited Kohn–Sham determinants, we considered all determinants, beginning from the highest contribution, which are needed to represent the full CIS wave function by 95 % (or until a maximum of 10 configurations was reached). The resulting orthogonalized plane waves were then renormalized.

Within the manifold of the electronic ground and first excited bound anion state as well as 96,000 discretized continuum states, the electronic Schrödinger equation 4 was solved along the classical trajectories utilizing Adams’s method as implemented in the integrate.ode class of the scipy Python module [85], including the nonadiabatic and diabatic couplings between the considered states. The time-dependent electronic state coefficients obtained this way were subsequently used to evaluate the hopping probabilities at each nuclear time step. In order to prevent an unphysical loss of zero-point energy due to electron ejection, all hops were suppressed in which the free electron’s kinetic energy $E_{\textrm{kin}}$ was larger than the vibrational excess energy $\hbar \omega _{27}+E^{\textrm{a}}_{\textrm{ZP}}-E^{\textrm{n}}_{\textrm{ZP}}$, where $\hbar \omega _{27}$ is the energy of one vibrational quantum of mode 27 and $E^{{\mathrm {a/n}}}_{\textrm{ZP}}$ are the zero-point energies of anion (a) and neutral (n).

For the decay of electronically unstable anion states, a gradual population loss modeled from free wave packet dispersion as described in Sect. 2 (cf. Equation 13) was taken into account.

Since our interest lies in the dynamics of the anion, the trajectories were not propagated further in the neutral state in the event of a hop into the continuum. This allows us to improve the hopping statistics by assigning to each nuclear trajectory a number of individual “hoppers,” for which the stochastic hopping process is performed individually in each time step. For each instance of a successful hop, the actual number of hoppers, $n_\text {sub}(t)$, is reduced by one. Starting with $n_\text {sub}(t_0) = 1000$, hopping is then checked $n_\text {sub}(t)$ times at a given time step, and the hopping procedure is terminated when the number of hoppers has reached zero.

4 Results

In the experiment, electron detachment of the 1-nitropropane anion can be achieved by vibrational excitation of specific normal modes, in particular, the infrared-intensive lowest lying C–H stretching mode ($\nu _{27}$), which involves the C–H bond closest to the NO$_2$ group [36, 38]. To account for this experimental situation, we modeled the vibrational excitation by sampling initial coordinates and velocities from a phase space distribution as outlined in the Computational section. In addition, it has to be taken into account that 1-nitropropane exhibits several energetically close-lying conformers that can be obtained by rotation around the C–N and the first C–C axis. We considered the energetically lowest three of them, which are shown in Fig. 2 with their respective state populations according to a Boltzmann distribution at a temperature of 250 K (taken from Ref. [38]) with a partition function only taking into account these three conformers. In the energetically lowest conformer, the dihedral angle $\angle$CCCN amounts to 60$^\circ$. The second conformer is linear in the sense that $\angle$CCCN is 180$^\circ$ and the third conformer differs from the first by a NO$_2$ group rotated by 60 degrees around the CN axis. The dynamics simulation was performed with a total of 100 trajectories, with the initial conditions sampled from these three different conformers, amounting to 78 trajectories for the first one, 14 for the second and 8 for the third. Using Eq. 14 and setting $\upsilon = 1$ for the C-H stretching mode closest to the NO$_2$ group ($\nu _{27}$ in the case of the most stable conformer at the $\omega$B97XD/6-31+G** + 2s2p level of theory) and $\upsilon = 0$ for all other modes then results in initial conditions corresponding to the experimental realization by Weber et al. [38].

Fig. 2

Optimized geometries of the 1-nitropropane anion (at the $\omega$B97XD/6-31+G** + 2s2p level of theory) in its most stable conformers. The relative energies with respect to the global minimum and the Boltzmann populations at 250 K are given

Fig. 3

a Experimental electron kinetic energy distribution of the 1-nitropropane anion after excitation of the energetically lowest C-H stretching mode at 2922 cm$^{-1}$ (reproduced from Ref. [38]), b simulated electron kinetic energy distribution of all hopping events after excitation of the lowest C–H stretching mode and propagation for 1 ps (orange histogram), c time-dependent population of all bound anion states (dark green) and exponential fit with a time constant of $\tau = 1105$ fs (light green)

In the dynamics simulation, autodetachment events are observed and analyzed in terms of the kinetic energy of the ejected electrons, which for a nonadiabatic hop is obtained as $E_{\textrm{kin}}=\frac{\hbar ^2{\textbf{k}}_i^2}{2m_e}$ from the k-vectors of the continuum state occupied in the hop, or, in the adiabatic case of an electronically unbound anionic state, from the energy difference between the anionic and the neutral state. The histogram obtained for all hopping events observed during the simulation represents the kinetic energy distribution and is shown in Fig. 3b. The signal exhibits its highest values between 0.01 and 0.02 eV, both decaying toward 0 eV as well as toward higher energies. Above 0.25 eV, most hopping events are suppressed by the condition of zero-point energy conservation. These theoretical findings share several important characteristics with the experimentally observed spectrum of Weber $et\ al.$ [38] which is reproduced in Fig. 3a, such as the position of the signal maximum and the decreasing signal strength for higher and lower energies. However, the simulated signal is broader than the experimental one, approaching zero only at higher energies.

Furthermore, the ultrafast autodetachment timescale can be estimated from the dynamics simulation. From the population of bound anionic states, a fast exponential decay with a time constant of 1105 fs is revealed, as shown in Fig. 3c, with approximately 40% of the total population remaining in a bound anion state after the simulation time of 1000 fs.

Fig. 4

a Pyramidalization angle of the NO$_2$ group at all hopping events into a free electron state for the whole ensemble of trajectories (green histogram) and for the whole ensemble at all time steps (blue curve), b and c larger and smaller of the two N–O atomic distances at all hopping events (orange histograms) and at all time steps (blue curves) of the whole ensemble, respectively

In addition to revealing the timescale, our simulations offer the possibility to analyze the molecular geometries for characteristic changes connected with detachment transitions. Two distinct parameters can be identified in which the structure differs from the average at all instants of time, as shown in Fig. 4. In a), the pyramidalization of the NO$_2$ group, which is quantified as the angle between the O–N–O plane and the C–N axis, at the hopping events and at all time steps is shown in green and blue, respectively. It can be seen that throughout the dynamics, the pyramidalization angle is centered around 30$^\circ$, which is also found as the equilibrium angle in the anionic ground state. At the instances where surface hops occur this angle is significantly reduced, indicating a planarization of the NO$_2$ group similar to the case of neutral nitropropane. This is in accordance with the expectation that for very weakly bound anions, the extra electron is so far away from the molecular core that the molecule assumes a geometry more resemblant of the neutral equilibrium geometry. Furthermore, the N–O atomic distances are shifted to smaller values, as shown in Fig. 4b and c for both the larger and smaller of the two bonds. In blue, the distribution for all time steps is provided, exhibiting an average value of 1.33 Å and 1.28 Å. At the hopping events, these values are shortened to $\sim$1.26 Å and $\sim$1.22 Å for the longer and shorter of the N–O distances, respectively.

Fig. 5

C–H distance for the initially excited mode $\nu _{27}$ a at the hopping events, b for all time steps of the dynamics (blue) as well as the initial distribution at $t=0$ (green), c kinetic energy as a function of time for mode $\nu _{27}$ averaged over the ensemble of trajectories

Interestingly, the initially excited C–H stretching mode $\nu _{27}$ is not directly responsible for enhancing the detachment efficiency, as can be seen from the histograms in Fig. 5a and b showing the distribution of the associated C–H distance both at the hopping events and at all time steps. The distributions are rather similar, which indicates that the C–H bond distance is no relevant factor in the detachment process. Instead, the main role of mode $\nu _{27}$ is to absorb the energy provided by the infrared irradiation, which is facilitated by its quite large oscillator strength (152 $\frac{\textrm{km}}{\textrm{mol}}$ at the DFT/$\omega$B97XD/6-31+G** + 2s2p level of theory). The absorbed energy is rapidly redistributed as can be inferred from the quickly decreasing kinetic energy of the mode shown in Fig. 5c, dropping below 50 % of the initial value after less than 100 fs. No single modes were found that were excited particularly strongly in the energy redistribution process; instead, the excess energy deposited in the C-H stretching mode is quickly dispersed over all vibrational degrees of freedom in accordance with the interpretation of Weber et al. [38] for their experimental findings. This process eventually allows the molecule to assume also geometries with a planar NO$_2$ moiety. Consequently, since in this way the structure of the neutral species is approached, the electron binding energy of the anion gets strongly reduced, leading to an increased autodetachment rate.

In our simulation, we also included the anion’s first excited electronic state, which is of dipole-bound character. Although at a given instant of time, only a small population of this state is reached, which never exceeds $\sim$3% and, on average, amounts to about 0.4% throughout the whole simulation time, it is striking that it is the source of about 56% of the detachment events. This is because the equilibrium structure of the dipole-bound state closely resembles that of the neutral molecule, and the two lie very close in energy, thus increasing the tendency of the system to ionize. In addition, the dipole-bound state gets populated predominantly when the nonadiabatic couplings to the anionic ground state are large, which also is the case in regions of low VDE. Therefore, the dipole-bound state can be regarded as a mediator strongly enhancing autodetachment.

5 Conclusion

We have simulated the autodetachment dynamics of the 1-nitropropane anion following vibrational excitation of a CH stretching mode using our recently developed surface hopping approach augmented by the inclusion of discretized continuum states and their coupling to the bound states. Our simulations allow us to infer some general features of the autodetachment dynamics and unravel its mechanism. The additional vibrational energy upon excitation is quickly redistributed throughout the molecule, enabling access to regions on the potential energy surface where the vertical detachment energy (VDE) of the anion is sufficiently small that autodetachment takes place. The geometric structures at which transitions into the detachment continuum mostly occur are characterized mainly by reduced N–O bond lengths and planarization of the NO$_2$ group, thus bearing similarity to the equilibrium structure of the neutral nitropropane. This is strongly mediated by the dipole-bound first excited state, which is transiently populated in these regions, followed by rapid population decay of the bound anion. The timescale on which these processes take place was found to be $\sim$1 ps. Overall, we were able to gain detailed insight into the intricate dynamics leading to autodetachment and could therefore not only reproduce characteristic features of the available experimental data, but also extend the understanding of the ultrafast deactivation processes after vibrational excitation of the 1-nitropropane anion.

Appendix

In this section, we present calculations of bound-to-continuum couplings and detachment probabilities on an exactly solvable one-dimensional model system, namely the attractive one-dimensional delta potential. The Hamiltonian is given by

$$\begin{aligned} {\hat{H}}=-\frac{\hbar ^2}{2m}\frac{d^2}{dx^2} + v_0 \delta (x). \end{aligned}$$

(15)

For $v_0<0$, this system bears a single bound state,

$$\begin{aligned} \psi _b(x)=\sqrt{\kappa }\,\textrm{e}^{-\kappa |x|}, \end{aligned}$$

(16)

of energy $E_b=-\frac{\hbar ^2\kappa ^2}{2m}$ $\left( \kappa =-mv_0/\hbar ^2\right)$, and has scattering states of the form

$$\begin{aligned} \psi _{sc}(x)=\frac{{\mathcal {N}}(k,\kappa ,L)}{\sqrt{\kappa ^2+k^2}}\left( k\cos kx - \kappa \sin k|x|\right) . \end{aligned}$$

(17)

When the continuum is discretized within a spatial box of length L, the normalization constant ${\mathcal {N}}$ is determined such that $\int _{-L}^L |\psi _{sc}(x)|^2dx=1$.

For the approximative solution, we construct the bound state as a linear combination of Gaussian functions,

$$\begin{aligned} {\tilde{\psi }}_b(x)=\sum _i C_i \,b_i(x), \end{aligned}$$

(18)

with

$$\begin{aligned} b_i(x)= \left( \frac{2\alpha _i}{\pi }\right) ^{1/4} \textrm{e}^{-\alpha _i x^2}, \end{aligned}$$

(19)

where the $C_i$ and the approximate energy ${\tilde{E}}_b$ are obtained variationally by constructing the Hamiltonian matrix in the orthogonalized basis of the Gaussian functions and determining eigenvalues and eigenvectors. The approximate scattering states are constructed as

$$\begin{aligned} {\tilde{\psi }}_{sc}=\frac{1}{\sqrt{1-S^2}} \left( \tilde{{\mathcal {N}}}(k,L)\cos kx - S {\tilde{\psi }}_b(x)\right) , \end{aligned}$$

(20)

with ${\tilde{\mathcal{N}}}$ as the box normalization constant for the cosine function and S as the overlap integral between the normalized cosine and ${\tilde{\psi }}_b$, resulting in the scattering states being orthogonal to the bound state. The exact bound state energy is shown in Fig. 6a) together with two approximate energies obtained employing (i) a basis set of 8 Gaussian functions with exponents of $\alpha _i=$0.005, 0.01, 0.1, 0.5, 1.0, 5.0, 10.0 and 50.0 (in the following termed "small basis," all values in $a_0^{-2}$) and (ii) another one with the same functions plus two additional Gaussians with exponents of 0.0002 and 0.001 $a_0^{-2}$ (henceforth termed “large basis”). Evidently, for the large basis (solid red line) the energy closely follows the exact value (blue dots), while for the small basis (dashed red line), a markedly larger value is obtained which even crosses the line of zero energy. This behavior is due to the inability of the small basis set to properly reproduce the exact wave function for small potential strengths ($v_0$ values approaching zero), since the latter becomes more spatially diffuse than the fixed-width Gaussian functions contained in the basis set allow. This fact is further quantified by the overlap integral ${\tilde{S}}$ between the exact and the approximate wave functions depicted in Fig. 6b), which shows that for the large basis set, the overlap is close to unity for almost all values of $v_0$ except the very smallest ones, while for the small basis set, notable deviations already occur for intermediate potential strengths of about 0.04 $E_Ha_0$.

Fig. 6

a Bound state energies for the attractive delta potential as a function of the potential strength $v_0$. Blue dots: exact energy, red lines: approximate energies (solid: large basis, dashed: small basis, for the basis definition see text). b Overlap integral between the exact and the approximate wave functions (solid line: large basis, dashed line: small basis)

Both for the exact and the approximate wave functions, we calculate a nonadiabatic coupling with respect to changes in the parameter $v_0$:

Fig. 7

Comparison of exact and approximate bound-to-continuum couplings and transition probabilities for the attractive delta potential. Left column: data for a continuum state of $E_{kin}=1$ meV, right column: data for a continuum state of $E_{kin}=100$ meV. a)/b): Diabatic couplings $H_{dia}$ for large (solid line) and small (dashed line) basis set. c)/d): Nonadiabatic couplings ${\tilde{d}}$ for large (solid orange) and small (dashed orange) basis set as well as exact value d (blue). e)/f): Coupling strengths calculated from (blue) the exact nonadiabatic coupling ($V=\hbar {\dot{v}}_0d$), (orange) approximate nonadiabatic coupling ($V=\hbar {\dot{v}}_0{\tilde{d}}$), (green) diabatic coupling $V=H_{dia}$, (red dots) total approximate coupling ($V=\textrm{i}\hbar {\dot{v}}_0{\tilde{d}}+H_{dia}$). The approximate values are computed for the large basis set. g)/h) Same as e)/f) for the small basis set

$$\begin{aligned} d&=\left\langle \psi _{sc}|\frac{d}{dv_0}|\psi _b\right\rangle \end{aligned}$$

(21)

$$\begin{aligned}&=\frac{2m}{\hbar ^2}\left\langle \psi _{sc}|\,|x|\,|\psi _b\right\rangle \end{aligned}$$

(22)

$$\begin{aligned} {\tilde{d}}&=\left\langle {\tilde{\psi }}_{sc}|\frac{d}{dv_0}|{\tilde{\psi }}_b\right\rangle \end{aligned}$$

(23)

$$\begin{aligned}&=\sum _i \frac{dC_i}{dv_0}\left( \tilde{{\mathcal {N}}}\langle \cos kx | b_i\rangle - S \sum _j C_j\langle b_j | b_i\rangle \right) . \end{aligned}$$

(24)

In addition, for the approximate solutions also nondiagonal matrix elements of the Hamiltonian (diabatic couplings) are considered:

$$\begin{aligned} H_{dia}&=\left\langle {\tilde{\psi }}_{sc}|{\hat{H}}|{\tilde{\psi }}_b\right\rangle \end{aligned}$$

(25)

$$\begin{aligned}&=\left( \frac{\hbar ^2k^2}{2m} -{\tilde{E}}_b\right) S+v_0\tilde{{\mathcal {N}}}{\tilde{\psi }}_b(0) \end{aligned}$$

(26)

In Fig. 7, panels a–d, the values of these couplings are presented for two different continuum states: one bearing an energy of 0.0085 $E_H$ (1 meV, left column) and the other one corresponding to 0.085 $E_H$ (100 meV, right column). In addition, in panels e–h, the transition probability following from these couplings is exemplified under the assumption that only the bound and a single discretized continuum state of the aforementioned energy are present. For a general coupling V, this probability is calculated as

$$\begin{aligned} P = \frac{|V|^2}{|V|^2+(E_b-E_{sc})^2}. \end{aligned}$$

(27)

In the calculation, a discretization parameter of $\Delta k=10^{-5} a_0^{-1}$ is employed, which determines the normalization constants of the discretized continuum states. The effective nonadiabatic coupling is calculated as $V_{na}=\hbar \,{\dot{v}}_0 \,d$, and for the "velocity" ${\dot{v}}_0$, a value of $6\cdot 10^{-4} \frac{E_H^2a_0}{\hbar }$ is assumed. This corresponds roughly to the maximal change of $v_0$ when the bound state energy oscillates with a frequency corresponding to a typical molecular vibration of 1000 cm$^{-1}$ around an energy of $-0.5$ eV with an amplitude of 0.49 eV, i.e., within an energy range of roughly 1 eV.

The following main observations can be made from inspecting Fig. 7: The diabatic couplings shown in (a) and (b) are larger for the approximate wave functions obtained with the small basis set (dashed green lines) than for those employing the large basis set (solid green lines), since the former deviate more from the exact solution. Vice versa, the nonadiabatic couplings presented in (c) and (d) are smaller for the small basis set (dashed orange lines), while for the large basis set (solid orange lines) they are comparable to the exact couplings (blue lines). Consequently, the approximate transition probabilities including both nonadiabatic and diabatic couplings are for the large basis set (e) and (f), red dots) dominated by the nonadiabatic contribution (orange) while the diabatic contribution (green) is almost negligible. At the same time, they are close to the exact probabilities (blue) due to the bound states being good approximations to the exact adiabatic states. Still persisting deviations can be attributed to the fact that approximate continuum state functions are used. Notice that in order to increase visual comparability of the probabilities changing over several orders of magnitude with varying $v_0$, their decadic logarithm is actually shown in Fig. 7.

For the transition probabilities obtained with the small basis set (Fig. 7g/h), the situation changes. Now the approximate wave functions are less adiabatic in character for small potential strengths. Therefore, notable deviations from the nonadiabatic couplings obtained in the exactly adiabatic picture arise. At the same time, the relevance of the diabatic couplings increases, such that, e.g., for the case depicted in panel h) it is the diabatic coupling (green line) that dominates the total transition probability (red dots). Additionally, for small continuum state energies, the occurrence of energy crossings leads to marked maxima in the transition probabilities (cf. the peak in panel g).

Overall, it can be clearly seen that the different types of couplings are complementary to each other: If one of them is large, the other tends to be smaller and vice versa. The correct quantum mechanical treatment of our system requires both of them to be present; however, cases can be found in which one or the other is dominant.