Ultrafast optical techniques allow us to study ultrafast molecular dynamics involving both nuclear and electronic motion. To support interpretation, theoretical approaches are needed that can describe both the nuclear and electron dynamics. Hence, we revisit and expand our ansatz for the coupled description of the nuclear and electron dynamics in molecular systems (NEMol). In this purely quantum mechanical ansatz, the quantum-dynamical description of the nuclear motion is combined with the calculation of the electron dynamics in the eigenfunction basis. The NEMol ansatz is applied to simulate the coupled dynamics of the molecule NO_{2} in the vicinity of a conical intersection (CoIn) with a special focus on the coherent electron dynamics induced by the non-adiabatic coupling. Furthermore, we aim to control the dynamics of the system when passing the CoIn. The control scheme relies on the carrier envelope phase of a few-cycle IR pulse. The laser pulse influences both the movement of the nuclei and the electrons during the population transfer through the CoIn.

## I. INTRODUCTION

The continuous development of attosecond laser pulses enables spectroscopic techniques, which allow for the time resolved investigations of ultrafast photo-initiated processes in atoms, molecules, and solids. Nowadays, it is possible to study electronic correlation and ultrafast molecular dynamics through pump–probe experiments.^{1–7} Within these experiments, attosecond, broad-band pulses are used to generate electron wavepackets in highly excited states of molecules, leading to the discovery of effects such as electron localization in diatomic molecules^{3,8} and, later, of purely electronic charge migration in biologically relevant molecules.^{5–7} To explain and interpret the observations of these experiments, theoretical approaches are needed that can describe the dynamics of electrons in molecules. Most approaches use time-dependent analogs of well-established quantum-chemical methods such as time-dependent Hartree–Fock theory (TD-HF)^{9} or time-dependent density-functional theory (TD-DFT).^{10} Furthermore, time-dependent post-Hartree–Fock methods such as time-dependent configuration-interaction (TD-CI),^{11,12} time-dependent coupled-cluster (TD-CC),^{13,14} and multi-configuration time-dependent Hartree–Fock^{15} are available for the correlated description of electron dynamics in molecular systems. In other theoretical approaches, the electronic wavefunction is propagated directly in time with the help of Green’s function^{16} or in the basis of molecular orbitals.^{17} All these theories focus on the evolution of the electronic subsystem driven by electronic correlation^{18,19} and predict long-lived coherences. The neglect of the nuclear motion is justified by the assumption that the dynamics of the electrons is much faster than that of the heavier nuclei. This results in charge migration, an oscillatory motion of electron density with frequencies defined by the energy gaps among the states populated with the initial laser pulse. If the states of the superposition are close together, the electron dynamics becomes slow, and therefore, the nuclear motion can no longer be neglected. However, as shown in numerous theoretical works,^{8,20–25} nuclear motion, in general, causes decoherence in molecular systems and should not be neglected in any cases. This decoherence causes the electronic wavepackets to exist only for short time scales.^{24} For small systems such as $H2+$ or $D2+$, a full quantum treatment of the coupled electron and nuclear dynamics is possible.^{26} Beyond these three particle problems, there are computationally very demanding methods available based on a multi-configurational ansatz^{27} or on the coupled description of nuclear and electronic flux.^{28,29} Further techniques are based on the coupled propagation of the nuclear and electronic wavefunctions on a single time-dependent potential energy surface.^{30–33} However, for larger molecular systems, the main techniques used are mixed quantum classical representations.^{34–37} For example, the electron dynamics is described using TD-DFT and the nuclear motion is considered using an Ehrenfest approach.^{34,35} However, these methods do not reflect the quantum nature of the nuclei, which becomes important for ultrashort pulse excitation and non-adiabatic transitions.

In this paper, we want to revisit and expand an ansatz for the coupled description of the nuclear and electron dynamics in molecular systems^{8,38,39} (shortened NEMol) developed in our group. It is based on electronic structure calculations and nuclear quantum dynamics. In its initial formulation, the electronic wavefunctions are represented as Slater determinants and propagated in the eigenstate basis. The coupling of the nuclear motion to the electron motion is incorporated explicitly through the nuclear wavepacket motion as well as through a coherence term with contributions from the nuclear and electronic wavefunctions. Compared to similar approaches,^{30–33} the feedback of the electron motion to the nuclear dynamics is less directly introduced by simulating the nuclear dynamics on coupled potential energy surfaces (PES’s). The central equation of the original NEMol ansatz^{8,38,39} relates the dynamics of the coupled one-electron density to the temporal evolution of the expected value of the nuclear positions. In the first part of this work, we want to generalize the NEMol ansatz by extending beyond this single geometry approximation. Therefore, we introduce the NEMol-grid in order to represent the electron dynamics at multiple points on the grid used for the nuclear wavepacket propagation. In the limit, the NEMol-grid is equal to the grid representing the nuclear wavepacket, but in practice, we choose a coarser one. By means of a simple approximation, it is possible to obtain a condensed representation of time-dependent electron density in the one-electron-two-orbital (1e-2o) picture.

In the second part, we want to explore the potential of our NEMol ansatz. For this purpose, we consider a situation that can generate coherent electron dynamics in excited states of molecules even without a laser pulse present. Such a scenario occurs in the vicinity of a conical intersection (CoIn).^{40–44} For this ubiquitous but nevertheless extraordinary points in a molecular system, the adiabatic separation between nuclear and electronic motion breaks down^{40,45,46} and the electronic states involved become degenerate. Beside the creation of funnels for radiationless electronic transitions, a coherent electron wavepacket is created whose dynamics approaches the time scale of the nuclear dynamics. All these properties of CoIn’s are determined by the shape and size of the non-adiabatic coupling elements (NAC’s) and the topography of the vicinity. As a realistic molecular system that provides such a situation, we have chosen the NO_{2} molecule. After excitation into the first excited state, a CoIn enables an ultrafast non-adiabatic transition back to the ground state within less than 100 fs. This fast relaxation and the photophysics of NO_{2}, in general, have been widely explored both theoretically^{47–59} and experimentally.^{56,60–67} Beside the free relaxation of NO_{2}, we also studied the influence on the coupled electron dynamics when applying a few-cycle IR laser pulse in the vicinity of the CoIn. The variation of the carrier envelope phase *ϕ* (CEP) of such a few-cycle pulse offers the possibility to steer electrons and nuclei.^{56,57,59,64,68–76} Similar to previous studies,^{56,57,59,64} we apply this CEP-control-scheme to NO_{2} and evaluate the CEP-dependence of the resulting coupled nuclear and electron dynamics.

