To obtain the numerically exact rate constant for the same model, we use hierarchical equations of motion (HEOM) to simulate the population dynamics and obtain the VSC-modified rate constant, with the details provided in Sec. VII of the Supplementary Material. The HEOM simulation requires a linear systembath coupling Hamiltonian. To this end, we follow the previous work and assume that the dipole operator is linear, 𝜇( R) = R. As a result, the light-matter coupling term in Eq. 1 (for a single molecule case) is simplified as 𝜔 c qc 𝜆 • 𝜇( R) = 𝜔 c 𝜆q c R. Further, we follow Ref. 39 by defining the normalized light-matter coupling strength as below,