Driving chemical reactions with polariton condensates

When molecular transitions strongly couple to photon modes, they form hybrid light-matter modes called polaritons. Collective vibrational strong coupling is a promising avenue for control of chemistry, but this can be deterred by the large number of quasi-degenerate dark modes. The macroscopic occupation of a single polariton mode by excitations, as observed in Bose-Einstein condensation, offers promise for overcoming this issue. Here we theoretically investigate the effect of vibrational polariton condensation on the kinetics of electron transfer processes. Compared with excitation with infrared laser sources, the vibrational polariton condensate changes the reaction yield significantly at room temperature due to additional channels with reduced activation barriers resulting from the large accumulation of energy in the lower polariton, and the many modes available for energy redistribution during the reaction. Our results offer tantalizing opportunities to use condensates for driving chemical reactions, kinetically bypassing usual constraints of fast intramolecular vibrational redistribution in condensed phase.

L ight and matter couple strongly when a large number of molecules are placed within optical cavities that confine light [1][2][3] . As a result, hybrid light-matter excitations called polaritons form when a collective molecular transition and a photon mode coherently exchange energy faster than the individual decay from each component. Light-matter strong coupling (SC) opens up a new path to modify material properties by controlling their electromagnetic environment 4,5 . For instance, vibrational strong coupling (VSC), where an infrared cavity mode couples to an ensemble of localized molecular vibrations in a film or solution, influences chemical reactivity even without external pumping 6,7 . However, the microscopic mechanism for modification of molecular processes through hybridization with light is still poorly understood [8][9][10] , since it could be limited by the presence of a large number of quasi-degenerate dark modes that do not possess any photonic character and are likely to behave similarly to uncoupled molecules 10 .
A Bose-Einstein condensate of polaritons 11 offers a solution to this problem since the macroscopic occupation of polaritonic states enhances the effects from SC. In the last decade, Bose-Einstein condensation has been demonstrated in several organic exciton-polariton systems at room temperature [12][13][14][15] . Recently, organic polariton condensates were used to build polariton transistors 16 , and theoretical predictions suggest they may also modify incoherent charge transport 17 . Interestingly, the consequences of polariton condensation on chemical reactivity have not been addressed in the literature prior to the present study.
Ideas of using Bose-Einstein condensates of molecules in chemistry have been previously proposed, but they require ultracold temperatures due to the large mass of the condensing entities 18,19 . The low effective mass that polaritons inherit from their photonic component, along with the large binding energy of Frenkel excitons, enables room-temperature condensation 20 . The partly photonic character of polaritons also offers additional benefits such as delocalization and remote action for manipulating chemistry 21 .
Here, we investigate the effect of polariton condensation on electron transfer (ET) (Fig. 1). While the reaction yield under infrared laser excitation, without SC, already differs from that under thermal equilibrium conditions 22,23 , polariton condensation amplifies this difference by changing the activation barrier for the forward and backward reactions unevenly, tilting the equilibrium towards either reactants or products.

