Electronic decoherence following photoionization: full quantum-dynamical treatment of the influence of nuclear motion

TL;DR

This paper presents a comprehensive quantum-dynamical approach to understand how nuclear motion causes electronic decoherence after photoionization, revealing that slow vibrational modes involving the molecular skeleton rapidly destroy electronic coherence within femtoseconds.

Contribution

It introduces a full quantum-dynamical adiabatic framework to model nuclear motion and its impact on electronic decoherence post-photoionization, considering multiple nuclear degrees of freedom.

Findings

Electronic decoherence occurs within a few femtoseconds.

Slow vibrational modes involving the molecular skeleton rapidly destroy coherence.

Fast hydrogen vibrations do not significantly affect short-time electronic coherence.

Abstract

Photoionization using attosecond pulses can lead to the formation of coherent superpositions of the electronic states of the parent ion. However, ultrafast electron ejection triggers not only electronic but also nuclear dynamics---leading to electronic decoherence, which is typically neglected on time scales up to tens of femtoseconds. We propose a full quantum-dynamical treatment of nuclear motion in an adiabatic framework, where nuclear wavepackets move on adiabatic potential energy surfaces expanded up to second order at the Franck-Condon point. We show that electronic decoherence is caused by the interplay of a large number of nuclear degrees of freedom and by the relative topology of the potential energy surfaces. Application to $\mathrm{H_2O}$, paraxylene, and phenylalanine shows that an initially coherent state evolves to an electronically mixed state within just a few femtoseconds. In these examples the fast vibrations involving hydrogen atoms do not affect electronic coherence at short times. Conversely, vibrational modes involving the whole molecular skeleton, which are slow in the ground electronic state, quickly destroy it upon photoionization.

Equations

$$H^{(\mu )}(Q_1,\cdots ,Q_f)=T+\Delta E^{(\mu )}+\sum _{i=1}^f{\kappa }_i^{(\mu )}{Q}_i+\frac{1}{2}\sum _{i,j=1}^f{\gamma }_{ij}^{(\mu )}{Q}_i{Q}_j.$$

$$\Psi (r,Q,t)=\sum _{\mu =1}^n{c}_{\mu }{\chi }_{\mu }(Q,t){\varphi }_{\mu }(r;Q),$$

$$⟨\hat{O}(r)⟩$$

$$=\sum _{\mu \nu }\int \mathrm{d}Q{c}_{\mu }^{\ast }{\chi }_{\mu }^{\ast }(Q,t)O_{\mu \nu }(Q)c_{\nu }{\chi }_{\nu }(Q,t)$$

$$\approx \sum _{\mu \nu }{O}_{\mu \nu }{\rho }_{\nu \mu },$$

$${\rho }_{\nu \mu }(t)=c_{\nu }{c}_{\mu }^{\ast }\int \mathrm{d}Q{\chi }_0^{\ast }(Q,0){\mathrm{e}}^{i{H}^{(\mu )}t}{\mathrm{e}}^{-i{H}^{(\nu )}t}{\chi }_0(Q,0).$$

$${\tilde{\rho }}_{\nu \mu }=c_{\nu }^{\ast }{c}_{\mu }{\mathrm{e}}^{i\left(\Delta E^{(\mu )}-\Delta E^{(\nu )}\right)t}.$$

$${\Psi }_0=\prod _{i=1}^f{\chi }_{0,i}|0⟩,{\chi }_{0,i}={\mathcal{N}}_i{\mathrm{e}}^{-{\left(Q_i-Q_{0,i}\right)}^2/2}$$

$${\rho }_{\nu \mu }(t)=c_{\nu }{c}_{\mu }^{\ast }\prod _{i=1}^f\int \mathrm{d}{Q}_i{\chi }_{0,i}^{\ast }(Q_i,0){\mathrm{e}}^{i{H}_i^{(\mu )}t}{\mathrm{e}}^{-i{H}_i^{(\nu )}t}{\chi }_{0,i}(Q_i,0)=c_{\nu }{c}_{\mu }^{\ast }\prod _{i=1}^f{\rho }_{\nu \mu }^{(i)}(t).$$

$$\Psi (\mathbf{Q},0)={\chi }_{0,1}(Q_1,0){\chi }_{0,2}(Q_2,0){\chi }_{0,3}(Q_3,0)\frac{1}{\sqrt{4}}\left(|1⟩+|2⟩+|3⟩+|4⟩\right)$$

