In vibrational strong coupling (VSC), molecular vibrations strongly interact with the modes of an optical cavity to form hybrid light–matter states known as vibrational polaritons. Experiments show that the kinetics of thermally activated chemical reactions can be modified by VSC. Transition-state theory, which assumes that internal thermalization is fast compared to reactive transitions, has been unable to explain the observed findings. Here, we carry out kinetic simulations to understand how dissipative processes, namely, those introduced by VSC to the chemical system, affect reactions where internal thermalization and reactive transitions occur on similar timescales. Using the Marcus–Levich–Jortner type of electron transfer as a model reaction, we show that such dissipation can change reactivity by accelerating internal thermalization, thereby suppressing nonequilibrium effects that occur in the reaction outside the cavity. This phenomenon is attributed mainly to cavity decay (i.e., photon leakage), but a supporting role is played by the relaxation between polaritons and dark states. When nonequilibrium effects are already suppressed in the bare reaction (the reactive species are essentially at internal thermal equilibrium throughout the reaction), we find that reactivity does not change significantly under VSC. Connections are made between our results and experimental observations.

## I. INTRODUCTION

Strong light–matter interaction, or simply strong coupling (SC), occurs when an optical cavity mode and a material excitation coherently exchange energy faster than either species decay.^{1–4} This interaction results in light–matter eigenstates called polaritons. In addition to inheriting both photonic and material properties, the hybrid states have different energies than their constituents.

These features make polaritons a promising platform for modifying the physicochemical properties of molecular systems.^{5,6} Typically, SC in such systems is achieved by the interaction between a cavity mode and *N* ≫ 1 quasidegenerate excitations in a molecular ensemble. This collective SC produces two polariton states and *N* − 1 dark states. At light–matter resonance, the so-called lower and upper polaritons are redshifted and blueshifted, respectively, from the bare molecular (photonic) excitations by the collective light–matter coupling strength $gN$, where *g* is the single-molecule light–matter coupling strength and $\Omega =2gN$ is the so-called Rabi splitting. In contrast, the dark states remain unchanged in energy. Nevertheless, the polaritons can dissipate into the dark states,^{7–13} enriching the relaxation dynamics of the SC regime.

In light of the considerable impact that polariton formation has on molecular excited states^{14–17} and of the photonic character of polaritons, it is not surprising that SC is emerging as a versatile tool for manipulating photochemistry. It has been shown that, in organic materials, polaritons can lower the activation barrier of spin-conversion processes,^{18–22} enhance conductivity,^{23–27} change the dynamics of photoinduced charge transfer,^{28–33} and alter the photoisomerization yield.^{34–39} In some situations, dissipation associated with SC plays a key role in the modification of photochemical processes. Notably, cavity decay (photon leakage from the cavity) can enable the suppression of photodegradation by diverting population from polaritons to the electronic ground state instead of to the product.^{40–42} The effect of cavity decay on photochemical reactivity has also been demonstrated theoretically for photodissociation^{43,44} and photoisomerization.^{38,39} Another contributor to the suppression of photodegradation is the relaxation between polaritons and dark states.^{40,42} This dissipative coupling can even mediate polariton-assisted remote energy transfer.^{45–49} Moreover, relaxation from polaritons to dark states can be harnessed to realize remote control of photoisomerization.^{50}

Much less is understood about the ability to modify thermally activated reactions using SC^{51,52} of molecular vibrations, commonly known as vibrational SC (VSC). Recent experiments^{53} demonstrate that reactions can be enhanced or suppressed by VSC without external pumping (e.g., laser excitation). Reactions affected by VSC include organic substitution,^{54–56} organic rearrangement,^{57} and hydrolysis.^{58} However, the intriguing VSC reaction kinetics still lacks a theoretical underpinning. Transition-state theory, the most commonly used reaction-rate theory, has been unsuccessful at explaining the VSC reactions.^{59–62} While a model of nonadiabatic electron transfer under VSC captures some of the observed trends,^{63} additional connections to experiments have yet to be made. Common to the attempted theoretical approaches is the following assumption: internal thermalization (i.e., within the states of each chemical species) is much faster than reactive transitions (i.e., between the states of different chemical species). Said differently, states of each chemical species are assumed to remain at internal thermal equilibrium throughout the reaction.

Nonequilibrium effects may therefore be relevant in thermally activated reactions modified by VSC. In the context of adiabatic reactions, recent findings^{64} reveal that non-Markovian dynamics along the cavity-mode coordinate can induce trapping of population in a high-energy region of the potential energy surface, preventing internal thermalization and promoting backward reactive transitions. Alternatively, it would be interesting to explore how VSC affects reactions with significant nonequilibrium effects outside the cavity. For instance, some organic reactions involve the formation of vibrationally hot intermediate species that react before they fully thermalize.^{65–67} This may happen in desilylation reactions, which feature reactive intermediates^{68} and appear in studies reporting suppression of product formation^{54,55,69} and rise in product selectivity^{55} using VSC. If nonequilibrium effects are indeed relevant to the VSC reactions, then another intriguing topic is the role of dissipation, which is what enables thermal equilibrium to be reached.

Here, we carry out a kinetic study of how the dissipative processes that VSC introduces to the chemical system—specifically, cavity decay plus incoherent energy exchange among polaritons and dark states—affect thermally activated reactions with significant nonequilibrium effects. The class of reactions we consider is electron transfer, which is modeled using Marcus–Levich–Jortner (MLJ) theory.^{70–72} We examine reactions where vibrational relaxation and reactive transitions occur on similar timescales. For the case of VSC, we use a generalization of the MLJ model to calculate rates of reactive transitions. Rates associated with internal thermalization are calculated using a standard treatment^{9} of polariton relaxation. By comparing reaction kinetics under VSC with those in the bare case and under weak light–matter coupling, we show that dissipation associated with VSC can alter reaction kinetics by suppressing nonequilibrium effects. It follows that the presence of nonequilibrium effects in the bare reaction can be important in determining the influence of VSC on reactivity. Before concluding this paper, we discuss our results in the context of the aforementioned experiments.

