Efficient Semiclassical Evaluation of Electronic Coherences in Polyatomic Molecules

Exposing a molecule to intense light pulses may bring this molecule to a nonstationary quantum state, thus launching correlated dynamics of electronic and nuclear subsystems. Although much had been achieved in the understanding of fundamental physics behind the electron-nuclear interactions and dynamics, accurate numerical simulations of light-induced processes taking place in polyatomic molecules remain a formidable challenge. Here, we review a recently developed theoretical approach for evaluating electronic coherences in molecules, in which the ultrafast electronic dynamics is coupled to nuclear motion. The presented technique, which combines accurate ab initio on-the-fly simulations of electronic structure with efficient semiclassical procedure to compute the dynamics of nuclear wave packets, is not only computationally efficient, but also can help shed light on the underlying physical mechanisms of decoherence and revival of the electronic coherences driven by nuclear rearrangement.


Introduction
Photoinduced molecular processes play a key role in physics, chemistry, and biology. In nature, light triggers a large variety of chemical reactions such as those involved in photosynthesis, vision, and the formation of vitamin D, but also can cause the radiation damage of biomolecules and photolysis. [1] The interaction of light with matter forms the basis of important technological applications such as solar cells in which photoinduced charge transfer and light harvesting are essential. [2] In all of these processes, molecules capture the energy of light and transform it into other forms, including electric energy, chemical energy, or heat. On a microscopic level, this energy conversion process is the result of a correlated motion of nuclei and electrons following photoionization or photoexcitation of a molecule. Therefore, understanding fundamental principles behind the light-induced electron and nuclear dynamics in molecules is of crucial importance if one wants to comprehend the diversity of natural phenomena.
Experimental studies of the ultrafast motion of electrons and nuclei in molecules are, however, extremely challenging because the molecules vibrate on a time scale of tens to hundreds of femtoseconds (1 fs = 10 -15 s) which corresponds to the atomic movement with the speed ≈ 1 km/s, while the motion of electrons takes place on an even faster, attosecond time scale (1 as = 10 -18 s). Clearly, experimental observations of such ultrafast molecular dynamics require advanced techniques with very high temporal and spatial resolution.
Fortunately, developments of coherent light sources in the beginning of the twenty-first century enabled the generation of subfemtosecond laser pulses with precisely controlled parameters. [3] Intense ultrashort laser pulses advanced the field of molecular physics by providing the researchers with a unique tool capable of initiating and tracing dynamics of both nuclei and electrons in a molecule with atomic spatial resolution and in real time. This progress made it possible to study the fundamental concepts such as electron-nuclear correlation, charge transfer and migration, electronic coherence as well as its decoherence and revivals due to nuclear motion.
In particular, recent experiments utilizing a sequence of an isolated attosecond pulse and an intense few-cycle infrared pulse have demonstrated the possibility to resolve the correlated elec-where is the time-dependent nuclear wave packet propagated on the PES associated to the I-th electronic state , and and denote electronic and nuclear coordinates, respectively. The electronic states, together with the corresponding PESs , are obtained by solving the stationary Schrödinger equation for the electronic subsystem for all possible nuclear configurations of a molecule. A major challenge in this approach is to describe the quantum mechanical behavior of nuclear wave packets in many dimensions and in the presence of nonadiabatic effects caused by couplings betweeen different electronic states.
The multi-configurational time-dependent Hartree (MCTDH) method is one of the most accurate and, at the same time, efficient approaches to solving the quantum nuclear dynamics taking place on several coupled electronic states. [16] MCTDH has been applied successfully to describe the nonadiabatic molecular dynamics [17] as well as the coherent electron-nuclear motion. [18,19] Although MCTDH takes into account all quantum effects, such as tunneling and nonadiabatic transitions, this rigorous technique scales exponentially with the dimensionality and, in addition, requires the expensive construction of global PESs.
To avoid the precomputation of globally fitted surfaces, various «direct dynamics» approaches have been proposed that calculate the PESs only along trajectories guiding a Gaussian basis, in which the evolving wave packet is expanded. The «on-the-fly» evaluation of the electronic structure guarantees that only the relevant regions of the configuration space are sampled. Whereas the variational multi-configurational Gaussians (vMCG) [20][21][22][23] are closest in spirit to MCTDH, many variations exist, including ab initio multiple spawning, [24] various Gaussian basis methods, [25][26][27] and more approximate semiclassical [28][29][30] and mixed quantumclassical [31][32][33] approaches.

