In this work, we use the numerically exact hierarchical equations of motion (HEOM) approach to propagate the quantum dynamics of this VSC model. For a practical calculation, truncation has to be made upon the number of matter states, restricting the dynamics in a relatively low-energy subspace while ensuring numerical accuracy. Here, we use the lowest F = 10 vibrational eigenstates to construct the matter Hilbert subspace (see Sec. II C.) For the model Hamiltonian in Eq. 18, the quantum subsystem is considered as ĤS + Ĥren , projected in the Hilbert subspace spanned by {|ν 0 , ...|ν 9 } (vibrational eigenstates of ĤS ), whereas ĥeff B is treated as bath DOFs and propagated implicitly based on the HEOM formalism. The theoretical details of the HEOM approach we used in this work are provided in Sec. I-A of the Supplementary Material. For HEOM propagation, there are several convergence control parameters, including (1) the number of bath terms obtained from the decomposition of the bare environment time-correlation function (TCF), (2) the time step for integration, (3) the depth of the EOMs (or the number of tiers), and (4) on-the-fly filtering error tolerance. We have carefully checked all of the above convergence parameters. More specifically, we use the fourth order Runge-Kutta (RK-4) integrator with a time step of 0.025 fs, together with the on-the-fly filtering algorithm with an error tolerance of 1×10 -7 . More details about the bath TCF decomposition schemes and numerical calculations using HEOM are provided in Sec. I-B and I-C of the Supplementary Material.