Cavity frequency-dependent theory for vibrational polariton chemistry

Recent experiments demonstrate the control of chemical reactivities by coupling molecules inside an optical microcavity. In contrast, transition state theory predicts no change of the reaction barrier height during this process. Here, we present a theoretical explanation of the cavity modification of the ground state reactivity in the vibrational strong coupling (VSC) regime in polariton chemistry. Our theoretical results suggest that the VSC kinetics modification is originated from the non-Markovian dynamics of the cavity radiation mode that couples to the molecule, leading to the dynamical caging effect of the reaction coordinate and the suppression of reaction rate constant for a specific range of photon frequency close to the barrier frequency. We use a simple analytical non-Markovian rate theory to describe a single molecular system coupled to a cavity mode. We demonstrate the accuracy of the rate theory by performing direct numerical calculations of the transmission coefficients with the same model of the molecule-cavity hybrid system. Our simulations and analytical theory provide a plausible explanation of the photon frequency dependent modification of the chemical reactivities in the VSC polariton chemistry.

T R +T r is the kinetic energy operator, whereT R andT r represent the kinetic energy operator for nuclei and for electrons, respectively, andV (x) is the potential operator that describes the Coulombic interactions among electrons and nuclei.
The cavity photon field Hamiltonian under the single-mode assumption is expressed aŝ where ω c is the frequency of the mode in the cavity,â † andâ are the photonic creation and annihilation operators, andq c = h/2ω c (â † +â) andp c = i hω c /2(â † −â) are the photonic coordinate and momentum operators, respectively. Choosing the Coulomb gauge, ∇ ·Â = 0, the vector potential becomes purely transverseÂ =Â ⊥ . Under the long-wavelength approximation, where A 0 = h/2ω c ε 0 Vê, with V as the quantization volume inside the cavity, ε 0 as the permittivity, andê is the unit vector of the field polarization.
We further introduce the Power-Zienau-Woolley (PZW) gauge transformation operator 1,2 aŝ The PZW transformation operator can also be expressed asÛ = exp − ī h 2ω c /hμA 0qc = exp − ī h ( j z jÂ x j ) . Recall that a momentum boost operatorÛ p = e − ī h p 0q displacesp by the amount of p 0 , such thatÛ pÔ (p)Û † p =Ô(p + p 0 ). Hence,Û is a boost operator for both the photonic momentump c by the amount of 2ω c /hμA 0 , as well as for the matter momentump j by the amount of z jÂ . UsingÛ † to boost the matter momentum, one can show that henceĤ C can be obtained 3 by a momentum boost with the amount of −z jÂ forp j , then addinĝ H ph .
The QED Hamiltonian under the dipole gauge (the "d · E" form 1,4 ) can be obtained by performing the PZW transformation onĤ C as followŝ where we have used Supplementary Eq. 5 to expressĤ C , and the last three terms of the above equation are the results ofÛĤ phÛ † . Usingq c andp c , one can instead show that because the PZW operator boosts the photonic momentump c by 2ω c /hμA 0 . The term ωc h (μA 0 ) 2 is commonly referred to as the dipole self-energy (DSE).
The Pauli-Fierz (PF) QED Hamiltonian 5-7 can be obtained by using a unitary transformation DenotingÛ φ = exp iφâ †â = e −Â (with φ = π 2 ), henceÂ = −iφâ †â . Using the BCH identity, we have Similarly, we have e −iφâ †ââ † e iφâ †â = e −iφâ † . Choosing φ = π 2 results inÛ † φâÛ φ → iâ andÛ † φâ †Û φ → −iâ † . Using these results, and applyingÛ φ onĤ D , we have the PF Hamiltonian as followŝ Note that we have used the fact thatÛ φĤMÛ † φ =Ĥ M , i.e.,Û φ does not contain any matter DOFs. Hence, the role ofÛ φ is to switchp c andq c , and for a photon field, they are inter-changeable due to the pure harmonic nature of the quantized field. The PF Hamiltonian has the advantage as a pure real Hamiltonian and the photonic DOF can be viewed 5,7 and computationally treated 8,9 as "nuclear coordinates". Projecting the above Hamiltonian in the ground electronic state of the molecule |Ψ g (which is obtained by solvingĤ el |Ψ g = E(R)|Ψ g , we obtain the model Hamiltonian, √h ωc ) 2 , in Eq. 1 of the main text, which is depicted in Fig. 1b of the main text.