## II. COUPLED NUCLEAR AND ELECTRON DYNAMICS (NEMol)

In the original NEMol ansatz,^{8,38,39} the coupled one-electron density *ρ*(*r*, *t*; ⟨*R*⟩(*t*)) is defined according to Eq. (1). For convenience, the detailed derivation of this equation can be found in the Appendix adapted to the current notation

with

The first summation consists of the state specific electronic density *ρ*_{jj}(*r*, *t*; ⟨*R*⟩(*t*)) weighted with the corresponding time-dependent population *A*_{jj}(*t*). The second summation defines the coherent contribution to the coupled electron density and consists of the time-dependent overlap *A*_{jk}(*t*), the one-electron transition density *ρ*_{jk}(*r*, *t*; ⟨*R*⟩(*t*)), and its pure electronic phase defined by the energy difference Δ*E*_{jk} between the electronic states involved. All quantities related to the electronic wavefunction are calculated for one nuclear geometry per time step, which is defined by the time-dependent expected value of the position ⟨*R*⟩(*t*) (for definition, see the Appendix). As long as we are focusing on situations with quite localized wavepackets and/or one-dimensional systems,^{8,38,39} this approximation works quite well. However, in order to treat higher dimensional systems and more complex processes, we want to generalize the NEMol ansatz in this work. To extend the ansatz, the integration over the full nuclear coordinate space is split up in segments to improve the resolution of the spatial dependence of the electronic phase term. For this purpose, a second grid, the NEMol-grid, is introduced. The resulting modified NEMol ansatz is described using an exemplary system with two nuclear coordinates *c*_{1} and *c*_{2}. The complete two-dimensional coordinate space is split up into *M* × *L* segments defined by their boundaries *m*_{min}, *m*_{max} and *l*_{min}, *l*_{max}. For each of these segments *ml*, the population terms $\alpha jjml(t)$ and the overlap terms $\alpha jkml(t)$ are calculated,

The sum of these segment terms results in the corresponding total population and overlap,

At the center *R*_{ml} of each segment, the state specific electronic densities, the one-electron transition densities, and the eigenenergies are determined, and with these values, the coupled one-electron density for each segment *ρ*_{ml}(*r*, *t*; *R*_{ml}) is calculated,

with

It should be noted that for each segment, the Δ*E*_{jk} values and the electron densities are no longer dependent on ⟨*R*⟩(*t*). In contrast to the original NEMol ansatz, now many Δ*E* values are simultaneously contributing to the overall electron dynamics. They are addressed whenever the nuclear wavepacket is located there. To obtain the total coupled electron density, the individual contributions of each segment are summed up,

This total coupled electron density *ρ*(*r*, *t*; *R*) describes the electron dynamics coupled to multiple grid points on which the nuclear wavepacket is represented.

A second aspect that we would like to introduce is a further simplification. For clarity reasons, it is here formulated in terms of the original NEMol ansatz. We now consider a system of two electronic states described by their electronic wavefunctions *φ*_{1} and *φ*_{2}. In the simplest case, the wavefunctions of both states are described by two Slater determinants, which only differ in the occupation of one spin orbital *θ*. Now, the coupled total electron density can be simplified by expressing the densities and transition densities using the spin orbitals,

The summation at the beginning includes the densities of all equally occupied orbitals and is followed by the densities of the remaining two orbitals *θ*_{1} and *θ*_{2} weighted with the populations *A*_{11}(*t*) and *A*_{22}(*t*). The coherent part contains the product of the orbitals *θ*_{1} and *θ*_{2}. Within this simplification, it is now possible to neglect the contributions of the equally occupied orbitals in order to study the coupled electron dynamics in a one-electron-two-orbital (1e-2o) picture. Under the above mentioned approximation, this 1e-2o picture is a possibility to examine the coherent part of the electron dynamics in a very condensed way. This simplification can also be made in combination with the NEMol-gird.

## III. NO_{2} COUPLED DYNAMICS

We apply our extended NEMol approach to the non-adiabatic dynamics of NO_{2}. In this molecule, a CoIn [depicted in Fig. 1(b)] between the *D*_{1} and the *D*_{0} state enables radiationless relaxation. The ultrafast non-adiabatic transition takes less than 100 fs and has been widely explored both theoretically^{47–59} and experimentally.^{56,63–67} First, we analyze the relaxation itself, and next, we apply a few-cycle IR laser pulse to control the dynamics in the vicinity of the CoIn, similar to previous studies.^{56,57,59,64} With our NEMol ansatz, we can study its influence on the motion of the nuclei and the electrons.

