Search for long lasting electronic coherence using on-the-fly ab initio semiclassical dynamics

Using a combination of high-level ab initio electronic structure methods with efficient on-the-fly semiclassical evaluation of nuclear dynamics, we performed a massive scan of small polyatomic molecules searching for a long lasting oscillatory dynamics of the electron density triggered by the outer-valence ionization. We observed that in most of the studied molecules, the sudden removal of an electron from the system either does not lead to the appearance of the electronic coherence, or the created coherences become damped by the nuclear rearrangement on a time scale of a few femtoseconds. However, we report several so far unexplored molecules with the electronic coherences lasting up to 10 fs which can be good candidates for experimental studies. In addition, we present the full-dimensional simulations of the electronic coherences coupled to nuclear motion in several molecules which were studied previously only in the fixed nuclei approximation.


I. INTRODUCTION
The recent progress in laser technologies 1,2 made it possible to study properties of matter with unprecedented resolution. Ultrashort intense laser pulses revolutionized the field of atomic and molecular physics 3 providing the scientific community with a unique tool to initiate and trace dynamics of both electrons and nuclei of a molecule in real time and with atomic spatial resolution. 4,5 Exposing a molecule to a short light pulse may bring the system to a nonstationary quantum state, thus launching a coupled dynamics of electronic and nuclear wave packets. Coherent superposition of multiple electronic states often results in ultrafast evolution of molecular observables, such as the electron density. Although driven by purely electronic effects, 6,7 the dynamics of the electron density is strongly coupled to the nuclear motion in the system. It was demonstrated in numerous studies 8-14 that the slow nuclear rearrangement has a dramatic impact on the electronic dynamics which leads to the fast decoherence of electronic oscillations on a time scale of just a few femtoseconds.
Capturing the interplay between nuclear rearrangement and the ultrafast electron motion requires a concerted description of the electron-nuclear dynamics. One of the most powerful approaches allowing the accurate description of this truly molecular quantum dynamics is the multi-configurational time-dependent Hartree (MCTDH) method. 15,16 MCTDH has been used successfully for simulating nonadiabatic wave packet dynamics in molecules 17 and recently applied for the description of electronic coherence. 12,14 Although this rigorous technique, on the one hand, makes it possible to take into a) Electronic mail: nik.v.golubev@gmail.com account all the quantum effects, such as tunneling and nonadiabatic transitions, on the other hand, it suffers from an exponential scaling problem and also requires the costly construction of global potential energy surfaces (PESs).
To overcome the constraints of techniques utilizing the precalculated form of PESs, methods evaluating the electronic structure "on the fly" have been developed. These "direct dynamics" approaches calculate the PESs only along trajectories, thus sampling only the relevant regions of the configuration space and avoiding the precomputation of globally fitted surfaces. Direct dynamics techniques range from fully quantum methods, such as the variational multi-configurational Gaussians (vMCG), 18,19 ab initio multiple spawning, 20 coupled coherent states, 21 multi-configurational Ehrenfest, 22 and Gaussian dephasing representation, 23 to more approximate mixed quantum-classical and semiclassical approaches including, e.g., the surface hopping, 24 Ehrenfest dynamics, 25 and Herman-Kluk propagator, 26 together with its extensions. 27,28 Both the approximate Ehrenfest-based schemes [8][9][10][11] and the numerically exact vMCG method 13,29 have been used recently for computing the influence of nuclear structure and dynamics on the electronic coherences.
Despite the success of both MCTDH and trajectoryguided direct dynamics techniques in accurate description of the electronic dynamics coupled to nuclear motion, the above mentioned methods are still rather expensive computationally and thus only a few relatively small molecules have been studied so far. The very limited number of molecules analyzed to date is insufficient to make clear conclusions on how the molecular structure and the presence of specific functional groups influence the time scale of electronic coherence. At the same time, the possibility to measure experimentally the ultrafast electron motion in a molecule depends crucially on the number of oscillations which the electron density has time to perform before the electronic coherence is destroyed by the nuclear rearrangement. 30,31 Therefore, computational preselection of molecules having the desired properties, including the long lasting electronic coherence, is essential for the successful experimental studies. 32 Interestingly, the high-level quantum and multitrajectory methods are somewhat redundant for the treatment of electronic coherence in polyatomic molecules. In most of the previously studied systems, the electronic coherence becomes suppressed by the slow nuclear rearrangement within the first ten femtoseconds of dynamics. This very short time scale typically implies that the nuclei remain very close to their original positions. Accordingly, the simulations of ultrafast electronic dynamics under the influence of nuclear motion require the capturing of only the first initial instant of molecular rearrangement. The latter makes the usage of expensive high-level techniques, initially designed for simulating long lasting nuclear dynamics, both unnecessary and impractical for systems, in which one is interested primarily in the evolution of the electronic subsystem.
Recently, 33 it was demonstrated that a simple semiclassical approach, 34 in which the propagating nuclear wave packet is approximated by a single Gaussian function, can compete in accuracy with the full-dimensional quantum techniques for calculating the electronic coherence. In this semiclassical approach, the center of the Gaussian follows classical Hamilton's equations of motion while the width and the phase of the wave packet are propagated using the local harmonic approximation of the PES. Because the width of the Gaussian evolves in time, the approach was termed the thawed Gaussian approximation (TGA) in the literature. [34][35][36] Importantly, the TGA gives the exact solution of the time-dependent Schrödinger equation in the limit when the propagation time approaches zero and thus is particularly suited for the treatment of processes taking place on short time scales.
Here, we use on-the-fly ab initio semiclassical TGA to perform a massive scan of small polyatomic molecules searching for the long lasting electronic coherence. We concentrate on studying the electronic dynamics triggered by ionization of a molecule because the real-time tracing of a positive charge is currently more affordable from the experimental point of view. 30,31,37,38 Furthermore, we preselect only those molecules which demonstrate special features in their ionic spectra: namely (I) the presence of a large energy gap between the low-energy valence ionic states and the remaining ones, which can be important for the experimental possibility to create the desired superposition of only a limited number of electronic states, and (II) the presence of a strong correlation between valence electrons of the neutral molecule, the socalled hole-mixing phenomenon, 39 which guarantees the existence of non-trivial dynamics of the electron density resulting from the created superposition.
We use a free online database PubChem 40 maintained by the National Institute of Health as a source of molecules in our search procedure. PubChem contains millions of molecular structures with additional information, such as physical properties, toxicity data, molecular identifiers etc. The database has an efficient search interface to easily access molecules with the desired properties and allows one to download the chemical structures in many standard formats. We limit the search procedure to consider the subset of small neutral molecules composed of C, H, O, and N atoms. Since thousands of molecules fulfill the aforementioned conditions, we randomly select about 250 molecules for the subsequent simulations. Therefore, the major goal of this study is not to test all the molecules in the database but to propose an efficient procedure to scan a large variety of molecules and to find systems best suited for the experimental investigations.
In addition to searching for long lasting electronic dynamics in so far not investigated molecules, we perform calculations of electronic coherence in several experimentally interesting systems which were previously studied employing the fixed nuclei approximation, namely 2-propyn-1-ol, 41 2phenylethyl-N,N-dimethylamine (PENNA), 42 3-buten-N,N-dimethylamine (BUNNA), 43 and 3-methylen-4penten-N,N-dimethylamine (MePeNNA). 43 We show that the electronic coherence in these molecules becomes suppressed by the nuclear rearrangement within first few femtoseconds of the dynamics which is still enough to observe at least one clear oscillation of the electron density along a molecular chain.
The remainder of the paper is organized as follows. In Sec. II the theoretical basis of the methods used to compute electronic coherences is presented. In Sec. III, we discuss the methodology and computational details of the approaches used in this study. The results of our simulations, including the discovered molecules with interesting properties, are presented in Sec. IV. In Sec. V, we summarize the results and conclude.

