, and the values are also provided in the supplementary information (SI). Note that as this is an anharmonic system, transitions beyond the nearest To comprehend the underlying reaction mechanisms, we find it helpful to calculate the reaction rate k for an incremental N mol v , the number of the vibrational eigenstates considered in the simulation. The black line with solid circles in Fig. shows the reaction rate k as a function of N mol v . When considering only the lowest two levels (N mol v = 2), the reaction rate has a finite value, k = 2 × 10 -6 fs -1 , which stems from the tunneling process, as illustrated by Path 1 in Fig. . When the third level is included, the solvent can induce the vibrational transition between the states |v mol 1 ⟩ and |v mol 2 ⟩, opening up a cotunneling pathway, Path 2 , as exhibited in Fig. . Interestingly, we observe that the introduction of Path 2 leads to a decreased reaction rate, which could be due to destructive interference with Path 1 . The addition of the fourth state |v mol 3 ⟩ introduces another cotunneling path, 3 as shown in Fig. (a) and leads to a net increase in the reaction rate. Notably, the rate for N mol v = 4 is already close to the converged result obtained with N mol v > 6. Higher vibrational states above the barrier contribute minimally to the reaction, due to small transition probabilities (from the vibrational states |v mol 0 ⟩ and |v mol 1 ⟩) and negligibly low population probabilities. This can be seen from Fig. , where we display the population dynamics of different vibrational states.