Time-dependent Observable Properties
To study the electron motion in a given molecule, one must compute an observable property. Utilizing the BOH form of the molecular wavefunction, Eq.
tron-nuclear dynamics during the dissociative ionization of H 2 and D 2 molecules, [4] as well as to measure similar processes in other more complicated diatomic molecules. [5,6] In their pioneering experiments, Calegari and co-workers observed ultrafast charge migration in amino acids phenylalanine [7,8] and tryptophan [9] by measuring the dependence of the yield of doubly charged ions on the delay between the ionizing pump and doubly ionizing probe pulses. Kraus et al. [10] reconstructed and controlled the attosecond charge migration in ionized iodoacetylene by analyzing the high-harmonic generation spectrum emitted after irradiating the molecule with strong infrared pulses of different wavelengths (800 and 1300 nm). Ve ry recently, the decoherence and revival of ultrafast charge oscillations driven by nonadiabatic dynamics in neutral silane have been observed using X-ray attosecond transient absorption spectroscopy (ATAS). [11] Although the advances in experimental techniques have generated a lot of information about the dynamics and properties of molecules, a reliable interpretation of the experimental results often requires sophisticated theoretical simulations. As a consequence, the role of theory in the design and interpretation of new experiments has become much more important. [12] In this short review, we discuss how the electronic coherence and decoherence can be computed and analyzed with a recently developed methodology that combines accurate «on-the-fly» ab initio electronic structure calculations with efficient semiclassical dynamics of nuclear wave packets. We show that this efficient single-trajectory method estimates the electronic coherence time in polyatomic molecules accurately by comparing its results with the results of full-dimensional quantum calculations. Next, we explain how this method was, due to its computational efficiency, employed to scan a large number of molecules searching for new systems with long-lasting electron density oscillations. Then we describe an extension, which helped to clarify the physical mechanism of decoherence and of the revival of ultrafast charge oscillations in silane. Finally, we discuss several improvements that could increase the accuracy of such on-the-fly semiclassical simulations.

Theoretical Description of Electron-nuclear Dynamics in Molecules
By exposing a molecule to intense light pulses, one creates a superposition of several electronic and vibrational quantum states. To describe the time evolution of such an initial superposition, one must solve the time-dependent Schrödinger equation. Although nearly hundred years have passed since the discovery of this equation, such simulations remain extremely difficult for systems with more than a few degrees of freedom. To treat systems containing many interacting particles, one is, therefore, forced to introduce approximations.
A common starting point of many methods for simulating the dynamics of molecules is the Born-Oppenheimer-Huang (BOH) expansion of the molecular wavefunction. [13,14] This expansion exploits the fact that atomic nuclei are much heavier than electrons, which makes it possible to separate the fast electron motion from the typically much slower nuclear motion and to represent molecular dynamics as the dynamics of nuclei moving on the potential energy surfaces (PESs) formed by the electrons in a particular electronic state. Importantly, the BOH approach is formally exact if the set of electronic states used in the expansion of the molecular wavefunction is complete. It is only when this expansion is truncated that an approximation is made. [15] In the BOH picture, the electron dynamics results from the coherent superposition of multiple stationary electronic states: (1) difference between the involved electronic states I and J at the stationary geometry , and the second term is the reduced action quantifying the difference in frequency due to the divergent motion of wave packets on PESs I and J.
The modulation of electronic coherence by nuclear dynamics can be understood easily from the analytical expressions (6) and (7). In particular, the increasing phase-space distance between the two nuclear wave packets causes the decay of coherence and, at the same time, affects the frequency of electronic oscillations. If the positions of nuclei were fixed, the time scale of electronic dynamics would be determined by the energy gap between the two electronic states. Due to the nuclear motion, however, the constant frequency is modulated by the signed phase-space area between the two nuclear trajectories (see ref. [40] for more details).

