Theoretical model

The model QED Hamiltonian used in this work is expressed as^31–33

Fig. 1

Vibrational strong coupling (VSC) regime in polariton chemistry.

a Schematic representation of a molecule placed inside an optical cavity. b Proton-coupled electron transfer reaction of a Shin–Metiu model. Ground-state potential energy surface (PES) of the molecule as a function of the mass-weighted proton coordinate $$ \sqrt{M}R$$ (in atomic unit) for the Shin–Metiu molecular model system. The ground-state electronic density at two different nuclear configurations (at the donor and acceptor minima) are illustrated in the insets. c Ground-state permanent dipole (solid red line) as a function of the mass-weighted proton coordinate $$ \sqrt{M}R$$. d Cavity Born–Oppenheimer (CBO) surface along photonic coordinates q_c and mass-weighted reaction coordinate $$ \sqrt{M}R$$, with the white dash line representing the minimum energy path at the resonant frequency ℏω₀ = ℏω_c = 0.1706 eV and with a coupling strength η = 0.047. e A zoom-in to the reactant well of the CBO surface at the resonant frequency ℏω_c = 0.1706 eV and η = 0.376. The arrows in d and e represent the directions of two polariton normal modes. f Schematic diagram showing the Rabi splitting ℏΩ_R due to the light-matter coupling between photon-dressed vibronic-Fock states, $$ ∣{\nu }_0,1⟩$$ (photonic excitation) and $$ ∣{\nu }_1,0⟩$$ (vibrational excitation).

Vibrational polariton Rabi splitting

At the equilibrium position of the reactant R₀, one can approximate the permanent dipole as $$ \mu \left(R\right)\approx {\mu }_0+{\mu }_0^{\prime }\left(R-R_0\right)$$, where μ₀ = μ(R₀) and $$ {\mu }_0^{\prime }=\frac{\partial \mu \left(R\right)}{\partial R}{\mid }_{R_0}$$. The light-matter interaction term in $$ \widehat{H}$$ (Eq. (1)) at R₀ becomes^15,17$$ \sqrt{\frac{2{\omega }_{\mathrm{c}}}{\hslash }}{\widehat{q}}_{\mathrm{c}}\chi \cdot \mu \left(R_0\right)=\sqrt{\frac{\hslash }{2M{\omega }_0}}\chi \cdot {\mu }_0^{\prime }\left({\widehat{a}}^{†}+\widehat{a}\right)\left({\widehat{b}}^{†}+\widehat{b}\right)+\sqrt{\frac{2{\omega }_{\mathrm{c}}}{\hslash }}{\widehat{q}}_{\mathrm{c}}\chi \left({\mu }_0-{\mu }_0^{\prime }{R}_0\right)$$, where $$ {\omega }_0=\frac{\partial E^2\left(R\right)}{\partial R^2}{\mid }_{R_0}$$ is the vibrational frequency at the equilibrium nuclear configuration R₀, M is the effective mass of the nuclear vibration, $$ {\widehat{b}}^{†}$$ and $$ \widehat{b}$$ are the creation and annihilation operators for the nuclear vibration associated with the coordinate R. At the resonant condition of ω_c = ω₀, the photon-vibration interaction couples photon-dressed vibronic-Fock states $$ ∣{\nu }_0,1⟩$$ (photonic excitation) and $$ ∣{\nu }_1,0⟩$$ (vibrational excitation), inducing a Rabi splitting ℏΩ_R as follows^15,17

Fig. 2

Decrease in the rate constant as increasing light-matter couplings.

a Infrared absorption spectrum by changing the normalized light-matter coupling strength η (see Eq (2)). b The transmission coefficient κ (under the limit t → t_p) under various light-matter coupling strength (indicated by Ω_R) at the resonant frequency ℏω_c = ℏω₀ = 0.1706 eV. c The "effective change" of the Gibbs free energy barrier Δ(ΔG^‡) with respect to the coupling strength Ω_R at 300 K. d Time-dependent transmission coefficient κ(t) at various light-matter coupling strengths at the resonant frequency ℏω_c = 0.1706 eV. e, f Cavity Born–Oppenheimer surfaces $$ {\widehat{V}}_{\mathrm{CBO}}\left(R,q_{\mathrm{c}}\right)$$ at η = 0.047 and η = 0.376, respectively, with representative reactive trajectories indicated with black solid lines.