$$\Psi (\mathbf{Q},0)=\prod _{i=1}^{48}{\chi }_{0,i}(Q_i,0)\frac{1}{\sqrt{2}}\left(|1⟩+|2⟩\right).$$

Full text

I Introduction

The advent of attosecond pulses allows studying electronic correlation and ultrafast molecular dynamics through pump-probe experiments with unprecedented time resolution Hentschel et al. (2001); Rozzi et al. (2013); Calegari et al. (2014). Attosecond, broad-band pulses can be used to generate coherent superpositions of cationic states Calegari et al. (2014); Pabst et al. (2011); Timmers et al. (2014); Kraus et al. (2013). The ionization triggers electronic and nuclear dynamics, finally leading to electronic decoherence. Decoherence is defined as the process where an initially pure state evolves to a statistical ensemble Zhu et al. (2004); Vacher et al. (2014). Electronic coherences are assumed to be responsible for a number of chemical processes, for example the high quantum efficiency of the energy conversion in photosynthesis, which is a matter of intensive current debate Romero et al. (2014).

Long-lived coherences are predicted by theories that focus on the evolution of the electronic subsystem, driven by electronic correlation Kuleff and Cederbaum (2007); Golubev and Kuleff (2015). Nuclear motion is neglected, because the electrons move much faster than the heavier nuclei. This results in charge migration, an oscillatory motion of hole and electron density with frequencies defined by the energy gaps among the cationic states populated in the ionization process. In Ref. Calegari et al. (2014), ultrafast dynamics in polyatomic molecules on a time scale shorter than the vibrational response were attributed to charge migration.

In recent theoretical work, the quantum nature of the nuclei was approximately taken into account Vacher et al. (2015a, b); Jenkins et al. (2016a, b). The authors sampled nuclear geometries within the width of the Gaussian wave packet of the nuclear ground state. This leads to a superposition of coherent oscillations with different frequencies that average out within a few femtoseconds. The cancellation of the oscillations is due to the energy gap between the cationic states at the respective nuclear geometries. This approach takes into account the spatial delocalization of the nuclear wave packet. However, it does not consider the time evolution of the nuclear wave packet on the different potential energy surfaces. In the context of high-order harmonic generation (HHG) in molecules, the propagation of vibrational wave packets on adiabatic potential energy surfaces was studied in Ref. Patchkovskii (2009). It was pointed out that the spatial overlap of vibrational wave packets in different electronic states depends on the topology of the potential energy surfaces.

In this paper, we present a model for an ab-initio, full quantum-mechanical treatment of nuclear motion on short time scales. Nuclei move on adiabatic potential energy surfaces that are expanded up to second order around the Franck-Condon point, i.e., the equilibrium geometry of the ground state. Bilinear mode-mode couplings are included. Due to the different values of the mode-mode couplings on each potential energy surface, they cannot be all removed simultaneously by rotation to another set of normal modes, and the model must be solved numerically. Using the multi-configuration time-dependent Hartree method (MCTDH) Worth et al. ; Beck et al. (2000); Meyer et al. (1990) for the wave packet propagation, this allows us to numerically study the electronic coherence in large molecules considering all internal degrees of freedom at short times, as long as the Taylor expansion of the potential energy surfaces remains valid. We observe electronic decoherence on a femtosecond time scale that can be attributed to the interplay of various vibrational modes in a molecule. Contrary to the widely held belief that fast modes are responsible for decoherence at short times Lépine et al. (2014), it is the small displacement of many slow modes that is ultimately responsible for the loss of coherence. In Sec. II we introduce the model and argue that, generally, electronic decoherence is a multi-modal effect. This is illustrated in Sec. III by application to $\mathrm{H_{2}O}$, paraxylene Vacher et al. (2015b), and phenylalanine Calegari et al. (2014).

II Model for the dynamics following photoionization

We study the dynamics following photoionization in an electronically adiabatic framework. In the Born-Oppenheimer approximation, the nuclei move on adiabatic potential energy surfaces defined by the spectrum of the electronic Hamiltonian at the respective nuclear geometries. We consider the situation where a coherent attosecond pulse is used to excite a molecule above the ionization threshold. As the pulse is wide in the frequency domain, it is likely that two or more cationic states are populated in a coherent way Pabst et al. (2011). We assume vertical excitation, where the ground state nuclear wave packet is placed on the potential energy surfaces of the cationic electronic states. Furthermore, we neglect the photoelectron, so that the cation is initially in a pure electronic state.