Supplementary Note 2. Details of the Normal Mode Analysis.
Here we provide the detailed derivations of the normal mode frequencies used in the GH rate theory.

The polariton Hamiltonian is defined asĤ
where M = 1836 a.u. for the nuclear reaction coordinate (proton coordinate).
The equilibrium point (R 0 , q 0 ) at the reactant side of H pl is located at , H pl and E(R) share the same equilibrium position along the R direction. At R 0 , the curvature of the reactant well is ∂E 2 (R) and ω 0 is the bottom of the well frequency of the reactant. We further denote the molecular dipole and its derivative at R 0 as µ 0 ≡ µ(R 0 ) and µ 0 ≡ ∂µ(R) ∂R R 0 , respectively. At R = R 0 and q c = q 0 , the Hessian matrix element The other terms are also straight-forward to evaluate, resulting in where C 0 = 2ωc Mh χ · µ 0 . Note that making an approximation 10 of the dipole operator as µ(R) ≈ µ 0 + µ 0 (R − R 0 ) in H pl gives the same results of H, although it is not necessary to make such an approximation.
The normal mode frequencies Ω ± are obtained by diagonalizing Supplementary Eq. 17, resulting in Similarly, the saddle point (R ‡ , q ‡ ) of H pl is located at ∂R R ‡ = 0, i.e., H pl and E(R) share the same saddle point (as well as same minimum as shown before). At R ‡ , the curvature of the reaction barrier is and ω b is the barrier frequency. We further denote the molecular dipole and its derivative at R ‡ as µ ‡ ≡ µ(R ‡ ) and µ ‡ ≡ ∂µ(R) ∂R R ‡ , respectively. At R = R ‡ and q c = q ‡ , the The other terms are also straight-forward to evaluate, resulting in where C ‡ = 2ωc Mh χ · µ ‡ . The normal mode frequenciesω are obtained by diagonalizing Supplementary Eq. 22, resulting in

Supplementary Note 3. Model Molecular Hamiltonian
The potential energy surface (PES) and permanent dipole moment are taken from a Shin-Metiu model, 11 which is illustrated in Fig. 1 of the main text. The Shin-Metiu model is an one dimensional molecular system that describes a proton-coupled electron transfer reaction between a donor and an acceptor ion. The model consists of a transferring proton with a mass of m p and charge z p , an electron, and two fixed ions (a donor and an acceptor ion, with the charge of z D and z A , respectively). The molecular Hamiltonian isĤ M =P 2 2M +Ĥ el +Ĥ vib , where M is the mass of the nuclei (proton in this model),Ĥ el =T r +V eN +V NN is the electronic Hamiltonian, wherê T r =p 2 r /2m e represents the kinetic energy operator of the electron with mass m e ,V eN describes the interaction between the electron and the three nuclei, which is written as a modified Coulomb where r is the position of the electron and e = 1 a.u. is the fundamental charge, R is the position of the proton, while R D and R A are the positions of the donor and acceptor ion, respectively. R c is a parameter that controls the strength of the modified Coulomb potentials. The nucleus-nucleus interaction potential V NN that describes the Coulomb repulsion between the proton and the static ions takes the form of The parameters in the molecular HamiltonianĤ M used in this work is tabulated in Supplementary The vibrational HamiltonianĤ vib is modeled by the Caldeira-Leggett 12 system-bath Hamilto- describes the interaction between the solvent mode R s and a dissipative bath, where R k represents the k th bath mode with a conjugate momentum P k and a mass M k = M s . The coupling constant c k and the frequency ω k is characterized by an ohmic spectral density with a characteristic phonon frequency ω p and a friction constant ζ. The vibrational frequencies and the coupling coefficients are sampled based on the following expressions 13 where ω N = ωp N (1 − e −ωm/ωp ). The index k refers to the k th bath mode, with M k = M = 1836 a.u., and N represents the total number of total bath modes, and the maximum bath frequency is ω m .
In this work, we use N = 80 bath modes and we choose the characteristic phonon bath frequency ω p = ω 0 = 170.6 meV, the maximum bath frequency ω m = 5ω p , and the friction constant, ζ = 49.6 meV (400 cm −1 ).

