For the model we considered, the molecular phonon bath is described by the Drude-Lorentz spectral density, so that its time-correlation function (TCF) decomposition is computed using the Padé spectral decomposition (PSD) scheme. Here, we use the [N -1/N ] scheme with 2 low-temperature correction terms. For HEOM propagation, we use the fourth order Runge-Kutta (RK-4) integrator with a time step of 0.005 fs, together with the on-the-fly filtering algorithm with an error tolerance of 1 × 10 -6 . The number of tiers is set as 20. The convergence of the calculation is carefully checked with the above parameters. Also, a factorizable initial fulldensity matrix is applied, which is the same as Eq. 39.