Reaction rate constant

The VSC polariton chemical kinetics can be viewed as a barrier crossing process on the cavity Born–Oppenheimer surface (CBO)^17,31,35 $$ V_{\mathrm{CBO}}\left(R,q_{\mathrm{c}}\right)=E\left(R\right)+\frac{1}{2}{\omega }_{\mathrm{c}}^2\left(\right)q_{\mathrm{c}}+\sqrt{\frac{2}{\hslash {\omega }_{\mathrm{c}}^3}}\chi \cdot \mu \left(R\right){\left)\right)}^2$$, which is a function of both q_c and R. Note that the correct QED description in Eq. (1) includes the dipole self-energy (DSE) (χ⋅μ(R))²/ℏω_c (see “Methods”). Without this term, one would get artificial changes of the barrier height and predicts a significant modification of the polariton potential energy barrier¹⁷ (see Supplementary Fig. 3b). Since we are interested in the VSC regime, the cavity mode has a similar range of frequency as the molecular vibrations, meaning that q_c evolve at a similar time scale as R. Based upon this consideration, we decide to follow the previous work^17–19 to treat both nuclear and photonic DOF classically. The electronic DOF is considered fully quantum mechanically, described by the adiabatic electronically ground-state wavefunction $$ \left.∣\right){\mathrm{\Psi }}_g\left(R\right)\left.⟩\right)$$.

It is formally rigorous to express the rate constant as the TST rate k_TST and the transmission coefficient κ as follows

To obtain a more intuitive understanding of how VSC light-matter interactions influence κ, let us consider a simplified model, $$ \widehat{H}-{\widehat{H}}_{\mathrm{vib}}$$ which only has two DOFs {R, q_c} such that we can obtain an analytic expression of the rate as k = k_TST ⋅ κ_GH. The transmission coefficient κ_GH (under the limit t → t_p) can be obtained from the Grote–Hynes (GH) theory^29,39–43. The TST rate is $$ k_{\mathrm{TST}}=\frac{{\omega }_0}{2\pi }{e}^{-\beta E_{\mathrm{b}}}$$, where E_b = E(R_‡) − E(R₀) is the potential energy barrier height measured from the bottom of the well R₀ to the top of the barrier R_‡ (see Fig. 1b), and ω₀ is the vibrational frequency of the reactant at R = R₀, and $$ \beta ={\left(k_{\mathrm{B}}T\right)}^{-1}$$. When explicitly considering the DSE, E_b remains invariant as changing the light-matter coupling strength or the photon frequency (see Eq. (1)), explaining why one can not observe any effects from a simple TST analysis¹⁸. The total rate constant k in the GH theory can be obtained using the multidimensional TST^40,44–46 (see “Methods”) as follows

Central hypothesis

With the above analysis, we conjecture that the cavity radiation mode inside the optical cavity is effectively acting as a “solvent” degree of freedom (DOF) that is coupled to the molecular reaction coordinate R, such that the presence of photonic coordinate enhance the recrossing of the reaction coordinate and reduces the transmission coefficients. A similar phenomenon is commonly referred to as the “dynamical caging” regime in simple organic reactions^30,47,48 and enzymatic catalysis^49–51, which have been successfully explained by the GH theory. Due to the low frequency of the photonic cavity mode (which is in the same range of the vibrational frequencies), we treat both R and q_c as the classical DOFs^17–19, and use the GH theory to explore the role of the cavity mode on reaction dynamics.

Decreasing κ as increasing Ω_R

