First, we provide the numerical details of the following calculations. To be consistent with the previous studies, we adopt the same parameters of the potential energy surface (PES) in Eq. ( ), i.e. E b = 2250cm -1 and x 0 = 44.4 a.u. , which yields a barrier frequency ω b = 1000cm -1 when the molecular mass M = 1 a.u. . For the molecular dissipative bath, the characteristic frequency Ω i = 200cm -1 and a lowfriction strength with the reorganization energy λ i = 0.05ω b Ω i are used. The characteristic frequency for the cavity bath remains fixed at Ω c = 1000cm -1 and we assume a cavity lifetime of τ c = 2J c (ω c ) ω c (1-e -β ωc ) = 200fs. Thus, the reorganization energy λ c is related to the cavity lifetime τ c and frequency ω c as given by In this work, we adopt the Padé decomposition scheme for the exponential expansion of the time correlation function in Eq. ( ). The environment is thermalized at the ambient temperature T = 300K. A Padé pole number P = 3 is found to be sufficient for obtaining the converged results throughout the simulations at this temperature. The convergence is meticulously verified concerning the hierarchal tier L α , the maximal bond dimension, and the time step. For each dissipative bath mode, the range of the occupation number n α p is defined by 0 ≤ n α,p ≤ L α . For simplicity, we take all L α to be the same and L = 10 assures the convergence for the parameters given above.