II. THEORETICAL BACKGROUND
The starting point of our simulations is a closed-shell neutral molecule in its ground electronic and vibrational state. Applying an intense laser pulse with an appropriate energy range, one can kick out an electron from the molecule and produce an ion. Explicit modeling of the ionization process requires simultaneous accounting for both the long-range continuum wavefunction of the leaving electron and the intricate correlated short-range structure of the bound electrons which is currently beyond the reach for all but the smallest systems. [44][45][46] A practical way to describe the ionization process in realistic molecules is to employ the sudden approximation, 47 i.e., assuming an instant removal of an electron from the system. In this case, the ionic state is prepared by projecting the electronic state of the neutral molecule onto the ionic subspace of the system. Furthermore, we em-ploy the Franck-Condon approximation assuming that the nuclear wavefunction of the neutral ground state remains unchanged during its transfer to the ionic subspace.
After ionization, a single nuclear Gaussian wave packet on each involved ionic surface is propagated independently from the others. In our simulations, we concentrate only on those molecules in which the energy gaps between the ionic states remain large enough to justify the negligibility of the nonadiabatic population transfer. 48 This condition is controlled by tracing the energies of the corresponding ionic states along each of the propagated trajectories. Although this scheme inevitably leads to the omission of potentially interesting molecules, it allows us to guarantee the applicability of the TGA and thus high accuracy of the reported results.
In Sec. II A, we review briefly the electronic structure methods which we use for computing the ionic states. In Sec. II B, we describe the sudden ionization approximation and the hole-mixing mechanism. Section II C discusses a general procedure for computing observables. In Sec. II D, we present the equations of motion for the evolution of a Gaussian wave packet within the TGA. The theory part ends with Sec. II E, where we describe in detail a simple phase-space approach allowing us to analyze the electronic coherence in physically intuitive terms.

