For our application, the number vectors in the Krylov subspace remains more or less constant for all basis set sizes and from 2 to 25 degrees of freedom, but it depends on the time step. For 1 or 5 fs time step, the number of vectors are about 20 or 50, respectively. This is clearly an advantage over Chebyshev scheme for which, the number of vectors depends on the full Hamiltonian spectral range which increases with the basis set size. Furthermore, for the SIL scheme and by construction, the wavepacket norm is conserved and the energy is virtually conserved. The basis sets are defined as follows. (i) A real Fourier series for the first mode, associated to the coordinate ϕ in a [-π, π] range. The quadrature grid is defined by equidistant points and the weights are constant. (ii) A scaled Harmonic oscillator basis set, HO, for the other modes. The not-normalized expression of the j th basis function of the i th mode, q i , is b i j (q i ) = H j-1 (x(q i )) • Exp(-1/2 • x(q i ) 2 ) where x(q i ) = s HO i • (q i -q HO 0,i ). Each parameter s HO i is adapted to the curvature and the mass of the 1D harmonic oscillator model associated to the mode q i . For the second mode of the system, s HO 2 = 0.9592722 while for all the bath modes s HO i = 1.0. The parameters q HO 0,i are defined, usually, at the corresponding equilibrium positions of the potential. However, since they are different for the two diabatic potentials, q HO 0,i are defined in the middle of the two minima. From the expressions of the diabatic potentials (Eqs. ( ) and ( )), q HO 0,2 = -κ/(2ω) for second mode of the system and q HO 0,i = -c i /(2ω i ) for the bath modes. The size of the basis sets, nb i λi , has been checked carefully with respect to the convergence of the diabatic populations. Therefore, for the direct-product calculations, nb 1 , nb 2 and the nb bath (for all bath modes) are, respectively, 256, 30 and 8. With those values, the diabatic and adiabatic populations up to 250 fs, the total energy and the other properties of interest are well converged. This means that, when the nb i values are increased, respectively, to 300, 40 and 12, the analysed properties remain unchanged (see Table ).