Figure 2 presents the influence of increasing light-matter coupling η (thereby increasing Ω_R) on the reaction transmission coefficient κ with the model Hamiltonian presented in Eq. (1). Figure 2a presents the IR spectrum computed based on the quantum light-matter interaction (Eq. (13) in “Methods”). The numerically exact Rabi splitting ℏΩ_R is slightly deviated from 2ℏω_c ⋅ η (as indicated by Eq. (2)) due to the linear approximation ($$ \mu \left(R\right)\approx {\mu }_0+{\mu }_0^{\prime }R$$) used in Eq. (2) (see the Supplementary Fig. 2). Figure 2b presents the transmission coefficient κ obtained from direct numerical simulations (Eq. (4) under the t → t_p limit) as well as from the GH theory (solid lines) κ_GH (by solving Eq. (12) in “Methods”). The GH theory quantitatively agrees with the results from the direct numerical simulations. With an increasing Rabi splitting Ω_R, the transmission coefficient κ decreased by almost one order of magnitude, whereas the TST rate k_TST remains unchanged (due to the unchanged barrier height in the PF QED Hamiltonian). These numerical results corroborate our hypothesis that the suppression of chemical rate originates from κ, which closely resembles the experimental result (e.g., Fig. 3D in ref. ²²).

Figure 2c presents another interesting result in this work. For the PF Hamiltonian description that explicitly includes the DSE term, there is no change in k_TST because there is no change of potential energy barrier (see Fig. 1d) nor free energy barrier¹⁸. The only change in the rate comes from κ. However, one can back out the “effective change” of the free energy barrier height due to the changing κ. To this end, we use the Eyring rate equation (see “Methods”) to convert the change of rate from κ into an effective Δ(ΔG^‡). The 4 times decrease in κ presented in Fig. 2b results in ~4 kJ/mol change in “effective” Δ(ΔG^‡) in Fig. 2c at ~700 cm⁻¹ of Ω_R. We emphasize that this is not the real change of the free energy barrier height, but rather an “effective” change of ΔG^‡ according to the change of κ based on our theoretical analysis. Interestingly, the experimentally measured results of Δ(ΔG^‡) (Fig. 3C in ref. ²², for example) closely resemble our theoretical finding in Fig. 2c, with the key difference that our theoretical results suggest that these are not the actual free energy barrier changes, but entirely due to the change of κ, i.e., kinetics. Note that if one hypothesizes that an unknown mechanism to force the upper or lower vibrational polariton states to be a gateway of VSC polaritonic chemical reaction⁵², then the activation energy change should shift linearly¹⁸ with Ω_R. The experimental results, on the other hand, demonstrate a non-linearity of reaction barrier²². Our theory indicates a non-linear increase of the “effective” Δ(ΔG^‡) as increasing Ω_R due to the change of κ, closely resembles the experimental discoveries (Fig. 3C, D in ref. ²²).

Figure 2d presents the time-dependent simulation of the transmission coefficient κ(t) defined in Eq. (4). With an increasing light-matter coupling hence a larger Ω_R, the plateau value of κ(t) keeps decreasing, and at the same time, κ(t) becomes more oscillatory. This is a typical behavior of the reaction dynamics in the solvent caging regime⁵³. As the coupling between q_c and R increase, the non-Markovian dynamics of q_c can significantly influence the recrossing dynamics of the reaction coordinate R, from the “non-adiabatic” limit of a weak coupling regime to the “dynamic caging” of a strong coupling regime^39,53.

To clearly demonstrate the difference between these two regimes, we further present the Cavity BO surface $$ V_{\mathrm{CBO}}=H-\frac{P^2}{2M}-\frac{p_{\mathrm{c}}^2}{2}-H_{\mathrm{vib}}=E\left(R\right)+\frac{1}{2}{\omega }_{\mathrm{c}}^2\left(\right)q_{\mathrm{c}}+\sqrt{\frac{2}{\hslash {\omega }_{\mathrm{c}}^3}}\chi \cdot \mu \left(R\right){\left)\right)}^2$$ along R and q_c in panel (e) and (f), with a representative reactive trajectory on top (black solid curve). Figure 2e presents a typical non-adiabatic case of the GH theory. When the instantaneous friction is weak ($$ \frac{\mid {\mathcal{C}}_{\mathrm{‡}}\mid }{{\omega }_{\mathrm{c}}}\ll {\omega }_{\mathrm{b}}$$), the GH theory becomes a model of non-equilibrium solvation, where the friction from the photonic coordinate q_c does not severely impede the transitions⁵³. In this case, the transmission coefficient remains close to the case without the cavity (black curve in Fig. 2d), and the reactive trajectory crosses the barrier without much influence from q_c. Figure 2f presents a typical “dynamical caging” regime of the GH theory, where the instantaneous friction from q_c to R is strong ($$ \frac{\mid {\mathcal{C}}_{\mathrm{‡}}\mid }{{\omega }_{\mathrm{c}}}\gg {\omega }_{\mathrm{b}}$$), such that the reaction coordinate R becomes trapped in a narrow “solvent cage” on the barrier top⁵³. At longer times, the bath relaxations of $$ {\widehat{H}}_{\mathrm{vib}}$$ allow the R to move away from the barrier top, but at shorter times, the reaction coordinate R oscillates within the cavity-induced “solvent” cage⁵⁴. The trajectory recrosses the dividing surface (R_‡ = 0) many times, resulting in oscillations of κ(t) at a short time and with a small plateau value of κ(t) at t_p (see red curve in Fig. 2d). Similar dynamical caging effects from the solvent have been extensively studied in simple organic reactions (S_N1 and S_N2)^30,47,48 and enzymatic reactions^49–51, where the solvent dynamics significantly influences the reaction rate constant^{39,40,53,55,56}. Here, the cavity photonic coordinate q_c acts like a “solvent coordinate”, and for strong couplings between q_c and R, the system exhibits the dynamical caging effect which effectively slows down the reaction rate constant. This is our theoretical explanation for the observed suppression of the rate constant for VSC polariton chemical reactions^20,21,23,24.