A. Electronic structure approaches
The centerpiece of our simulations is an accurate description of electronically excited states. As was explained in the introduction, we concentrate on studying the electronic dynamics in singly ionized molecules. From the computational point of view, the key difference in the description of neutral and ionized systems is the number of unpaired electrons in the ground electronic state. While most of the neutral molecules have an even number of electrons, thus forming a closed-shell ground electronic configuration, the removal of an electron during the ionization leads to the appearance of an open-shell species. Accurate description of the electronically excited states of an open-shell molecule is often case specific 49 and thus is not well suited for the automatic search procedure.
The latter difficulty can be efficiently circumvented by utilizing electronic structure approaches designed to describe (N − 1)-electron system starting from the N -electron reference ground state. To this end, several computational schemes, such as the outer-valence Green's function (OVGF) method, 47,50 two-particle-onehole Tamm-Dancoff approximation (2ph-TDA), 50 a socalled partial third-order (P3) 51 and its renormalized variants 52 schemes, as well as the algebraic diagrammatic construction (ADC) 53,54 and the equation-ofmotion coupled-cluster for ionization potential (EOM-CC-IP) 55-57 techniques, have been developed in last few decades. Among these methods, ADC and EOM-CC-IP have been proven to be both highly accurate and computationally efficient 58-60 which make them appropriate techniques for studying electronically excited states of ionized molecules.
The main weakness of the ADC and EOM-CC-IP techniques is the computation of nuclear gradients and Hessians required for the propagation of the thawed Gaussian wave packets. While nuclear gradients are available in EOM-CC-IP method at a large computational cost, 61 the procedure for computing analytically nuclear gradients in ADC and Hessians in both ADC and EOM-CC-IP has not been reported yet. For that reason, we utilize the time-dependent density-functional theory (TDDFT) 62,63 for a fast initial pre-screening of the electronic coherences created by ionization. The TDDFT is used only for a preliminary scan of large number of molecules and finding promising candidates, which are then re-analyzed using high-level ab initio ADC and EOM-CC-IP methods.

B. Sudden ionization and hole-mixing mechanism
Within the sudden approximation, the ionization event, i.e., the removal of an electron from a molecule, takes place on an infinitely short time scale. In practice, the time needed for the remaining electrons to respond to a sudden ionization was found to be about 50 attoseconds. 64 It was shown 64 that this time is universal, i.e., it does not depend on the particular system, and as such appears as the time scale of the electron correlation. Therefore, experimentally the sudden ionization can be achieved by applying an intense high-energy pulse with a duration shorter than the electron correlation time.
Mathematically, the removal of an electron from a system can by modeled by applying the annihilation operatorĉ i to the ground neutral state where index i refers to a particular molecular orbital of the neutral system and |Ψ N −1 i (0) is the initial state of the obtained ion at time t = 0. The initial state |Ψ N −1 i (0) is, in general, not an eigenstate of the cationic Hamiltonian and thus will evolve in time. The ultrafast multielectron dynamics triggered by the ionization can manifest as the time-evolution of the created hole along a molecular chain. Driven by purely electronic effects, this mechanism was termed "charge migration" 6,7 to distinguish it from a more common charge transfer driven by nuclei. 65,66 To understand the subsequent dynamics, let us expand the initial electronic wavefunction in a basis of ionic eigenstates |Ψ N −1 I : where a I are the expansion coefficients. Expanding the cationic eigenstates in a configurational interaction se-ries, (3) where c I k and c I kla amplitudes are known from ab initio calculations, one can solve the system of linear equations for the unknown coefficients a I in Eq. (2). To illustrate the procedure, we consider an idealized model situation when only two strongly correlated electrons occupying orbitals i and j of a neutral system contribute to the resulting ionic states. In this case, the electron correlation leads to the appearance of hole-mixing between corresponding orbitals in the ionic states 39 where, due to the orthonormality of the states, the two real coefficients c 1 and c 2 , representing hole-mixing amplitudes, satisfy the relation c 2 1 + c 2 2 = 1. The electronic wave packet formed from these two states at time t = 0 reads which can be expressed in terms of the ionization out of the ground state as Comparing Eqs. (1) and (6), it is seen that the ionization from the orbital i corresponds to the conditions a I c 1 + a J c 2 = 1 and a I c 2 − a J c 1 = 0 for the expansion coefficients and hole-mixing amplitudes. Solving this system of equations, one finds the expansion coefficients a I = c 1 and a J = c 2 . Similarly, if we ionize an electron from the orbital j, the initial state |Ψ N −1 j (0) will be superposition of the ionic states I and J with weights a I = c 2 and a J = −c 1 . Importantly, in the absence of electron correlation, i.e. when each of the ionic states is formed from a single hole configuration only, the hole will stay in the orbital in which it has been initially created and does not migrate. Therefore, many-electron effects play a central role in ultrafast electron dynamics triggered by sudden ionization of a molecule. 6