Generally speaking, the equilibrium geometry of the ground state does not correspond to an equilibrium geometry of the excited states. After photoionization, the system enters a nonstationary state, and electronic and nuclear dynamics set in. In the adiabatic approximation, the nuclear wave packets move independently of each other on their respective potential energy surfaces.

II.1 Adiabatic Hamiltonian

Finding the shape of the adiabatic potential energy surfaces is the bottleneck towards establishing an exact solution of the problem outlined in the previous paragraph. The computational effort for determining the potential energy surfaces grows exponentially with the number of degrees of freedom. For small times, the potential energy surfaces can be approximated by Taylor polynomials up to second order, as proposed in the vibronic-coupling Hamiltonian Worth and Cederbaum (2004). This corresponds to vibrations of the nuclei around their equilibrium configuration Patchkovskii (2009). In this work, we neglect nonadiabatic couplings, which will play a role only if the cationic electronic states feature conical intersections or avoided crossings close to the Franck-Condon region Timmers et al. (2014). The shorter and thus spectrally broader the ionization pulse, the less important nonadiabatic effects should become. The Hamiltonian of the cationic electronic state $\mu$ in a molecule with $f$ nuclear degrees of freedom is then given by

[TABLE]

We use atomic units and mass- and frequency-weighted normal-mode coordinates $Q_{i}$ throughout this paper. $T=\sum\limits_{i=1}^{f}-\frac{\omega_{i}}{2}\partial_{Q_{i}}^{2}$ refers to the kinetic energy, where the ground-state normal-mode frequencies $\omega_{i}$ are employed. The energy difference to the lowest cationic state is denoted by $\Delta E^{(\mu)}$. The coefficients $\kappa_{i}^{(\mu)},\gamma_{ii}^{(\mu)},\gamma_{ij}^{(\mu)}$ are obtained by central-difference approximation for small displacements from the Franck-Condon point along the normal modes. They describe the gradient, curvature, and mode-mode coupling, respectively.

For propagating the nuclear wave packet, we use the multi-configuration time-dependent Hartree method (MCTDH) Meyer et al. (1990) in its multi-set formulation Fang and Guo (1994); Worth et al. (1996). The nuclear wave packet in each electronic state is described by an independent product of time-dependent single particle functions, thus optimally treating the independent nuclear evolution on the different cationic potential energy surfaces.

II.2 Electronic Decoherence

To quantify electronic decoherence, the coherence among the $n$ cationic states is calculated from the reduced density matrix $\rho$ of the electronic subsystem. Starting from the Born-Huang ansatz Born and Huang (1954), the decomposition of the total wave function in terms of nuclear and adiabatic electronic wavefunctions $\chi,\phi$, respectively, reads

[TABLE]

where $r$ refers to electronic coordinates. The expectation value of an electronic observable $\hat{O}(r)$ can be obtained via

[TABLE]

where for the last line we make the assumption $O_{\mu\nu}(Q)=O_{\mu\nu}$. The electronic reduced density matrix elements in Eq. (5) correspond to the nuclear wave packet overlaps:

[TABLE]

In the single-geometry model Vacher et al. (2015b); Jenkins et al. (2016b), the nuclear functions reduce to delta peaks, $\chi_{0}(Q)=c_{0}\delta(Q-Q_{0})$. Neglecting nuclear motion, Eq. (6) reduces to the long-lived electronic coherences predicted in charge-migration theories Calegari et al. (2014); Raab et al. (1999); Breidbach and Cederbaum (2003); Golubev and Kuleff (2015)

[TABLE]

Evolution of this single geometry along a trajectory as a consequence of an averaged force corresponds to Ehrenfest dynamics, also resulting in a perpetually pure electronic subsystem without decoherence. A full quantum-mechanical treatment of the nuclei, however, reveals that electronic decoherence is caused by the fast spread of nuclear wave packets along all degrees of freedom. The associated speed is determined by the relative topology of the potential energy surfaces, as shown in Fig. 1. Consider an initial state formed by a product of Gaussians and the electronic ground state $\ket{0}$

[TABLE]

with appropriate prefactors $\mathcal{N}_{i}$ ensuring normalization. This ground-state wave packet is placed on the cationic surfaces following the ionization by a coherent pulse. Neglecting the mode-mode couplings $\gamma_{ij}^{(\mu)},i\neq j$, the Hamiltonian allows one to factorize the wave packet along the coordinates $Q_{i}$ for all times, where $H^{(\mu)}=\sum_{i=1}^{f}H_{i}^{(\mu)}$. From Eq. (6),

