Resonance theory of vibrational polariton chemistry at the normal incidence

3.1 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) | ν R ′ j 〉 , which then relax to the vibrational ground state manifold on the product side (right well), | ν R j 〉 . 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 ₁ ≪ k ₂, k ₃. 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 ≈ N e − β ℏ ω 0 , such that 1 ≪ n _ex ≪ N. For example, when N = 10¹², βℏω ₀ ≈ 5, n _ex ∼ 10⁹. 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 ₁ will occur [30], [31], [32], [34], [35], [39]. This has been extensively discussed in recent theoretical work [39], [47]. If we explicitly assume that k ₁ ≪ k ₂, k ₃, then | G 〉 → k 1 | ν L ′ j 〉 is the rate limiting step, and the population of intermediate states will reach a steady-state behavior. As such, due to the steady-state approximation, the overall rate constant is [47]

where k ₀ 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 ₀ accounts for. Note that Eq. (10) assumes that the pure cavity effect k _VSC can be added with k ₀, 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], [56], [57] to the Hamiltonian in Eq. (1) and obtaining an effective Hamiltonian, where now the cavity modes { q ̂ k } and the photon loss bath modes { x ̂ k , ζ } (described by H ̂ loss ) are transformed into effective photonic normal mode coordinates { x ̃ ̂ k , ζ } , which are collectively coupled to the system DOF through the following term,

where S ̂ ≡ ∑ j = 1 N μ ( R ̂ j ) ⋅ 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, F ̂ k = cos ϕ k ⋅ ∑ ζ c ̃ k , ζ x ̃ ̂ k , ζ is the stochastic force exerted by the k-th effective bath, { x ̃ ̂ k , ζ } are the normal modes of { q ̂ k , x ̂ k , ζ } , and the coupling constants c ̃ k , ζ 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 F ̂ 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 S ̂ , and the transition is mediated by the effective photon bath operators F ̂ k with their spectral densities J _eff(ω _k , ω) in Eq. (12). Using FGR to estimate the transition with the frequency ω = ω ₀ (the | G 〉 → k 1 | ν L ′ j 〉 transition), 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,

where D denotes the dimension of the in-plane direction in a FP cavity. The collective Jaynes-Cummings-type [50] coupling strength g N 2 (without cavity frequency dependence) is defined as

and the 1/N factor accounts for the normalized rate constant per molecule. Furthermore,

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, ω ₀ ≈ 1200 cm⁻¹ and room temperature 1/β = k _B T ≈ 200 cm⁻¹, such that βℏω ₀ ≫ 1. Finally, P k represents the thermal weighting factor for accessing the cavity mode ω _k , with