C. Charge migration analysis
To analyze and visualize the electron motion in a system, an observable property must be computed. Starting from the Born-Huang representation 67 of the molecular wavefunction where a I (t) is, in general, a time-dependent complex amplitude of the electronic state I, χ I (R, t) represents the normalized time-dependent nuclear wave packet propagated on the I-th PES, and r and R denote electronic and nuclear coordinates, respectively, the expectation value of a general operatorÔ(r, R) can be written as Here, O IJ (R) are matrix elements ofÔ(r, R) operator between electronic states I and J If both the operatorÔ(r, R) and the electronic states Ψ I (r, R) depend on nuclear coordinates R only weakly, Eq. (8) for the time-dependent expectation value can be reduced to where represent the populations of electronic states when I = J and the electronic coherences 12-14 when I = J. A convenient observable quantity to describe the electron motion triggered by the ionization of a molecule is the hole density, defined as the difference between the stationary electron density of the neutral system before ionization and the time-dependent electron density of the cation 6 : whereρ(r) is the one-body electron density operator. Explicit equations for constructing and analyzing the hole density in Eq. (12) within the ADC method can be found, e.g., in Refs. 39,41 .

D. Thawed Gaussian approximation
As one can see from the general Eq. (10) and a specific example in Eq. (12), the time dependence of the expectation value of an electronic operator is determined exclusively by the time dependence of the electronic populations and coherences χ IJ (t). Evaluation of the quantities χ IJ (t) via Eq. (11) requires, in turn, the knowledge of nuclear wave packets χ I (R, t) for all the involved states at each moment of time. Here, we perform the wave packet propagation with the TGA, 34 one of the simplest approaches to semiclassical wave packet dynamics in a general potential. We neglect the nonadiabatic transitions in our simulations and thus the populations of electronic states a I remain constant in time.
The Gaussian wave packet considered in TGA is given in the position representation as where R I t and P I t are the phase-space coordinates of the center of the wave packet, A I t is a complex symmetric width matrix with a positive-definite imaginary part, and γ I t is a complex number whose real part is a dynamical phase and imaginary part ensures the normalization at all times.
In the TGA, the wave packet is propagated in the effective time-dependent potential given by the local harmonic approximation (LHA) of the true potential V I (R) around the center R I t of the wave packet at time t.
Inserting the wave packet ansatz (13) and the effective potential (14) into the time-dependent Schrödinger equation, we obtain an equivalent system of ordinary differential equations for the evolution of wave packet parameters 34 where m is the real symmetric mass matrix and L I t denotes the Lagrangian Numerical solution of the system of differential equations (15) is straightforward, although additional transformations can be performed in order to make the integration of Eqs. (15c) and (15d) more stable (see, e.g., Ref. 68 ). Since the TGA requires only a single classical trajectory, along with the corresponding Hessians, it is particularly suitable for an on-the-fly implementation, where the local properties of the PES are obtained as-needed from an ab initio electronic structure calculations. In practice, Eqs. (15) can be solved in three stages: (I) the first two equations [Eqs. (15a) and (15b)] define a classical trajectory of a particle moving on the corresponding PES, (II) the Hessians are computed in parallel for every position in the classical path and (III) the reconstruction of the width and the phase of the wave packet is performed using Eqs. (15c) and (15d).