## II. THEORY

### A. Hamiltonian

Consider *N* identical molecules inside an optical cavity. Each molecule can undergo a series of nonadiabatic electron transfer reactions. Such reactions occur via transitions between diabatic electronic states, each representing a different reactive species (e.g., reactant, intermediate, and product). In the spirit of MLJ theory, the electronic states experience vibronic coupling to high-frequency intramolecular vibrational modes and low-frequency solvent vibrational modes,^{72} i.e., electron transfer is coupled to high- and low-frequency vibrations. Both types of modes are hence reactive degrees of freedom. For simplicity, suppose that each molecule has only one high-frequency reactive mode (hereafter conveniently referred to as vibrational mode) and assume that VSC takes place between this mode and a single cavity mode. Already, we can see that the role of VSC is to modify the dynamics of a reactive mode.

The Hamiltonian describing the electron-transfer molecules under VSC is

The first term,

characterizes the “system,” the subspace whose evolution we are interested in. Comprising the system are the electronic, vibrational, and cavity degrees of freedom. Electronic states are given by

where |*ϕ*^{(i)}⟩ has energy *E*_{ϕ} and represents the molecule *i* belonging to reactive species *ϕ*. Next,

characterizes the vibrational modes, where the mode belonging to molecule *i* is represented by the creation (annihilation) operator $ai\u2020$ (*a*_{i}). Vibronic coupling involving the high-frequency vibrational modes is described by

For any molecule of reactive species *ϕ*, *λ*_{ϕ} is the dimensionless displacement of the corresponding vibrational mode. The following two terms of *H*_{S} [Eq. (2)] are the cavity Hamiltonian and the light–matter interaction:

respectively. The cavity mode has frequency *ω*_{c} and is represented by the creation (annihilation) operator $a0\u2020$ (*a*_{0}). In writing Eq. (7), we have assumed that molecules of all reactive species undergo VSC. Presumably, this condition has been realized in studies showing that organic desilylation^{55} and TPPK-acetone cyclization^{58} can be modified by VSC, namely that involving a vibrational mode which is optically bright and nearly identical in energy for both reactant and product species. For cases where only one of reactant and product experiences VSC, theoretical modeling^{63,73} predicts interesting effects such as autocatalysis;^{63} we do not focus here on situations where different reactive species couple unequally to the cavity. Finally, the electron-transfer interaction reads

where *J*_{ϕφ} = *J*_{φϕ} is the diabatic coupling between reactive species *ϕ* and *φ*. The remaining degrees of freedom, which constitute the “bath,” are captured in *H*_{B}. The bath includes the low-frequency solvent modes that participate in electron transfer, as well as modes that induce decay and incoherent energy exchange among the polaritons and dark states. The system–bath coupling responsible for these relaxation processes is characterized by *H*_{S−B}. For the sake of brevity, we will not explicitly define *H*_{B} and *H*_{S−B} here.

The eigenstates of the system, treating the electron-transfer interaction (*V*_{ET}) as perturbative, can be found by diagonalizing the zeroth-order Hamiltonian

where *H*_{CV} = *H*_{v} + *H*_{c} + *H*_{c−v} governs the subsystem that includes cavity and vibrational modes. Two polariton and *N* − 1 dark modes make up the eigenmodes of *H*_{CV}. Let operators $\alpha q=\u2211i=0Ncqiai$ for *q* = ±, 2, …, *N* represent the eigenmodes. The polariton modes have frequencies $\omega \xb1=12\omega c+\omega v\xb1(\omega c\u2212\omega v)2+4g2N$, respectively, and are given by

where $\theta =12tan\u22121[2gN/(\omega c\u2212\omega v)]$. The dark modes all have frequency *ω*_{v}, and each one is represented by *α*_{q} for *q* = 2, …, *N*. Due to the degeneracy of the dark states, there exist multiple ways to express the dark modes in terms of the bare modes (see Refs. 63 and 74). We now define multiparticle states for various system degrees of freedom. Define $|\varphi \u2261\varphi 1,\varphi 2,\u2026,\varphi N$ as the multiparticle electronic state where molecule *i* belongs to reactive species *ϕ*_{i}. For the cavity-vibrational modes, we define |**m**⟩ ≡ |*m*_{+}, *m*_{−}, *m*_{2}, …, *m*_{N}⟩, which is an eigenstate of *H*_{CV} with *m*_{q} excitations in mode *q* and with energy $\u2211q=\xb1,2,\u2026,Nmq\u210f\omega q$. Moving to this multiparticle representation and carrying out some additional rearrangements, we can write the zeroth-order system Hamiltonian [Eq. (9)] in the diagonal form,

The eigenstates

are expressed as products of an electronic state and a displaced cavity-vibrational state. The latter is given by

where $Dq(\lambda )=exp\lambda \alpha q\u2020\u2212\lambda *\alpha q$ is the displacement operator for eigenmode *q*. The displacement of mode *q*, when the system is in electronic state |** ϕ**⟩, is

where

is the contribution due to the vibronic coupling in molecule *i* (when this molecule belongs to reactive species *ϕ*_{i}). From Eqs. (14) and (15), it is evident that VSC redistributes the displacements of the bare modes among the polariton and dark modes. How each eigenmode is displaced depends on its frequency and its overlap with the bare modes [Eq. (15)]. Returning to the electronic-cavity-vibrational eigenstates |** ϕ**;

**m**⟩, the corresponding energies are

### B. Kinetic model