On-the-fly ab initio Simulations of Electronic Coherences in Polyatomic Molecules
Let us now consider the quantum dynamics triggered by ionization of a molecule with an ultrashort laser pulse. We assume that the duration of the applied pulse is shorter than the electron correlation, i.e. the ionization event occurs on an infinitely short time scale, which allows us to employ the so-called sudden approximation. [46,47] In this approximation, the initial ionic state is obtained by projecting the electronic wavefunction resulting from the removal of an electron from the ground neutral state of a molecule onto its ionic subspace.
In general, the created initial state is not a stationary state of the cationic Hamiltonian and, therefore, will evolve in time. The ionization thus triggers ultrafast multielectron dynamics, which can be reflected in the propagation of the created hole along a molecular chain. This mechanism, driven exclusively by electrons, has been called "charge migration", [48,49] to avoid confusion with a more frequent "charge transfer" driven by nuclear motion. [50,51] A showcase example where the ultrafast ionization initiates the charge migration oscillations is the propiolic acid (HC 2 COOH). This molecule provides an excellent example for testing the semiclassical TGA because the electronic coherences in this system were already calculated with the exact quantum MCTDH method. [19] While dynamically exact, the MCTDH simulations employed an approximate, vibronic-coupling (VC) Hamiltonian model [15] to represent the PESs of the involved electronic states. Ab initio thirdorder algebraic diagrammatic construction [ADC(3)] method [52] with the standard double-zeta plus polarization (DZP) [53] basis set was used to compute the corresponding PESs. Assuming the ionization out of the highest occupied molecular orbital (HOMO), a coherent superposition of the first and the third cationic states of the molecule was created at the initial moment of time according to the corresponding hole-mixing weights [54] and employing the sudden ionization mechanism.
The exact quantum calculations confirm the important effects that nuclear motion has on the electronic coherence-the electronic oscillations are completely suppressed after 15 fs (see the red solid line in Fig. 1). [19] To estimate the importance of nonadiabatic dynamics, we repeated the quantum calculation on a reduced, "diabatic" version of the VC model, where the diabatic PESs are uncoupled. In this case, the TGA is exact and its results are equal to those of MCTDH method (yellow solid line in Fig.  1). The small deviations between the full and reduced versions of VC model are due to the nonadiabatic effects, but start to appear only at later times.
We also performed semiclassical simulations using adiabatic PESs obtained by diagonalizing the diabatic VC model (blue dashed line in Fig. 1). Although in this case the adiabatic form of the Hamiltonian does not provide additional insights about the dynamics of nuclear wave packets, it plays a role of an intermediate step between computations performed using the precomputed