E. Semiclassical description of decoherence
The simple Gaussian form of the nuclear wave packet (13) imposed by the TGA makes it possible to perform the integration step in Eq. (11) for the electronic coherences analytically (see, e.g., Ref. 69 ) where we introduced the notation δΛ := Λ J − (Λ I ) * for the difference of tensors Λ J and (Λ I ) * , I and J denote the corresponding electronic states, and Λ I can be the matrix A I t , vector or scalar Evolution of the width of the wave packet within TGA allows one to account for stretching and compression of the Gaussian along the propagated trajectory. Although this additional flexibility of the wave packet provides, in general, a more accurate solution of the Schrödinger equation, it also complicates Eq. (17) and thus the analysis of the electronic coherences. Let us assume, for the sake of simplicity, that the widths of the Gaussians propagated in the electronic states I and J remain fixed along trajectories, i.e., A I t = A J t = iΓ, where Γ is a real positive-definite symmetric matrix. The latter makes it possible to further simplify Eq. (17) as (see Ref. 33 ) where is the phase-space distance between the centers of the two Gaussian wave packets in mass-and frequency-scaled coordinates and momentã R I and the phase term ∆S(t) is defined as The integral of Laplacian on each surface can be rewritten via Legendre transformation as Here, we used the relation H I (R I t , P I t ) = H I (R 0 , P 0 ) because the Hamiltonian is a constant of motion, and H I (R 0 , P 0 ) = V I (R 0 ) because P 0 = 0. Then, we can write the phase ∆S(t) as where is the energy gap between the electronic states I and J at R 0 and is the reduced action equal to the signed area within the closed curve C consisting of the two classical paths connecting (R 0 , P 0 ) with (R I t , P I t ) or (R J t , P J t ), and a straight line connecting (R I t , P I t ) and (R J t , P J t ). 33 The analytical expression (20) permits a simple semiclassical interpretation of the effect of nuclear dynamics on electronic coherence. The decay of coherence takes place due to the increase of distance between the nuclear wave packets in phase space, Eq. (21). Within the frozen Gaussian approximation the magnitude of the electronic coherence can be interpreted as a product of the wave-packet overlap in coordinate space, i.e., their spatial overlap, and the overlap of the wave packets in momentum space, which is referred to as dephasing. 13 The third mechanism responsible for the decoherence, namely the change of populations of electronic states, is not taken into account in our simulations. Besides the influence on the absolute value of the electronic coherence, the diverging nuclear trajectories affect the frequency of electronic oscillations. In the absence of nuclear motion, the time scale of electronic dynamics is defined by the energy difference between the corresponding electronic states, Eq. (26). Due to the nuclear motion, the static frequency is modified by the area S red (t) between two nuclear trajectories in the phase space, Eq. (27).
Importantly, ignoring the mode-mode mixing (Duschinsky rotation), 70 the term responsible for the decay of the electronic coherence in Eq. (20) can be decomposed into individual coordinate and momentum contributions associated with each of the vibrational modes of the system where the multiplication is made over all the vibrational modes, and [k] denotes k-th component of the corresponding vector. As one can see from Eq. (28), the loss of wave packets overlap in coordinate or momentum space within a single vibrational mode will lead to the overall decoherence (see also Ref. 12 ).

III. METHODOLOGY AND COMPUTATIONAL DETAILS
The screening procedure consists of the following four steps: (I) the molecule of interest is downloaded from the PubChem database; (II) the ionic spectrum of the molecule at a fixed reference geometry is computed utilizing the non-Dyson version of the ADC scheme 54,71 at the third order of perturbation theory [nD-ADC(3)]; (III) if the molecule has an appropriate ionic spectrum, namely both sufficient energy gaps between ionic states and strong hole-mixing, an initial prescreening of the electronic coherences is performed using the TGA combined with the TDDFT 62,63 ; (IV) the best candidates showing a long lasting coherence are verified by the TGA using the EOM-CC-IP 55,57 with single and double excitations (EOM-CCSD-IP).
Since the molecules in the PubChem database are collected from different sources, it is hard to estimate the reliability of molecular geometries presented in the database.
For this reason, we downloaded molecules in the simplified molecular-input line-entry system (SMILES) format, 72 which describes only the atoms connectivity in a molecule. Preliminary molecular geometries were constructed using OpenBabel 3.0.0 73 with the MMFF94 force field. 74 Ground-state geometries of the neutral molecules were optimized using the density functional theory (DFT) 75 at the wB97XD/6-311++G(d,p) 76 level. The optimization was performed with the Gaussian 16 package. 77 The noncorrelated reference Hartree-Fock (HF) orbitals for subsequent ADC computations were obtained using GAMESS-UK 7.0 package 78 with the standard double-zeta plus polarization (DZP) 79 basis set. The active space for every molecule was chosen to contain all the orbitals except 1s core states of all heavy atoms. All one-hole (1h) and two-hole-one-particle (2h1p) electronic configurations were taken into account in the calculations of ionic states and in the follow-up time-dependent density analysis.
The TDDFT prescreening was performed at the wB97XD/6-311++G(d,p) level of theory. The classical nuclear trajectories were propagated with the velocity Verlet algorithm, 80 where the time step was set at 0.25 fs for a total simulation time of 25 fs. The same classical propagation procedure with PESs computed using the EOM-CCSD-IP method with DZP basis set was applied to verify the TDDFT results for molecules that showed the long lasting electronic coherence. The on-the-fly evaluation of the PESs was done with the Q-Chem package. 81 The time-dependent nuclear wave packets were reconstructed from the classical trajectories using the single-Hessian variant of the TGA method. 82 The initial wave packets were prepared from the Hessian of the reference ground state of the neutral molecule and distributed between the involved electronic states according to the corresponding hole-mixing weights. The Hessians of the excited states were computed numerically at geometries corresponding to the vertical ionization from the reference ground state.