The origin of the resonant effect

Figure 3a presents the transmission coefficient κ (when t → t_p) as a function of the photon frequency ω_c with three normalized coupling constant η (defined in Eq. (2)). The results are obtained from the GH theory (solid line) as well as the direct numerical simulation of Eq. (4) (filled circles). One can clearly see a resonant behavior of κ when changing the photon frequency, agreeing with the analytical result (Eq. (7)) of a simpler model. These findings in Fig. 3a closely resemble recent experimental results of desilylation reaction (Fig. 3A in ref. ²¹, Fig. 3B in ref. ²⁰), aldehyde/ketone Prins cyclization (Fig. 3 in ref. ²⁴), and enzymatic reaction in pepsin (Fig. 3C in ref. ²³). Note that under a relatively small light-matter coupling η = 0.047 (green), the resonant frequency that gives a minimal κ is close to ω_b, which is also close to the reactant equilibrium frequency of the reactant ω₀ in our SM model. For the parameter regime η < 0.1 (not entering into the USC), we find that the resonant condition (based on Eq. (7)) is close to ω_b.

Fig. 3

Resonant effect in vibrational strong coupling regime of polariton chemistry.

a Transmission coefficient κ as a function of the photon frequency at three different values of the coupling strength η. b–d The cavity Born–Oppenheimer surfaces V_CBO(q_c, R) under the normalized coupling strength η = 0.094 (corresponding to the blue solid line in panel (a)) at the photon frequency b ℏω_c = 2.5 meV, c 80 meV, and d 1.0 eV, with the representative reactive trajectories indicated with black solid curves.

Note that experimentally, one often plot the cavity frequency-dependent reaction kinetics against the absorption curve of vibrational polariton. With our theoretical understanding and model calculations, we conclude that these two resonant behavior have two different origins and resonant frequencies. The resonant condition observed in the IR spectrum for Rabi Splitting requires ω_c = ω₀, whereas the resonant effects for a minimum of the rate constant require ω_c ≈ ω_b. However, it is possible for a given molecular system which has ω₀ ≈ ω_b. For example, in a theoretical work (at the level of MP2 perturbation theory) by Merkel and co-workers⁵⁷, a well-studied S_N2 reaction (CH₃F+ H⁻ → CH₄+F⁻) has a ω_b = 975.5 cm⁻¹, which is close to one ground-state vibrational frequency ω₀ = 978.7 cm⁻¹. In fact, this reaction could be also an ideal one subject to future investigations of VSC modifications of reactivities. On the other hand, there are also cases where ω₀ and ω_b are different. For example, in an S_N2 reaction involving a Si–C bond cleavage in 1-phenyl-2-trimethylsilylacetylene, we find (using the geometries reported in ref. ⁵⁸ at the same level of electronic structure theory) that the computed imaginary barrier frequency to be ω_b ≈ 74 cm⁻¹, whereas the Si–C stretching frequency in the reactant well⁵⁸ is ω₀ ≈ 860 cm⁻¹.

