where γ ν = 200 cm -1 and λ ν = 0.1ω b γ ν /2, according to the previous work in Ref. 8. The frequency of bath modes ω j and the coupling strength c j are sampled based on a numerical procedure that discretizes the J ν (ω) discussed in Ref. 22 (or in general, the algorithm described in Ref. 23), with 300 bath modes and the maximum bath frequency ω m = 2000 cm -1 . Further, we treat cavity mode qc and its loss bath Ĥc together, described by Eq. 5, characterized by J eff (ω) in Eq. 6. As such, the cavity mode and the lossy environment are described by the normal modes { Xj }. We sample the frequency Ωj and the coupling strength Cj according to J eff (ω) in Eq. 6 with an equally spaced sampling strategy. The bath frequencies Ωj are around the peak value ω c . The initial conditions of the classical DOFs {p i , x i } (phonon vibrational modes in Ĥν , Eq. 3) and { Pj , Xj } (photonic normal modes in Ĥeff , Eq. 5) are sampled based on their Wigner distributions. Computational details are provided in the Supporting Information, Sec. III. We use the symmetrical quasi-classical (SQC) window function approach (γ-SQC), the linearized semiclassical spin mapping approach (spin-LSC), Ehrenfest dynamics, and the globalflux surface hopping (GFSH) approaches as the MQC method to simulate the vibrational quantum dynamics influenced by the cavity. We have used these methods to simulate exciton-polariton dynamics in a model Shin-Metiu potential coupled to the cavity. A brief outline of each method is also provided in the Supporting Information, Sec. III. The reaction rate constant is obtained by fitting the population dynamics of |ν R ⟩ state. The rate constant from the HEOM simulations is obtained using the procedure detailed in Ref. 8 or Ref. 9, with details provided in the Supporting Information, Sec. II.