The constant mean-field (CMF) integration scheme with an initial stepsize of 0.1 fs and an error tolerance of 10 -8 is used for all the MCTDH calculations. The MCTDH coefficients (A vector) and SPFs are propagated by short iterative Lanczos (SIL) algorithm (integration order 20 and error tolerance of 10 -8 ) and Bulirsch-Stoer (BS) (integration order 8 and error tolerance of 10 -8 ) integrators, respectively. The multi-set formalism is adopted in the two-state MCTDH dynamics, and two successive modes are grouped together, creating five combined modes. In the ML-MCTDH calculations, a variable mean-field (VMF) scheme is employed, where the equations of motion are solved by the Runge-Kutta integrator of 8 th order with an accuracy of 10 -8 and initial stepsize of 10 - 4 fs. Moreover, we have conducted systematic convergence tests by varying the number of single-particle functions (SPFs), n, and the grid size of the primitive basis functions, N. The tree structure of the ML-MCTDH wavefunction can be found in section 5 of the SI.