Direct Dynamics Methods for Computing Electronic Coherence
Bearpark, Robb and co-workers [34][35][36][37] pioneered the use of trajectory-based techniques for treating the electronic dynamics coupled to nuclear motion. These authors propagated multiple trajectories representing an initially delocalized wave packet using the mean-field Ehrenfest approximation. Although Ehrenfest dynamics captures decoherence due to the destructive interference of coherent oscillations with different frequencies along different trajectories, it neglects completely the decoherence caused by the decay of overlaps of the nuclear wave packets moving on different PESs. Avoiding the mean-field approximation by allowing the nuclear wave packets launched on different PESs to evolve independently, the decay of nuclear overlaps was captured and the electronic coherence oscillations induced by the ionization of several molecules were analyzed with the direct dynamics versions of vMCG scheme. [38,39] Recently, a surprisingly simple semiclassical approach, in which the nuclear wave packet is approximated by a single timedependent Gaussian function, was used to evaluate the electronic coherence with an accuracy comparable to the full-dimensional quantum calculation. [40] In this approach, originally proposed by Heller, [41] the position and momentum of the Gaussian's center follow classical Hamilton's equations of motion, whereas the width and phase of the Gaussian are propagated using a time-dependent quadratic potential given by the local harmonic approximation of the original PES. Because the width of the Gaussian evolves in time, Heller's method has been called the "thawed Gaussian approximation" (TGA). [42][43][44] Expressed in position representation, the Gaussian wavepacket is where is the population of the corresponding electronic state, and are the phase-space coordinates of the center of the wave packet, is a complex symmetric width matrix with a positivedefinite imaginary part, and the complex number contains the dynamical phase (in its real part) and ensures the normalization at all times (by its imaginary part).
The Gaussian form of the nuclear wavefunction (5) permits an analytical evaluation of the integration in Eq. (4) for the electronic coherences. [45] To simplify the resulting expression, we disallow the time dependence of the widths of the Gaussians evolving in electronic states I and J (this is often called the "frozen Gaussian approximation"), thus obtaining: [40] Here, and are the spatial and momentum separations, respectively, between the centers of the two Gaussians at time t in mass-and frequencyscaled phase-space coordinates and , and is the classical action, where the first term is reponsible for electronic oscillations with the frequency given by the energy moiety with a period of 3.8 fs, corresponding to the energy spacing between the involved cationic states. The charge is eventually trapped by the nuclear motion and distributed along the molecular chain, but only after performing several clearly visible oscillations.
We would like to point out also the obvious limitations of the presented methodology. It is well known that the TGA cannot capture tunnelling and wave packet splitting, but these effects are typically unimportant at the ultrashort time scale on which decoherence occurs. In the present context, the main weakness of the TGA is its inability to take nonadiabatic transitions into account. To justify the negligibility of nonadiabatic effects in our simulations, we focused on molecules with large energy gaps between the ionic states. If conical intersections appear near the point of excitation, the subsequent wave packet dynamics can change dramatically (see, e.g., ref. [56]) and more refined methods, capable of describing nonadiabatic dynamics, must be used. Nevertheless, the ultrafast electron dynamics studied in this work occurs on subfemtosecond time scale, which usually implies that the nuclei do not have enough time to move far from their original positions and reach regions of the configurational space where the nonadiabatic couplings are strong. Unless the excitation or ionization places the PESs and those performed with the exact PESs obtained on the fly along the trajectories. When the wave packet is propagated with the exact Hamiltonian computed on the fly, it can capture anharmonicity effects and even reach regions inaccessible to the model VC Hamiltonian. As a result, the electronic oscillations evaluated with our on-the-fly calculations (green dash-dotted line in Fig. 1) have a similar period, but decay slightly faster than the electronic oscillations computed with the VC model. Because the nonadiabatic effects are here much smaller than the effect of approximating the potential energy surface, the semiclassical on-the-fly result is probably more accurate than the "exact" quantum result based on the VC model.
Simulations of electronic coherences with the semiclassical TGA require the propagation of a single nuclear Gaussian wave packet on each involved PES, while the wave packets moving in different states are independent from each other. The computational cost of the employed scheme is comparable to the conventional ab initio molecular dynamics (in particular, the costs are the same if the width of the Gaussian is fixed). Therefore, in this setting, the TGA not only demonstrates an accuracy comparable to full-dimensional quantum methods, but can also be very computationally efficient.
Encouraged by this observation, we applied the TGA to perform a massive scan of small polyatomic molecules and searched for systems exhibiting long-lasting electronic oscillations. We used a free online database PubChem maintained by the National Institute of Health to analyse the correlated electron-nuclear dynamics in about 250 randomly selected small molecules composed of C, H, O, and N atoms.
Our simulations demonstrate (see ref. [55] for more details) that the electronic coherences created after the ultrafast ionization are suppressed by the nuclear rearrangement within a few femtoseconds in most of the studied molecules. However, we were able to identify several molecules with electronic coherences lasting for as long as 10 fs, which are, therefore, promising candidates for future experiments. [55] Figure 2 shows electronic coherences and the dynamics of a positive charge generated after the sudden ionization of the HOMO electron of the but-3-ynal molecule. As one can see, the hole, initially localized around the carbon triple bond of the molecule, starts to oscillate between the triple bond and the aldehyde  [55] wave packet directly at a conical intersection, the TGA becomes exact in the limit of short propagation time and is, therefore, well suited for treating ultrafast processes.