Results
Bose-Einstein condensation of vibrational polaritons. Bose-Einstein condensation of vibrational polaritons has not been experimentally achieved yet; however, as we shall argue, there are compelling reasons to believe that they are presently within reach. Most theoretical investigations on polariton condensation with organic microcavities involve systems under electronic strong coupling (ESC) 24,25 ; polariton condensation under VSC requires a separate treatment due to the difference in energy scales and the involved relaxation pathways 26 . While typical bare exciton energies range from 2-3 eV with Rabi split-ting~200 meV under ESC, the bare frequency of vibrations is 100-300 meV with Rabi splitting~20−40 meV under VSC. Since the Rabi splitting is of the order of k B T at room temperature, thermal effects are crucial for vibrational polaritons. Under ESC, polariton relaxation is assisted by high-frequency intramolecular vibrations 27 , whereas, under VSC, low-frequency solvent modes play a key role in this process 28,29 , similar to what happens in THz phonon Fröhlich condensation in biomolecules 30,31 .
We model the polariton system as a set of N vibrational modes (â vib;j ), with frequency ω vib , strongly coupled to a single photon mode (â ph ) with frequency ω ph . In the Hamiltonian of the system, we have applied the rotating wave approximation. Upon diagonalizing this Hamiltonian, we get normal modes: lower and upper polaritons, and N − 1 dark mode with frequencies ω − , ω + , and ω k D , respectively: where Ω ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi is the Rabi splitting and Δ = ω ph − ω vib the detuning between cavity and molecular vibrations. To model polariton population dynamics, we use Boltzmann rate equations where the polariton system is weakly coupled to a low-frequency solvent bath, which enables scattering between modes 32 . These rate equations also account for final-state stimulation, where n i is the population, γ i is the decay rate and P i is the external pumping rate of the i th mode. The scattering rate from mode j to i, W ij , satisfies detailed balance: W ij =W ji ¼ e Àβðϵ i Àϵ j Þ .
Here, β = 1/(k B T), k B is the Boltzmann constant, T is the temperature, and ϵ i = ℏω i where ω i is the frequency of the i th mode. The decay from different modes is γ i ¼ jc i vib j 2 Γ # þ jc i ph j 2 κ, where jc i vib j 2 and jc i ph j 2 are the molecular and photon fraction, respectively, Γ ↓ is the decay rate of the molecular vibrations, and κ is the cavity leakage rate. Two factors play a determining role in the condensation threshold: (i) the rate of scattering between polariton and dark modes relative to losses from the system, i.e., the rate of thermalization and (ii) the abundance of modes close in energy to the condensing mode 33 . For all calculations, we assume fast thermalization. As mentioned in (ii), the presence of many modes close to the lower polariton would deter condensation by distributing the energy pumped into the system among all these modes. Thus, the energetic proximity between the dark state manifold, which has a large density of states (DOS), and the lower polariton pose one of the biggest challenges for polariton condensation under VSC.
The distribution of excitations between the polariton and dark modes is shown in Fig. 2 for different detunings and we observe a condensation transition at ℏΔ ≈ −1.5k B T (see Supplementary Note 2 for details). Above the condensation threshold, a large fraction of excitations resides in the lower polariton Àβ_ðΩÀΔÞ=2 . The average population per molecule at the condensation threshold n ¼ P th =NΓ # is a good measure of the feasibility of vibrational polariton condensation. For instance, demanding population inversion, n > 0:5, would be experimentally difficult to achieve in general. In Fig. 3, we plot n for different light-matter coupling strengths, 2_g ffiffiffiffi N p , and detunings, ℏΔ. Here, we numerically obtain P th as the pumping rate when 10% of the excitations are in the lower polariton. The threshold obtained this way closely corresponds with the theoretical condition for condensation n k D > n solvent where, n k D ¼ 1 NÀ1 ∑ N k¼2 n k D is the average occupation of a dark mode, and n solvent (E) is the Bose-Einstein population of a solvent mode with energy E at room temperature T room 33 . The energy difference between the lower polariton and the dark state reservoir ℏ(Ω − Δ)/2 determines the condensation threshold. From Fig. 3 we see that vibrational polariton condensation is feasible for water even at room temperature for up to zero detuning.
Our model does not include disorder; as a result, all dark modes are degenerate at frequency ω vib , but in experimental systems, inhomogeneous broadening of transitions can lead to a nonzero density of dark states even at the bottom of the lower polariton branch 34 . This fact will affect the condensation threshold and should be considered in the future while looking for experimental systems that can demonstrate vibrational polariton condensation. Stimulating the lower polariton directly by shining a resonant laser on it 16 or using a Raman scattering scheme 35 can help overcome this issue by dynamically lowering the condensation threshold.
Our system consists of N molecules placed inside an optical cavity supporting a single photon mode with a bosonic operator a ph and frequency ω ph . The molecules can be in the reactant R or product P electronic state; for the i th molecule, these states are denoted by R i and P i , respectively. Each electronic state is dressed with a high-frequency intramolecular vibrational mode with a bosonic operatorâ x;i and frequency ω vib where x = R, P; this mode couples to the photon mode. The equilibrium geometry of this vibrational mode depends on the electronic state  ffiffiffi ffi N p =k B T < 2), the threshold for condensation is high, n ) 1, and polariton condensation is difficult to achieve experimentally. The threshold for condensation is much lower, n < 0:1, for the lighter colored (yellow, orange) regions. In our plot above only the upper polariton is pumped and we use cavity leakage rate Þ and S is the Huang-Rhys factor. Apart from the intramolecular vibrations, an effective lowfrequency solvent mode surrounding each molecule facilitates ET. It is treated classically, with q S,i and p S,i being its position and momentum.
The HamiltonianĤ for the full system is a generalization of Eq. (1) to account for the chemical reaction, whereĤ 0 describes the photon (Ĥ ph ), intramolecular vibrations and solvent modes of the ith molecule (Ĥ x;i ), and light-matter couplings (V x;i ). The diabatic couplingV react is a perturbation that couples R and P electronic states with coupling strength J RP .
We have taken the dipole moment to be zero when the vibrational coordinate is in its equilibrium position in both R and P electronic states. Relaxing this assumption will add terms of the form c x ðâ ph þâ ph Þ x i x i to the Hamiltonian and will not affect the reaction rates calculated.
where ΔG is the free-energy difference of each individual molecule reaction. We construct potential energy surfaces (PES) by parametrically diagonalizingĤ 0 as a function of the solvent coordinate q S,i . The with H 0 and correspond to the number of R and P molecules, respectively. While dynamics underĤ 0 conserves N R , N P , the effect ofV react is to induce reactive transitions that modify those quantities while keeping N = N R + N P constant. We assign the label 1 ≤ i ≤ N R to R molecules, and N R + 1 ≤ i ≤ N to P molecules. We also reorganize the intramolecular vibrations into a single bright mode, that possesses the correct symmetry to couple with light and N − 1 dark modes (D k ), labeled by an additional index 2 ≤ k ≤ N, which do not couple with light. The dark modes are orthogonal to the bright mode Unless mentioned otherwise, the number of R and P molecules is N R and N P , respectively, and for brevity, we will drop (N R , N P ) dependence in the subscript. The bright and photon modes combine to form the upper polariton (UP)â þ , and lower polariton (LP)â À , modes: with mixing angle, where Ω ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi is the Rabi splitting, and Δ = ω ph − ω vib the detuning between the cavity and molecular vibrations. The eigenstates ofĤ 0 are the dark, upper, and lower polariton modes with frequencies given in Eq. (2).
According to the MLJ theory, the rate constant for ET outside of an optical cavity depends on properties of the intramolecular and solvent modes [40][41][42] . Under laser driving, this rate constant is, where Here, P n ðnÞ is the Poisson distribution with average mode population n, λ S is the solvent reorganization energy, E z f is the activation energy, and j njn þ f 0 j 2 is the Franck-Condon (FC) factor, where n j i and n þ f 0 are the intramolecular initial and final states, respectively. P n ðnÞ has been taken to correspond to the ideal laser-driven-damped harmonic oscillator, leading to a coherent state in the vibrational mode. The presence of anharmonic couplings would lead to intramolecular vibrational energy redistribution (IVR) 43 , reducing the value of P n ðnÞ for high-lying Fock states. However, as we shall see below, even under these ideal circumstances, the condensate can outcompete the laser-driven situation in terms of reactivity. We thus expect the benefits of the condensate to be enhanced when IVR processes are taken into account.
Apart from vibrations within the reacting molecule, under VSC, the ET process also depends on vibrations in all other molecules and the photon mode, and can be represented by, Here and hereafter, the primed and unprimed quantities refer to electronic states with (N R , N P ) and (N R − 1, N P + 1) reactantproduct distributions, respectively. The symmetry of the lightmatter coupling allows us to use the dark state basis introduced in 44 and 38 to reduce the number of modes involved in the reaction from N + 1 to three, Here, the c th molecule is reacting, while D x,c and D 0 x;c are dark modes highly localized in it, with corresponding operatorsâ ðR;cÞ D andâ ðP;cÞ 0 D (see Supplementary Note 1). We perform all our calculations in this subsection using parameters from point A in Fig. 3 but while pumping the lower ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-022-29290-9 polariton. Here, ℏΔ = −0.5k B T, 2_g ffiffiffiffi N p ¼ 0:72k B T, k B T = 0.1389ℏω vib (T = 298 K when ℏω vib = 185 meV) and N = 10 7 ; we choose pumping rate P − = 0.08NΓ ↓ , which leads to average mode populations N + = 0.064, N − = 1.94 × 10 4 and N D = 0.078 under symmetric coupling g R = g P = g. Here, 2.4% of all excitations reside in the lower polariton. To compare the reaction rates under polariton condensation and outside the cavity under pumping, we take n ¼ 0:08 in Eq. (13). Under condensation, the initial vibrational state of the system can be described by ρ ¼ ∑ n þ ;n À ;n D Pðn þ ; n À ; n D Þ n þ ; n À ; n D n þ ; n À ; n D , where the entries in n þ ; n À ; n D label number of quanta in the UP, LP, and D R,c modes, respectively. The results from the rate equations Eq. (3) provide us only with the average steady-state mode populations, N + , N − , and N D , and not the distribution P(n + , n − , n D ). For simplicity, we assume the semiclassical approximation Pðn þ ; n À ; n D Þ % δ n þ ;0 P th N D ðn D Þδ n À ;N À , where P th N D ðnÞ is the thermal distribution with average population N D . This approximation is reasonable for populations The product vibrational states are ν þ ; ν À ; ν D 0 .
We assume that cavity leakage and the rate of scattering between modes is much faster than the rate of the chemical reaction. For a cavity with~100 ps lifetime and ET reactions with 1/k R→P~1 0 6 − 10 2 ps 45 , this assumption is valid. Therefore, if the populations of polariton modes change during the course of the reaction, they quickly reach a steady state before the next molecule reacts. Similarly, we also assume that the polariton and dark mode populations reach a steady state before the backward reaction takes place while computing the rate constant k cond P!R . Generalizing the cavity MLJ theory presented in ref. 38 , we calculate the rate constant for the forward reaction under polariton condensation, where The FC factor jF f ;n ν þ ;ν À ;ν D j 2 ¼ 0; N À ; njν þ ; ν À ; ν D 0 2 , and activation energy E f ;nz ν þ ;ν À ;ν D play an important role in determining the rate constant.
While many methods have been developed for computing multimode FC factors [46][47][48] , the focus has been on increasing the number of modes while keeping their occupation small. The current problem, however, offers a new technical challenge: the large occupation of LP makes the aforementioned methods computationally expensive. Instead, we draw inspiration from previous work that employs generating functions 47 and combine those techniques with the powerful Lagrange-Bürmann formula 49 to obtain analytical expressions for the required three-dimensional FC factors (see details in Supplementary Section S3).
The activation energies for the various channels of reactivity are, When condensation takes place, the number of quanta in the lower polariton N −~1 0 5 is so large that the summation in k cond R!P ðnÞ becomes difficult to estimate. To simplify the computation and gain intuition, we group channels into sets with the same change in the total number of intramolecular vibrational quanta f = ν + + ν − + ν D − N − − n upon ET. The closeness in energy between PES with the same f, and hence similar activation barriers, is the rationale for this grouping. k cond R!P ðnÞ then goes from a free summation over three indices ν + , ν − and ν D into a summation over four indices f, ν + , ν − and ν D with the constraint To understand the qualitative difference between reactions under polariton condensation and external pumping without SC, in Fig. 4a, b we plot the PESs (not to scale) showing the forward reaction under symmetric light-matter coupling and zero detuning. The yellow (black) parabolas in Fig. 4a, b represent PESs for a molecule in the electronic state R j i ( P j i) and vibrational state 2 j i ( 2 þ f 0 ) in Fig. 4a and 0; N À ; 2 ( 0; N À ; 2 þ f 0 ) in Fig. 4b. The red parabolas in Fig. 4b are additional final PESs provided by the vibrational polariton condensate (hereafter referred to solely as "condensate") that account for all other final vibrational states ν þ ; ν À ; ν D 0 .
The net rate of ET is, where k z R!P and k z P!R (z = IR, cond) are the rate constants for the forward and backward reactions, respectively, which are themselves functions of N R and N P when g R ≠ g P . We find the steady-state solution N ss R from this equation and compute the reaction yield N ss P =N. The difference in yield between the condensate and bare case is particularly large when λ S ≪ ℏω vib < |ΔG| (see Fig. 4c, d for symmetric coupling g R = g P ). To understand the underlying reason, we define the dominant channel f min as the one with minimum activation barrier outside of the cavity.
Setting the derivative in Eq. (24) equal to zero and taking into account the discrete nature of f, we find the dominant channel, . When λ S ≪ ℏω vib , |ΔG|, this channel contributes most to the rate constant because  24)) will be in the normal regime when λ S ≪ ℏω vib , |ΔG|.
In Fig. 4d, we see periodic yield modification in ΔG with period~ℏω vib , which decays for large ΔG/ℏω vib due to concomitant decline in FC factor for large changes in the number of vibrational quanta between the initial and final states. Outside of the cavity, we only see the first fringe (Fig. 4c). To observe the full periodic structure in yield, we need a large occupation of higher vibrational states which requires very large pumping rates outside of the cavity. However, under polariton condensation, the macroscopic population of the lower polariton enables these interesting periodic features to be observed at room temperature with experimentally attainable pumping rates. Additionally, polariton condensation not only modifies the reaction yield under symmetric light-matter coupling strengths, as seen in Fig. 4, it also changes the yield when the reactant and product asymmetrically couple with light (Fig. 5). Potential energy surfaces (not to scale) and reaction yield. a, c, e are results for a laser-driven system without light-matter strong coupling (SC) and b, d, f are for the same system under SC and 2.4% of the population in the lower polariton (condensation). All these plots are for symmetric lightmatter coupling g R = g P . a, b For a clearer qualitative picture, we plot the PESs under zero detuning Δ = 0. Initial (yellow) and final (black) PESs for a molecule undergoing the forward reaction with solvent coordinate q S,c . While the energy separation between black PESs is ℏω vib , the condensate provides many additional final PESs (red, separated by ℏΩ/2 at resonance). c Reaction yield N ss P =N at temperature k B T = 0.1389ℏω vib (T = 298 K when ℏω vib = 185 meV), Huang-Rhys factor S = 3.5, and average occupation of the intramolecular vibrational mode n ¼ 0:08. d Reaction yield N ss 08NΓ ↓ , N = 10 7 , temperature and Huang-Rhys factor are the same as c. The contributions of the red PESs through the condensate provide a broader tunability of reaction yields with respect to ΔG than under laser driving without SC. Notice that originally endergonic (exergonic) reactions in the absence of optical pumping can become exergonic (endergonic) under the featured nonequilibrium conditions. e, f A cross-section of plot (c, d) when λ S = 10 −2 ℏω vib . The pink shaded regions correspond to cases where the dominant forward (backward) channel is in the inverted (normal) regime; the opposite is true for the green shaded regions. The condensate amplifies the forward (backward) reaction in the pink (green) shaded regions.
Condensation provides many additional channels for the forward and backward reactions (separated by~ℏΩ/2, see red curves in Fig. 4b, showing only the forward channels at resonance Δ = 0) due to transfer of quanta between the polariton and dark modes during the reaction. Condensation speeds up the reaction when the dominant channel is either in the inverted or in the normal regime (Fig. 6b). This is the case because there are additional channels with energy both higher (benefiting the inverted regime) and lower (benefiting the normal regime) than the dominant channel (see red curves in Fig. 4b). Apart from reduced activation energy, the additional channels provided by the condensate also have large enough FC factors to affect the rate constant. Reaction channels that involve changes in the number of quanta in the LP during the reaction have significantly larger FC factors (~10 20 times) under condensation N − = 0.1N than without any pumping N − = 0 (Fig. S1 in the supplementary information compares them). Changes in the rate constant as a function of λ S (Fig. 6a) and ΔG (Fig. 6b) are large for small λ S / ℏω vib and when ΔG/ℏω vib = n/2 where n is an odd integer since activation energy effects are large for these set of parameters.