Having obtained the system eigenstates and energies, we proceed to discuss the kinetic model for simulating thermally activated nonadiabatic electron transfer under VSC. We are interested in reactions where nonequilibrium effects are significant outside the cavity. To focus on how the reactions are impacted by VSC, including dissipation brought about by VSC, we limit ourselves to the case of *N* = 2 and neglect multiply excited cavity-vibrational states; the state-space truncation is a good approximation in our simulations (Sec. III), where we choose parameters such that (1) virtually all initial population resides in the zero- and first-excitation manifolds of the cavity-vibrational subspace and (2) transitions within these lower manifolds are much faster than transitions into higher manifolds. With these simplifications, we keep the dynamics tractable enough for physical interpretation while still working in the regime of collective VSC (i.e., *N* > 1). In Sec. IV, we discuss the case of many-molecule VSC.

We now present the kinetic model. The evolution of the system is governed by the master equation

The quantity *p*_{(ϕ;m)} represents the population of |** ϕ**;

**m**⟩, and

*k*(

**′;**

*ϕ***m**′|

**;**

*ϕ***m**) denotes the rate of transition |

**;**

*ϕ***m**⟩ → |

**′;**

*ϕ***m**′⟩. Processes of three types are included in the kinetic model: reactive transitions, cavity-vibrational loss and gain, and relaxation among polaritons and the dark state (there is only one dark state for

*N*= 2). The latter two sets of processes occur within the same electronic state and contribute to internal thermalization.

A reactive transition occurs between states |** ϕ**;

**m**⟩ and |

**′;**

*ϕ***m**′⟩ if the initial state is converted to the final state by the electron-transfer reaction of molecule

*i*from reactive species

*φ*to reactive species

*φ*′ ≠

*φ*(i.e., |

**⟩ ≠ |**

*ϕ***′⟩,**

*ϕ**ϕ*

_{i}=

*φ*,

*ϕ*

_{i}′ =

*φ*′, and

*ϕ*

_{j}=

*ϕ*

_{j}′ for

*j*≠

*i*). The corresponding transition rate, derived in analogy to the MLJ electron-transfer rate,

^{75}is

where $\lambda s(\phi \phi \u2032)=\lambda s(\phi \u2032\u2061\phi )$ is the low-frequency reorganization energy for the reaction between species *φ* and *φ*′, *k*_{B} is the Boltzmann constant, and *T* is the temperature. The rate depends on a generalized Franck–Condon factor $m\u0303(\varphi \u2032)\u2032|m\u0303(\varphi )$ for the cavity-vibrational states, where

We have relabeled the single dark mode of the two-molecule system as *q* = d. The undisplaced cavity-vibrational state |*m*_{q}⟩ is the single-particle eigenstate of *H*_{CV} with *m*_{q} excitations in mode *q*. Matrix elements of the displacement operator $Dq(\lambda )$ with respect to the undisplaced cavity-vibrational states are evaluated according to the property.^{76}

where $Lnk(x)$ is an associated Laguerre polynomial; for $mq\u2032$ < $mq\u2032$, the matrix element is evaluated using the same expression except with *m*_{q}′ ↔ *m*_{q} and *λ* → −*λ*^{*}.

We next discuss cavity-vibrational loss and gain. First, we consider loss. For the optical cavities used in previous experiments of VSC reactions, bare cavity states decay via the leakage of photons out into free space. In contrast, bare vibrational states decay through dissipation into solvent and intermolecular modes. In VSC, the cavity-vibrational states can decay through a combination of cavity and vibrational channels. The decay of $|\varphi ,e^q$, where $e^q$ is the unit vector representing a single excitation of mode *q* = ±, d, has the rate^{9,18}

The bare cavity decay rate is *κ*, and the bare vibrational decay rate is *γ*. For nonzero temperatures, the reverse process, gain, can occur. The corresponding rate is determined by detailed balance: $k\varphi ;e^q|\varphi ;0=k(\varphi ;0|\varphi ;e^q)exp(\u2212\u210f\omega q/kBT)$.

Finally, we have relaxation among polaritons and the dark state. Processes of this type originate from anharmonic coupling between vibrational states and low-frequency molecular modes of the local environment.^{9,11,12,77} The relaxation from $|\varphi ,e^q$ to $|\varphi ,e^q\u2032$ (where *q*, *q*′ = ±, d and *q* ≠ *q*′) has the rate^{9,18}

where Θ(*ω*) is the Heaviside step function, $n\xaf(\omega )\u2009=\u2009exp(\u210f\omega /kBT)\u22121\u22121$ is the Bose–Einstein distribution function, and $J(\omega )$ is the spectral density of the anharmonically coupled low-frequency bath modes. For the simulations in Sec. III, we assume an Ohmic spectral density:^{9} $J(\omega )=\eta \omega \u2061exp\u2212(\omega /\omega cut)2$, where *η* is a dimensionless parameter representing the anharmonic system–bath interaction and *ω*_{cut} is the cutoff frequency of the low-frequency bath modes.

Now that we have presented the various processes captured by our kinetic model, we note, for completeness, that *k*(** ϕ**′;

**m**′|

**;**

*ϕ***m**) = 0 for any transition |

**;**

*ϕ***m**⟩ → |

**′;**

*ϕ***m**′⟩ not induced by these processes, i.e., not described by Eqs. (18), (21), and (22).

To conclude this section, we summarize how VSC influences reactive transitions and internal thermalization within our kinetic model. We do so in the context of reactions with nonequilibrium effects, which arise when internal thermalization is not fast compared to reactive transitions. First, recall that VSC alters energies and redistributes the vibronic coupling of localized bare modes among the delocalized eigenmodes (Sec. II A). In view of this, Eq. (18) indicates that VSC affects the rate of reactive transitions by changing activation energy—which is related to the energies of the initial and final states—and the Franck–Condon factor between initial and final states. These rate changes can lead to modified reactivity when reactive species are at internal thermal equilibrium and when they are not. On the other hand, the additional relaxation channels created by VSC, i.e., cavity loss and gain for polaritons [Eq. (21)] and relaxation among polaritons and the dark state [Eq. (22)], help thermalize the vibrational states and only change populations when reactive species are not at internal thermal equilibrium. Thus, if creating these additional relaxation channels is the dominant effect of VSC, then nonequilibrium effects will be suppressed, as we will see in Sec. III.