Revealing the Physical Mechanism of Decoherence
As we discussed already in Sec. 2.2, the semiclassical Gaussian method makes it possible to reveal the physical mechanism of decoherence beyond the simple but vague justification by the nuclear motion. Here, we apply an extension of the technique developed in ref. [40] to interpret the recent experimental measurements of the decoherence and revival of the ultrafast charge oscillations in neutral silane. [11] The experiment, performed by the group of H. J. Wörner at ETH Zurich, utilizes an intense infrared 5.2-fs laser pulse to excite the molecule and an isolated sub-200-as X-ray pulse to probe the induced dynamics using the ATAS technique. [11] It has been shown [57] that the major ingredient required to simulate the ATAS is the electronic coherences. We therefore concentrate here on analyzing the electronic coherences between the excited states of a neutral silane created after the interaction of the molecule with the pump pulse.
Ab initio computations of the valence excited electronic states of silane demonstrated [11] that the corresponding PESs have both nonadiabatic and Jahn-Teller interactions which cannot be taken into account within the discussed semiclassical scheme. Therefore, fully quantum MCTDH method with a second-order VC Hamiltonian was used to reproduce the experimentally measured electronic coherences (see the Supplemental Information of ref. [11] for details). The simulations revealed that the experimentally measured signal results from the superposition of a pair of electronically excited states, which we shall denote simply as I and J. Even though the results of MCTDH calculations agree very well with the experimental measurements [11] and it is clear that the nuclear dynamics is responsible for the modulation of the coherence, the internal mechanism of the process is hard to deduce from sophisticated numerical quantum methods such as MCTDH.
To understand the mechanism of dephasing and revival of the electronic coherence observed in the experiment and reproduced by the accurate MCTDH simulations, we performed a theoretical analysis using a simple Eq. (6) for the overlap of Gaussian wave packets. Because our scheme cannot take the nonadiabatic transitions into account, we do not run the semiclassical propagation explicitly as we did in ref. [40] Instead, we analyze the fully quantum MCTDH wave packet dynamics using the formalism presented in Sec. 2.2.
To establish the correspondence between the quantum wave packet and its semiclassical representation, we set the position and momentum of the Gaussian by computing the expectation values of the corresponding operators for the quantum wave packet evolving in the electronic state I: Extracting in the same way the position and momentum of the nuclear wave packet evolving in the electronic state J, we can use Eqs. (6) and (7) to analyze the electronic coherences. We will focus on the magnitude of the electronic coherence, i.e., on the quantity , and will not discuss the phase . Figure 3 illustrates the MCTDH-calculated electronic coherence (panel A in Fig. 3) and the semiclassical analysis which elucidates the origin of the decay and revival of the electronic coherence between states I and J. In order for coherence between a pair of the electronic states to exist, the nuclear wave packets evolving in these states must overlap in both coordinate and momentum (8) space (see panel B in Fig. 3 for the spatial and momentum separations between the wave packets, and panel C in Fig. 3 for the corresponding contributions to the coherence). The comparison of the normalized magnitude of the electronic coherence (panel D in Fig. 3) demonstrates a good agreement between fully quantum MCTDH simulations and the semiclassical results. The semiclassical coherence slightly overestimates the quantum one because the evolving wave packets do not conserve their Gaussian form during their propagation on the corresponding PESs.

Conclusions and Outlook
We described a recently proposed efficient and simple onthe-fly semiclassical method for computing electronic coherences in polyatomic molecules. We also reviewed its applications exploring the effects that nuclear motion has on the electronic coherence in the propiolic acid and on the decoherence and revival of the coherence in silane. On one hand, the remarkable computational efficiency of the method made it possible to scan a large number of polyatomic molecules and search for systems with long-lasting oscillatory electronic dynamics. On the other hand, the simplicity of the semiclassical method (in comparison with sophisticated wave packet grid-based quantum methods) permitted us to disentangle several physical mechanisms contributing to decoherence and to the revival of coherence.  [11] .
The semiclassical TGA used in this paper has a plenty of room for improvement. Semiclassical on-the-fly simulations can be made more accurate with various extensions of the TGA, which can propagate not only a Gaussian wave packet but also a Gaussian multiplied with a linear [58,59] or general [43,60] polynomial. If combined with thermo-field dynamics, the TGA can account for finite temperature effects with almost no computational overhead. [45] Yet, because of its efficiency, the simple approach reviewed here should make it possible to analyze coherence and decoherence in larger, biologically relevant molecules. Since the TGA can treat molecules with hundreds of atoms, the method could be used to preselect molecules suitable for further experiments studying attosecond electron dynamics and contribute to the continuing discussions on the role of quantum coherence in biology. [61][62][63]