Within this framework, as depicted schematically in Fig. , the operators p c and x c in the Hamiltonian (see Sec. II A) for the cavity radiation mode are parameterized as the classical trajectories p t c and x t c , respectively. Meanwhile, all other DoFs remain in the quantum-mechanical domain, and the quantum dynamics are described by the density operator ρ(p t c ,x t c ,t). We start with a factorized initial state between the molecular vibrational mode and the baths, ρ(0) = ρ m (0) ⊗ ∏ α ρ α , with each bath being in its thermal equilibrium at temperature T . To obtain the molecular reaction dynamics, we adopt the HEOM method. A key idea in the HEOM method involves the exponential expansion of the bath correlation function, defined as