[TABLE]

The oscillating $\rho_{\nu\mu}^{(i)}(t)$ will dephase, because, in general, the frequencies defined by the gaps among the cationic electronic states are $Q$-dependent. Essentially, coherence is lost because the product of a large number of factors ranging between 0 and 1 tends to 0. The process is speeded up if at least one mode contributes a factor close to 0. This corresponds to a vanishing spatial overlap caused by differing gradients in the potential energy surfaces. The frequency of the mode itself in the ground state does not determine its influence on decoherence.

We calculate the purity $\mathrm{Tr}(\rho^{2})$ to monitor the evolution of an initially pure state into a mixed state. The former yields $\mathrm{Tr}(\rho^{2})=1$, the latter, in an $n$-state system with equal weights $\frac{1}{n}$ for each state, $\mathrm{Tr}(\rho^{2})=1/n$. This approach is convenient especially if more than two cationic states are considered; furthermore, it is representation independent. Note that if the decay of the off-diagonal matrix element can be expressed as $\rho_{\mu\nu}=c\mathrm{e}^{-\gamma t}$, the purity decays twice as fast; for an equally weighted two-level system, $\mathrm{Tr}(\rho^{2})=\frac{1}{2}+2|c|^{2}\mathrm{e}^{-2\gamma t}$.

III Results and discussion

In the following, we apply our model to $\mathrm{H_{2}O}$, paraxylene, and phenylalanine cations. For $\mathrm{H_{2}O}$ and paraxylene, the potential energy surfaces were calculated at the multi-configurational self-consistent field (MCSCF) level of theory restricting the active orbitals of the ion to occupied orbitals in a Hartree-Fock calculation of the neutral system and using the GAMESS software package Schmidt et al. (1993). This level of accuracy was found to be sufficient in previous work Li et al. (2013a, 2014). For phenylalanine, we approximate the potential energy surfaces of the cation using a single configuration, using Koopmans’ theorem. We compare one-dimensional simulations, where only one normal mode at a time is considered, to full-dimensional simulations. By this, we demonstrate that decoherence is an effect that can be attributed to the interplay of a large number of vibrational modes, rather than to a small number of fast modes.

III.1 $\mathrm{H_{2}O}$

$\mathrm{H_{2}O}$ was chosen as a model molecule for its simplicity. The coefficients for the adiabatic potential energy surfaces were calculated for four cationic states using seven electrons in four active orbitals. The energy spacings with respect to the cation ground state are 0.06 eV, 0.27 eV, 0.87 eV, respectively.

Figure 2 shows the evolution of the coherence in the $\mathrm{H_{2}O}$ molecule. The results for one-dimensional simulations with one normal mode each are shown together with the simulation of the full three-dimensional system.

We observe that the initially pure state

[TABLE]

evolves towards a mixture in the electronic subsystem. A constant value is reached already after 1 fs. The individual modes are identified with the asymmetric stretch ($3756\,\mathrm{cm^{-1}}$), the symmetric stretch ($3657\,\mathrm{cm^{-1}}$) and the bending motion ($1595\,\mathrm{cm^{-1}}$). While the fast modes alone maintain electronic coherence, the slow mode is mainly responsible for the overall loss of coherence. The rate of the decoherence, and the suppression of recurrence, is attributed to the dephasing of the oscillations along single modes. Even in $\mathrm{H_{2}O}$ with only three degrees of freedom, the loss of coherence is thus seen to be a multi-modal effect with a strong participation of the slower bending motion.

The full quantum-mechanical treatment of nuclear motion reveals its high influence on electronic decoherence. In Ref. Li et al. (2013b), long-lived oscillations were observed in ionized liquid water and were attributed to coherent hole dynamics including the two lowest-lying cationic states. While the chemical environment of a single water molecule does not directly translate to liquid water or water clusters, the loss of coherence can be expected to occur even faster if more modes are included. Restriction to two states in $\mathrm{H_{2}O}$ does not change the decoherence time scale.

III.2 Paraxylene

