, where Δ𝑡 is the integration timestep and ∆𝑟 𝑖𝑗 and ∆𝑎 𝑖𝑗 are the relative distance and acceleration for the pair of particles 𝑖 and 𝑗. 𝜂 is the timestep parameter. For a harmonic oscillator system with a reasonably high value (taken from Ref. ) of 𝜂 = 2 -3 , the timestep relation reads Δ𝑡 ≈ = 0.1 fs for a typical nuclear vibration of 3000 cm -1 . Hence, a simulation of 10 ps with a timestep of 0.1 fs yields 10 5 integration steps. Multiple computational methods have been developed in the last few decades to simulate excited state dynamics beyond the BO approximation. Some of the most well-known are multiconfiguration time-dependent Hartree (MCTDH) and its multilayer variant ML-MCTDH, trajectory surface hopping in its multiple flavors, mean-field Ehrenfest (MFE), multiconfigurational Ehrenfest (MCE), mixed quantum-classical Liouville equation (QCLE), ab initio multiple spawning (AIMS) and methods based on the exact factorization of the electron-nuclear wave function.