Resonant catalysis of thermally activated chemical reactions with vibrational polaritons

Interaction between light and matter results in new quantum states whose energetics can modify chemical kinetics. In the regime of ensemble vibrational strong coupling (VSC), a macroscopic number \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N$$\end{document}N of molecular transitions couple to each resonant cavity mode, yielding two hybrid light–matter (polariton) modes and a reservoir of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N-1$$\end{document}N−1 dark states whose chemical dynamics are essentially those of the bare molecules. This fact is seemingly in opposition to the recently reported modification of thermally activated ground electronic state reactions under VSC. Here we provide a VSC Marcus–Levich–Jortner electron transfer model that potentially addresses this paradox: although entropy favors the transit through dark-state channels, the chemical kinetics can be dictated by a few polaritonic channels with smaller activation energies. The effects of catalytic VSC are maximal at light–matter resonance, in agreement with experimental observations.

Here, we show the relation between reactant and product harmonic oscillator operators. Let us consider the vibrational Hamiltonians for the single-molecule reactant and product electronic states (we omit label (i) for simplicity hereafter),Ĥ where m is the reduced mass of the mode, ω A is the frequency of the mode in each electronic state (A = R, P), d P is the difference between nuclear equilibrium configurations, ∆E is the energy difference between the electronic states, andp andx are the momentum and position operators for the described mode; therefore, the harmonic oscillator potential energy surface for P is a displaced-distorted version of that for R. The creation and annihilation operators are defined in terms of position and momentum (d R = 0), conversely, the position-momentum representation is written in terms of the creation and annihilation operators aŝ whered P = 2m/ d P ; therefore, the reactant operators are written in terms of product ones aŝ These transformations can be written in terms of a squeezing and a displacement operator [1]: with actions given byŜ † P (r)â † PŜ P (r) =â † P cosh r −â P sinh r, for Supplementary Note 2.
Here, we introduce the initial and final many-body vibronic states. The rate to calculate corresponds to the stoichiometric process where N is the number of molecules in the product electronic state P, and M − N is the number of molecules in the reactant electronic state R, such that M is the total number of molecules in the reaction vessel. Assigning labels to each molecule, without loss of generality, the transformation of the N + 1-th molecule can be written in the form which reduces to The charge transfer is ruled by the adiabatic couplingĴ = J RP ; then, the matrix element that describes the process of our focus is with many-body vibronic states given by where |Φ Y X is an eigenfunction of a vibrational Hamiltonian of the form In Supplementary equation (19), we have used the notation introduced in equation (2), are the Hamiltonians of the upper/lower and k-th dark modes, respectively, all with creation and annihilation operators as defined in equation (3). Therefore, the matrix element corresponding to the transition becomes Here, we derive the tridimensional Franck-Condon factor in equation (12) of the main text. The non-vanishing overlaps between the vibrational ground state of the reactants and an arbitrary vibrational excitation with quantum numbers {v + , v − , v D } on the products can be written in terms of creation operators as These operators acting in the UP and LP can be written as linear combinations of the operators acting on the electromagnetic mode and the bright mode [equation (3)], i.e., The binomial theorem yields Since [â 0 ,â B(N ) ] = 0, the only non-vanishing terms are those with m = n = 0, otherwise the overlap in the photonic mode would be between non-displaced states with different excitations; therefore, Moreover, the creation operators acting on the bright and dark modes can be expressed as linear combinations of operators acting on the N -th molecule and the bright mode that excludes it [equation (6)], i.e., By expanding the binomials as before, and discarding the terms that excite the B(N − 1) mode, we arrive at Acting the creation operator on the N -th mode allows us to write (27) Therefore, the square of the Franck-Condon factor in equation (21) is Here, we discuss the integration of the rate law.

Chemical Master Equation.
The chemical master equation for the reaction in equation (7) is given by where Pr(n, t | m, 0) is the conditional probability to observe n molecules of the donor at time t given that there were m at t = 0, and a(n) = nk V SC R→P (n) is the propensity function [2]. . ( This probability density function can be used to determine the average number of donor molecules at a given time: Taking the time derivative of this average yields equation (18). However, for the number of molecules considered, M = 10 7 , this calculation becomes intractable; therefore, we resort to the strategy described in the Materials and Methods section of the main manuscript and corroborate its validity with the stochastic simulation algorithm [2].
Stochastic Simulation Algorithm (SSA). For the decomposition reaction in equation (18), we can define as the conditioned probability density function for the time of the next reaction (τ ) given that the number of donor molecules left is M − n at t. This function enables the construction of an exact numerical realization of the reaction with the following algorithm: 1. Initialize the system at N R (0) = M .
2. With the system in state N R (t) = M − n(t), evaluate a(N R ).
3. Generate a value for τ = − ln(r)/a(N R ), where R is a uniformly distributed random number.

5.
Register N R (t) as needed. Return to 2 or else end the simulation.
In Supplementary Table 1, we show the correlation (r 2 ) between the reaction times calculated according to the meanfield finite difference approach described in the manuscript and the reaction times corresponding to the same step in the reaction with populations obtained from the mean of 100 trajectories computed with the SSA algorithm. Since these correlations are very close to the unity, we conclude that the compared methods are numerically equivalent [3]. These observations are consistent with a recent study that shows that mean-field theories provide good descriptions for polaritonic systems involving a large number of molecules [4].