and Z is the partition function such that ∑ k P 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 ∑ k f ( k ) → ∫ d k D ( Δ k ‖ ) D f ( k ) , 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 P ( ω ) = e − β ℏ ω / Z (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 N 2 = N g c 2 (cf. Eq. (8)). When the dipole orientations are fully isotropic, ∑ j = 1 N ⁡ cos 2 φ j = N 〈 cos 2 ⁡ φ 〉 = N / 3 .

For 2D FP cavities, similarly, one has (cf. Eq. (13))

where g N 2 is defined in Eq. (14), and the 2D DOS weighted by cos² ϕ _k is defined as

whereas the standard 2D-DOS is defined as

Note that g 2 D ′ ( ω ) = g 2 D ( ω ) / 2 (see the proof in Section VI of the Supplementary Material). Since there is only a 1/2 factor difference between g 2 D ′ ( ω ) and g _2D(ω), which does not influence the shape of the rate profiles, we will regard g 2 D ′ ( ω ) as g _2D(ω) in the following discussions.

We further define the accumulated spectral function A ( ω ) as follows, (cf. Eqs. (17) and (19))

and k VSC D in Eq. (13) can then be written as

3.2 The resonance effect at the normal incidence

Next, we work to provide an analytic expression of A ( ω ) 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 ( H ̂ 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 A ( ω 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, Z = ∑ k e − β ℏ ω k = ∫ d ω g 1 D ( ω ) e − β ℏ ω ≈ 2 e − β ℏ ω c / ( c Δ k ‖ ) , so 1 / Z cancels the e − β ℏ ω c and the 2/(cΔk _‖) factor that arises from the integral. This leads to an approximate analytic expression of A ( ω 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 A ( ω 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.

Figure 3 presents the cavity dispersion relation of ω _k (θ) (see Eq. (2)) in panels (a) and (d), the 1D DOS g _1D(ω) (see Eq. (24)) in panels (b) and (e), and the 1D accumulated spectral function A ( ω ) (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 = ω ₀ at θ = 0, A ( ω ) is maximized at ω _c = ω ₀ and accordingly, the rate constant will also maximize based on the FGR expression (Eq. (13)). In the detuned case (ω _c ≠ ω ₀ or |θ| > 0) in panels (d)–(f), the intensity of A ( ω ) still peaks at θ = 0, but the value of A ( ω 0 ) diminishes at the “resonance condition” ω _k = ω ₀ (for generating Rabi splitting).

Figure 3:

Dispersion relation, DOS, and the accumulated spectral function A ( ω ) for a 1D cavity. (a) The cavity dispersion relation (Eq. (2)), (b) the schematic 1D DOS g _1D(ω) (Eq. (24)) and (c) the 1D accumulated spectral function A ( ω ) (Eq. (22), evaluated using Eq. (25)) for the normal incidence case ω _c = ω ₀, where the resonance condition is reached at θ = 0. (d)–(f) corresponds to the red-detuned case (oblique incidence), with ω _c = 0.85 ω ₀, whose resonance condition is reached at θ ≈ 32°. The cavity lifetime is taken as τ _c = 200 fs.

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 = ω ₀ for k _‖ = 0. Thus, the VSC-modified rate constant occurs only when ω _c = ω ₀. This is because there is a van-Hove-type singularity [58] in the 1D DOS, g _1D(ω), which manifests itself as the 1 / ω 2 − ω c 2 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 D be the effective lateral size of a given mode (which is not the cavity length), and L be the mirror distance (along the k _⊥ direction in Figure 1), so the effective quantization volume (per mode) V = L ⋅ D 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 H ̂ loss term, which was assumed to be identical for all cavity modes ω _k , being independent of k _‖).

An estimation for D is provided as follows. As mentioned before, the typical values for the VSC experiments are N ≈ 10⁶ ∼ 10¹² [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 / V ≈ 1 0 20 cm⁻³ [59]. Using V = L ⋅ D 2 , and the typical value for the mirror distance L = 1 μ m, we have D ≈ 1 0 − 1 ∼ 100 μ m (or 10² ∼ 10⁵ 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 D , one can also estimate the range of D / c ≈ 1 ∼ 100 fs. For example, when D ∼ 300 n m , D / c ∼ 1 0 − 15 s − 1 = 1 f s . Note that D 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.,

where Γ₀₁ is the rate for the |0 _k ⟩ → |1 _k ⟩ photonic Fock states transition due to thermal excitation, and Γ₁₀ = 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 P k with P eff ( ω k ) , defined as follows

where τ _‖ (defined in Eq. (28)) is k _‖-dependent. This can also be viewed as putting a τ c − 1 / τ c − 1 + τ ‖ − 1 correction factor to P k in Eq. (16), where τ _‖ explicitly depends on k _‖ (Eq. (28)). Further, the partition function is also modified as Z → Z eff = ∑ k τ c − 1 e − β ℏ ω k / τ c − 1 + τ ‖ − 1 . As expected, when k _‖ = 0, τ ‖ − 1 = 0 , one should have P eff ( ω k ) → P k = e − β ℏ ω k / Z . 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 P eff ( ω k ) in Eq. (30), the accumulated spectral function A ( ω 0 ) in Eq. (22) for 2D cavity is modified as

where we used P eff ( ω ) in Eq. (31) and τ _‖(ω) in Eq. (28), and we further define the following weighting factor

This F ( ω ) 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 = ω ₀, fulfilling the normal incidence condition. Note that the correction factor τ c − 1 / τ c − 1 + τ ‖ − 1 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 P eff ( ω k ) → P k (as τ _‖ → ∞ when k _‖ = 0).

Figure 4 presents the cavity dispersion relation of ω _k (θ) (see Eq. (2)) in panels (a) and (d), the weighting factor F ( ω ) (see Eq. (32)) in panels (b) and (e), and the 2D accumulated spectral function A ( ω ) (see Eq. (22)) for the 2D cavity case in panels (c) and (f). Figure 4(b) shows the numerical behavior of the weighting factor F ( ω ) under different D / c values (see Eq. (28)), among which D / c = 1000  fs, 10 fs, 1 fs, and 0.1 fs, corresponding to D = 3 × 1 0 5  nm, 3 × 10³ nm, 300 nm, and 30 nm, respectively. All are within the reasonable range of D 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 D / c decreases. Note that c / D is usually a very large quantity, so that when the incident angle θ is slightly larger, Γ 10 ′ ≫ Γ 10 becomes dominant.

Figure 4:

Same as Figure 3, but with a 2D FP cavity. (a) Same as Figure 3(a). (b) The weighting factor F ( ω ) (see Eq. (32)) under different D / c values, where D / c = 0.1 fs corresponds to D ≈ 0.03  µm, D / c = 1000  fs corresponds to D ≈ 300 μ m. (c) The accumulated spectral function A ( ω ̃ ) (see Eq. (22)) for the normal incidence case ω _c = ω ₀, where the resonance condition is reached at θ = 0. (d)–(f) corresponds to the red-detuned case (oblique incidence), with ω _c = 0.85 ω ₀, whose resonance condition is reached at θ ≈ 32°. The cavity lifetime is taken as τ _c = 200 fs.

Figure 4(c) presents the behavior of the accumulated spectral function A ( ω ) , which is calculated by evaluating Eq. (22) numerically using trapezoidal integration within the region of ω _c ≤ ω _k ≤ 5ω _c using 4 × 10⁶ grid points, where numerical convergence is carefully checked. One can see that A ( ω ) peaks at ω > ω _c when D / c is large, and gradually moves to ω = ω _c when D / c decreases and Γ 10 ′ dominates the behavior for k _‖ > 0. Additionally, compared to Figure 3(c), here A ( ω ) tails towards the higher energy. This is because the weighting factor F ( ω ) is not truly singular at ω _c. The smaller the D / c value, the sharper A ( ω ) will be. When taking the limit of D / c → 0 , Figure 4(c) reduces back to Figure 3(c). On the other hand, when D / c → ∞ , there will be no loss in the in-plane direction, corresponding to a much wider A ( ω ) (see the dark blue curve in Figure 4(c–f)), which is only bounded by e^−βℏω . Under this condition, A ( ω ) still peaks at a particular frequency, but with ω _c > ω ₀. Figure 4(d–f) corresponds to the red-detuned case under oblique incidence, where ω _c = 0.85ω ₀.

With the above analysis, we have theoretically justified why the VSC-modified chemical kinetics only occurs at the normal incidence when ω _c = ω ₀ for a 2D FP cavity, which agrees with experimental observations [1], [11], [12], [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 A ( ω ) 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 ω ₀ > ω _c there will always be a mode available that satisfies ω _k = ω ₀, 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 = ω ₀ 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.

3.3 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 S ̂ = ∑ j μ ( R ̂ j ) . Evaluating Eq. (23) using Eqs. (25) and (26) leads to

where we have explicitly approximated n ( ω 0 ) ≈ e − β ℏ ω 0 (cf. Eq. (16)). As a special case of Eq. (33), when ω _c = ω ₀, Eq. (33) becomes

where Ω R = 2 N g c ⋅ ω 0 . The cavity quality factor is often defined as Q = τ c − 1 ω 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 = τ c − 1 ≈ 53  cm⁻¹ of the cavity transmission spectra). If the cavity frequency is ω _c = ω ₀ = 1200 cm⁻¹, 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 ≈ N e − β ℏ ω 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² φ⟩ = 1/3.

3.4 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 ω ₀ (for example, Eq. (33)), and a broadening function (assumed to be a Gaussian), expressed as follows [47]

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(ω ₀) 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 q ̂ c λ ⋅ μ ( R ̂ ) = ω c λ q ̂ c 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 R / P P denotes the ratio of equilibrium population between the reactant and product, see Section VII of the Supplementary Material. The time derivative P ̇ R ( 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 ₀ (outside the cavity rate constant) obtained from the HEOM simulation and report k/k ₀ = 1 + k _VSC/k ₀.

We report the numerical value of k/k ₀ as a function of the cavity frequency ω _c. For the rate constant predicted by FGR, we only report the value of k/k ₀ = 1 + k _VSC/k ₀ (see Eq. (10)), where k _VSC is evaluated using Eq. (36), and the variance defined in Eq. (37b) is estimated as σ = 30.74 cm⁻¹ for the model parameters we used. See Supplementary Material, Section VII for details. And we directly use the numerical result of k ₀ 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 = ω ₀, 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 < ω ₀ where ω _c = 1000 cm⁻¹, and there is no apparent η _c dependence of k(t), indicating that the coupling to the cavity has no effect. Figure 5(c) presents the value of k/k ₀ 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 = ω ₀ = 1190 cm⁻¹. 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.

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 = ω ₀ = 1190 cm⁻¹). (b) The flux-side correlation functions are calculated by HEOM but off-resonance (with ω _c = 1000 cm⁻¹). (c) The profile of the resonant VSC rate constant k/k ₀ 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.

Next, we explicitly consider going beyond the single-mode limit. For the 1D FP cavity, k VSC 1 D 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 C = 8 π ( c Δ k ‖ ) 2 Z eff , and τ ‖ ( ω ) = ω D / c ω 2 − ω c 2 (c.f. Eq. (28)). Note that this expression also peaks at ω _c = ω ₀ (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 ₀ using k _VSC expression in Eq. (40), corresponding to the case of many modes inside a 2D FP cavity. Here, we choose D / c = 3.33  fs, corresponding to D = 1 μ m. This should be viewed as the typical value of D , which is the effective lateral size of a given mode. Results obtained with a range of other choices of D are provided in the Supplementary Material, Section VIII, all of which show a sharp peak at ω _c ≈ ω ₀. Note that the broadening factor (Eq. (36)) was not included for k VSC 2 D 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 A ( ω ) , and the convergence is carefully checked. One can observe that the resonance peak is still centered around ω _c = ω ₀ with minor red-shift, which demonstrates the normal incidence condition. The resonance peak is asymmetric due to the asymmetry of A ( ω ) (see Figure 4(c–f)). Moreover, the rate profile tails toward the lower energy regions, which is the opposite of the trend in A ( ω ) (see Figure 4). Compared to the single mode version of the theory, considering many modes in a 2D FP cavity predicts that the “action spectrum” of the VSC-modified rate constant has an asymmetric behavior around ω _c = ω ₀, with a longer tail when ω _c < ω ₀. This is an interesting prediction from the current theory in Eq. (40). In recent VSC experiments by Simpkins [8], it does seem that the ω _c < ω ₀ side has a longer tail than the ω _c > ω ₀ side of the action spectrum (k _VSC vs ω _c plot, see Figure 3(a) of Ref. [8]). However, this seemingly asymmetrical rate constant profile in Ref. [8] could be caused by a lack of more experimental data points for a blue-tuned cavity (ω _c > ω ₀) due to the experimental difficulty of obtaining such measurements. More experimental data are required to definitively test this trend. Note that in the Simpkins experiment [8] the rate constant was resonantly suppressed. Recent quantum dynamics simulations [39] suggest that by resonantly coupling the cavity mode to a spectator mode (which in turn couples to the reaction coordinate), the rate constant can be suppressed by the cavity. Future work is needed to investigate such a resonance suppression effect.

Figure 6:

FGR rate profiles of k/k ₀ as a function of ω _c. (a) FGR rate profiles for the single mode case (or the many modes case inside a 1D FP cavity) calculated using Eq. (36) (same as the solid lines in Figure 5(c)). (b) FGR rate profiles k VSC 2 D for many mode cases inside a 2D FP cavity calculated using Eq. (40). Here, we use D / c = 3.33  fs, which corresponds to D = 1 μ m. Note that both k VSC 1 D and k VSC 2 D are rescaled by a factor of 0.4 to be consistent with Figure 5.