We now apply the model to paraxylene, a benzene molecule with two methyl groups on opposite sides of the benzene ring. Paraxylene was included in the aforementioned study Vacher et al. (2015b), which makes it suitable to compare the conceptually different approaches. With 18 atoms, paraxylene has 48 internal degrees of freedom. We chose to project out the two methyl group rotations, which otherwise should have been treated as hindered rotors with a periodicity of $60^{\circ}$ Vendrell et al. (2002). Two additional very low frequency modes were left out of the model to avoid numerical problems. Leaving modes out of the model can only slow down decoherence, not accelerate it, thus leaving our conclusions unchanged.

Following Ref. Vacher et al. (2015b), we include the first two cationic states of paraxylene, separated by $0.51\,\mathrm{eV}$ at the Franck-Condon point. The cationic potential energy surfaces were determined at the MCSCF level of theory using five electrons in three active orbitals.

Initially, the system is prepared in an equally weighted superposition of the two cationic states,

[TABLE]

The high number of degrees of freedom restricts us to one time-dependent single particle function (SPF) per degree of freedom and state, which is however fully flexible and expanded in a harmonic DVR grid Beck et al. (2000). For selected states, more SPF were included to estimate their contributions.

Figure 3 shows the evolution of the initially pure state to a mixed state on the timescale of $2-3\,\mathrm{fs}$. The one-dimensional simulations, represented by the blue dashed curves, show that coherence is maintained at the single-mode level for most of the modes. Recurrence is observed after the associated oscillation period. The overall coherence, as outlined in Sec. II, can be described as the product of the one-dimensional coherences, if mode-mode couplings are neglected. Due to the dephasing of the oscillations along the different normal modes, the overall coherence drops to zero and does not recur. This shows that the loss of coherence can be explained from the high number of degrees of freedom. There is no single mode that is responsible alone for the fast decoherence. The breathing mode $(1330\,\mathrm{cm^{-1}})$ and the two Kekule modes $(1777\,\mathrm{cm^{-1}},1831\,\mathrm{cm^{-1}})$ show recurrence of the purity on a femtosecond timescale, if treated individually.

In Fig. 4, we study the influence of various selected modes on the decoherence. The ten fastest normal modes correspond to C-H vibrations. Electronic coherence is maintained if only this subset of modes is considered. The 34 slowest modes, with their relatively long recurrence times, lead instead to rapid loss of coherence. It is thus not the speed of the nuclear motion that defines the influence on electronic decoherence, but rather the topology of the potential energy surfaces relative to each other, as illustrated in Fig. 1. Note that, within our adiabatic model, the decoherence rate is independent of the weights assigned to the initial states, as shown in Fig. 5.

In a typical MCTDH calculation, the results are taken to be converged if the population of the highest SPF is below $10^{-3}$, or if the effect of leaving the highest SPF out is below $10^{-5}$ (see Ref. Beck et al. (2000)). This level of accuracy is certainly not feasible with a large number of degrees of freedom. However, we can approximate by the following reasoning: There are many pairs of vibrational modes without discernible mode-mode coupling. Along these modes, the nuclear wave packet is very well described using only one SPF. More SPFs are needed to capture the dynamics with respect to coupled modes. The one-dimensional calculations are of course exact with one SPF per surface.

Figure 6 shows the effect of including a second SPF for a subset of twelve coupled modes $(|\gamma_{ij}^{(\mu)}|>10^{-3})$. The results are comparable to the ones obtained by using only one SPF per mode. It is noticeable that the recurrence at $8\,\mathrm{fs}$ is suppressed if two SPF are used in this group of modes and that decoherence becomes slightly faster.

In Ref. Vacher et al. (2015b), electron dynamics following photoionization was studied in paraxylene. Decoherence was attributed to the width of the nuclear wave packet, because the energy gap and hence the frequency of $\rho_{12}$ changes over its extension. After averaging over all nuclear coordinates, coherence is lost without recurrence, in the case of paraxylene with a half-life time of $t_{1/2}=4\,\mathrm{fs}$. This approach does not allow for transfer of density between the different nuclear geometries within the wave packet and decoherence becomes a consequence of dephasing by averaging the coherent matrix elements over many different geometries. It follows from Eq. (6) if the kinetic energy operator is neglected and the nuclear wavepacket is replaced by a geometry distribution $p(Q)$. Note that non-adiabatic effects are thereby excluded. Although this model provides a physically reasonable time scale for decoherence in multidimensional systems, it is unclear when and how the approximation breaks down and how to systematically improve it. In contrast, in the model proposed here, the whole nuclear wave packet evolves fully quantum-mechanically and the treatment of decoherence, and potentially of time-resolved electronic spetroscopies, becomes exact at short times. By comparing, within our model, the decoherence rate with and without including kinetic energy, as shown in Fig. 7, we find that neglecting the kinetic energy operator in the quantum-dynamical calculation slightly overestimates electronic decoherence and misses the persisting linear component seen in the long-term behavior.