## III. SIMULATIONS

In this section, we perform kinetic simulations of thermally activated nonadiabatic electron transfer under VSC. We simulate a set of representative reactions whose reactive transitions are comparable in timescale to internal thermalization. As discussed in Sec. I, we are interested in understanding how such reactions are affected by VSC-induced dissipative processes. To this end, we choose reaction parameters that highlight the impact of cavity decay, as well as relaxation among polaritons and the dark state. In addition, parameters are chosen such that the dynamics of states with more than one cavity-vibrational excitation are insignificant compared to the dynamics of states with one or zero of such excitation (Sec. II B). Parameters specific to each reaction are listed in Table I. Throughout the simulations, we use $\u210f\omega v=2000\u2009cm\u22121$ and assume a resonant cavity mode (i.e., *ω*_{c} = *ω*_{v}). We also fix the temperature at *T* = 298 K. Unless otherwise stated, we use the following parameters to characterize internal thermalization: *γ* = 0.01 ps^{−1}, *κ* = 1 ps^{−1}, *η* = 0.001, and *ω*_{cut} = 0.1*ω*_{v}. We note that the chosen cavity decay rate (*κ*) is much faster than the chosen vibrational decay rate (*γ*) and that both rates are similar to those found in VSC experiments.^{11,12,49} In calculations of reactions under VSC, we choose $g=(0.03\omega v)/2$, which yields a Rabi splitting of Ω = 0.06*ω*_{v}.

Reaction . | Figures . | Reaction type . | Parameters^{a}
. |
---|---|---|---|

1 | 1 | A → B | E_{B} = −0.6, λ_{B} = 1.5, J_{AB} = 0.01, $\lambda s(AB)=0.08$ |

2 | 2 | A → B | E_{B} = 0.95, λ_{B} = 1, J_{AB} = 0.002, $\lambda s(AB)=0.05$ |

3 | 3 | A → B → C | E_{B} = −1.05, E_{C} = −1.35, λ_{B} = 1.5, λ_{C} = 4.5, J_{AB} = 0.0003, J_{BC} = 0.02, $\lambda s(AB)=0.05$, $\lambda s(BC)=0.3$ |

Reaction . | Figures . | Reaction type . | Parameters^{a}
. |
---|---|---|---|

1 | 1 | A → B | E_{B} = −0.6, λ_{B} = 1.5, J_{AB} = 0.01, $\lambda s(AB)=0.08$ |

2 | 2 | A → B | E_{B} = 0.95, λ_{B} = 1, J_{AB} = 0.002, $\lambda s(AB)=0.05$ |

3 | 3 | A → B → C | E_{B} = −1.05, E_{C} = −1.35, λ_{B} = 1.5, λ_{C} = 4.5, J_{AB} = 0.0003, J_{BC} = 0.02, $\lambda s(AB)=0.05$, $\lambda s(BC)=0.3$ |

^{a}

For all reactions, *E*_{A} = 0 and *λ*_{A} = 0. All *J*_{ϕφ} not listed here are equal to 0. All parameters, except $\lambda \varphi $ (dimensionless), have units $\u210f\omega v$.

For comparison, we also simulate the reactions in the bare case (see Appendix A for kinetic model) and in the case of weak light–matter coupling (see Appendix B for the kinetic model and its derivation). In the latter, we set the light–matter coupling to 1% of the value used for VSC. For all simulations in the weak-coupling regime, cavity decay is fast compared to the overall reaction dynamics (which can have a different timescale than the reactive transitions; see Sec. IV for further discussion). Therefore, the weak light–matter coupling effectively induces relaxation between vibrational and cavity states ( Appendix B).^{78,79}