Discussion
Our result is a first step towards understanding the effect of Bose-Einstein condensation of polaritons on chemical reactivity. We demonstrate this effect using a simple electron transfer model (MLJ) with molecular vibrations strongly coupled to light. In particular, we show that one can counteract the massive degeneracy of dark modes and enhance polaritonic effects by having a macroscopic occupation of the lower polariton mode i.e., Bose-Einstein condensation. Our results indicate that the latter is feasible for experimentally realizable pump powers and Rabi splittings, despite the close proximity in the energy of the dark state manifold with ℏΩ~k B T. These results can guide the choice of suitable materials for condensation under VSC. While laser driving without SC modifies the reaction yield, this change is amplified by the condensate, due to the availability of many additional reactive channels that differ in energy by~ℏΩ/2 rather than~ℏω vib . For a wide range of parameters, we find that this leads to a periodic dependence of reaction yield as a function of ΔG (with period~ℏω vib ), rendering a set of originally endergonic reactions exergonic, and vice versa. These effects are substantially weaker under laser driving and highlight both the energetic (availability of additional channels with lower activation energy) and entropic (redistribution of vibrational energy from the condensate into the polariton and dark modes upon reaction) advantages of exploiting polariton condensates for reactivity. To summarize, vibrational polariton condensation offers a novel strategy to accumulate energy into a well-defined normal mode, a holy-grail in the field of vibrational dynamics that has been historically hindered by IVR. Its successful demonstration could revive hopes of "mode selective chemistry" 50 , beyond electron transfer processes. In future work, it will be interesting to explore how the studied phenomena generalize to molecular polariton condensates in different spectral ranges. b a d c Fig. 5 Reaction yield for asymmetric light-matter coupling. a, c The yield of the reaction when only the product (reactant) weakly couples with light. b, d Analogous plots under strong coupling 2g P ffiffiffi ffi N p ¼ 0:1ω vib , g R = 0 (g P = 0, 2g R ffiffiffi ffi N p ¼ 0:1ω vib ). We use parameters Δ = −0.0695ω vib , k B T = 0.1389ℏω vib (T = 298K when ℏω vib = 185 meV), S = 3.5, P − = 0.08NΓ ↓ and N = 10 7 . We assume the same scattering parameters W ij and decay rates Γ ↓ , κ as in

Data availability
Datasets generated by our code are available by email upon request to the authors.

Code availability
Computational scripts used to generate the plots in the present article are available by email upon request to the authors.
Received: 29 October 2021; Accepted: 9 February 2022; Fig. 6 Rate constant. Ratio of the rate constants inside k cond R!P and outside k IR R!P of the cavity under laser excitation with Δ = −0.0695ω vib , k B T = 0.1389ℏω vib (T = 298 K when ℏω vib = 185 meV), S = 3.5, 2g ffiffiffi ffi N p ¼ 0:1ω vib , P − = 0.08NΓ ↓ and N = 10 7 for symmetric coupling g R = g P = g. a Relative rate constant as a function of reorganization energy, λ S , with ΔG = −3.33ℏω vib (the pink curve; here, the dominant channel lies in the inverted regime) and ΔG = −3.73ℏω vib (the green curve; here, the dominant channel lies in the normal regime) and b as a function of ΔG with λ S = 10 −2 ℏω vib . Here, the pink region corresponds to the dominant channel in the inverted regime and the green to the normal regime.