, driving the evolution of the nuclear wavefunction via a nuclear TDSE . The electron-nuclear coupling operator, in the electronic evolution equation , as it is explicitly dependent on the nuclear wavefunction χ(R, t), and via the time-dependent vector potential, A ν (R, t), also contains implicitly the electronic wavefunction Φ(r, t; R). Quantum molecular dynamics simulations performed in this work are based either on the Born-Huang representation or on the exact factorization. More specifically, vibronic wavepacket dynamics is carried out using the electronic diabatic representation, rather than the adiabatic representation, where essentially in Eq. ( ), the BO PESs are replaced by the diabatic PESs and the nonadiabatic operator is replaced by the diabatic potential coupling. The transformation from diabatic-to-adiabatic basis allows us to reconstruct information along the dynamics in the adiabatic representation, for a direct comparison with trajectory-based results. Trajectory-based calculations are performed in the adiabatic basis, and even when working with the exact factorization, the electronic conditional amplitude is represented as a Born-Huang-like expansion as discussed above. More details on the trajectory-based algorithms will be given in Section III.