When increasing the coupling strength to the USC regime (0.1 < η < 1.0), the resonant frequency is significantly red-shifted from ω_b. For example, when η = 0.188 (red curve), the resonant condition for reaching a minimum value of κ is 25 meV. Nevertheless, in the range of 10 meV < ℏω_c < 100 meV, κ remains a very low value around 0.2, similar to the value at ω_c = ω_b. This red-shift of resonant frequency at which the rate constant is most significantly reduced has not been observed experimentally. Our theory predicts that if VSC experiments can reach the ultra-strong coupling regime, then the resonant frequency will be significantly shifted.

The origin of this resonant behavior in VSC chemical reaction rate constant (as indicated in Eq. (7)) can also be intuitively understood by examining representative trajectories (black solid curves) on the cavity BO potential energy surfaces presented in Fig. 3b–d, with the black solid lines, indicate representative trajectories. At a very low frequency ℏω_c = 2.5 meV shown in Fig. 3b, the photon coordinate essentially remains frozen compared to the dynamics of the reaction coordinate R during the course of the reaction. As a result, under this frozen solvent limit, the transmission coefficient remains close to the no-coupling scenario. At ℏω_c = 80 meV in Fig. 3c, with $$ \frac{\mid {\mathcal{C}}_{\mathrm{‡}}\mid }{{\omega }_{\mathrm{c}}}\gg {\omega }_{\mathrm{b}}$$, the light-matter interactions lead to the dynamical caging of the reaction coordinate at the barrier top, leading to significant decrease in the transmission coefficient κ_GH. When the photon frequency is further increased (ℏω_c = 1 eV), the reactant and the product wells become separated with a narrow channel as shown in Fig. 3d. At such a high photon frequency, the reactive channel connecting the reactant and product becomes extremely narrow⁵⁹ (much narrower than the usual dynamical caging scenario depicted in Fig. 3c or Fig. 2f), such that the reactive trajectories almost follow a straight path and is no longer caged near the dividing surface. As opposed to the dynamical caging regime, the transmission coefficient in Fig. 3d is less suppressed than the minimum κ when the photon frequency is near ω_b. Similar behavior of the reaction dynamics is also observed for the USC regime (η = 0.188 in Fig. 3), where the results are provided in Supplementary Fig. 6. Therefore, the suppression of the chemical kinetics through the dynamical caging effect by the photon mode is highly sensitive to the photon frequency, proving a plausible mechanism for explaining the resonant behavior^21,23,24 of the reaction rate constant in VSC polariton chemistry.

Pauli–Fierz QED Hamiltonian

The minimal coupling QED Hamiltonian in the Coulomb gauge (the “p ⋅ A” form) is expressed as

Grote–Hynes rate theory

In multidimensional transition state theory, the reactant to product rate constant is given as^40,44–46