IV. RESULTS AND DISCUSSION
We have applied our search methodology to find polyatomic molecules in which the sudden ionization leads to the appearance of a long lasting electronic coherence. Although various classes of small organic molecules are present in the PubChem database, the hole-mixing is found to appear mostly in molecules containing an electron donor site, such as unsaturated carbon atoms or aromatic rings, and an electron acceptor site, e.g. aldehyde or amine functional groups. The presence of a strong hole-mixing in these types of molecules has been pointed out in previous works on this subject. 43,83 The results section is divided into two parts. In Sec. IV A, we first present molecules found in the Pub-Chem database which demonstrate the long lasting electronic coherences as well as clear charge migration oscillations throughout the molecular structure. In Sec. IV B, we present the simulations of the electronic coherences coupled to nuclear motion in molecules which were previously studied only at a fixed nuclear geometry.
A. Molecules with long lasting electronic coherence

But-3-ynal
Let us start with the but-3-ynal molecule which is composed of a chain of four carbon atoms with an alkyne group at one end and an aldehyde group at the opposite end. Interestingly, the oxygen and hydrogen atoms of the aldehyde group are displaced in the opposite directions from the molecular plane formed by the carbon atoms making the molecule asymmetric.
The valence molecular orbitals of the neutral but-3ynal and the ionization spectrum resulting from the removal of an electron from these orbitals are shown in panels (a) and (b) of Fig. 1, respectively. As one can see, the three lowest ionic states of the molecule are a mixture of one-hole contributions of the valence molecular orbitals. In particular, a strong hole-mixing is seen between the first and the third ionic states: an electron missing in the highest occupied molecular orbital (HOMO) [blue sticks in Fig. 1(b)] and an electron missing in the HOMO-2 [green sticks in Fig. 1(b)]. The electron density of the HOMO is localized primarily around the alkyne group, while the HOMO-2 is localized in the vicinity of the aldehyde group. The HOMO-1 [shown by orange sticks in Fig. 1(b)] is less correlated with the other orbitals and forms an ionic state lying between the pair of states corresponding to the ionization from the HOMO and HOMO-2. A sudden removal of an electron either from HOMO or from HOMO-2 will create an electronic wave packet, which will initiate charge migration oscillations between the carbon triple bond and the aldehyde group with a period of about 3.8 fs, determined by the energy gap between the first and the third cationic states.
We simulated the ionization from the HOMO of but-3-ynal populating the first and the third electronic states in proportion ≈81%/19%, respectively, according to the hole-mixing weights. Since the contribution of the HOMO orbital to the second ionic state is marginal, we do not include this state in the initial superposition. The phase between the electronic states is chosen in a way to localize the initial charge in the HOMO. The evolution of electronic coherence created by such a superposition is shown in Fig. 1(c). As one can see, the oscillations of the coherence gradually dephase within 10 fs due to the coupling to the nuclear motion. The period of the oscillations is slightly different from that predicted with the static nuclei. This is partially due to the difference between ADC and EOM-CC-IP electronic energies (see Sec. III for details) but also reflects the influence of nuclear motion on the time scale of electronic oscillations [see Eqs. (25) and (27), and also Ref. 33 ].
The hole density Q(z, t) computed along the molecular axis of the but-3-ynal is shown in Fig. 1(d). The quantity Q(z, t) was obtained by integrating Q(r, t), Eq. (12), over the x and y components, perpendicular to z, while the axis z was chosen to pass through the longest spatial extension of the molecule, thus called "molecular axis". As one can see from Fig. 1(d), the charge performs a couple of clear oscillations before being trapped by the nuclear motion and distributed along the molecular chain.
To ensure that the nonadiabatic effects do not lead to a significant population transfer between the neighboring electronic states, we monitored the energies of the three lowest states along the trajectories propagated in the first and the third states. As one can see from Fig. 2, the energy gaps between the states along each of the trajectories are sufficiently large to neglect the nonadiabatic transitions. Note that, in general, the population transfer can occur also between energetically distant adiabatic electronic states due to the momentum term in the Hamiltonian (see, e.g., Ref. 84 for more details). However, in our scheme the initial wave packets have zero momentum at the beginning of the propagation and the propagation time is short enough to assume the negligible influence of the nonadiabatic effects caused by the momentum operator.