The nuclear dynamics is performed on the two-dimensional adiabatic potential energy surfaces of the *D*_{1} and the *D*_{0} state shown in Fig. 1. The coordinates spanning the PES’s are the gradient difference and derivative coupling vectors defining the branching space of the *D*_{1}/*D*_{0}-CoIn depicted in Fig. 1(b). These two vectors correspond to the bending angle *α* and the asymmetric stretching coordinate *b*, defined as half the difference between the two NO distances. The last internal degree of freedom, the symmetric stretch coordinate, is kept constant at the value of the optimized *D*_{1}/*D*_{0}-CoIn (1.267 Å). As shown by Richter *et al.*,^{57} the population dynamics obtained within this two-dimensional coordinate space is in very good agreement with the full dimensional simulations performed by Arasaki *et al.*^{56} We performed our dynamics simulations in the adiabatic representation, and the corresponding NAC’s between *D*_{1} and *D*_{0} are shown in Fig. 1(c). It should be mentioned that in previous studies,^{56,57,59} the simulations were performed in the diabatic representation, and therefore, small deviations may occur due to the limitation of the grid spacing. Further information about the simulation setup can be found in Sec. II of the supplementary material.

In order to calculate the coupled electron density according to Eq. (7), we define a NEMol-grid of 15 × 13 points that are equally distributed between 1.34 and 2.86 rad in the *α*-coordinate and between −0.33 and 0.33 Å in the *b*-coordinate. The necessary population- and overlap-terms are calculated for equal-spaced segments around these grid points. To cover the entire PES, the segments for the boundary grid points are larger. The transformation of the full wavepacket onto the NEMol-grid, the overlap terms, and the resulting coherence terms are visualized in Fig. S6 (free propagation) and Fig. S10 (propagation with laser pulse) of the supplementary material. The two active orbitals that are required to describe the NEMol-dynamics in the one-electron-two-orbital (1e-2o) picture are shown in Fig. 2 at the optimized CoIn. The non-binding orbital *n*_{N} with contributions at the nitrogen atom is associated with the *D*_{1} state, and the non-binding orbital *n*_{O} located only at the oxygen atoms is attributed to the *D*_{0} state. The energy difference Δ*E* between the *D*_{0} and *D*_{1} state for each grid point is shown in Fig. S4 of the supplementary material.

### A. Free dynamics of NO_{2}

To initiate the dynamics simulation in the *D*_{1} state, we assumed a delta pulse excitation. The temporal evolution of the population of both states is shown in the upper panel of Fig. 3, and the dynamics of the nuclear wavepackets integrated over the *α*-coordinate and the *b*-coordinate is depicted in Fig. S5 for both surfaces. The nuclear wavepacket started in *D*_{1} reaches the vicinity of the CoIn after ∼7 fs for the first time. While passing the coupling region in the time interval from 7 to 15 fs, the population of the electronic ground state increases to over 60%. The part of the nuclear wavepacket remaining in the *D*_{1} state reaches its turning point around 15 fs and then propagates backwards. This leads to a second passage through the CoIn area and an increase in the population of the *D*_{0} state around 22 fs. The nuclear wavepacket evolving on the lower adiabatic surface re-encounters the CoIn region later at around 30 fs. During this third passage, a substantial part of the population is transferred back into the excited state. After 35 fs, the wavepacket is delocalized on both surfaces and the population is nearly equal in both states. Toward the end of the simulation at around 50 fs, a fourth passage occurs. The wavepacket remains symmetrical with respect to the *b*-coordinate for the whole simulation time. For the wavepacket on the lower PES [see right sight of Fig. S5(b) of the supplementary material], the formation of a nodal structure for *b* = 0.0 Å is clearly visible, which is a signature of destructive self-interference due to the geometric phase effect.^{45,46,77}

In the lower part of Fig. 3, snapshots of the electron density in the 1e-2o picture are shown. For a better visualization, also the difference in density with respect to *t* = 0 fs is depicted. The molecule is orientated in such a way that the molecular plane is equivalent to the *yz*-plane and the center of mass defines the origin of the laboratory frame. Therefore, the internal *α*-coordinate points to the same direction as the *y*-coordinate and the internal *b*-coordinate is associated with the *z*-coordinate. The orientation of the molecule is shown in the upper right corner of Fig. 4. In correspondence to the non-adiabatic transition from the *D*_{1} state to the *D*_{0} state, the main feature of the electron dynamics is the loss of density at the nitrogen atoms and the corresponding gain of density at the oxygen atoms. In addition, the change in the electron density attributed to the motion of the nuclei (Born–Oppenheimer part) is present. Due to the high symmetry of NO_{2}, the electron density is mirror-symmetrical with respect to the *xy*-plane, which is equivalent to the symmetric behavior of the nuclear wavepacket with respect to the *b*-coordinate.

To analyze the electron dynamics, we calculated the dipole moment of the electron density within the 1e-2o picture. In the upper panel of Fig. 4, the temporal evolution of its three components is shown; for the molecular orientation, see the upper right corner of Fig. 4. To distinguish the Born–Oppenheimer part of the dynamics from the coherent electron dynamics, the density was calculated once with the coherent part included and once without. For both quantities, the respective dipole moments were determined as well as their difference, hereinafter labeled as Δ 1e-2o, and are shown in the lower panel of Fig. 4. The active orbitals do not change along the *x*-coordinate, and thus, the 1e-2o-*x*-component of the dipole moment stays zero and is excluded from further discussions. The 1e-2o-*y*-component shows the largest values and the strongest changes over time. Its evolution follows the dynamics of the population. In the initial 20 fs, the first passage through the CoIn region occurs and, simultaneously, the value of the 1e-2o-*y*-component changes from 0.3 to −0.3 a.u. The zero crossing occurs at 10 fs. For later times, when dephasing and partial recurrence of the nuclear wavepackets become important, the *y*-component approaches zero at about 40 fs and becomes negative thereafter again. These main features disappear for the Δ 1e-2o-*y*-component (lower panel Fig. 4), and only fast oscillations with one order of magnitude smaller amplitudes are left. The largest amplitudes are observed around 10, 30, and 50 fs. These amplitudes coincide with the passages of the wavepacket through the CoIn region. The large difference between the 1e-2o and the Δ 1e-2o value means that the dynamics of the *y*-component is dominated by the nuclear motion. That is understandable since the *y*-coordinate is aligned along the main direction of dynamics (*α*-coordinate), which mediates the non-adiabatic transition. The temporal evolution of the 1e-2o-*z*-component is an order of magnitude smaller and almost identical to its Δ value. The dynamics of the *z*-component is not dominated by the nuclear motion but solely induced by the coherent electron dynamics. Therefore, we can use the *y*- and the *z*-component to distinguish between the two contributions of the coupled electron dynamics. As the Δ values of both components lie amplitude-wise in the same region and show a similar pattern, they are suitable to monitor the coherent electron dynamics in the system. Overall, the nuclear motion has a much larger impact on the dipole moment than the coherent electron dynamics.

By applying the Fourier transform to the temporal evolution of the dipole moments, the corresponding frequencies are determined. Beside the Δ 1e-2o- and 1e-2o-components, the dipole moment calculated with the full density was also used. The resulting spectra for the *y*- and the *z*-component for all three cases are shown in Fig. 5. The spectra are all normalized to one individually. The relative magnitude between all quantities can be estimated from Fig. 4. All frequencies with an intensity larger than 0.1 are listed in Tables S3 and S4 of the supplementary material.

The Δ 1e-2o spectra (Fig. 5, blue), reflecting the coherent electron dynamics, cover the largest frequency range from 0.2 to 2.3 eV for both components, whereas the energy differences Δ*E* (0.0–1.0 eV, see Fig. S4) in the vicinity of the CoIn, which enter in the coherent part of the electronic wavepacket, are smaller. These discrepancies can be rationalized when taking a closer look at the coherence term [see Eq. (1)]. Two of the factors in the product contribute to the overall phase, the nuclear overlap and the electronic phase term containing the Δ*E* values. The phase of the overlap term relates to the difference in momentum of the nuclear wavepackets involved. In our test system NO_{2}, the wavepacket on *D*_{1} approaches the CoIn with a high momentum, larger than the Δ*E* gaps near the CoIn. In other words, the coherent dynamics of the electronic wavepacket is in the NO_{2} case also significantly influenced by the phase-differences of the nuclear wavepackets moving on different potentials. This correlation is illustrated in Fig. S7 of the supplementary material for two individual NEMol-grid points. The frequencies for the 1e-2o-components (Fig. 5, red) are dominated by the slower nuclear dynamics (Born–Oppenheimer part) giving rise to the strong peaks below 0.5 eV. Simultaneously, high energy parts lose intensity. This effect is stronger for the *y*-component, whereas for the *z*-component, the initial pattern is still recognizable. This behavior is further increased for the full density (Fig. 5, green). For both components, some peaks appear in all three cases, especially in the energy range between 0.5 and 0.75 eV. They can be attributed to the coherent electron dynamics and may also be experimentally observable.

Further information can be gained by extracting the time when these frequencies occur. This allows us to connect them to a specific movement in the system. Therefore, we performed short-time Fourier transform spectra for the Δ 1e-2o-*y* and Δ 1e-2o-*z* components using a Gaussian windowing function with a width of 180 data points corresponding to a time of 18.14 fs. The resulting two spectrograms are shown in Fig. 6. The Δ 1e-2o-*y* spectrogram (left) shows two main pairs of signals around 10 fs (0.5–1.7 eV) and 50 fs (0.7–1.7 eV), which correspond to the first and the fourth passage of the wavepacket through the CoIn region. The signals are most pronounced at the first passage and significantly attenuated at the fourth passage. There are considerably weaker peaks observable at 25 and 30 fs, which can be attributed to the second and the third passage. In addition, the Δ 1e-2o-*z* spectrogram (right) shows two main signals. The first one appears around 10 fs (first passage through the CoIn) and covers a frequency range from 0.5 to 1.7 eV. The third passage around 30 fs can be attributed to the second signal, which extends over low-frequency components (0.1–1.0 eV) and has a lower intensity. Again, considerably weaker peaks can be found around 20–25 fs (second passage) and after 50 fs (fourth passage). Thus, each passage of the nuclear wavepacket through the CoIn region induces coherent electron dynamics, although not to the same extent for both components. The coherent dynamics is only short-lived for 5–7 fs, and the intensity of its signal decreases with time. The highest intensities are observed for the first transition when the localized initial nuclear wavepacket hits the CoIn. The subsequent dephasing and branching of the nuclear wavepacket blur the electronic coherence. In summary, we observe a short but recurring appearance of the coherent electron dynamics that is modulated by the nuclear wavepacket motion. In the following, we focus on the first passage (10 fs) for applying a few-cycle IR pulse to influence the coupled dynamics of NO_{2} since here the largest electronic coherence in the field-free case exists.

### B. Dynamics in the presence of a few-cycle IR pulse

Again, a delta pulse excitation is used to initiate the dynamics. With the appropriate time delay, a few-cycle IR laser pulse is applied to influence the first passage through the CoIn and thereby the subsequent coupled dynamics. The used few-cycle pulse has a Gaussian shape and is defined as

with

with the central frequency *ω*_{0}, the time zero *t*_{0}, the maximal field amplitude *E*_{max}, the full width half maximum (FWHM), and the carrier envelope phase *ϕ* (CEP). The time zero *t*_{0} of the pulse, defining the position of its maximum, was chosen to match the time window when the wavepacket is located near the CoIn (*t*_{0} = 10 fs). For this time, the nuclear wavepacket is still very localized and the electronic coherence is maximal. The central frequency *ω*_{0} is chosen to be resonant with the actual energy gap Δ*E* = 0.76 eV between the electronic states. The remaining three pulse parameters, the field amplitude *E*_{max}, the full width half maximum (FWHM), and the CEP *ϕ*, are set to *E*_{max} = 0.103 GV cm^{−1} (which corresponds to a maximum intensity of 1.4 × 10^{13} W cm^{−2}), FWHM = 8 fs, and *ϕ* = 0*π*. In comparison with the pulse parameters used by Richter *et al.*,^{57,59} all values are quite similar. Only our intensity is lower to stay in the range where the influence of the CEP pulse is mainly determined by the interplay of the non-adiabatic transition and the light induced electronic coherence.^{76} By this, we also ensure to stay below or at the threshold of ionization. The light–matter interaction is treated within the dipole approximation (for details, see Sec. I of the supplementary material). We assume that the electric component of the pulse is optimally aligned with the transition dipole moment. The absolute value of the TDM is used, which is shown in Fig. S3(a) of the supplementary material. As stated by Richter *et al.*,^{57} already a moderate molecular alignment distribution is sufficient to observe the effect of such a control pulse.

The evolution of the adiabatic populations influenced by the few-cycle IR-field is shown in the upper panel of Fig. 7. The related nuclear wavepacket dynamics on both surfaces integrated over the *α*-coordinate and the *b*-coordinate is depicted in Fig. S8 of the supplementary material. During the first transition through the CoIn region (7–15 fs), a 50:50 population of both states is created. The interaction with the light pulse is reflected in the small wriggles around 10 fs. The subsequent dynamics is comparable to the field-free case up to 30 fs. Thereafter, no clear passage through the CoIn region is observable. Thus, the IR pulse induces a change in the nuclear dynamics, which persists beyond the pulse duration. As an important consequence, the nuclear motion becomes asymmetric with respect to the *b*-coordinate, and the nuclear wavepacket even loses its nodal structure [compare both Figs. S5(b) and S8(b) of the supplementary material], which was also observed by Richter *et al.*^{57} This asymmetry leads to the partly deviations from the CoIn region after 30 fs. On the lower panel of Fig. 7, snapshots of the electron density in the 1e-2o picture are shown. Again, the difference in density with respect to *t* = 0 fs is depicted. The main features in the dynamics are quite similar to the field-free case. However, like for the nuclear motion, the dynamics of the electron density becomes asymmetric with respect to the *xy*-plane, i.e., the *b*-coordinate. This asymmetry persists after the laser pulse is no longer active (for example, see the snapshots at 30.0 fs). The oscillation of the electron density from the right to the left oxygen is most prominently observable for the snapshots at 7.6 and 9.6 fs.

The oscillations of the electron density are again recorded by the three dipole moment components, shown in the upper panel of Fig. 8. The coherent part of electron dynamics is visualized by the Δ 1e-2o dipole moment components for the *y*- and *z*-coordinates in the lower panel. Again, the 1e-2o-*x*-component stays zero for the whole simulation time. As the few-cycle IR pulse induces the asymmetry mainly along the *b*-coordinate, the overall temporal evolution of the 1e-2o-*y*-and the Δ 1e-2o-*y* components is similar to the field-free case. The 1e-2o-*z*-component experiences the main changes. During the pulse, strong and fast oscillations are observed with an amplitude nearly thirty times larger than for the field-free case. The oscillations stay up to ten times larger after the pulse. The superimposed slow oscillation with a period of about 20 fs can be assigned to the asymmetry in the nuclear motion. It does not appear for the Δ 1e-2o-*z* component reflecting solely the coherent electron dynamics. By breaking the symmetry of the nuclear motion with the laser pulse, the electronic coherence induced in the NO_{2} molecule is significantly larger. Again, it is observable mainly in the *z*-component and *b*-coordinate. During the light pulse, it is now the coherent electron dynamics that is responsible for the largest changes in the dipole moment.

The corresponding frequencies for the Δ 1e-2o-components, the 1e-2o-components, and the dipole moment calculated with the full density are again determined by Fourier transform. Their spectra are shown in Fig. 9. All frequencies with an intensity larger than 0.1 are listed in Tables S5 and S6 of the supplementary material. In both Δ 1e-2o spectra, frequencies up to 4.0 eV appear, which are higher compared to the field-free case. As expected, the main peaks of the Δ 1e-2o-*y* spectra [Fig. 9(a), blue] are in the same energy region as in the field-free case and only the Δ 1e-2o *z*-spectrum [Fig. 9(b), blue dotted line] shows differences. The main peaks are shifted to higher energies by roughly 0.7 eV. The laser pulse injects energy (0.76 eV) into the system, which influences the momentum of the nuclear wavepacket and thereby the phase of the overlap term [Eq. (1)], which subsequently leads to higher frequencies observed in the coherent electron dynamics. The correlation between the phase of the overlap term, the electronic phase, and the laser pulse is illustrated in Fig. S11 of the supplementary material for two individual grid points. The frequencies for the *y*-component determined with the 1e-2o-density [Fig. 9(a), red] and the full-density [Fig. 9(a), green] exhibit the same behavior as in the field-free case. The high energy parts significantly lose intensity since the slower nuclear dynamics (Born–Oppenheimer part) dominates this signal. The dominance of the oscillating dipole moment originating from the coherent electron dynamics shows up in the nearly identical spectra for 1e-2o-*z* [Fig. 9(b), red] and Δ 1e-2o-*z* [Fig. 9(b), blue]. For the *z*-spectra of the full-density [Fig. 9(b), green], the high energy parts lose some intensity but still more high energy contributions survive compared to the field-free case.

The results of the short-time Fourier transform for the Δ 1e-2o-*y* and the Δ 1e-2o-*z* dipole moment component using a Gaussian windowing function with a width of 180 data points corresponding to a time of 18.14 fs are shown in Fig. 10. Both spectrograms show a dominant signal that is attributed to the first passage through the CoIn region. The observable electron dynamics is significantly strengthened by the simultaneous light pulse interaction. In the case of the Δ 1e-2o-*y* spectrogram (left), some new features between 10 and 30 fs appear. Due to the symmetry breaking of the nuclear motion by the laser pulse, signals with very low frequencies and an extended signal around 1.0 eV appear. For the more affected Δ 1e-2o-*z* component, only one dominant peak is observed. In summary, the presence of a few-cycle IR pulse modifies the coupled dynamics by breaking the symmetry of the nuclear motion and changing the temporal evolution of the population. Both factors lead to a significant increase in electronic coherence in the molecule especially along the *z*-coordinate (laboratory frame) and the *b*-coordinate (internal frame).

## IV. WAVEFORM CONTROL OF MOLECULAR DYNAMICS

In the last part, we investigate the controllability of the nuclear and electron dynamics by the variation of the CEP *ϕ* of a few-cycle IR laser pulse. As shown in the literature,^{68,70–72,74–76} the CEP control scheme offers the possibility to steer electrons and nuclei not only in the ionization process but also during the passage through a CoIn. The few-cycle IR pulse builds up a coherent electronic and nuclear wavepacket with a well-defined phase-relationship controllable by the CEP. In the vicinity of a CoIn, also the non-trivial geometric phase (Pancharatnam–Berry phase) is introduced.^{45,46,77,78} The interplay of both phase-terms lead to an interference process when the CoIn is passed. The interference (constructive or destructive) can be manipulated by the CEP.

### A. Control of the nuclear dynamics

As a first step, we focus on the controllability of the nuclear dynamics. Therefore, we define control objectives that are directly accessible via the nuclear wavepacket and use the population *P*_{D0}(*t*, *ϕ*) of the *D*_{0} ground state as reference

One objective is the CEP efficiency Γ(*t*),^{76} which is calculated as the difference of the maximum and the minimum population *P*_{D0}(*t*, *ϕ*) for each time step,

For its maximum value, the population of the target state shows the highest CEP-dependence and consequently the highest degree of controllability with respect to the population transfer. The light pulse amplifies the coherent electron dynamics in the system by breaking the symmetry with respect to the asymmetric stretching coordinate *b*, as shown in Sec. III B. Therefore, the second objective is the CEP-dependent asymmetry parameter *AN*(*t*, *ϕ*), quantifying the CEP induced asymmetry in the nuclear motion with respect to the coordinate *b*,

where $PD0L(t,\varphi )$ and $PD0R(t,\varphi )$ are defined as follows:

In the spirit of the efficiency Γ(*t*), a maximal asymmetry *AN*_{max}(*t*) is calculated as

For its maximum, the motion of the nuclear wavepacket shows the highest asymmetry and controllability. Its CEP dependence is illustrated in Fig. 11.

The temporal evolution of Γ(*t*) and the CEP dependent population *P*_{D0}(*t*, *ϕ*) at three selected times are shown in Fig. 12. The CEP efficiency (blue line) reaches its global maximum (13%) nearly simultaneously with the peak intensity (*t*_{0} = 10 fs) of the laser pulse (gray area). The increase in Γ(*t*) is slightly delayed, and the subsequent decrease to 3% occurs in two steps. After the laser pulse, approximately at 15 fs, Γ(*t*) has a finite oscillating value with a maximum of about 5% around 20 fs, which indicates the second passage through the CoIn region. The later passages through the CoIn region at 30 fs and after 40 fs can roughly be seen in the increase in Γ(*t*). The deviation (violet curve) of the mean population (averaged over all CEPs) from the population in the field-free case is significant, especially during the IR pulse and after 30 fs. As discussed with respect to Fig. 7, the induced asymmetry leads to a partial missing of the CoIn region after 30 fs, which is almost independent of the CEP chosen. The CEP-dependence of the population *P*_{D0}(*t*, *ϕ*) [see Fig. 12(b)] is recorded for three selected times marked as vertical lines in 12(a). For better visualization, the mean difference is used here and, unless otherwise stated, in all following respective figures. The first line at 15 fs (green) matches the end of the laser pulse. The second (red line) and the third point (yellow line) correspond to the second and fourth passages through the CoIn region. For all three times, *P*_{D0}(*t*, *ϕ*) shows a sinusoidal oscillation with a periodicity of approximately *π*. For interference, a periodicity of 2*π* should emerge. Thus, the observed *π* dependence of the population is an indication that it is mostly due to the temporal asymmetry of the few-cycle laser pulse.^{69,76}

An analog analysis is performed for the asymmetry of the nuclear motion along the stretching coordinate *b* and shown in Fig. 13. The maximal asymmetry *AN*_{max}(*t*) shows its global maximum around 8 fs. As it is defined with respect to the population in *D*_{0} alone, the values for the early times (in the beginning of the laser pulse) are overestimated compared to the actual population in the *D*_{0} state. Nevertheless, we can deduce that *AN*_{max}(*t*) follows the envelope of the laser pulse. The subsequent peaks between 15 and 20 fs, at 30 fs, and between 42 and 48 fs correspond to the passages through the CoIn region. The decreasing height of the maxima reflects again the delocalization of the nuclear wavepacket with time. The CEP-dependence of the asymmetry of the nuclear motion [see Fig. 13(b)] *AN*(*t*, *ϕ*) is recorded for the same times as previously selected for the CEP-dependent populations *P*_{D0}(*t*, *ϕ*). It should be mentioned that the entire value of *AN*(*t*, *ϕ*) is shown here and not the mean difference. The asymmetry in the nuclear motion along the coordinate *b* shows a sinusoidal oscillation, now with a periodicity of 2*π* for all three times, which is typical for interference. This means that for the two quantities *P*_{D0}(*t*, *ϕ*) and *AN*(*t*, *ϕ*), we observe a different CEP-dependence. In other words, there are two different mechanisms active in the system, which can be projected out by using different observables.

In addition, we calculated the temporal evolution of Γ(*t*) and *AN*_{max}(*t*), as well as the CEP-dependence of *P*_{D0}(*t*, *ϕ*) and *AN*(*t*, *ϕ*) using the *y*-component and the *z*-component of the TDM. Since the results are quite similar, for the ones obtained with the absolute value of the TDM, the orientation of the molecule with respect to electric field of the pulse should not play a major role. For more details, see Sec. IV of the supplementary material.

### B. Control of the electron dynamics

As shown in Sec. III B, the laser pulse is creating a coherent electronic superposition in the vicinity of the CoIn. Therefore, we also examined the influence of the CEP variation on the electron density. The first control objective is the CEP-dependent asymmetry parameter *AE*(*t*, *ϕ*) of the 1e-2o-density *ρ*(*r*, *t*, *ϕ*),

with the probabilities *N*^{L}(*t*, *ϕ*) and *N*^{R}(*t*, *ϕ*) to find the electron on the left or the right side of the molecule given by

The maximal asymmetry of the electron density *AE*_{max}(*t*) is calculated as follows:

For its maximum, the electron dynamics shows the highest CEP-dependence and thus the highest controllability. The temporal evolution of *AE*_{max}(*t*) and the CEP-dependent asymmetry of the electron density *AE*(*t*, *ϕ*) at three selected times are shown in Fig. 14. The maximal asymmetry *AE*_{max}(*t*) is highest during the laser pulse (gray area). It decreases within 8 fs and becomes smaller by a factor of ten. However, during this time period, two peaks at 12 and 15 fs can be recognized. Afterward, the maximal asymmetry oscillates between nearly zero and 0.125 until the end of the simulation time. Comparing the maximal asymmetry of the electron density *AE*_{max}(*t*) with the one of the nuclei [*AN*_{max}(*t*)], faster oscillations are observed. To further analyze the response of the electron density [see Fig. 14(b)], *AE*(*t*, *ϕ*) is recorded for three selected points in time marked as vertical lines in 14(a). The first line at 10 fs (green) corresponds to the main peak of *AE*_{max}(*t*) and is taken at the maximum of the pulse. The second point (red line) is taken at 15 fs when the laser pulse is approximately over. The last point in time (yellow line) is at 40 fs. At all three times, *AE*(*t*, *ϕ*) shows a sinusoidal oscillation with a periodicity of approximately 2*π* and a decreasing amplitude with time. The asymmetry of the electron density thus has the same periodicity as the nuclear asymmetry *AN*(*t*, *ϕ*), which is as previously mentioned typical for an interference process.

As already discussed in Sec. III B, the response of the dipole moment to the applied laser field is an observable directly connected to the electron motion. In this case, the 1e-2o-*y*- and the 1e-2o-*z*-component are of interest. Their maximal CEP-dependence *γ*_{y}(*t*) and *γ*_{z}(*t*) are evaluated as the difference of the maximum and the minimum value of 1e-2o-*y*-DM(*t*, *ϕ*) and 1e-2o-*z*-DM(*t*, *ϕ*), respectively, for each time step. The maximal CEP-dependence *γ*_{y}(*t*) is depicted as a function of time in Fig. 15(a) and its related component 1e-2o-*y* in Fig. 15(b) at three selected times.

The maximal CEP-dependence *γ*_{y}(*t*) like all other objectives shows its maximum simultaneously with the maximum of the IR pulse. In this period, the shape of the *γ*_{y}(*t*) curve is similar to the Γ(*t*) curve [see Fig. 12(a)], only the decrease with decaying pulse intensity is even more asymmetric. After the pulse in the time window from 20 to 40 fs, the CEP-dependence oscillates. Again, the oscillations are significantly faster than for the nuclear objectives. The CEP-dependence of the 1e-2o-*y*-component is recorded in Fig. 15(b) for the same three selected times as for *AE*(*t*, *ϕ*). It shows a sinusoidal oscillation with a periodicity of approximately *π* and a decreasing amplitude with later times. Thus, the component shows the same periodicity as Γ(*t*) even with the same phase.

The temporal evolution of the maximal CEP-dependence *γ*_{z}(*t*) and its 1e-2o-*z*-component as a function of the CEP are shown in Fig. 16. The maximal CEP-dependence *γ*_{z}(*t*) is significantly larger than *γ*_{y}(*t*) in consistency with our finding in Sec. III B that the *z*-component reacts more strongly to the laser pulse. The overall shape of *γ*_{z}(*t*) is quite similar to the temporal evolution of *AE*_{max}(*t*) [see Fig. 14(a)], and the 1e-2o-*z*-component shows the same periodicity of 2*π* as *AE*(*t*, *ϕ*). The only difference is a phase shift of *π*.

In summary, two different responses on the CEP variation are present in the nuclear and electron dynamics. Both asymmetry parameters *AN*(*t*, *ϕ*) and *AE*(*t*, *ϕ*) as well as the 1e-2o-*z*-component of the dipole moment provide a distinction between left and right within the molecular plane (*yz*-plane). The associated 2*π* periodicity is typical for an interference process. Γ(*t*) and the 1e-2o-*y*-component of the dipole moment are directly sensitive to the main direction of motion along the *α*-coordinate and the *y*-coordinate. The motion in this direction mediates the non-adiabatic transfer between the *D*_{1} and *D*_{O} states. For these cases, the CEP-dependence shows a *π* periodicity, arising from the temporal asymmetry of the few-cycle pulse itself.^{69,76} Both mechanisms are present for the nuclear as well as for the electron dynamics and can be detected depending on the chosen observable.

## V. CONCLUSION

In this paper, we expand our ansatz for the description of the coupled nuclear and electron dynamics in molecular systems^{8,38,39} (NEMol). We applied our method to the photoinduced ultrafast dynamics in NO_{2}, which is dominated by a CoIn. We observe the appearance of a coherent electronic wavepacket at each passage of the CoIn. The coherence is not strong and only short lived due to the high symmetry of the molecule, which cancels out the individual contributions.^{79} Beside the field-free relaxation, we also studied the influence of a few-cycle IR laser pulse applied in the vicinity of the CoIn. The induced symmetry breaking significantly enhances the degree of coherence and its life time. Inspired by previous works,^{56,57,73,75,76} we varied the carrier envelope phase *ϕ* (CEP) of the IR pulse to control the movement of electrons and nuclei during the passage through the CoIn.

In the first part, we generalized our NEMol ansatz. The principle advantage of this ansatz is based on the combination of highly developed quantum-chemical methods with the accurate description of the nuclear quantum dynamics. In the original ansatz,^{8,38,39} an expression for the time-dependent electronic wavepacket is formulated, where the electronic part of the total wavefunction is propagated in the electronic eigenstate basis. Its dynamics is extracted from the nuclear wavepacket propagation on coupled potential energy surfaces by introducing the parametric dependence on the time-dependent expected value of position ⟨*R*⟩(*t*). By extending the NEMol ansatz with a grid representation, it is possible to couple the electron dynamics to multiple grid points on which the nuclear wavepacket is represented. Through a simple approximation, we were able to condense the coupled dynamics of the one-electron excitation process in the density of one active electron (1e-2o-picture). In the second part, we compared the coupled nuclear and electron dynamics of NO_{2} with and without an IR pulse present when the system reaches the CoIn for the first time. Using the NEMol ansatz, we characterized the coherent electron dynamics by analyzing the temporal evolution of the induced dipole moment. The observed frequencies of the coherent electron dynamics cover a range up to 2.3 eV. These high values originate from the nuclear overlap term and the electronic phase term. In NO_{2}, the phase contribution of the nuclear overlap term is high and therefore provides a significant contribution to the electron dynamics. The applied few-cycle IR laser pulse generated an asymmetric movement of the nuclear and electronic wavepackets, which is vital for the controllability at the CoIn. The induced oscillating dipole reflects an enhanced buildup of the coherent electron dynamics by the laser pulse, which survives for several 10 fs. In the last part, the CEP of the IR pulse was varied to influence both the nuclear dynamics and the electron dynamics. The CEP-dependent effect lives considerably longer than the pulse in all investigated observables. Depending on the chosen observable, a *π* or 2*π* periodicity can be found indicating two mechanisms, one based on an interference process (2*π*) and the other induced by the temporal asymmetry of the few-cycle pulse itself (*π*). Both periodicities are observed for the nuclear and the electron dynamics. In each case, they can be projected out by using different observables.

We demonstrated the potential of our NEMol ansatz to describe the coupled nuclear and electron dynamics in molecular systems beyond diatomics. In NO_{2}, we followed the dynamics in the excited state dominated by fast changing wavepacket interference effects. The ansatz is expandable to simulate the induced coherent electron dynamics in the excitation process itself as well as the higher-dimensional molecular system as long as the underlying nuclear dynamics can be treated quantum mechanically. Two electron processes could be realized by using pair densities.

## SUPPLEMENTARY MATERIAL

See the supplementary material for the details of the wavepacket simulation setup, the underlying quantum chemical data of NO_{2}, and additional figures and tables for the NEMol-dynamics. A section contains the results for the CEP-control obtained with the *y*-component and the *z*-component of the TDM. Animations of the coupled electron density in the 1e-2o picture for the free propagation and in the presence of a few-cycle IR laser pulse with a CEP of 0.0*π* and 1.4*π* are also added in the supplementary material.

## AUTHORS’ CONTRIBUTIONS

T.S. performed all calculations. T.S. and R.d.V.-R. analyzed the results and contributed equally to the final version of the manuscript.

## ACKNOWLEDGMENTS

The authors gratefully acknowledge the DFG Normalverfahren and the Munich Center of Advanced Photonics (MAP).

There are no conflicts to declare.

## DATA AVAILABILITY

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

### APPENDIX: THE ORIGINAL NEMol ANSATZ

The following detailed formulation of the NEMol ansatz^{8,38,39} is given here in the improved notation. The total molecular wavefunction Ψ_{tot}(*r*, *R*, *t*) is setup as the sum over the electronic states, with *χ*(*R*, *t*) being the nuclear wavefunctions, *φ*(*r*, *t*; *R*) being the electronic wavefunctions, *R* and *r* being the nuclear and electronic coordinates, respectively, and *t* being the time,

Applying the Born–Oppenheimer approximation, the uncoupled electronic wavefunctions *φ*_{i} are hereby parametrically depending on the nuclear coordinates *R* and define a multi-dimensional vector *φ*_{tot}. The total nuclear wavefunction *χ*_{tot} also represents a multi-dimensional vector, spanned by the coupled wavefunctions *χ*_{i}. For details on how the temporal evolution of the nuclear wavefunctions *χ*_{i} on coupled potential energy surfaces (PES’s) is determined, see Sec. I of the supplementary material. Multiplying Ψ_{tot}(*r*, *R*, *t*) from the left with *χ*_{tot} and subsequently integrating over the nuclear coordinates results in an expression of the coupled total electronic wavefunction,^{8,38,39}

with

The coupled total electronic wavefunction is parametrically depending on the time-dependent expected value of the position ⟨*R*⟩(*t*). In other words, Φ_{tot} is evaluated at one single nuclear geometry, which changes with time. The individual components Φ_{j} are defined by the following equation:

with

The first part depends on the population *A*_{jj} of the respective state *j*, while all others summands include the nuclear overlap term *A*_{jk}, which specifies the degree of coherence induced between the two states *j* and *k*. The population and coherence of the electronic states as well as the influence of all coupling terms are already determined by the nuclear quantum-dynamics simulation. If the coupling between the electronic states is weak, the nuclear wavefunctions propagate independently and the coherence term becomes zero. In this case, the coupled electronic wavefunctions Φ_{j} in Eq. (A4) become equivalent to the uncoupled electronic wavefunction *φ*_{j}. Standard quantum-chemical calculations at the ⟨*R*⟩(*t*) structure yield the real-valued wavefunctions *φ*_{j}(*r*; ⟨*R*⟩(*t*)) of the relevant electronic states and their eigenenergies. The temporal evolution of *φ*_{j}(*r*, *t*; ⟨*R*⟩(*t*)) is determined by the deformation of the electronic structure induced by the nuclear motion (Born–Oppenheimer part) and an oscillation through phase space defined by a pure electronic phase,^{8,38,39}

The phase term *ξ*_{j}(*t*) depends on the eigenenergies *E*_{j}(⟨*R*⟩(*t*)) and has to be calculated recursively,

This recursive evaluation is necessary to retain the memory of the progressing electronic phase. Thereby, the propagation velocity of the phase in the complex plane changes smoothly in time while the nuclear wavepacket propagates. Using the coupled total electronic wavefunction Φ_{tot}(*r*, *t*; ⟨*R*⟩(*t*)), the associated electron density *ρ*(*r*, *t*; ⟨*R*⟩(*t*)) can be determined by multiplying Φ_{tot}(*r*, *t*; ⟨*R*⟩(*t*)) from the left with *φ*_{tot} and subsequently integrating over *N* − 1 electronic coordinates (with *N* being the total number of electrons),

with

The first summation consists of the state specific electronic density *ρ*_{jj}(*r*, *t*; ⟨*R*⟩(*t*)) weighted with the corresponding time-dependent population *A*_{jj}(*t*). The dynamics of these contributions to the coupled electron density is determined by the temporal evolution of the nuclear wavepacket, i.e., its expected value of the position ⟨*R*⟩(*t*). The second summation defines the coherent contribution to the coupled electron density and consists of the time-dependent overlap *A*_{jk}(*t*), the one-electron transition density *ρ*_{jk}(*r*, *t*; ⟨*R*⟩(*t*)), and its pure electronic phase defined by the energy difference Δ*E*_{jk} between the involved electronic states.

## REFERENCES

_{2}and CO

^{+}

_{2}: A three-mode nuclear dynamics investigation

^{2}A

_{1}and $A\u0303$

^{2}B

_{2}states of NO

_{2}on the negative ion photoelectron spectra of NO

_{2}

^{−}

_{2}on the $X\u0303$

^{2}A′/A

^{2}A′ conical intersection

^{2}A

_{1}/A

^{2}B

_{2}electronic states of NO

_{2}based on new ab initio potential energy surfaces

^{2}A′/A

^{2}A′ conical intersection

_{2}: Global potential energy surfaces of the ground (1

^{2}A

_{1}) and the first excited (1

^{2}B

_{2}) electronic states

_{2}close to the bottom of the X

^{2}A

_{1}-A

^{2}B

_{2}conical intersection

_{2}conical intersection

_{2}in the second absorption band: Ab initio and quantum dynamics calculations

_{2}

_{2}after excitation to the A-band

_{2}from time-resolved photoelectron emission

_{2}studied by time-resolved imaging

_{2}

_{2}

_{2}

_{2}probed by homodyne high-harmonic spectroscopy

_{2}electronic relaxation

_{2}