Considering a simplified (classical) model of the molecule-cavity hybrid system, $$ H-H_{\mathrm{vib}}=\frac{P^2}{2M}+E\left(R\right)+\frac{p_{\mathrm{c}}^2}{2}+\frac{1}{2}{\omega }_{\mathrm{c}}^2{\left(q_{\mathrm{c}}+\sqrt{\frac{2}{\hslash {\omega }_{\mathrm{c}}^3}}\chi \cdot \mu \left(R\right)\right)}^2$$ which only contains two DOFs {q_c, R} (and are viewed as classical DOFs), the normal-mode frequencies at R₀ are $$ {\mathrm{\Omega }}_{\pm }^2=\frac{1}{2}\left({\omega }_0^2+\frac{{\mathcal{C}}_0^2}{{\omega }_{\mathrm{c}}^2}+{\omega }_{\mathrm{c}}^2\right)\pm \frac{1}{2}\sqrt{{\left({\omega }_0^2+\frac{{\mathcal{C}}_0^2}{{\omega }_{\mathrm{c}}^2}+{\omega }_{\mathrm{c}}^2\right)}^2-4{\omega }_{\mathrm{c}}^2{\omega }_0^2}$$, where $$ {\mathcal{C}}_0=\sqrt{\frac{2{\omega }_{\mathrm{c}}}{M\hslash }}\chi \cdot {\mu }_0^{\prime }$$. The normal-mode frequencies at R_‡ are $$ {\mathrm{\Omega }}_{\pm }^{‡2}=\frac{1}{2}\left(-{\omega }_{\mathrm{b}}^2+\frac{{\mathcal{C}}_{‡}^2}{{\omega }_{\mathrm{c}}^2}+{\omega }_{\mathrm{c}}^2\right)\pm \frac{1}{2}\sqrt{{\left(-{\omega }_{\mathrm{b}}^2+\frac{{\mathcal{C}}_{‡}^2}{{\omega }_{\mathrm{c}}^2}+{\omega }_{\mathrm{c}}^2\right)}^2+4{\omega }_{\mathrm{c}}^2{\omega }_{\mathrm{b}}^2}$$, where $$ {\mathcal{C}}_{‡}=\sqrt{\frac{2{\omega }_{\mathrm{c}}}{M\hslash }}\chi \cdot {\mu }_{‡}^{\prime }$$. Details of the derivation of these normal-mode frequencies are provided in Supplementary Note 2. Using these normal-mode frequencies, the rate constant in Eq. (10) for the $$ \widehat{H}-{\widehat{H}}_{\mathrm{vib}}$$ model is expressed as $$ k=\frac{1}{2\pi }\frac{{\mathrm{\Omega }}_{+}{\mathrm{\Omega }}_{-}}{{\mathrm{\Omega }}_{\mathrm{+}}^{‡}}{e}^{-\beta E_{\mathrm{b}}}$$, where E_b is the energy barrier. Using the fact that $$ {\left({\mathrm{\Omega }}_{+}^{‡}{\mathrm{\Omega }}_{-}^{‡}\right)}^2=-{\omega }_{\mathrm{b}}^2{\omega }_{\mathrm{c}}^2$$ and $$ {\left({\mathrm{\Omega }}_{+}{\mathrm{\Omega }}_{-}\right)}^2={\omega }_{\mathrm{0}}^2{\omega }_{\mathrm{c}}^2$$ (see the general proof in ref. ⁴²), the rate constant can be further expressed as follows

When considering the phonon bath $$ {\widehat{H}}_{\mathrm{vib}}$$ under the Markovian limit (while consider q_c as the non-Markovian coordinate), λ can be obtained by solving the following equation

Model molecular Hamiltonian

The potential energy surface (PES) and permanent dipole moment are taken from a SM model⁶⁶, which is illustrated in Fig. 1. The SM model is a one-dimensional molecular system that describes a proton-coupled electron transfer reaction between a donor and an acceptor ion. The model consists of a transferring proton, an electron, and two fixed ions. The molecular Hamiltonian is $$ {\widehat{H}}_{\mathrm{M}}=\frac{{\widehat{P}}^2}{2M}+{\widehat{H}}_{\mathrm{el}}+{\widehat{H}}_{\mathrm{vib}}$$, where M is the mass of the nuclei (proton in the SM model), $$ {\widehat{H}}_{\mathrm{el}}={\widehat{T}}_r+{\widehat{V}}_{\mathrm{eN}}+{\widehat{V}}_{\mathrm{NN}}$$ is the electronic Hamiltonian, where $$ {\widehat{T}}_r={\widehat{p}}_{\mathrm{r}}^2/2m_{\mathrm{e}}$$ represents the kinetic energy operator of the electron with mass m_e, $$ {\widehat{V}}_{\mathrm{eN}}$$ describes the interaction between the electron and the three nuclei, and $$ {\widehat{V}}_{\mathrm{NN}}$$ that describes the Coulomb repulsion between the proton and the static ions. The resulting PES $$ E\left(R\right)=⟨{\mathrm{\Psi }}_g\left(R\right)\mid \left({\widehat{H}}_{\mathrm{M}}-{\widehat{T}}_R\right)\mid {\mathrm{\Psi }}_g\left(R\right)⟩$$ and the permanent dipole moment $$ \mu \left(R\right)=⟨{\mathrm{\Psi }}_g\left(R\right)\mid \widehat{\mu }\mid {\mathrm{\Psi }}_g\left(R\right)⟩$$ are shown in Figs. 1b and 1c, respectively. The details of this model, as well as the numerical procedure to obtain the ground-state potential and dipole are provided in Supplementary Note 3 and Supplementary Note 5.a.