Supplementary Note 4. Derivation of the κ GH from the Generalized Langevin Equation
Here, we provide an alternative derivation of the GH theory based upon the generalized Langevin equation [14][15][16] (GLE). We begin by expressing the total Hamiltonian in Eq. 1 of the main text as follows We further make the approximation of the permanent dipole by expanding it at the saddle point due to the fact that µ(R ‡ ) = 0 and R ‡ = 0 for the model system used in this work. In this work we obtain µ ‡ = −1.887 a.u. by fitting µ(R) in the vicinity of R ‡ (see Supplementary Fig. 2). In addition, we expanding the potential energy of the molecule at R = R ‡ as where we used the results that ∂E(R) ∂R R ‡ = 0 and R ‡ = 0.
To simply the derivation we use the mass-weighted coordinates and momentum, such that where in the second line we introduced and c k is defined in Supplementary Eq. 27-28. The Hamilton's equations of motion of To solve the above coupled equations, we perform the Laplace transform of function The Laplace transform ofẌ,Ẍ k andq c are given as Using the above expressions in Supplementary Eq. 36-37, and plugging them back into the Laplace transformed Eq.36-37, we haveX The Laplace transformed Supplementary Eq. 35 can be expressed as where we have used Supplementary Eq. 41-42, Supplementary Eq. 38, as well as the fact that We further define the Laplace transformed frictionsξ c (γ) (cavity photon friction) andξ p (γ) (phonon bath friction) as well as the random noiseF c (γ) (cavity photon noise) andF p (γ) (phonon bath noise) as followsξ Using the above definitions in Supplementary Eq. 43, we have We take the inverse Laplace transform of Supplementary Eq. 48 as follows where we denote the inverse Laplace transform as Each term in the above equation are inverse Laplace transformed as follows where in the second and third line we have used L −1 [γX(γ)](τ ) =Ẋ(t). Further, using the convolution theorem, we have Using Supplementary Eq. 49-56, we arrive at the following generalized Langevin equation (GLE) for the mass weighted coordinate X as follows where F c (t) and F p (t) are the total random noise associated with the cavity mode and phonon bath, respectively, each satisfies fluctuation-dissipation theorem as F c (0)F c (t) = M k B T ξ c (t) and In the original coordinate R, the GLE is expressed as For the simplified model system ofĤ −Ĥ vib , the GLE is expressed as Using the above generalized Langevin equation in Supplementary Eq. 58, we can obtain a rate constant expression (transmission coefficient), which is commonly referred to as the Grote- Using the fact that F p (t) = F c (t) = 0, taking the ensemble average of the above equation over the initial distribution R(0), and recognize that have the following Grote-Hynes equation for the frequency λ as follows whereξ c (γ) andξ p (γ) are expressed in Supplementary Eq. 44 and Supplementary Eq. 45, respectively. The transmission coefficient in the GH rate theory (by using the positive solution of λ in Supplementary Eq. 61) is given as The same result can also be derived using the normal-mode analysis of the Hamiltonian as shown Using Supplementary Eq. 63 and numerically solving Supplementary Eq. 61 gives the GH transmission coefficients of the model molecule-cavity hybrid system.
Alternatively, one can take the Markovian limit of the phonon bathĤ vib . Note that for the model system used in this work, J(ω) = ζωe −ω/ωp (Supplementary Eq. 27), we take the Markovian limit ω p → ∞ (ω p larger than all physical frequency of the system), which leads to The Laplace transformed kernel under the Markovian limit is theñ using the above Markovian limit, Supplementary Eq. 61 is simplified into following equation of λ which facilitate the procedure of solving λ numerically.
Finally, for the simplified model system withĤ −Ĥ vib , the GLE is expressed in Supplementary Eq. 59, and the GH equation of λ is which is the same as The analytical solution of the above equation is Taking the positive solution, the analytical κ GH has the following expression is the slope of the dipole moment on the dividing surface R ‡ . This is the result presented in Eq. 5 of the main text.
whereV eN (r, R) andV NN (R) are defined in Eq. 24 and Eq. 25, respectively. Note that when computing the adiabatic electronic energies, the nuclear position R is viewed as a parameter.
Further, r i |T r |r j is given analytically 18,19 as follows Directly diagonalizing the matrix of r i |Ĥ el |r j at a given nuclear position R in this grid basis gives the adiabatic electronic statesĤ el |Ψ g (R) = E(R)|Ψ g (R) . In this work, we only focus on the electronic ground state where c i (R) = r i |Ψ g (R) is the expansion coefficient obtained by diagonalizing the matrix of r i |Ĥ el |r j . The ground state permanent dipole momentμ = R +r is computed as The numerical results of E(R) and µ(R) are presented in Fig. 1b  Supplementary Fig. 3 presents the adiabatic electronic potential E α (R) (bold curves) by solvinĝ H el |Ψ α = E α (R)|Ψ α , as well as photon dressed potentials E α (R) + (n + 1 2 )hω c (thin curves). One can clearly see that the the ground adiabatic surface (black thick curve) is well separated from the first excited adiabatic electronic surface (blue thick curve) by at least 3.5 eV, and the photon dressed state |Ψ g , 30 is at the same energy of |Ψ S1 , 0 . However, they are not directly coupled to each other from the light matter interactions. Thus the presence of electronic excited states has a negligible influence on the polariton absorption spectrum and the reaction dynamics.
We emphasize that it is crucial to consider the dipole self-energy (DSE) term ωc h (μA 0 ) 2 in Supplementary Eq. 10. Ignoring this term 20 will cause an artificial change of the barrier height on the Cavity Born-Oppenheimer Surface defined as Supplementary Fig. 3(b), resulting in a large error of the rate constant.