Pent-4-enal
The pent-4-enal molecule is structurally similar to but-3-ynal but is composed of a chain of five carbon atoms of the I-th cationic state along the nuclear trajectory propagated on the J-th PES of the but-3-ynal molecule. The sufficient energy gaps between the involved states along both trajectories justify the negligible influence of the nonadiabatic effects on the dynamics.
with an alkene group at one end and an aldehyde group at the opposite end. Due to the lack of symmetry the molecule belongs to the C 1 point group.
The valence molecular orbitals of the neutral pent-4enal and the ionization spectrum are shown in panels (a) and (b) of Fig. 3, respectively. Similarly to the but-3ynal, the strong correlation between valence electrons of the molecule leads to the appearance of the hole-mixing between the ionic states. As one can see, the two lowest ionic states are a mixture of the one-hole contributions of the HOMO and HOMO-1, while the third and the fourth states are composed from the HOMO-2 and HOMO-3.
We simulated the ionization from the HOMO populating the first and the second electronic states in proportion ≈87%/13%, respectively, according to the hole-mixing weights. The evolution of electronic coherence created by such a superposition is shown in Fig. 3(c). Similarly to the but-3-ynal, the electronic coherence dephases within 10 fs due to the nuclear rearrangement. However, since the energy gap between the involved states in the case of pent-4-enal is smaller than that of the but-3-ynal, the electronic coherence performs a smaller number of oscillations within the same time. As for the ionization from the HOMO-2 and HOMO-3, the small energy gap between the third and the fourth ionic states makes the oscillations of the electron density so slow that the decoherence takes place faster than the period of a single oscillation.
The hole density Q(z, t) computed along the molecular axis of the pent-4-enal is shown in Fig. 3(d). As one can see, the charge oscillates between the double bond and the aldehyde group of the molecule before being trapped by the nuclear motion.

2,5-dihydrofuran and 3-pyrroline
Interesting examples of a strong electron correlation between an unsaturated carbon bonds and an electron acceptor atom can be found in the cyclic molecules 2,5-dihydrofuran and 3-pyrroline. Both structures contain a five-atom ring with a heteroatom (oxygen in 2,5dihydrofuran and nitrogen in 3-pyrroline) opposite to a carbon double bond. The 2,5-dihydrofuran molecule belongs to the C 2v symmetry, while the additional hydrogen atom in 3-pyrroline reduces its symmetry to the C s point group.
The valence molecular orbitals of the neutral molecules and the ionization spectra of 2,5-dihydrofuran (Fig. 4) and 3-pyrroline (Fig. 5) are shown in panels (a) and (b), respectively. As one can see, in both cases the first two cationic states consist of a strong mixture of the HOMO and HOMO-1 in similar proportions. The HOMOs of both molecules are localized in the double bond while the HOMOs-1 are localized in the vicinity of the heteroatom.
The electronic coherences appearing after the removal of the HOMO electron from 2,5-dihydrofuran and 3pyrroline are shown in panels (c) of Figs. 4 and 5, respec-tively. In both molecules the coherence performs several oscillations before being destroyed by the nuclear rearrangement. The hole densities Q(z, t) computed along the molecular axes passing through the midpoint of the double bond and the heteroatom are shown in panels (d) of Figs. 4 and 5. In both molecules the charge initially localized in the double bond oscillates from one side of the system to the other with a period of about 3 fs before being trapped and distributed along the molecule.
The 2,5-dihydrofuran and 3-pyrroline are interesting candidates for experimental studies involving X-ray transient absorption measurements 85 because of the different heteroatoms contained in the molecules. By taking advantage of element-specific core-to-valence transitions induced by X-ray radiation, one can trace the dynamics of electron density with atomic spatial resolution. Depending on the available laser setup, the energy window corresponding to the absorption by the O or N atom in 2,5-dihydrofuran and 3-pyrroline, respectively, can be exploited.