III.3 Phenylalanine

In Fig. 8, we show the evolution of electronic coherence in phenylalanine, which was studied in Ref. Calegari et al. (2014) for an initial state arising from the photoionization cross-sections. In our analysis, initially, the cation ground state and the first excited state are in a coherent superposition. Again, we use one SPF per surface and mode, and we compare one-dimensional simulations to the propagation including all 63 degrees of freedom. Most mode-mode couplings have a value $|\gamma_{ij}^{(\mu)}|<10^{-3}$. We find that the initially pure state evolves to a mixed state within about $1\,\mathrm{fs}$, if all modes are considered. Recurrence is suppressed, as the oscillations dephase along the different modes. While most modes preserve electronic coherence over the period of investigation, we observe a stronger participation of relatively slow vibrational modes with ground-state frequencies of $1816\,\mathrm{cm^{-1}},1809\,\mathrm{cm^{-1}}$ , and $1340\,\mathrm{cm^{-1}}$, respectively.

IV Conclusion and Outlook

It is often assumed that, if nuclear motion plays a role in electronic decoherence, it is the fast modes that should be considered Calegari et al. (2014). In an organic molecule, these are the C-H vibrations with a typical frequency of about $3000\,\mathrm{cm^{-1}}$, corresponding to a period of vibration of $11\,\mathrm{fs}$. This is typically slower than calculated charge migration times and is therefore not considered to have an effect. Our study, however, indicates that it is not the fast, but the interplay of the slow modes (in the sense of their inverse vibrational frequencies in the neutral ground state) that causes decoherence on a femtosecond timescale. In Fig. 4, we show that, in paraxylene, if only the 10 C-H vibrations are considered, the system stays coherent within at least $6\,\mathrm{fs}$. This is because the potential energy surfaces for the fast modes turn out to be only vertically displaced from each other. As shown schematically in Fig. 1 and discussed in Sec. II.2, in this case, the nuclear wave packets move synchronously on the potential energy surfaces, and their spatial overlap does not change much. In the slower modes, the potential energy surfaces can also be horizontally displaced from each other, and the spatial overlap is reduced.

We presented an ab-initio model for electronic decoherence following photoionization that takes the quantum nature of the nuclei into account and allows for a full-dimensional treatment of the molecule. Combined with a proper treatment of the preparation step and with the consideration of the relevant electronic observables it can be used to interpret and predict the outcome of current experiments in molecular attoscience Calegari et al. (2014). The probe step can be excluded from the analysis, as the loss of electronic coherence would be apparent in any probe technique sensitive to electronic structure. We showed that electronic decoherence can be explained by considering the topologies of different potential energy surfaces. It is the interplay of a large number of vibrational modes that are causing decoherence in the electronic density matrix, not a set of e.g. fast C-H vibrations. Our results suggest that in molecular systems, purely electronic dynamics that may be described in terms of a coherent electronic wave packet, exists only for sub-femtosecond time scales, and nuclear motion cannot be neglected. With the approach shown in this paper, one can calculate the time-dependent density matrix of the electronic subsystem and from it the time-dependent expectation values of any observable that depends only on the electronic subsystem via $\langle A\rangle=\mathrm{Tr}(A\rho),$ for example the hole density. Observables of the form $O_{\mu\nu}(Q)$ that depend on the electronic states and the nuclear geometry can be easily computed as well since the full nuclear wave packets are available.

The model is currently limited to short time scales. In the future, we plan to explore quantum-classical approaches that allow us to propagate the nuclei for longer times. However, we note that for describing the short-time electronic response with nuclear-induced decoherence, this model (with possible non-adiabatic extensions) is all what is required, as decoherence sets in before the nuclear wave packets abandon the quadratically expanded region of the potential energy surfaces. The exact quantum-mechanical calculations will serve as a reference for future development of these methods.

Acknowledgements.

This work has been supported by the excellence cluster The Hamburg Centre for Ultrafast Imaging - Structure, Dynamics and Control of Matter at the Atomic Scale of the Deutsche Forschungsgemeinschaft.