## 2 Methods  
### Wavepacket Dynamics  
Ultrafast experiments follow the evolution of a time-dependent wavepacket, \(|\Psi(\mathbf{r},\mathbf{R},t)\rangle\), where \(\mathbf{r}\) and \(\mathbf{R}\) are the electronic and nuclear coordinates, and \(t\) is time. The wavepacket is the solution to the time-dependent Schrodinger equation,  
\[t\hbar\frac{d}{dt}|\Psi(\mathbf{r},\mathbf{R},t)\rangle=\hat{H}|\Psi(\mathbf{ r},\mathbf{R},t)\rangle, \tag{1}\]  
where the Hamiltonian \(\hat{H}\) can be decomposed into a field-free molecular Hamiltonian \(\hat{H}_{0}\) and an interaction term \(\hat{H}_{\text{int}}\) that describes the interaction with the pump pulse (and, later, the probe pulse). Most commonly, the pump will be an optical laser pulse, so that the light-matter interaction is given in the dipole approximation by \(\hat{H}_{\text{int}}(t)\approx-de(t)\), where \(\epsilon(t)\) is the time-dependent electromagnetic field and \(d\) the electric dipole operator. The wavepacket excited by the pump pulse will have certain qualitites imprinted by the pump laser, including a bandwidth, phase-properties, and overall energy distribution. There are many strategies for solving the time-dependent Schrodinger equation in Equation (1) above. One possible ansatz expands the wavepacket in the complete basis of orthonormal eigenstates \(|\Psi_{j}\rangle\) of the molecular Hamiltonian \(\hat{H}_{0}\) (see e.g., Ref. [36]),  
\[|\Psi(t)\rangle=\sum_{j}c_{j}(t)e^{-iE_{j}t/\hbar}|\Psi_{j}\rangle, \tag{2}\]  
where \(|\Psi(t)\rangle\) is the excited wavepacket, \(|\Psi_{j}\rangle\) are orthonormal rovibronic wave functions and \(E_{j}\) are the corresponding energies. According to first-order perturbation theory, the time-dependent coefficients \(c_{j}(t)\) become,  
\[c_{j}(t)=iD_{js}\left[\frac{e}{\hbar}\int_{-\infty}^{t}dte^{i\omega_{js}t} \epsilon(t^{\prime})\right], \tag{3}\]  
where \(D_{js}\) are the combined dipole transition moments and the Franck-Condon factors for excitation from the initial state \(s\) to the final state \(j\). The expression in the square brackets corresponds to the energy- and the time-dependent complex excitation function with \(e\) the charge of an electron, the angular frequency \(\omega_{js}=\left(E_{j}-E_{s}\right)/\hbar\) and \(\epsilon(t^{\prime})\) the excitation field [28; 37; 38]. The trouble with this solution is that it requires that the eigenstates of the system are calculated first, which renders the approach practical only in molecules with a small number of degrees of freedom, such as the diatomic molecule D\({}_{2}\) in the present article. In the general case, one must propagate the wavepacket numerically while circumventing as much as possible the effect of the exponential scaling of the required grid (basis) with the number of degrees of freedom. One approach is to expand the molecular wavepacket in a basis of non-stationary Gaussian wavefunctions (coherent states). Methods that take this path include the ab-initio multiconfigurational Ehrenfest method (AI-MCE) [39; 40] used in this article, the closely related ab-initio multiple spawning (AIMS) [41], the hybrid ab-initio multiple cloning (AIMC) method [42], the variational multiconfigurational Gaussians (v-MCG) [43] method, and the coupled coherent-states (CCS) method [44; 45; 46]. All of these methods trace their roots to seminal semi-classical work by [PERSON] [47].