b. Polariton Eigenspectrum and Absorption Spectrum. The vibrational polaritonic
Hamiltonian is defined as followŝ is the ground state adiabatic surface, µ(R) = Ψ g (R)|μ|Ψ g (R) is the ground state permanent dipole obtained from Eq. 75. In this Hamiltonian, we treat the nuclear DOF quantum mechanically. We numerically obtain the vibrational polariton eigenvalue E ν and eigenvector {|Φ ν } by solving the following eigenequation using the R-grid-Fock basis {|R j ⊗ |n ≡ |R j , n}, where {|R j } are the nuclear grid points, |{n } are vacuum Fock states, i.e., the eigenstates of the photonic HamiltonianĤ ph = (â †â + 1 2 )hω c . The matrix elements of the vibrational polaritonic HamiltonianĤ vpl is given as where the matrix elements of the nuclear kinetic energy operator R i |T r |R j = R i |P 2 2M |R j are evaluated using the Fourier grid method 18,19 as follows Diagonalizing to obtain the polariton eigenvalue E ν , as well as the polariton eigenvector where {b ν n,j (R) = R j , n|Φ ν } are the eigenvectors ofĤ vpl .
Supplementary Fig. 4 presents the Rabi splitting calculated by directly diagonalizingĤ vpl (black dots) compared to the simple analytical expression (red lines)hΩ ) momentum operator, respectively, whereâ andâ are the photon creation and annihilation operators. Under the dipole gauge, the matter interact with the quantized radiation mode of the cavity by displacing the photonic coordinate with the amount of characterizes the coupling strength between the molecule and the cavity (see Method). In this study, we have also explicitly assumed that the dipole moment is always aligned with the cavity polarization direction. Vibrational Polariton Rabi Splitting. At the equilibrium position of the reactant R 0 , one can approximate the permanent dipole as µ(R) ⇡ µ 0 + µ 0 0 (R R 0 ), where µ 0 = µ(R 0 ) and µ 0 0 = @µ(R) @R | R 0 . The light-matter interaction term inĤ (Eq. 1) at R 0 becomes [15,17] is the vibrational frequency at the equilibrium nuclear configuration R 0 , M is the e↵ective mass of the nuclear vibration,b † andb are the creation and annihilation operators for the nuclear vibration associated with the coordinate R. At the resonance condition of ! c = ! 0 , the photon-vibration coupling induces a Rabi splitting~⌦ R as follows [15,17] where the normalized coupling strength ⌘ p c = i p~! c /2(â † â) are the photonic coordinate and momentum operator, respectively, whereâ † andâ are the photon creation and annihilation operators. Under the dipole gauge, the matter interact with the quantized radiation mode of the cavity by displacing the photonic coordinate with the amount of characterizes the coupling strength between the molecule and the cavity (see Method). In this study, we have also explicitly assumed that the dipole moment is always aligned with the cavity polarization direction.
Vibrational Polariton Rabi Splitting. At the equilibrium position of the reactant R 0 , one can approximate the permanent dipole as µ(R) ⇡ µ 0 + µ 0 0 (R R 0 ), where µ 0 = µ(R 0 ) and µ 0 0 = @µ(R) @R | R 0 . The light-matter interaction term inĤ (Eq. 1) at R 0 becomes [15,17] is the vibrational frequency at the equilibrium nuclear configuration R 0 , M is the e↵ective mass of the nuclear vibration,b † andb are the creation and annihilation operators for the nuclear vibration associated with the coordinate R. At the resonance condition of ! c = ! 0 , the photon-vibration coupling induces a Rabi splitting~⌦ R as follows [15,17] where the normalized coupling strength ⌘ where the normalized coupling strength ⌘ To compute the absorption spectrum of the molecule-cavity hybrid system, we use a simple analytical scheme, 21 and a phenomenological width parameter ε = 21.95 cm −1 to account for the broadening of the absorption spectrum observed in the recent experiments. 22 The absorption cross section σ(E) as a function of excitation energy E is expressed 21,23 as follows where E ν is the ν th vibrational polaritonic eigenenergy obtained by solving Supplementary Eq. 77, and E 0 is ground vibrational polaritonic eigenenergy ofĤ vpl , and c is the speed of the light. The transition matrix element is computed as follows where µ(R) is expressed in Supplementary Eq. 75, i and j are indices for electronic grid points r i and nuclear grid points R j , respectively, and n is the index of the vacuum's Fock state. In addition, is the expansion coefficients of the adiabatic electronic states obtained by diagonalizinĝ H el (R j ) = k,l r l |Ĥ el (R j )|r k |r l r k | (matrix elements provided in Supplementary Eq. 72) at the particular nuclear configuration R j .
c. Parameterization of the Potential Energy Surface. To facilitate the numerical simulation of the reactive recrossing dynamics, the potential energy E(R) and the permanent dipole µ(R) are obtained from FGH quantum calculation are fitted to analytical functions. For E(R), we used a cosine series to fit the potential energy surface (PES) Similarly, we fitted µ(R) with the following sine series where N = 6 sine functions to fit the DM with the parameters detailed in Supplementary Table 84.
Supplementary Table 3: Parameters used to fit the ground-state dipole µ(R).
where the friction kernel and the random force for the phonon bath are To simplify our theoretical analysis of the transmission coefficient, We assume the Markovian approximation under the limit of ω p → ∞ (ω p is larger than all physical frequency of the system), leading to the friction kernel as follows Under this limit, the and the random force correlation function becomes and the memory kernel in GLE is evaluated as where the first equality is from the convolution theorem of Laplace transform, . The GLE in Supplementary Eq. 85 under the Markovian limit is thus expressed as the following Langevin equation where the random force F p (t) satisfies F p (0)F p (t) = M k B T ξ p (t) = 2k B T ζδ(t),the force acting on the nuclear reaction coordinate is , and the classical equation of motion for the photonic DOF is To solve the above coupled equations numerically, we have used a modified velocity Verlet integrator 24,25 for the molecular DOF R and a velocity Verlet integrator for the photonic DOF q c as follows where the function A(t) is 24,25 and ξ(t) and θ(t) are Gaussian random numbers generated at time t.
Alternatively, we also directly simulate the full Hamiltonian H(R, q c , {R k }) in Supplementary Eq. 29 with the explicit time-dependent propagation of the phonon bath {R k }. A simple velocity Verlet algorithm is used to solving the following equations of motioṅ The transmission coefficient is numerically calculated from the flux-side correlation function formalism 26-28 as follows