Resonance theory of vibrational polariton chemistry at the normal incidence

: We present a theory that explains the resonance effect of the vibrational strong coupling (VSC) modiﬁed reaction rate constant at the normal incidence of a Fabry–Pérot (FP) cavity. This analytic theory is based on a mechanistic hypothesis that cavity modes promote the transition from the ground state to the vibrational excited state of the reactant, which is the rate-limiting step of the reaction. This mechanism for a single molecule coupled to a single-mode cavity has been conﬁrmed by numerically exact simulations in our recent work in [J. Chem. Phys. 159, 084104 (2023)]. Using Fermi’s golden rule (FGR), we formulate this rate constant for many molecules coupled to many cavity modes inside a FP microcavity. The theory provides a possible explanationfortheresonanceconditionoftheobservedVSC effect and a plausible explanation of why only at the normal incident angle there is the resonance effect, whereas, for an oblique incidence, there is no apparent VSC effect for the rate constant even though both cases generate Rabi splitting and forming polariton states. On the other hand, the current theory cannot explain the collective effect when a large number of molecules are collectively coupled to the cavity, and future work is required to build a complete microscopic theory to explain all observed phenomena in VSC.


Introduction
Recent experiments [1]- [6] have demonstrated that chemical reaction rate constants can be suppressed [1]- [4], [7]- [9] or enhanced [5], [6], [10] by resonantly coupling molecular vibrations to quantized radiation modes inside a Fabry-Pérot (FP) microcavity [11]- [13].This effect has the potential to selectively slow down competing reactions [3] or speed up a target reaction, thus achieving mode selectivity and offering a paradigm shift in chemistry.Despite extensive theoretical efforts [8], [14]- [41], the fundamental mechanism and theoretical understanding of the cavitymodified ground-state chemical kinetics remain elusive [14], [42]- [44].To the best of our knowledge, there is no unified theory that can explain all of the observed phenomena in the vibrational strong coupling (VSC) experiments [14], including (1) The resonance effect, which happens when the cavity frequency matches the bond vibrational frequency,  c =  0 , but also only happens when the in-plane photon momentum is k ‖ = 0 (the normal incidence), (2) The collective effect [1], [4], [5] which is the increase in the magnitude of VSC modification when increasing the number of molecules N (or concentration N∕), (3) The driving by thermal fluctuations without optical pumping [1], [3].(4) The isotropic disorder of the dipoles in the cavity, which is assumed in experiments with many molecules [14].
We aim to develop a microscopic theory to explain these observed VSC effects, especially focusing on understanding the resonance effect under normal incidence.Experimentally, only the resonance at normal incidence (k ‖ = 0) gives rise to VSC effects on the rate constant, while a red-detuned cavity that has a light-matter resonance at k ‖ > 0 (oblique incidence) does not give any VSC effect.This observation strongly suggests that forming Rabi splitting is not a sufficient condition for achieving the VSC-modified rate effect.Despite recent theoretical progress [18], [45], [46], the resonant condition under normal incidence remains an unresolved question.
In this work, we generalized our recently developed analytic Fermi's golden rule (FGR) rate theory of VSC in Ref. [47] by incorporating many molecules and many cavity modes for both 1D FP cavities (with only 1D in the in-plane direction) and 2D FP cavities (with 2D in the in-plane direction) cases.In particular, we evaluated the photonic mode density of states (DOS) inside a 1D FP cavity and found that it gives rise to a van-Hove-type singularity at k ‖ = 0; for a 2D FP cavity, it is found that due to the cavity modes with k ‖ > 0 propagating outside a given cavity mode extent area, the modified photon mode DOS still remains dominant around the bottom of the dispersion band where k ‖ = 0, which are the keys to account for the normal incidence condition of the VSC-modified chemical reaction rate constant.The current theory provides a possible explanation of the resonance condition for the observed VSC effect and provides a microscopic understanding of why only at the normal incident angle there is a resonance effect.

