where the diabatic potentials V 1 (x, y ) and V 2 (x, y ) are two identical two-dimensional harmonic oscillators with centers shifted with respect to one another in the x direction and displaced in energy. They are coupled through a linear potential λ(y ) in the y -direction. Figure shows the adiabatic PESs obtained from diagonalizing Eq. ( ) for the BMA model. The details of the models can be found in appendix A 1. For the quantum dynamics (QD) simulations, the time-dependent Schrödinger equation is solved on a grid in the diabatic basis using the split-operator method . The initial gaussian nuclear wavepacket is initialized in the excited where {χ l (x, y , t)} are the nuclear-coordinate and time-dependent coefficients in the BH expansion of the full state of the system and |χ(x, y , t)| 2 = l |χ l (x, y , t)| 2 is the total nuclear density. The trajectory-based (quantum-classical) analogous quantities are ensembleaverages