In addition, $$ {\widehat{H}}_{\mathrm{vib}}={\sum }_k\frac{P_k^2}{2M_k}+\frac{1}{2}{M}_k{\omega }_k^2{\left(R_k+\frac{c_k}{M{\omega }_k^2}\cdot R\right)}^2$$ is the vibrational system-bath Hamiltonian that describes the interactions between reaction coordinate R and other vibrational phonon modes in the molecule. The coupling constant c_k and the frequency ω_k is characterized by an ohmic spectral density $$ J\left(\omega \right)=\frac{\pi }{2}{\sum }_k\frac{c_k^2}{M_k{\omega }_k}\delta \left(\omega -{\omega }_k\right)=\zeta \omega e^{-\omega /{\omega }_{\mathrm{p}}}$$, with a characteristic phonon frequency ω_p and a friction constant ζ. In Table 1, we outline several key parameters in our model system, whereas the full details are provided in Supplementary Note 3.

Table 1

Key parameters of model.

ℏω0 (meV)	ℏωb (meV)	$$ {\mu }_0^{\prime }$$ (a.u.)	$$ {\mu }_{\mathrm{‡}}^{\prime }$$ (a.u.)
170.6	162.05	0.225	−1.887

Numerical simulation of κ(t)

All simulations were performed under T = 300 K by evolving the classical dynamics governed by H(R, q_c) in Eq. (1). Langevin dynamics is used to model the influence of H_vib on the light-matter hybrid system, whereas q_c is explicitly propagated in time and treated as a non-Markovian “solvent” DOF. The friction constant in the Langevin dynamics was chosen to be ζ = 400 cm⁻¹ according to the spectral density of the $$ {\widehat{H}}_{\mathrm{vib}}$$ (see details in Supplementary Note 3 and 4). The time step used in the simulation is dt = 4 a.u., which was carefully checked to produce stable integration for all simulations. From a long constraint MD trajectory on the dividing surface R_‡ = 0, the constrained configurations {q_c, R_‡} are sampled for every 270 fs along that constrained trajectory. A total of 100,000 trajectories are released from the dividing surface, with the initial velocities randomly sampled from the classical Maxwell–Boltzmann distribution. Each of the sampled configuration is propagated for 200 fs, which guaranteed that the flux-side correlation function would plateau. The flux-side correlation function in Eq. (4) is computed through the ensemble average. Details of the numerical simulation procedure are provided in Supplementary Note 5.d.

Effective Δ(ΔG^‡)

To account for the “effective change” of the Gibbs free energy barrier Δ(ΔG^‡) corresponding to the changes in κ, we consider the Eyring rate equation $$ k=\frac{k_{\mathrm{B}}T}{h}{e}^{-\frac{\mathrm{\Delta }{G}^{‡}}{k_{\mathrm{B}}T}}$$, and thus $$ \mathrm{\Delta }{G}^{‡}=-\frac{1}{\beta }\ln \left(2\pi \beta \cdot k\right)$$. With k = κ ⋅ k_TST, we can rewrite the above ΔG^‡ as $$ \mathrm{\Delta }{G}^{‡}=-\frac{1}{\beta }\ln \kappa -\frac{1}{\beta }\ln 2\pi \beta k_{\mathrm{TST}}$$. Because k_TST is a constant at any coupling strength and cavity frequency and is the same for bare molecular case, the effective Δ(ΔG^‡) solely depends on the change of κ. The change of free energy barrier compared to the bare molecular reaction (with κ₀ and $$ \mathrm{\Delta }{G}_0^{‡}$$) is then $$ \mathrm{\Delta }\left(\mathrm{\Delta }{G}^{‡}\right)=\mathrm{\Delta }{G}^{‡}-\mathrm{\Delta }{G}_0^{‡}=-\frac{1}{\beta }\ln \frac{\kappa }{{\kappa }_0}$$, which is used to compute the value presented in Fig. 2c.

Absorption spectrum

We employ a simple approach⁶⁷ to compute the absorption spectrum of the molecule-cavity hybrid system. The absorption cross section $$ \sigma \left(\mathcal{E}\right)$$ as a function of excitation energy $$ \mathcal{E}$$ is expressed^67,68 as follows