In all simulations, the initial population is a thermal distribution of reactant eigenstates. The distribution is determined by the Boltzmann probability function [i.e., *p*_{(ϕ;m)} ∝ exp(−*E*_{ϕ;m}/*k*_{B}*T*)] and normalized to 1. In accordance with our truncation of the cavity-vibrational subspace (Sec. II B), the initial distribution consists only of states with ≤1 excitation in this subspace.

After calculating the population evolution for each state in the system, the total population *N*_{φ}(*t*) of each reactive species *φ* at time *t* is computed according to

The projection operator

counts the number of molecules belonging to reactive species *φ*.

We first simulate a reaction where the vibrationally hot product undergoes the reverse reaction as fast as it decays [Fig. 1(a)]. The main reactive transition involves a 0 → 1 vibrational excitation upon going from the reactant to product. The vibrationally hot product either returns to the reactant via the reverse transition or relaxes to the ground state, with both processes occurring at similar rates. This scenario represents a nonequilibrium effect where the product does not fully thermalize before it undergoes the reverse reaction. Under VSC, the reaction is enhanced compared to the bare case [Fig. 1(c), blue solid line vs black dashed line]. The observed modification is consistent with the following mechanism: dissipative processes associated with VSC speed up internal thermalization and thereby force the product to decay instead of transforming back into the reactant. To test this mechanism, we carry out simulations for various values of the cavity decay rate *κ* [Fig. 1(d)] and anharmonic system–bath coupling parameter *η* [Fig. 1(e)]. The latter determines the rates of relaxation among polaritons and the dark state [Fig. 1(b), red dotted arrows]. As the cavity decay rate *κ* is lowered to zero, most of the reaction enhancement goes away [Fig. 1(e)]. In contrast, the modification largely survives as *η* is decreased. These trends reveal not only that the reaction is enhanced by accelerated internal thermalization but also that cavity decay plays the main role in speeding up the thermalization. Specifically, cavity decay of the polaritons [Fig. 1(b), yellow dashed arrows] enables the vibrationally hot product to cool much faster than it can undergo a reverse reactive transition. Relaxation from the dark state to polaritons [Fig. 1(b), red dotted arrows] plays a supporting role: it transfers the population to the (polariton) states that decay via cavity leakage.

The second reaction we investigate has a reactive transition that starts from a vibrational excited state and occurs on the same timescale as vibrational decay [Fig. 2(a)]. Although thermodynamically unfavorable, this reaction can serve as an instructional example for thermodynamically favorable situations where the most reactive channel involves a vibrational excited state in the reactant. Due to a large energy difference between reactant and product electronic states, there is only one reactive transition, which is accompanied by a 1 → 0 vibrational deexcitation. This transition happens at a rate comparable to the decay of the initial vibrational state. Since vibrational decay is extremely fast compared to its reverse process, the timescale of reactant internal thermalization is that of vibrational decay, i.e., that of the reactive transition. Intuitively speaking, the population of the reactive vibrational excited state does not reach its fully replenished (thermalized) value prior to subsequent reactive transitions. In the presence of VSC, the reaction experiences an enhancement [Fig. 2(c), blue solid line vs black dashed line]. By varying *κ* [Figs. 2(d)] and *η* [Fig. 2(e)], we again find that the reaction modification is attributed mostly to cavity decay and less so to the relaxation among polaritons and the dark state. Cavity decay facilitates rapid thermalization between the ground state and polaritons [Fig. 2(b), yellow dashed arrows], while the dark state quickly reaches thermal equilibrium through incoherent energy transfer with the polaritons [Fig. 2(b), blue dotted arrows]. As a result of these relaxation dynamics, the excited cavity-vibrational states are refilled with population before the next reactive transition takes place.

While the first two reactions are enhanced by VSC, suppression is found for the third reaction, which involves a vibrationally hot intermediate that reacts before it can fully thermalize [Fig. 3(a)].^{65–67} The first reactive transition is from reactant to intermediate and involves a 0 → 1 vibrational excitation. Next, the vibrationally hot intermediate can either make a reactive transition to the product or relax to the ground state, with similar rates for both processes. The former process is followed by vibrational decay, while the latter is followed by an intermediate-to-product reactive transition that is much slower than vibrational decay. Overall, a significant amount of intermediate population first makes a fast reactive transition and then cools, and the remaining intermediate population first cools and then makes a slow reactive transition. With VSC, the population kinetics of the reactant is unchanged [Fig. 3(c), blue solid line vs black dashed line] because the vibrationally hot intermediate, with or without VSC, is transformed into a different state much quicker than it can transition back to the reactant. However, there is greater accumulation of the intermediate [Fig. 3(d), blue solid line vs black dashed line] and suppressed formation of the product under VSC [Fig. 3(e), blue solid line vs black dashed line]. Carrying out the same analysis as done with the previous two reactions [see Figs. 3(f)–3(i)], we conclude that the rapid decay of polaritons via cavity leakage [Fig. 3(b), yellow dashed arrows], assisted by the relaxation from the dark state to polaritons [Fig. 3(b), green dotted arrows], causes the intermediate to reach thermal equilibrium before it can go through the faster reactive transition.

To obtain additional insight into nonequilibrium effects and how VSC alters reactivity, we study the reactions under weak light–matter coupling. Figures 1(c), 2(c), and 3(d)–3(e) show population kinetics for VSC (blue solid line), weak light–matter coupling (aquamarine dotted line), and the bare case (black dashed line). For all reactions, weak light–matter coupling leads to significantly modified kinetics. The direction of change is the same as that of VSC. This finding is consistent with the fact that, similar to VSC, weak light–matter coupling opens up relaxation between vibrations and the cavity and enables internal thermalization of the reactive species to be sped up by cavity decay. It also makes sense that the magnitude of change is less than that of VSC, since the vibration-cavity relaxation resulting from weak light–matter interaction is not as fast as cavity decay.^{78} Nevertheless, it is interesting that one only needs to enter the weak-coupling regime to manipulate reactivity. The same conclusion has been reached in a study of how cavity decay allows light–matter coupling to protect molecules from photodegradation.^{41}

The simulations discussed so far suggest that the VSC-induced modifications to the above reactions arise mainly from the suppression of nonequilibrium effects. To strengthen this claim, we simulate the same reactions except that the vibrational decay rate is made 100 times larger (i.e., *γ* = 1 ps^{−1} = *κ*). With this condition, we find that the VSC and bare kinetics (Fig. 4, blue solid line vs black dashed line) are either virtually identical or much more similar than in the case of slow vibrational decay. Thus, the changes in reactivity due to VSC are considerably reduced when, outside the cavity, internal thermalization is already rapid compared to reactive transitions.

While our focus has been on nonequilibrium effects and their suppression by VSC-related dissipation, we should note that the calculated reaction kinetics do exhibit some changes as a consequence of VSC modifying activation energies and redistributing vibronic coupling among the cavity-vibrational eigenmodes (Sec. II B). Because cavity decay has been identified as the major contributor to the altered kinetics in the above reactions, the impact of the minor contributors could be revealed by comparing simulations for VSC and no cavity decay (i.e., *κ* = 0 ps^{−1}) with those for the bare case. For the reaction of Fig. 1, the kinetics is noticeably different for the two conditions [Fig. 1(d), yellow solid line vs black dashed line]. It seems that the difference cannot be fully accounted for by the presence of relaxation among polaritons and dark state [Fig. 1(e)]. Hence, modified energy differences and Franck–Condon factors likely impart a small but appreciable influence on the reaction. This conclusion is corroborated by the slight enhancement observed for the VSC reaction when vibrational decay is made faster than all reactive transitions [Fig. 4(a), blue solid line vs black dashed line]; in this regime, dissipative processes introduced by VSC should have no effect. For further discussion of how changes in energetics and Franck–Condon factors due to VSC manifest as changes in reactivity, see Refs. 63 and 73 (see also Ref. 80, which discusses these topics but in the context of SC to electronic states).

## IV. CONNECTION TO EXPERIMENTS

A noteworthy difference between our simulations and the experiments of VSC reactions is the number *N* of molecules involved in VSC. In our simulations, *N* = 2. In the experiments, *N* is estimated to be 10^{6}–10^{12}.^{9,63,81}

For such large *N*, the particular reactions studied here are unlikely to be modified by VSC. To understand this statement and to identify reactions that would be significantly affected by many-molecule VSC, we present the following argument. Without loss of generality, consider a reaction featuring a single forward reactive transition, which involves a 0 → 1 vibrational excitation. This hypothetical reaction resembles that depicted in Fig. 1(a). Assuming the steady-state approximation (SSA) for the vibrationally excited product and ignoring vibrational gain (i.e., the reverse process of vibrational decay), the reaction rate can be written as *k*_{R→P} = *k*_{f}[*k*_{d}/(*k*_{r} + *k*_{d})], where *k*_{f} is the rate of the forward reactive transition, *k*_{r} is the rate of the reverse reactive transition, and *k*_{d} is the decay rate of the vibrationally excited product. The quantity in square brackets, *k*_{d}/(*k*_{r} + *k*_{d}), is the efficiency with which the vibrationally excited product decays rather than returning to the reactant. In the reaction with VSC, there are *N* + 1 forward reactive transitions, each involving the 0 → 1 excitation of a polariton or dark state; there are two and *N* − 1 such transitions, respectively. We first focus on the two polaritonic forward reactive transitions. Recall that the Franck–Condon factor for a 0-1 vibronic transition is proportional to the displacement between initial and final vibrational states [Eq. (20) with $mq\u2032$ = 1 and *m*_{q} = 0]. Since the overlap between a polariton mode and each bare vibrational mode is *O*(*N*^{−1/2}) [Eq. (10)], the displacements—and therefore the 0-1 Franck–Condon factors—associated with the polaritons are *O*(*N*^{−1/2}) times the corresponding bare values [Eqs. (14) and (15)]. Then, the polaritonic reactive transitions have rates that scale as *O*(*N*^{−1}) because the rate of a reactive transition depends on the square of a Franck–Condon factor [Eq. (18)]. Suppose that cavity decay is much faster than all other transitions from the polaritons, including reverse reactive transitions and relaxation to dark states. This scenario, which has been observed for cavity-coupled W(CO)_{6} in nonpolar solvents,^{82} results in a decay efficiency of 1 for product states with an excitation in a polariton mode. By applying the SSA to such states and neglecting cavity gain (i.e., the reverse process of cavity decay), we arrive at the following effective rate for the polaritonic reactive transitions: $kR\u2192P(p)=\epsilon kf/N$, where *ε* is a dimensionless constant that accounts for changes in the reaction rate caused by the polaritons having modified energies. In contrast, we assume that the reactive transitions involving the dark states do not afford changes in the reaction rate, given that the dark states have similar decay dynamics to the bare vibrations.^{7,9,11,12,77,82,83} Furthermore, relaxation from dark states to polaritons should be negligible: this dissipative process is mediated by local system–bath interactions (see Sec. II B) and has rates scaling as *O*(*N*^{−1}),^{9,18,47} corresponding to the probability of each molecule in either polariton [Eq. (10)]. Hence, in order for the reaction to be modified by VSC, the reaction rate ($kR\u2192P(p)$) due to the polaritonic reactive transitions should be, at least, comparable to the bare reaction rate (*k*_{R→P}), i.e.,

With arbitrary reactions in mind, criterion (25) says that for systems with experimentally relevant values of *N*, VSC can change the kinetics of thermally activated reactions if

the polaritons give rise to

*ε*≫ 1, which means that activation energies are significantly reduced compared to the bare case (this scenario is discussed in Ref. 63), and/orin the bare case, reverse reactive transitions are extremely fast compared to internal thermalization, i.e., populations are efficiently trapped in the reactants.

When at least one of these two conditions is satisfied, a large-*N* generalization of the kinetic model presented in Sec. II B should predict altered reactivity under VSC. Indeed, Ref. 63 theoretically demonstrates the ability for condition 1 to enable VSC catalysis of ground-state reactions when *N* = 10^{7}. While such a demonstration has not been carried out for condition 2, Ref. 84 [see Fig. 6(a) of that work] has shown that condition 2 enables electronic SC to boost the efficiency of harvesting triplet excitons into singlet polaritons, ultimately leading to enhanced photoluminescence compared to the bare organic material (in this case, the reverse reactive transition is the very rapid fission of singlets into triplets).

Despite the above and other differences between this theoretical work and the experiments, the suppression of nonequilibrium effects is a general phenomenon and could be relevant to the mechanism by which VSC modifies reaction kinetics in the latter. According to our calculations, the suppression of nonequilibrium effects can lead to significant increase or decrease in reactivity, depending on the specific properties of the reaction. Both types of modifications have been found experimentally.^{53} The resemblance does not stop there. Recall the experimental studies showing that the desilylation of PTA using TBAF is suppressed by VSC.^{54,69} This reaction involves an intermediate species.^{68} Based on kinetic measurements, the authors propose that VSC forces the reaction to go through a new pathway, which involves a different intermediate. Compared to the intermediate of the bare reaction, the intermediate in the VSC reaction forms more easily but reacts less readily. We would like to propose an alternative mechanism based on our results. In the bare case, there could be nonequilibrium effects that cause slow conversion of the reactant to intermediate (similar to the reaction of Fig. 2) and fast conversion of the intermediate to product (like in the reaction of Fig. 3; see Refs. 65–67). Under VSC, the intermediate would be the same as in the bare case, but nonequilibrium effects would be suppressed, thereby accelerating the reactant-to-intermediate transformation and decelerating the intermediate-to-product transformation. This mechanism is consistent with the experimentally proposed one. Future works should test the possibility of VSC-induced suppression of nonequilibrium effects in the desilylation and other reactions.

When evaluating this possibility, it is important to keep in mind that nonequilibrium effects can affect reactivity even if internal thermalization occurs on a much shorter timescale than the reaction. Whether or not nonequilibrium effects are present is determined by how fast internal thermalization is relative to reactive transitions, which may be faster than the (overall) reaction. These guiding principles are supported by our simulations. For the bare reactions shown in Figs. 1 and 3, internal thermalization occurs on a 100 ps timescale, while the reaction occurs on a 10 ns timescale [Figs. 1(a), 1(c), 3(a), and 3(c)]. Nevertheless, as already discussed in Sec. III, the ability of VSC—in particular, the associated relaxation processes—to change the reaction kinetics relies on having (in the bare case) reactive transitions with similar rates as internal thermalization. The above principles, when made specific to adiabatic reactions, can be stated as follows:^{85} regardless of the actual reaction rate or that predicted by transition-state theory, the influence that nonequilibrium effects have on reactivity depends on how the rate of internal thermalization compares to the reactive flux (i.e., the rate of product formation per unit population) at the transition-state barrier.

## V. CONCLUSIONS

In this work, we study how reactions with significant nonequilibrium effects are influenced by the dissipative channels that VSC introduces to the chemical system. By using the MLJ formalism of electron transfer as our reaction model, we present a kinetic framework that captures reactive transitions, vibrational decay, cavity decay, and relaxation among polaritons and dark states. We then simulate reactions where internal thermalization and reactive transitions occur on the same timescale. The considered reactions respectively exhibit three representative nonequilibrium effects: vibrationally hot product transforming back to reactant while thermalizing (Fig. 1), slow replenishing of reactive vibrational excited state after a reactive transition (Fig. 2), and vibrationally hot intermediate reacting before it can fully thermalize (Fig. 3). Under VSC, the first two reactions are enhanced, whereas the last reaction is suppressed. These modifications largely occur because nonequilibrium effects are suppressed by VSC-induced dissipative processes. The suppression of nonequilibrium effects is driven by cavity decay, which causes polaritons to thermalize much faster than reactive transitions. Relaxation between polaritons and dark states allows the latter to indirectly experience accelerated thermalization due to cavity decay. When we substantially increase the vibrational decay rate, i.e., (for the cavity-free case) make internal thermalization much faster than reactive transitions, the bare and VSC reactivities become much closer (Fig. 4). Finally, we discuss our work in the context of recent experiments,^{53} which demonstrate that thermally activated reactions can be enhanced or suppressed by the collective VSC of a large number of molecules. We identify types of reactions for which our theory would predict modified kinetics induced by many-molecule VSC. We also highlight resemblances between our results and experimental observations. These resemblances suggest that VSC-triggered suppression of nonequilibrium effects could play a role in the experiments.

## ACKNOWLEDGMENTS

The development of the presented models by M.D. and J.Y.-Z. was supported by the NSF Grant No. CAREER CHE 1654732. The comparison of these models with previous rate theories for thermally activated polariton chemistry by J.A.C.-G.-A. was supported by AFOSR (Award No. FA9550-18-1-0289). M.D. was also supported by a UCSD Roger Tsien Fellowship, while J.A.C.-G.-A. was also supported by a UC-MEXUS/CONACYT graduate fellowship (Reference No. 235273/472318). M.D. thanks Luis Martínez-Martínez, Sindhana Pannir-Sivajothi, and Kai Schwennicke for useful discussions.

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

### APPENDIX A: BARE CASE

For the kinetic model of reactions in the bare case (i.e., light–matter coupling strength *g* = 0), we write the zeroth-order system eigenstates [i.e., the eigenstates of *H*_{0}, Eq. (9)] in the basis of localized vibrational and cavity modes [see Eqs. (4) and (6), respectively] instead of the basis of the polariton and dark modes. In the chosen basis, the eigenstates are still of the form $|\varphi ;m\u2261|\varphi \u2297|m\u0303(\varphi )$ [Eq. (12)], except that the displaced cavity-vibrational states are given by

where |**m**⟩ ≡ |*m*_{0}, *m*_{1}, …, *m*_{N}⟩ represents the state where mode *i* = 0, 1, …, *N* has *m*_{i} excitations and $Di(\lambda )=exp\lambda ai\u2020\u2212\lambda *ai$ is the displacement operator for mode *i*. The zeroth-order energy of |** ϕ**;

**m**⟩ is

The bare system is evolved under the same assumptions (including *N* = 2) as in the case of VSC (Sec. II B). The master equation of the bare dynamics takes the general form of Eq. (17). For the reaction of molecule *i* = 1, 2 from reactive species *φ* to reactive species *φ*′ ≠ *φ*, the reactive-transition rate is calculated with Eq. (18), where the generalized Franck–Condon factors are evaluated using

for *j* = 2*δ*_{i1} + *δ*_{i2} ≠ *i* (here, *δ*_{ik} is the Kronecker delta) and the analog of Eq. (20) with *q* → *i*. The undisplaced single-particle state |*m*_{i}⟩ represents *m*_{i} excitations in mode *i*. The decay transitions have rates

The reverse rates are governed by detailed balance,

Among the singly excited vibrational states, there are no transitions of the same nature as the relaxation among polaritons and dark states, i.e., $k(\varphi ;e^i\u2032|\varphi ;e^i)=0$. As mentioned in Sec. II B, this incoherent energy exchange is caused by local system–bath interactions, which do not couple different local vibrational modes. All rates *k*(** ϕ**′;

**m**′|

**;**

*ϕ***m**) that have not been discussed in this section are taken to be 0.

### APPENDIX B: WEAK LIGHT–MATTER COUPLING

In this section, we derive the kinetic model that we use to simulate reactions where the light–matter interaction strength is $g=(3\xd710\u22124)\omega v/2$ and the cavity decay rate is *κ* = 1 ps^{−1} ≫ $gN(where\u2009\u2009N=2)$. Together, these parameters signify the regime of weak light–matter coupling. In this regime, the light–matter coupling can be treated as a perturbation to the system. Thus, the zeroth-order system eigenstates are those of the bare case ( Appendix A). However, the light–matter coupling does change the population dynamics of the system. Below, we follow a textbook approach (see Ref. 75, pp. 103–106 and 413) to show that, when the cavity and/or vibrational excitations decay on a timescale much shorter than that of the reaction dynamics, the effect of weak light–matter coupling is to induce incoherent energy exchange between the excitations. Essentially, the employed approach is a perturbative correction to the non-Hermitian dynamics of the bare case.

Formally, the dynamics of the bare system can be given by the quantum master equation,^{75}

where *ρ*(*t*) is the reduced density operator of the system and the superoperator $L0\u2212iR$ generates the bare dynamics. The Liouvillian superoperator $L0(\u22c5)=[H0,\u22c5]/\u210f$ generates the coherent (Hermitian) dynamics of the bare system. Here, *H*_{0} is the zeroth-order system Hamiltonian without light–matter interaction. The superoperator $R$ represents system–bath interaction and generates the incoherent (non-Hermitian) dynamics of the system, i.e., all reaction and relaxation dynamics of the bare case. If we act on Eq. (B1) with ⟨** ϕ**;

**m**| from the left and |

**;**

*ϕ***m**⟩ from the right, we recover the master equation governing the bare population dynamics (see Appendix A),

where, for convenience, we have defined $k(\varphi ;m)(out)\u2261\u2211(\varphi \u2032;m\u2032)\u2260(\varphi ;m)k(\varphi \u2032;m\u2032|\varphi ;m)$ as the outgoing rate of population transfer from state |** ϕ**;

**m**⟩. If we instead act on Eq. (B1) with ⟨

**;**

*ϕ***m**| from the left and |

**′;**

*ϕ***m**′⟩ ≠ |

**;**

*ϕ***m**⟩ from the right, we arrive at

where *ρ*_{(ϕ;m),(ϕ′;m′)} (*t*) = ⟨** ϕ**;

**m**|

*ρ*(

*t*)|

**′;**

*ϕ***m**′⟩ is a coherence. In writing Eq. (B3), we have assumed that each coherence is decoupled from populations and other coherences. We have also, for simplicity, neglected pure dephasing.

For the system under weak light–matter coupling, the light–matter interaction *H*_{c−v} [Eq. (7)] is added to the Hamiltonian (*H*_{0}) of the bare system. To account for this change in the quantum master equation, we simply add the superoperator $Lc\u2212v(\u22c5)=[Hc\u2212v,\u22c5]/\u210f$ to the superoperator ($L0\u2212iR$) that generates the bare dynamics. In doing so, we implicitly ignore the effect of light–matter coupling on system–bath interaction; this approximation should hold for sufficiently small values of light–matter interaction strength. Quantum master equation [Eq. (B1)] now reads

With this transformation, states with a single vibrational excitation (i.e., $|\varphi ,e^i$, *i* = 1, 2) evolve as

and states with a single cavity excitation (i.e., $|\varphi ,e^0$) evolve as

To arrive at Eqs (B5) and (B6), we have used the general relation $\rho (\varphi \u2032;m\u2032),(\varphi ;m)=\rho (\varphi ;m),(\varphi \u2032;m\u2032)*$. For states with no vibrational or cavity excitations, the equations of motion remain the same as in the bare case. From Eqs. (B5) and (B6)—in particular, the last line of each equation—it is clear that the evolution of vibrational and cavity populations now depends on the light–matter coherence $\rho (\varphi ;e^0),(\varphi ;e^i)$. In the presence of weak light–matter coupling, this light–matter coherence evolves according to

where *j* = 2*δ*_{i1} + *δ*_{i2} ≠ *i*. We have defined Δ ≡ *ω*_{c} − *ω*_{v}. Notice the appearance of $\rho (\varphi ;e^j),(\varphi ;e^i)$, a coherence involving different vibrational states, in the last line of Eq. (B7). In the weak-coupling regime, the evolution of this purely vibrational coherence is given by

It is apparent that the coherence $\rho (\varphi ;e^j),(\varphi ;e^i)$ does not directly couple to populations. In other words, the coupling of populations to the purely vibrational coherence is of higher order in *g*, the strength of the perturbative light–matter interaction. We are interested in the lowest-order correction to the population dynamics due to light–matter interaction. Hence, we neglect the contribution of $\rho (\varphi ;e^j),(\varphi ;e^i)$ in Eq. (B7). This truncation, followed by formal integration of Eq. (B7), leads to

where we have made the change in variable (*t* − *t*′) → *t*′. For our simulations, we are interested in how populations change over times much longer than cavity and vibrational decay. Since these two processes are included in $k(\varphi ;e^0)(out)$ and $k(\varphi ;e^i)(out)$, respectively, then the exponential term will decay on a timescale much shorter than the timescales of interest. We therefore make a Markovian approximation—i.e., the substitutions $p(\varphi ;e^i)(t\u2212t\u2032)\u2192p(\varphi ;e^i)(t)$ and $p(\varphi ;e^0)(t\u2212t\u2032)\u2192p(\varphi ;e^0)(t)$, as well as extending the integral to infinity—to obtain

where

Equation (B13) is in line with the standard expressions for the Purcell factor.^{1,86} Equations (B11) and (B12) reveal that the dynamics under weak light–matter coupling is identical to the bare dynamics [Eq. (B2)], except that vibrational mode *i* incoherently exchanges energy with the cavity at rate $\gamma \varphi i\u2032$. Thus, the master equation for the weak coupling regime is that of the bare case but with the following rates for relaxation between vibrational and cavity excitations:

We reiterate that this result rests on the condition that at least one of cavity and vibrational excitations decays on a timescale faster than that of the reaction dynamics. In all our simulations involving weak light–matter coupling, we only consider cases where this separation of timescales is satisfied, and so, we use the kinetic framework described by Eq. (B14) and the immediately preceding text.