B. Molecules from the literature
Let us turn to the analysis of the electron-nuclear dynamics in 2-propyn-1-ol, PENNA, BUNNA and MePeNNA molecules which were previously studied employing the frozen nuclei approximation (see Refs. [41][42][43]. Similarly to most of the molecules found in the present study, the strong electron correlation between molecular orbitals localized around donor and acceptor sites of the neutral 2-propyn-1-ol, PENNA, BUNNA and MePeNNA leads to the appearance of the hole-mixing in the ionic states. The ionization spectra of all the molecules are somewhat similar [41][42][43] and represent a few lines corresponding to the ionization of the outer-valence electrons separated by an energy gap from the remaining ionic states. To simplify the comparison of the electronic coherence for different molecules, here we assume that only two states with largest hole-mixing amplitudes become equally populated for each of the studied systems. Figure 6 shows the electronic coherence in the four considered molecules for various combinations of the electronic states. A convenient quantity allowing a comparison of the coherence times of different molecules is the purity function Tr[ρ(t) 2 ], where the electron density matrix ρ(t) is related to the matrix of nuclear overlaps, Eq. (11), by transposition: ρ IJ (t) = χ JI (t). Due to decoherence, the purity decays from the value Tr[ρ(0) 2 ] = 1 for the initially pure state to the value 1/n for the equally weighted mixture of n states. The bottom panel of Fig. 6 demonstrates that the electronic coherences in BUNNA, PENNA and MePeNNA are fully suppressed within about 5 fs (green, red and orange solid lines in Fig. 6, respectively). Due to the small energy gaps between the involved electronic states, the electron density has time to perform only a single oscillation in the case of PENNA and MePeNNA molecules (see the top panel of Fig. 6). Although the decoherence takes place on the similar time scale in the BUNNA molecule, the large energy gap between the states makes it possible to observe more oscillations of the electron density in this case. The electronic coherence in the 2-propyn-1-ol survives the nuclear motion for about 8 fs which makes it possible to observe a couple of oscillations of the electron density in this molecule.

V. CONCLUSIONS AND OUTLOOK
In this paper, we have performed a thorough theoretical scan of a large number of small polyatomic molecules searching for specific structural and dynamical properties which can be useful for the experimental measurements of ultrafast electronic dynamics and its coupling to the nuclear motion. In particular, we were interested to find molecules with long-lasting electronic oscillations, surviving the decoherence caused by the nuclear rearrangement for the longest possible time. We concentrated on studying the correlated many-electron dynamics induced by the outer-valence sudden ionization with short intense laser pulses. Using a combination of high-level electronic structure techniques with efficient on-the-fly semiclassical description of nuclear dynamics, we have analyzed in total about 250 molecules with potentially interesting properties. Our results show that in most of the studied molecules, the sudden ionization either does not lead to a superposition of states trough the hole-mixing mechanism, or the obtained electronic coherences become damped by the slow nuclear dynamics on a time scale of a few femtoseconds, which can make experimental measurements of laser-induced electron motion in these systems problematic. Yet, we have found a long lasting electronic coherences in several new molecules, which have thus become promising candidates for experimental studies. In addition, we have performed the full-dimensional simulations of electron-nuclear dynamics in 2-propyn-1ol, BUNNA, PENNA, and MePeNNA molecules which were previously studied only in the static nuclei approximation. We observed that the electronic coherence in the 2-propyn-1-ol lasts for about 8 fs before being captured by the nuclear motion, while in PENNA, BUNNA and MePeNNA the electron dynamics is fully suppressed within 5 fs.
We would like to emphasize again that we have concentrated here on studying only those molecules which demonstrated specific features in their ionic spectra. In particular, the dynamical simulations were carried out only for the molecules with a strong hole-mixing between the lowest ionic states. At the same time, there exist various other mechanisms of the correlation-driven ultrafast electron dynamics in molecules 7,39 which we do not investigate in the present study. In addition, we excluded molecules with the strong nonadiabatic effects which can not be accurately treated within the employed computational scheme. These systems, however, can be promising candidates for studying the appearance of electronic coherences in the vicinity of a conical intersection, 86 as well as for exploring the transfer of electronic coherence between electronic states due to the nonadiabatic processes. 87 Nonetheless, all of the reported molecules demonstrating a long lasting electronic coherence are readily available in the market and, thus, can be useful candidates for studying the ultrafast many-electron charge migration dynamics.
The semiclassical vertical-Hessians TGA used in this paper can be further improved by calculating Hessians along the propagated trajectory and thus take into account more complicated situations, e.g., dissociation of a molecule. Moreover, the improved versions of the TGA such as the extended thawed Gaussian approximation (ETGA), 88 which propagates a Gaussian wave packet multiplied by a general polynomial, or a so-called three thawed Gaussians approximation (3TGA), 89 benefiting from representation of the wave packet by multiple Gaussians, were recently reported which can make on-the-fly semiclassical simulations even more accurate.
Finally, we would like to point out that the efficient approach used in this work opens the door to the analysis of electron-nuclear dynamical processes in larger, biologically relevant systems. Being able to treat molecules with a few hundred atoms, the TGA technique combined with the appropriate electronic structure method can help shed light on the continuing debates on the role of quantum coherence in biology, 90-92 quickly preselect molecules suitable for further experimental investigations, and support theoretically recent experimental observations of attosecond electron dynamics in realistic molecular systems. We hope that our work will motivate such studies.