Model system
Let us consider N identical molecules coupled to many radiation modes inside a FP cavity, where Rj is the reaction coordinate for the jth molecule, V( Rj ) is the ground state potential for each reaction molecule (a double well potential for this paper as is typical for VSC simulations [22], [27], [34], [39]), and μ( Rj is the dipole operator associated with the ground electronic state of reaction coordinate Rj (electronic permanent dipole).In particular,  j,k is the angle between μ( Rj ) and the field polarization direction êk , where we only consider transverse electric (TE) polarization, such that μ( Rj ) ⋅ êk = ( Rj ) cos  j,k .A schematic illustration is provided in the top panel of Figure 1.
The photonic wavevector k (also referred to as the field propagation direction) has two components, one perpendicular to the cavity mirror k ⊥ , and the other coplanar with the cavity k ‖ .The FP cavity has the following dispersion relation, where c is the speed of light in vacuum, n c is the refractive index of the cavity, c∕n c is the speed of the light inside the cavity, and  is commonly referred to as the incident angle (where tan = k ‖ ∕k ⊥ ), which is the angle of the photonic mode wavevector k relative to the norm direction of the mirrors (see the top panel of Figure 1 for a schematic illustration).In most of the VSC experiments, n c ≈ 1.5 for the solution used inside the microcavity.Because n c ≈ 1, it will not influence the order of the magnitude of our discussion.
For simplicity, we explicitly drop n c throughout this paper.
Later, whenever we write c in an expression we should replace it with c∕n c , in principle.When k ‖ = 0 (or  = 0), the photon frequency is ( The cavity frequency  k in Eq. ( 1) is associated with the wavevector k, according to Eq. ( 2).Furthermore, qk , âk and â † k are the photonic field annihilation and creation operators for mode k, respectively.The light-matter coupling strength is where  0 is the effective permittivity inside the cavity and  is the cavity quantization volume.Each reaction coordinate Rj is coupled to its own local phonon bath described by Ĥ .
Each cavity mode qk couples to its independent bath {x k, }, accounting for the cavity loss.On the other hand, qk also couples to the dipole of each molecule μ( Rj ) with a relative angle  j,k .For 2D-FP cavities, we further define the angle between the dipole and the k ‖ plane as  j which varies from 0 to , and the angle between the field polarization and the projection of the dipole on the k ‖ plane as  j,k which varies from 0 to 2.It is easy to prove that cos  j,k = cos  j ⋅ cos  j,k .For 1D FP cavity,  j,k will reduce to  j since  j,k = 0. Details of the Hamiltonian and a schematic illustration of the orientations of the dipole operator and field polarization vector are provided in the Supplementary Material, Section I. Figure 1(a) presents a schematic illustration of the cavity dispersion relation in Eq. ( 2) (red dashed line).The molecular excitation dispersion (black dashed line) is insensitive to the incident angle and is a straight line, with energy ℏ 0 (see Eq. ( 7)).These two dispersions hybridize due to the light-matter interactions, generating polariton dispersions (the upper and lower branches with solid curves) with the color coding indicating the character of the states, with purely photonic (red), purely vibrational (black), and hybridized (yellow to orange). Figure 1(b) presents a schematic illustration of the typical cavity detuning dependence of the rate constant modifications, with the highest intensity of the modification arising at the frequency when  c =  0 (resonance condition at the normal incidence).
In this paper, we consider a reaction using a thermal barrier crossing model.Figure 2

state, |𝜈 ′
L ⟩ in the reactant well.Then, through the coupling between | ′ L ⟩ and | ′ R ⟩, a chemical reaction occurs.Finally, the vibrational excited state | ′ R ⟩ relaxes to the ground state of the product | R ⟩.The presence of the cavity mode qk explicitly enhances the transition | L ⟩ → | ′ L ⟩.The symmetric double-well model [39] is used to model the reaction, with details in Supplementary Material, Section II. Figure 2(b) shows the phonon spectral density J  () (blue) adapted from Ref. [39] as well as the effective spectral density J eff () (red) of the cavity and its associated environment that accounts for loss.Note that J eff () resembles the Brownian oscillator spectral density that centers at a particular frequency.When its peak frequency is in resonance with the quantum vibrational frequency  0 , J eff () could potentially accelerate the state-to-state quantum transitions | L ⟩ → | ′ L ⟩ (as indicated by the red arrows in Figure 2(a)).Note that the actual experimental system might not be able to be modeled as a simple symmetric double well potential as shown in Figure 2. Nevertheless, we do expect the mechanism obtained from investigating this simple model to be insightful and characteristic of the VSC problems.Recent quantum dynamics simulations [39] using models with asymmetrical double well potential, a much higher reaction barrier than ℏ 01 , or coupling cavity to the spectator mode (which in turn couple to the reaction coordinate) do show a sharp resonance modification of the rate constant.We anticipate that the current mechanistic explanation can also be used to explain these sharp resonance features.
Consider a simplified reaction mechanism outside the cavity as Note that this is the quantum description of the reaction based on quantized states, whereas the classical description is a barrier crossing along the reaction coordinate.These vibrational diabatic states can be directly obtained by computing the eigenspectum of V( R) and then diabatizing it.The dominant pathway enhanced by VSC effects is through the first excited states [47].The simplified mechanism for this reaction is that the thermal activation process causes the transition of  [47].
Considering many molecules, we focus on the single excitation subspace.This includes the ground state |G⟩ and N singly excited states | j ⟩ (where j ∈ [1, N] labels the molecules), defined as The vibrational transition dipole matrix element is which is identical for all molecules j.When measuring the absorption spectra of the molecule, the optical response shows a peak at the quantum vibrational frequency In the singly excited manifold, the light-matter coupling term, ∝ ∑ k, j qk ⊗ μ( Rj ) ⋅ êk in Eq. ( 1), will hybridize the bright excitonic state and the photon-dressed ground states, generating the polariton states [48], [49].When all dipoles are fully aligned with a given field polarization direction êk , such that μ( Rj ) ⋅ êk = ( Rj ), and under the resonance condition  k (k ‖ ) =  0 for this particular k, the light-matter hybridization generates the upper and lower polariton states, causing the Rabi splitting expressed as where Details of this standard analysis are provided in the Supplementary Material, Section III.
The formation of Rabi splitting/polariton states comes from a collective phenomenon, resulting in the well-known dependence of

√
N or equivalently √ N∕ for Ω R , which has been experimentally confirmed [4].It has been estimated that there are N ∼ 10 6 − 10 12 molecules effectively coupled to the cavity mode [14], [16], [52] for recent VSC experiments [1], [4], and Ω R ∼ 100 cm −1 when  0 ≈ 1000 cm −1 for typical VSC experiments [4], [5].Despite encouraging progress, what remains largely a mystery is how the collective light-matter coupling can induce the VSC modified rate constant [14].Another less investigated area [18], [45], [46] is why forming polaritons at a finite incident angle does not necessarily lead to the change of the VSC kinetics, and the only observed VSC effects occur at k ‖ = 0 (or  = 0).This strongly hints that forming polariton states is not a sufficient condition for VSC-modified effects, and polariton states might not be the best representation for explaining the VSC modification, because the polariton states are present in both normal and oblique incidences, and yet only the former case result in the VSC modified rate constant.
3 Theoretical results

Analytic rate constant expression
To provide a microscopic mechanism of VSC-modified reactions, we hypothesize that the cavity modes enhance the transition from ground states to a vibrationally excited state manifold of the reactant, leading to an enhancement of the steady-state population of both the delocalized states on the reactant side and the excited states manifold on the product side (right well) , which then relax to the vibrational ground state manifold on the product side (right well), For a single molecule strongly coupled to a single cavity mode, our numerical simulations [47] have confirmed the validity of this hypothesis.The proposed reaction mechanism is represented below among which k 1 ≪ k 2 , k 3 .Note that in the current work, we only consider the single excitation subspace (where one particular vibration is excited).In real experiments, many molecules could be simultaneously excited [13], with a number n ex ≈ Ne −ℏ 0 , such that 1 ≪ n ex ≪ N.For example, when N = 10 12 , ℏ 0 ≈ 5, n ex ∼ 10 9 .Future development is needed to fully account for such statistical distributions among molecules.When the molecular system is originally in the Kramers low friction regime (before the Kramers turnover [53], [54], or the so-called energy diffusion-limited regime), the cavity enhancement of the rate constant k 1 will occur [30]- [32], [34], [35], [39].This has been extensively discussed in recent theoretical work [39], [47].If we explicitly assume that where k 0 is the chemical reaction rate constant outside the cavity, and k VSC accounts for the pure cavity-induced effect.As this is a thermally activated reaction, there already exist some excited-state populations and transitions outside the cavity, which k 0 accounts for.Note that Eq. ( 10) assumes that the pure cavity effect k VSC can be added with k 0 , which is a fundamental assumption in the current theory.
To quantitatively express k VSC , we analyze the overall effect of the cavity and the photon loss environment on molecular systems by performing a normal mode transformation [55]- [57] to the Hamiltonian in Eq. ( 1) and obtaining an effective Hamiltonian, where now the cavity modes { qk } and the photon loss bath modes {x k, } (described by Ĥloss ) are transformed into effective photonic normal mode coordinates { xk, }, which are collectively coupled to the system DOF through the following term, (11) where  ≡ ∑ N j=1 ( Rj ) ⋅ cos  j is the collective system operator which does not depend on k, as we used the relation cos  j,k = cos  j ⋅ cos  j,k to separate the k-dependent and k-independent components.For simplicity, we have assumed that  j,k →  k is j-independent for the 2D cavity case, which means the molecular dipoles are distributed in a 2D plane that perpendicular to the mirrors (which is naturally true for the 1D cavity case).Further, k = cos  k ⋅ ∑  ck, xk, is the stochastic force exerted by the k-th effec- tive bath, { xk, } are the normal modes of { qk , xk, }, and the coupling constants ck, as well as bath frequencies ωk, are characterized by an effective spectral density, where  c is the cavity lifetime.Detailed derivation is provided in Supplementary Material, Section IV.The rate constant change k VSC in Eq. ( 10) originates from a purely cavity-induced effect, which promotes the transition from |G⟩ to the singly excited states manifold {| j ⟩}.
Note that this transition is mediated by the cavity operators k through the collective coupling between all molecules and the cavity modes, as is suggested by the light-matter coupling term in Eq. (11).We use FGR to estimate this transition rate constant.The coupling for this quantum transition is provided by , and the transition is mediated by the effective photon bath operators k with their spectral densities J eff ( k , ) in Eq. (12).Using FGR to estimate the transition with the frequency  =  0 (the |G⟩ , and assuming that the pathways are completely independent (i.e., no interference between pathways), we have the following expression for the overall reaction rate constant, ⋅ n( 0 ), (13) where D denotes the dimension of the in-plane direction in a FP cavity.The collective Jaynes-Cummings-type [50] coupling strength g 2 N (without cavity frequency dependence) is defined as and the 1∕N factor accounts for the normalized rate constant per molecule.Furthermore, and  is the partition function such that ∑ k  k = 1.Note that the same thermal average over different modes is also used in a recent study of electron transfer rate theory in Ref. [61].Detailed derivations are provided in Section V of the Supplementary Material.Under the continuous k ‖ limit, one can replace the sum in Eq. ( 13) with an integral as , where Δk ‖ is the spacing of the in-plane wavevector k ‖ (or the k-space lattice constant).See Section VI of the Supplementary Material for details.For 1D FP cavities, Eq. ( 13) becomes where () = e −ℏ ∕ (cf.Eq. ( 14)), and we have explicitly used cos  k = 1, and the 1D DOS is defined as Note that when all molecules are aligned with the cavity field polarization direction, such that cos  j = 1, g 2 N = N g 2 c (cf. Eq. ( 8)).When the dipole orientations are fully isotropic, For 2D FP cavities, similarly, one has (cf.Eq. ( 13)) where g 2 N is defined in Eq. ( 14), and the 2D DOS weighted by cos 2  k is defined as whereas the standard 2D-DOS is defined as Note that g ′ 2D () = g 2D ()∕2 (see the proof in Section VI of the Supplementary Material).Since there is only a 1∕2 factor difference between g ′ 2D () and g 2D (), which does not influence the shape of the rate profiles, we will regard g ′ 2D () as g 2D () in the following discussions.We further define the accumulated spectral function () as follows, (cf.Eqs. ( 17) and ( 19)) and k D VSC in Eq. ( 13) can then be written as

The resonance effect at the normal incidence
Next, we work to provide an analytic expression of () for the 1D and 2D FP cavities, which is one of the main theoretical results of this work.The 1D and 2D DOS defined in Eqs. ( 18) and ( 21) can be evaluated using the dispersion relation in Eq. ( 2).
For the one-dimensional FP cavity [46], if we ignore the influence of cavity loss ( Ĥloss in Eq. ( 1)), one can show that the DOS for the photonic modes (D = 1) is expressed as where Θ( −  c ) is the Heaviside step function.
Details of the derivations are provided in the Supplementary Material, Section VI.The DOS, g 1D (), in Eq. ( 24) has a singularity at  =  c , which is known as (the first type of) the van-Hove-type singularity [58].
Such a concentrated peak in g 1D () at  =  c has been numerically observed in Figure 1 of Ref. [46].We will turn to the case of including the effects of photon propagation in the in-plane direction later in this section.By using Eq. ( 24), we have the spectral function ( 0 ) in Eq. ( 22) for 1D FP cavities as follows, where  m → ∞ is the cutoff frequency.The integral in Eq. ( 25) gives a finite value despite the singularity in g 1D (), because only the contribution from  =  c survives.At the same time,  = ∑ k e −ℏ k = ∫ d g 1D ()e −ℏ ≈ 2e −ℏ c ∕(cΔk ‖ ), so 1∕ cancels the e −ℏ c and the 2∕(cΔk ‖ ) factor that arises from the integral.This leads to an approximate analytic expression of ( 0 ) for 1D FP cavity case as follows We have also numerically evaluated Eq. ( 25) and compared it with Eq. ( 26) for the VSC rates, presented in Figure S2 of the Supplementary Material, which shows a nearly identical behavior.The above theoretical results also suggest that for a 1D cavity, the commonly used single mode approximation [22], [27], [39] is indeed valid, because only the mode of frequency  c survives.Using the expression of ( 0 ) (Eq. ( 26)) in the rate constant expression of Eq. (23) and taking the limit of N = 1, one obtains the previous result of k VSC (see Eq. ( 35)) for a single molecule coupled to a single mode in Ref. [47].We should remind the reader that all of the existing VSC experiments were conducted with 2D cavities.
(see Eq. ( 24)) in panels (b) and (e), and the 1D accumulated spectral function () (see Eq. ( 22)) which is directly proportional to k VSC in panels (c) and (f).In panels (a)-(c), one can clearly see that under the normal incident condition  k =  0 at  = 0, () is maximized at  c =  0 and accordingly, the rate constant will also maximize based on the FGR expression (Eq.( 13)).In the detuned case ( c ≠  0 or || > 0) in panels (d)-(f), the intensity of () still peaks at  = 0, but the value of ( 0 ) diminishes at the "resonance condition"  k =  0 (for generating Rabi splitting).
This analysis also provides a possible explanation for the resonance effect at normal incidence (k ‖ = 0) for a 1D FP cavity.In Eq. ( 26), it is clear that the peak of this function is located at  c =  0 for k ‖ = 0. Thus, the VSC-modified rate constant occurs only when  c =  0 .This is because there is a van-Hove-type singularity [58] in the 1D DOS, g 1D (), which manifests itself as the 1∕ √  2 −  2 c term in Eq. ( 25), such that the integral only survives and gives a finite value at  =  c , and at  >  c , the integral becomes vanishingly small.
However, directly extending this simple consideration for the DOS cannot explain the normal incidence condition for a 2D FP cavity (even when only considering the TE polarization direction).This is because the 2D DOS g 2D () does not have any singularity.Specifically, the DOS for the photonic modes inside a 2D FP cavity is expressed as where Section VI of the Supplementary Material contain more details of this derivation.
For the 2D cavity case, one needs to consider beyond the simple DOS argument.Note that the photon loss associated with the lifetime  c only considers the loss in the k ⊥ direction.What we have not explicitly considered before was the photon traveling outside a mode area along the k ‖ direction.Let  be the effective lateral size of a given mode (which is not the cavity length), and  be the mirror distance (along the k ⊥ direction in Figure 1), so the effective quantization volume (per mode)  =  ⋅  2 .The mode lifetime can be estimated as which is propotional to k −1 ‖ when k ‖ ≪ k ⊥ , and the associated rate constant is Γ ′ 10 = 1∕ ‖ (which corresponding to the photon loss of |1 k ⟩ → |0 k ⟩).Note that  ‖ differs from the cavity lifetime  c introduced previously.Specifically,  ‖ accounts for thermal photon traveling outside a coupling area associated with a given mode  k in the in-plane direction,  c describes the loss channel only due to the escaping of the photon with a direction that is perpendicular to the mirror surface k ⊥ (and was introduced through the Ĥloss term, which was assumed to be identical for all cavity modes  k , being independent of k ‖ ).
An estimation for  is provided as follows.As mentioned before, the typical values for the VSC experiments are N ≈ 10 6 ∼ 10 12 [14], [16], [52], which is the effective number of molecules per mode (see estimations in Ref. [52]).
The effective density is estimated to be N∕ ≈ 10 20 cm −3 [59].Using  =  ⋅  2 , and the typical value for the mirror distance  = 1 μm, we have  ≈ 10 −1 ∼ 100 μm (or 10 2 ∼ 10 5 nm), which agrees with the numerical simulation in a FP cavity based on eigenfrequency analysis of the scalar Helmholtz equation [51].With the range of , one can also estimate the range of ∕c ≈ 1 ∼ 100 fs.For example, when  ∼ 300 nm, ∕c ∼ 10 −15 s −1 = 1 fs.Note that  is different than the typical length of the cavity in the in-plane direction (which is on the order of mm [6], [8]).For a photon traveling outside a particular cavity mode area, it is still within the cavity quantization area that contains many modes.On the other hand,  c usually varies from 100 fs [3] to 5 ps [60] in typical VSC experiments.
Note that the term e −ℏ k in Eq. ( 16) originates from the photon field thermal distribution, which can also be interpreted as the ratio between two photonic transition rate constants according to the detailed balance relation, i.e., e −ℏ k = Γ 01 ∕Γ 10 , (29) where Γ 01 is the rate for the |0 k ⟩ → |1 k ⟩ photonic Fock states transition due to thermal excitation, and Γ 10 = 1∕ c is the cavity loss rate along the k ⊥ direction (associated with |1 k ⟩ → |0 k ⟩), which was assumed to be identical for all k modes.Note that all of the above-mentioned excitation and decay processes are related to the thermally activated radiation (thermal photon), and not related to the pumping with an external radiation field.To account for the additional effect of photon propagating outside a given area associated with a specific mode  k , we modify the detailed balance relation (in Eq. ( 29)) by replacing the original  k with  eff ( k ), defined as follows where  ‖ (defined in Eq. ( 28)) is k ‖ -dependent.This can also be viewed as putting a  −1 c ∕ factor to  k in Eq. ( 16), where  ‖ explicitly depends on k ‖ (Eq.( 28)).Further, the partition function is also modi- Note that in Eq. ( 30), we have not considered the effect of photon leaving from the mode k ′ and re-entering into the mode k.This should be viewed as the limitation of the current theory.Future work is needed to consider this effect.
Using  eff ( k ) in Eq. ( 30), the accumulated spectral function ( 0 ) in Eq. ( 22) for 2D cavity is modified as where we used  eff () in Eq. ( 31) and  ‖ () in Eq. ( 28), and we further define the following weighting factor This  () takes a sharp maximum at k ‖ = 0 and decays quickly when k ‖ increases, because Γ ′ 10 increases quickly as k ‖ increases.This means that for a 2D cavity, as used in all existing VSC experiments, the VSC-modified rate constant is still maximized around  k (k ‖ = 0) =  c =  0 , fulfilling the normal incidence condition.Note that the correction factor can also be applied to the 1D FP cavity but does not introduce any difference in Figure 3, due to the van-Hove singularity in the DOS (see Eq. ( 24)) which dominates the entire integral, forcing Figure 4 presents the cavity dispersion relation of  k () (see Eq. ( 2)) in panels (a) and (d), the weighting factor  () (see Eq. ( 32)) in panels (b) and (e), and the 2D accumulated spectral function () (see Eq. ( 22)) for the 2D cavity case in panels (c) and (f). Figure 4(b) shows the numerical behavior of the weighting factor  () under different ∕c values (see Eq. ( 28)), among which ∕c = 1000 fs, 10 fs, 1 fs, and 0.1 fs, corresponding to  = 3 × 10 5 nm, 3 × 10 3 nm, 300 nm, and 30 nm, respectively.All are within the reasonable range of  values discussed previously.One can see that the maximal contribution still comes from k ‖ = 0, although no singularity is present.Moreover, the width becomes narrower as ∕c decreases.Note that c∕ is usually a very large quantity, so that when the incident angle  is slightly larger, Γ ′ 10 ≫ Γ 10 becomes dominant.
Figure 4(c) presents the behavior of the accumulated spectral function (), which is calculated by evaluating Eq. ( 22) numerically using trapezoidal integration within the region of  c ≤  k ≤ 5 c using 4 × 10 6 grid points, where numerical convergence is carefully checked.One can see that () peaks at  >  c when ∕c is large, and gradually moves to  =  c when ∕c decreases and Γ ′  where ∕c = 0.1 fs corresponds to  ≈ 0.03 μm, ∕c = 1000 fs corresponds to  ≈ 300 μm.(c) The accumulated spectral function ( ω) (see Eq. ( 22)) for the normal incidence case  c =  0 , where the resonance condition is reached at  = 0. (d)-(f) corresponds to the red-detuned case (oblique incidence), with  c = 0.85 0 , whose resonance condition is reached at  ≈ 32 • .The cavity lifetime is taken as  c = 200 fs.
is because the weighting factor  () is not truly singular at  c .The smaller the ∕c value, the sharper () will be.When taking the limit of ∕c → 0, Figure 4(c) reduces back to Figure 3(c).On the other hand, when ∕c → ∞, there will be no loss in the in-plane direction, corresponding to a much wider () (see the dark blue curve in Figure 4(c-f)), which is only bounded by e −ℏ .Under this condition, () still peaks at a particular frequency, but with  c >  0 .Figure 4(d-f) corresponds to the red-detuned case under oblique incidence, where  c = 0.85 0 .
With the above analysis, we have theoretically justified why the VSC-modified chemical kinetics only occurs at the normal incidence when  c =  0 for a 2D FP cavity, which agrees with experimental observations [1], [11]- [13].This is because even though there is no singularity in g 2D (), the photons propagating outside the mode area along the k ‖ direction force the 2D cavity spectra function () to peak at  =  c , forcing the normal incidence condition.
The condition for observing Rabi splitting (see Eq. ( 8)), on the other hand, is  k =  c √ 1 + tan 2  =  0 for any  ≥ 0. Although the modes with  > 0 barely contribute to k VSC , the mode density is finite (see Figure 3(e)) and for  0 >  c there will always be a mode available that satisfies  k =  0 , generating Rabi splitting at  > 0. As such, the theory provides a step forward towards understanding the fundamental difference between the condition for forming the Rabi splitting and that of the VSC resonance modification of the rate constant.This explains the experimentally observed resonance phenomena [11], [14] that occur only at  c =  0 at the normal incident angle when k ‖ = 0 (or  = 0), but not at a finite angle of  even though the resonance condition for generating Rabi splitting is fulfilled.

No apparent collective effect
For our discussion on collectivity, we begin by considering the FGR expression in Eq. (13).For simplicity, we just focus on the 1D cavity case, since for 2D cavity there is no apparent collective effect either.If all the molecules' dipoles are perfectly aligned with the cavity field polarization direction, then cos j = 1 for all molecules, j, and  = ∑ j ( Rj ).Eval- uating Eq. ( 23) using Eqs.(25) and (26) where we have explicitly approximated n( 0 ) ≈ e −ℏ 0 (cf.Eq. ( 16)).As a special case of Eq. ( 33), when  c =  0 , Eq. ( 33) where Ω R = 2 The cavity quality factor is often defined as Q =  −1 c  0 for the resonance condition.For the recent VSC experiment by Ebbesen [3], the typical values for these parameters are  c ≈ 100 fs (reading from a width of Γ c =  −1 c ≈ 53 cm −1 of the cavity transmission spectra).If the cavity frequency is  c =  0 = 1200 cm −1 , then the quality factor is Q ≈ 22.6.
However, for the current theory in Eq. ( 34), the overall rate constant would not explicitly depend on N (Eq.( 33)), meaning that only for the small N and strong coupling between molecules and the cavity mode there will be an appreciation amount of the cavity-modified effect.This is in contrast to the experimental observation of the collective effect and should be viewed as a major limitation of current theory.This limitation could be related to the fact that we have only considered the case of single excitation subspace in our theory, whereas in the experiments, a total of n ex ≈ Ne −ℏ 0 molecules could be simultaneously excited [13] due to the thermal statistics.Future work is needed for considering multiple excitations in n ex vibrations and the rate constant theory in this scenario.
When considering the disorder of the orientation between the dipole and the cavity field polarization direction, the FGR rate in Eq. ( 33) becomes upon statistical averaging of dipole orientations.For fully isotropically distributed dipoles, ⟨cos 2 ⟩ = 1/3.

Resonance behavior of k VSC
We want to demonstrate the numerical behavior of the current theory predicted by Eqs. ( 25) and (31).Because the current theory lacks the collective effect, we take the N = 1 limit and scale up the coupling strength between a single molecule and the cavity modes, as most previous work does [22], [23], [39].This leads to the expression of (cf.Eq. ( 33)) under the single mode limit (or under the 1D cavity case, see Eqs. ( 25) and ( 26)).When further considering the presence of homogeneous or inhomogeneous broadening of the molecular system, the FGR expression will be a convolution between the original FGR expression, which does not considering the broadening for the  0 (for example, Eq. ( 33)), and a broadening function (assumed to be a Gaussian), expressed as follows [47] k where where  is the variance of the Gaussian.
As expected, the k VSC expression in Eq. ( 33) should contain several characteristic physical constants, including the speed of light c in  c (see Eq. ( 3)) as it is related to light-matter interaction, Planck's constant ℏ in g c (see Eq. ( 8)) as it should be a quantum theory, and Boltzmann's constant k B in n( 0 ) as it is a thermally activated theory.We adopt a model system used in Ref. [39] to demonstrate the basic trend of k VSC predicted by the current theory.The schematic of the model is provided in Figure 2, whereas the details are provided in Supplementary Material, Section II.
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 Section VII of the Supplementary Material.The HEOM simulation requires a linear system-bath coupling Hamiltonian.To this end, we follow the previous work [22], [39] 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  qc R.
Further, we follow Ref. [39] by defining the normalized light-matter coupling strength as below, We use a similar range of  c as used in Ref. [39].
The forward rate constant from the HEOM simulation is obtained by evaluating [39], [47] where  eq ≡ P  ∕P  denotes the ratio of equilibrium population between the reactant and product, see Section VII of the Supplementary Material.The time derivative Ṗ (t) in Eq. ( 39) is evaluated numerically.For the symmetric double potential model considered in this work,  eq = 1.The limit t → t p represents that the dynamics have already entered the rate process regime (linear response regime) and t p represents the "plateau time" of the time-dependent rate which is equivalent to a flux-side time correlation function formalism.One can also view Eq. ( 39) as the flux-side correlation function that provides the time-dependent rate constant k(t), which captures both the initial transient dynamics (the oscillatory behaviors of k(t)) and the longer time rate process (plateau of k(t p )).For the FGR-based theory (Eq.( 35)), we use the value of the k 0 (outside the cavity rate constant) obtained from the HEOM simulation and report k∕k 0 = 1 + k VSC ∕k 0 .
We report the numerical value of k∕k 0 as a function of the cavity frequency  c .For the rate constant predicted by FGR, we only report the value of k∕k 0 = 1 + k VSC ∕k 0 (see Eq. ( 10)), where k VSC is evaluated using Eq. ( 36), and the variance defined in Eq. (37b) is estimated as  = 30.74cm −1 for the model parameters we used.See Supplementary Material, Section VII for details.And we directly use the numerical result of k 0 obtained from the HEOM simulation.Figure 5 presents the numerical simulations of the rate constant from HEOM as well as the FGR results.Figure 5(a) presents k(t) for the resonant case when  c =  0 , at various light-matter coupling strengths  c .One can see the plateau value of k(t) increases as  c increases.Figure 5(b) presents the case where  c <  0 where  c = 1000 cm −1 , and there is no apparent  c dependence of k(t), indicating that the coupling to the cavity has no effect.value of k∕k 0 from Eq. ( 36) (scaled by 0.4) as a function of  c , depicted by the thick solid lines.A range of light-matter coupling strength  c is explored.The FGR expression shows the sharp resonance behavior of the VSC-modified rate profile at  c =  0 = 1190 cm −1 .A similar sharp resonance has been observed in VSC experiments [1], [5], [6] and quantum dynamics simulations [39].Further, we provide the rate constant calculated from the numerically exact HEOM simulations (see Section VII of the Supplementary Material), depicted by the open circles with a thin guiding line.
Although the analytic FGR expression overestimates the rate constant by about two times, the overall agreement between the FGR expression and the HEOM numerical results is remarkable, across the range of  c and  c we explored.
Next, we explicitly consider going beyond the singlemode limit.For the 1D FP cavity, k 1D VSC reduces back to the single-mode approximation.For the 2D FP cavity, based on the expression in Eqs. ( 13) and ( 31), the VSC-modified rate constant is expressed as where Eq. ( 28)).Note that this expression also peaks at  c =  0 (as indicated in Figure 4(c)).In Eq. ( 40),  c is the lower limit of the integral with respect to d, as well as appearing explicitly in the expression of  ‖ .The result of this definite integral in Eq. ( 40) is not as simple as replacing  with  c as in the single-mode approximation (Eq.( 35)).
Figure 6 presents the FGR rates under different  c values. Figure 6(a) is the same as Figure 5(a), which corresponds to the single-mode case (or the many-mode case inside a 1D FP cavity).Figure 6(b) presents the estimated value of k∕k 0 using k VSC expression in Eq. ( 40), corresponding to the case of many modes inside a 2D FP cavity.Here, we choose ∕c = 3.33 fs, corresponding to  = 1 μm.This should be viewed as the typical value of , which is the effective lateral size of a given mode.Results obtained with a range of other choices of  are provided in the Supplementary Material, Section VIII, all of which show a sharp peak at  c ≈  0 .Note that the broadening factor (Eq. ( 36)) was not included for k 2D VSC for clarity, and one can in principle include it which will further broaden the width of the rate constant distribution.The numerical integration scheme is the same as the calculation of (), and the convergence is carefully checked.One can observe that the resonance peak is still centered around  c =  0 with minor red-shift, which demonstrates the normal incidence condition.The cavity; details of the quantum dynamics simulation results; effects of the ∕c values on the rate profiles for the 2D cavity case.

Figure 1 :
Figure 1: Top: Schematic illustration of the normal incidence condition for VSC-modified reactions.Bottom: (a) Schematic illustration of the dispersion relations of the cavity (red dashed line), the vibrational energy (gray dashed line), and hybrid polariton states (solid lines).(b) Schematic plot of reaction rate modification as a function of the cavity frequency  c .
(a) presents the first few vibrational states of the double well model, where | L ⟩ denotes the vibrational ground state of the reactant (left well), | ′ L ⟩ denotes the vibrationally excited state of the reactant, and similar for the product (right well).The red arrow represents the thermal activation process from the vibrational ground state, | L ⟩, to the vibrationally excited

Figure 2 :
Figure 2: Potential energy surface for the reaction model.The red arrows represent the thermal activation process from the vibrational ground state, | L ⟩, to the vibrationally excited state, | ′ L ⟩ in the reactant well (left side of the barrier).Through the coupling between | ′ L ⟩ and | ′ R ⟩, a chemical reaction occurs.Finally, the vibrational excited state | ′ R ⟩ relaxed to the ground state | R ⟩.(b) The effective spectral density J eff () (red curve), corresponds to the cavity and its associated loss, compared to the phonon spectral density J  () (blue).

( 15 )
is the Bose-Einstein distribution function, where  = 1∕(k B T) with k B as the Boltzmann constant and T as the temperature.For the typical parameters in VSC experiments,  0 ≈ 1200 cm −1 and room temperature 1∕ = k B T ≈ 200 cm −1 , such that ℏ 0 ≫ 1.Finally,  k represents the thermal weighting factor for accessing the cavity mode  k , with

Figure 5 :
Figure 5: Numerically exact simulation and the analytic FGR results of the rate constant.(a) The flux-side correlation functions computed by HEOM at resonance (with  c =  0 = 1190 cm −1 ).(b) The flux-side correlation functions are calculated by HEOM but off-resonance (with  c = 1000 cm −1 ).(c) The profile of the resonant VSC rate constant k∕k 0 as a function of  c with different light-matter coupling strengths,  c , obtained by FGR expressed in Eq. (36) (solid lines) and HEOM simulations (open circles with guiding thin lines), respectively.The cavity lifetime is set to be  c = 200 fs.

,
and  ‖ () = ∕ Then the reaction occurs through the diabatic couplings between | ′ L ⟩ and | ′ R ⟩, followed by a vibrational relaxation of the product state, | ′ R ⟩ → | R ⟩.The ratelimiting step for the entire process is k 1 , where k 2 ≫ k 1 such that the populations of both | ′ L ⟩ and | ′ R ⟩ reach a steady state (plateau in time), and from the steady-state approximation, the overall rate constant for the reaction is k 0 ≈ k 1 .This steady-state behavior of the | ′ L ⟩ and | ′ R ⟩ states has recently been verified by numerically exact quantum dynamics simulations