Thermodynamic coupling of reactions via few-molecule vibrational polaritons

Our theoretical model considers a single molecule placed inside a UV-vis cavity, such that the cavity frequency is off-resonant to any optically allowed transitions of the molecule. In this regime, the coupling between the cavity and the molecule is purely parametric through the molecule’s polarizability, with the vibration of the molecule causing a dispersive shift in the cavity resonance [49]. Due to better spatial overlap between the mode profile of the cavity and the molecule, we consider using a Fabry-Perót cavity. However, the formalism presented here is valid for other cavity types. We show that the cavity-molecule system, when pumped with a laser that is red-detuned from the cavity and the detuning is in the order of the molecule’s vibrational frequency (Figure 1b), yields an effective Hamiltonian resembling the vibrational polaritonic Hamiltonian. Importantly, the light-matter coupling can be tuned by varying the laser power.

We model the photon mode and the vibration of the molecule as harmonic oscillators with frequencies $\omega_{\text{cav}}$ and $\Omega_{\text{v}}$, and annihilation operators $a$ and $b$, respectively. The cavity has a decay rate of $\kappa$ and the vibrational mode has a decay rate of $\gamma$. In this work, the losses will be modeled using Lindblad master equations [58]. Since the polarizability, $\alpha$, of the molecule, to the leading order, depends on its vibrational displacement [53], $x_{v}={x_{\text{zpf},v}}(b^{\dagger}+b)$, the Hamiltonian for the cavity-molecule system is [49]

where $g_{0}={x_{\text{zpf},v}}\bigg{(}\omega_{\text{cav}}\frac{\partial\alpha}{\partial x_{v}}\frac{1}{\epsilon_{0}V_{\text{c}}}\bigg{)}$ is the vacuum cavity-molecule coupling, with $V_{\text{c}}$ and $\epsilon_{0}$ as the mode volume of the cavity and vacuum permittivity, respectively. This Hamiltonian formally resembles an optomechanical setup [49], where the displacement of a mechanical oscillator modulates the frequency of the cavity. The cavity then acts back on the oscillator through radiation pressure force, which is a function of the cavity’s photon occupation.

We drive the cavity with a laser red-detuned $(\omega_{\text{L}}=\omega_{\text{cav}}-\Delta$, $\Delta>0)$ from the cavity resonance. The full Hamiltonian for a laser mode with annihilation operator $l$ coupled to the cavity-molecule subsystem is given as

with the cavity-laser coupling $J=\sqrt{\frac{\kappa}{\tau_{\text{rt}_{\text{L}}}}}$. Here, $\tau_{\text{rt}_{\text{L}}}$ is related to the laser power $P=\frac{n_{\text{L}}\hbar\omega_{\text{L}}}{\tau_{\text{rt}_{\text{L}}}}$ with $n_{\text{L}}=\langle l^{\dagger}l\rangle$ being the mean photon number in the laser mode [59].

We make the rotating wave approximation (RWA) in the laser-cavity coupling and diagonalize the laser-cavity subsystem. The operators $\tilde{l}$ and $\tilde{a}$ represent the new ‘laser-like’ and ‘cavity-like’ normal modes of the laser-cavity subsystem with frequencies $\tilde{\omega}_{\text{L}}$ and $\tilde{\omega}_{\text{cav}}$, respectively. The Hamiltonian after making the approximation and change of basis is

where $\varphi=\frac{1}{2}\tan^{-1}\bigg{(}\frac{2J}{\omega_{\text{cav}}-\omega_{\text{L}}}\bigg{)}$ is the laser-cavity mixing angle [60].

Considering the red-detuned case for the cavity-laser detuning, $\Delta$, in the order of ${O}(\Omega_{\text{v}})$, we drop the off-resonant contributions in the cavity-molecule interaction term, simplifying $H_{\text{full}}$ to

We will later set the laser-cavity detuning $\Delta=\Omega_{\text{v}}$.

We define a composite laser-cavity photon mode with annihilation operator $\mathcal{A}_{\text{ph}}=\frac{\tilde{l}^{\dagger}\tilde{a}}{\sqrt{(\tilde{n}_{\text{L}}-\tilde{n}_{a})}}$ (Figure 1b), where $\tilde{n}_{\text{L}}=\langle\tilde{l}^{\dagger}\tilde{l}\rangle$ and $\tilde{n}_{a}=\langle\tilde{a}^{\dagger}\tilde{a}\rangle$ are the mean photon occupations in the ‘laser-like’ and ‘cavity-like’ normal modes, respectively. The Heisenberg equations of motion for the operators $\mathcal{A}_{\text{ph}}$ and $b$ in the mean-field approximation [61],

yield the effective Hamiltonian

when ${n}_{\text{L}}\gg\tilde{n}_{a}$ and $J\ll\Delta$ (see Supplementary note-1). Here we have suggestively defined $\omega_{\text{ph}}\equiv\Delta$. For an input laser drive with power $P$, the interaction strength between the composite photon and the molecular vibration transforms to $g_{\text{eff}}=\frac{g_{0}}{\Delta}\sqrt{\frac{P\kappa}{\hbar\omega_{\text{L}}}}$, consistent with the results obtained from the classical treatments of the laser mode [54, 52, 39, 59].

When $\omega_{\text{ph}}=\Omega_{\text{v}}$, $H_{\text{eff}}$ resembles a vibrational polaritonic Hamiltonian, where the composite photon mode is resonant with the vibrational DOF (Figure 1b) [5]. Here the coupling strength is tunable by changing the pumping power of the laser. This can, in principle, foster the SC regime when the coupling strength supersedes the decay processes in the system. To look at parameter sets yielding this regime and to compute spectra, we simulate the dynamics of the density matrix ($\rho$) of the system using Lindblad master equations [53, 58] given as

The last four terms on the right-hand side are the Lindblad-Kossakowski terms defined as $\mathcal{L}_{\mathcal{O}}[\rho]=2\mathcal{O}\rho\mathcal{O}^{\dagger}-[\mathcal{O}^{\dagger}\mathcal{O},\rho]$. Here, $\mathcal{L}_{\mathcal{A}}$ models the incoherent decay from the composite photon mode. The incoherent decay, thermal pumping, and pure dephasing of the vibrational mode by the environment at temperature $T$ are modeled by the $\mathcal{L}_{b}$, $\mathcal{L}_{b^{\dagger}}$, and $\mathcal{L}_{b^{\dagger}b}$ terms, respectively, where $n^{\text{th}}_{v}=(e^{\hbar\Omega_{\text{v}}/k_{\text{B}}T}-1)^{-1}$ is the Bose-Einstein distribution function at transition energy $\hbar\Omega_{\text{v}}$. Additionally, in the limit of large photon number in the laser ($n_{\text{L}}\gg 1$), assuming the photon occupation to be constant, and thus the laser mode to be non-lossy, the decay rate of the composite photon equals the cavity decay rate $\kappa$ (see Supplementary note-2).
The simulations have been performed using the QuTip package [62, 63] and the results are presented in Figure 2 for the molecule Rhodamine 6G [49], where $\hbar\Omega_{\text{v}}=0.17$ eV, $\omega_{\text{ph}}=\Omega_{\text{v}}$ ($0.17$ eV), $\gamma=\gamma_{\text{r}}/2+\gamma_{\text{pd}}=0.01\Omega_{\text{v}}$ ($1.7\times 10^{-3}$ eV), $\kappa=0.02\Omega_{\text{v}}$ ($3.4\times 10^{-3}$ eV), $g_{0}=1.5\times 10^{-3}\Omega_{\text{v}}$ ($2.6\times 10^{-4}$ eV). Here, $\gamma_{\text{r}}=10^{-4}\Omega_{\text{v}}$ ($1.7\times 10^{-5}$ eV) and $\gamma_{\text{pd}}$ are the rates for vibrational relaxation and pure dephasing, respectively [6] . The fluence of the lasers is chosen to be below $\sim 10$ MW$/\text{cm}^{2}$ [49, 64] with a beam area of $A=5$ $\mu$m². Figure 2a shows the effective light-matter coupling, $g_{\text{eff}}$, as a function of laser fluence, $P/A$, and single photon coupling strength, $g_{0}$. In Figure 2b, the vibrational spectrum of the molecule $S_{b}(\omega)=\text{Re}[\int_{0}^{\infty}e^{-i\omega t}\langle x_{\text{v}}(t)x_{\text{v}}(0)\rangle_{\text{ss}}\text{d}t]$, splits, demonstrating SC. We see the Rabi-splitting increases with laser power, thus giving us additional control over the light-matter coupling strength. Figure 2c shows the spectra with one, two, and four molecules for constant laser power. Figure 2d is the vibrational spectrum of the molecule as a function of the cavity-laser detuning. The avoided crossing at the detuning ($\Delta$) equal to the vibrational frequency ($\Omega_{\text{v}}$) demonstrates maximal hybridization between the photonic and matter DOF. Finally, Figure  2e and f show the emission spectra from the cavity at steady state (ss), $S_{\mathcal{A}}(\omega)=\omega^{4}\cdot\text{Re}[\int_{0}^{\infty}e^{-i\omega t}\langle\mathcal{A}_{\text{ph}}^{\dagger}(t)\mathcal{A}_{\text{ph}}(0)\rangle_{\text{ss}}\text{d}t]$ [65, 53, 66], also revealing the polariton peaks.

The matter component of the polariton modes is delocalized over many molecules under collective SC [67, 5]. This delocalization can be exploited more effectively with few-molecule polaritons, owing to reduced involvement of dark modes [31, 9], which remain parked essentially at the same energy as the original molecular transitions. In this work, we consider the molecular species undergoing electron-transfer reactions, modeled using Marcus-Levich-Jortner (MLJ) theory [68, 69, 70]. Our system consists of two reactive molecules A and B of different species placed inside an optomechanical cavity. Here molecule A features a spontaneous reaction (with negative free energy change, $\Delta G_{\text{A}}<0$), while molecule B features an endergonic reaction ($\Delta G_{\text{B}}>0$, with $\Delta G_{\text{B}}>k_{\text{B}}T$). We demonstrate thermodynamic coupling between the two molecular species via the composite photon mode, such that the spontaneous electron transfer in A can drive B to react. Schematically, electron transfer in A creates a vibrationally hot product (Figure 3a), which, outside the cavity, just decays to the product ground state. However, inside the cavity, in the timescale of the Rabi frequency, this excitation can be captured by the photon mode, which then can excite the reactant in B to its vibrational excited state. The electron transfer in B can then proceed spontaneously from the reactant’s excited state (Figure 3b). Notably, this scheme can also be generalized to other types of reactions.

Within the framework of MLJ theory, the molecules can exist in either of the two diabatic electronic states: $|\text{R}_{i}\rangle$ corresponding to the reactant and $|\text{P}_{i}\rangle$ corresponding to the product for molecule $i\in\{\text{A},\text{B}\}$ (Figure 3a). For molecule B, in this case the switching between $|\text{R}_{\text{B}}\rangle$ and $|\text{P}_{\text{B}}\rangle$ through electron transfer contributes to useful mechanical work, manifested in changes of nuclear configuration [71, 72, 73]. The electronic states for each molecule are dressed with a local high-frequency intramolecular vibrational coordinate represented by annihilation operator $a_{x,i}$, for $x\in\{\text{R},\text{P}\}$, and coupled to a low-frequency effective solvent mode treated classically with rescaled momentum and position as $p_{\text{S},i}$ and $q_{\text{S},i}$. respectively. We assume that the high-frequency modes of both species, being resonant, are the only ones that couple to the composite photon [74]. Upon reaction, the high-frequency mode undergoes a change in its equilibrium configuration according to, $a_{\text{R},i}=D^{\dagger}_{i}a_{\text{P},i}D_{i}$, where $D_{i}=\text{exp}[(a^{\dagger}_{\text{P},i}-a_{\text{P},i})\sqrt{S_{i}}]$ is the displacement operator, and $S_{i}$ is the Huang-Rhys factor [37]. The Hamiltonian describing the system is given as $H=H_{0}+V_{\text{react}}$, where

Here $H_{\text{ph}}=\hbar\omega_{\text{ph}}\big{(}\mathcal{A}_{\text{ph}}^{\dagger}\mathcal{A}_{\text{ph}}+\frac{1}{2}\big{)}$ is the bare Hamiltonian corresponding to the composite photon mode consisting of the laser and the cavity, $H_{x,i}$ represents the high-frequency mode and the solvent mode associated with molecule $x_{i}$,

with $d_{\text{S},i}$ and $\Omega_{\text{S},i}$ being the displacement and frequency along the solvent coordinate, respectively, and $\Delta G_{i}$, the free energy difference for the molecular species $i$. Additionally, $V_{x,i}=\hbar g_{x,i}(a_{x,i}\mathcal{A}_{\text{ph}}^{\dagger}+a^{\dagger}_{x,i}\mathcal{A}_{\text{ph}})$ is the effective coupling between the photonic and molecular DOF. For simplicity, we assume that the reaction involves a vibrational mode with nearly identical frequency and light-matter coupling strength for species A and B in both the reactant and product electronic states ($\Omega_{x,i}=\Omega_{y,j}\equiv\Omega_{\text{v}}$ and $g_{x,i}=g_{y,j}\equiv g$). Finally, the diabatic couplings between the electronic states $|\text{R}_{i}\rangle$ and $|\text{P}_{i}\rangle$ are given by $V_{\text{react}}$, where $J_{i}$ is the coupling strength.

We can solve $H_{0}$ parametrically as a function of the solvent coordinates to construct the potential energy surfaces(PES) (Figure  4a and b). Considering these diabatic couplings $V_{\text{react}}$ to be perturbative, $H_{0}$ can be diagonalized to obtain the two polariton modes, $a^{(\pm)}_{x_{\text{A}},y_{\text{B}}}$, with frequencies $\omega_{\pm}=\frac{1}{2}(\omega_{\text{ph}}+\Omega_{\text{v}}\pm\sqrt{(\omega_{\text{ph}}+\Omega_{\text{v}})^{2}+8g^{2}})$, and one dark mode, $a^{D}_{x_{\text{A}},y_{\text{B}}}$, with frequency $\Omega_{\text{v}}$, given as

such that, $c_{x,\text{A}}+c_{y,\text{B}}=0$ and $|c_{x,\text{A}}|^{2}+|c_{y,\text{B}}|^{2}=1$. Here, $\theta=\frac{1}{2}\tan^{-1}\bigg{(}\frac{2\sqrt{2}g}{\delta}\bigg{)}$ is the mixing angle, where $\delta=(\omega_{\text{ph}}-\Omega_{\text{v}})$ is the detuning between the composite photon mode and the molecule. Here we have chosen the composite photon mode to be resonant with the intramolecular vibration, i.e., $\delta=0$.

We now define multi-particle states $|\bm{\phi};\bm{\nu}_{\bm{\phi}}\rangle$ that span the Hilbert space of the system, where $|\bm{\phi}\rangle=|x_{\text{A}},y_{\text{B}}\rangle$, $x,y\in\{\text{R},\text{P}\}$ corresponds to the electronic DOF, and $|\bm{\nu}_{\bm{\phi}}\rangle=|\nu^{+}_{x_{\text{A}},y_{\text{B}}},\nu^{-}_{x_{\text{A}},y_{\text{B}}},\nu^{D}_{x_{\text{A}},y_{\text{B}}}\rangle$ to the cavity-vibrational mode of each electronic state $|\bm{\phi}\rangle$ [6]. To describe the reaction, we look at the population dynamics in the electronic states. The kinetic master equations governing the time evolution of the system are given as [6]

where $p_{(\bm{\phi};\bm{\nu}_{\bm{\phi}})}(t)$ represents the population in $|\bm{\phi};\bm{\nu}_{\bm{\phi}}\rangle$ state, and $k(\bm{\phi}^{\prime};\bm{\nu}^{\prime}_{\bm{\phi}^{\prime}}|\bm{\phi};\bm{\nu}_{\bm{\phi}})$ is the rate constant for population transfer from $|\bm{\phi};\bm{\nu}_{\bm{\phi}}\rangle$ to $|\bm{\phi}^{\prime};\bm{\nu}^{\prime}_{\bm{\phi}^{\prime}}\rangle$ due to processes like reactive transitions between the electronic states accompanied by solvent reorganization and decay through the cavity and vibrational DOF. The rate constant for the reactive transition at a temperature $T$ within the framework of MLJ theory is given as [74]

Here, $E_{\bm{\phi};\bm{\nu}_{\bm{\phi}}}=E_{x_{\text{A}}}+E_{y_{\text{B}}}+\hbar\big{[}\omega_{+}\big{(}\nu^{+}_{x_{\text{A}}y_{\text{B}}}+\frac{1}{2}\big{)}+\omega_{-}\big{(}\nu^{-}_{x_{\text{A}}y_{\text{B}}}+\frac{1}{2}\big{)}+\Omega_{\text{v}}\big{(}\nu^{D}_{x_{\text{A}}y_{\text{B}}}+\frac{1}{2}\big{)}\big{]}$ is the energy of the state $|\bm{\phi};\bm{\nu}_{\phi}\rangle=|x_{\text{A}},y_{\text{B}}\rangle\otimes|\nu^{+}_{x_{\text{A}},y_{\text{B}}},\nu^{-}_{x_{\text{A}},y_{\text{B}}},\nu^{D}_{x_{\text{A}},y_{\text{B}}}\rangle$, and $\lambda_{\text{S}}^{(\bm{\phi}\bm{\phi^{\prime}})}$ and $J_{\bm{\phi}\bm{\phi^{\prime}}}$ are the solvent reorganization energy and diabatic coupling, respectively, corresponding to the reacting species. Additionally, $\langle\bm{\nu^{\prime}}_{\bm{\phi}^{\prime}}|\bm{\nu}_{\bm{\phi}}\rangle$ represent the Franck-Condon factors for the hybrid photon-vibration states $|\bm{\nu^{\prime}}_{\bm{\phi}^{\prime}}\rangle$ and $|\bm{\nu}_{\bm{\phi}}\rangle$ corresponding to the electronic states $|\bm{\phi^{\prime}}\rangle$ and $|\bm{\phi}\rangle$, respectively. For the simulations in this work, the Franck-Condon factors have been computed numerically from eigenstates obtained using the standard discrete-variable representation (DVR) of Colbert and Miller [75].

The reactive transitions transfer populations across different electronic states, while the cavity and vibrational decays lead to dynamics within the same electronic state. In these simulations, since $k_{\text{B}}T\ll\hbar\Omega_{\text{v}}$, we restrict ourselves to the first excitation manifold in the photon-vibration DOF. With the bare vibrational decay rate for the intramolecular vibrations for molecule A and B as $\gamma_{\text{A}}$ and $\gamma_{\text{B}}$, respectively, and bare cavity decay rate as $\kappa$, we have

where $\bm{1}_{q,\bm{\phi}}$ represents a single excitation in the polaritons $(q=\pm$) or the dark ($q=D$) mode, and the $c_{qj}$’s correspond to the expansion coefficients of the excited eigenmode in terms of the cavity and vibrational modes [76].

Finally, the anharmonic couplings between the vibrational mode of interest and an other bath of low frequency modes leads to transitions between the polaritons and the dark mode, [6]

where $\Theta(\Omega)$ is the heavyside step function, $n^{\text{th}}(\Omega)$ is the Bose-Einstein distribution function at the transition energy $\hbar\Omega=\hbar(\Omega_{q^{\prime}}-\Omega_{q})$, and $\mathcal{J}(\Omega)$ is the spectral density of the low frequency modes. Assuming the spectral density to be Ohmic [77], we have $\mathcal{J}(\Omega)=\eta\Omega\exp[-(\Omega/\Omega_{\text{cut}})^{2}]$, where $\eta$ is a dimensionless parameter modeling the anharmonic system–bath interactions and $\Omega_{\text{cut}}$ is the cut-off frequency for the low-frequency modes.

The results of the simulations are presented in Figure 4. Here, $\hbar\Omega_{\text{R},i}=\hbar\Omega_{\text{P},i}=0.22$ eV (call $\hbar\Omega_{\text{v}}$), $g_{0}=2\times 10^{-3}\Omega_{\text{v}}$ ($4.4\times 10^{-4}$ eV), $P/A=6.4$ MW/cm², A=5 $\mu$m², $\gamma_{\text{A}}=\gamma_{\text{B}}=1\times 10^{-5}\Omega_{\text{v}}$ ($2.2\times 10^{-6}$ eV) [6], $\omega_{\text{cav}}=2.2$ eV, $\kappa=0.015\Omega_{\text{v}}$ ($3.3\times 10^{-3}$ eV), $\lambda_{\text{S}}^{\text{A}}=0.04\hbar\Omega_{\text{v}}$ ($8.8\times 10^{-3}$ eV), $\lambda_{\text{S}}^{\text{B}}=0.1\hbar\Omega_{\text{v}}$ ($2.2\times 10^{-2}$ eV), $\Delta G_{\text{A}}=-\hbar\Omega_{\text{v}}$ ($-0.22$ eV), $\Delta G_{\text{B}}=0.7\hbar\Omega_{\text{v}}$ ($0.15$ eV), $\hbar\Omega_{\text{cut}}=0.1\hbar\Omega_{\text{v}}$ ($2.2\times 10^{-2}$ eV), and $\eta=0.0001$ [6]. The decay rates have been chosen to be similar to those typically found in VSC experiments [78, 61, 79]. The diabatic couplings are chosen to be $\hbar J_{\text{A}}=\hbar J_{\text{B}}=0.005\hbar\Omega_{\text{v}}=(1.1\times 10^{-3}$ eV) and $T=298$ K. We start from the initial electronic state $|\text{R}_{\text{A}},\text{R}_{\text{B}};\bm{0}\rangle$. Independently, the reaction of molecule A is spontaneous due to its negative free energy change, $\Delta G_{\text{A}}<0$, while molecule B remains in its thermodynamically stable conformer $|R_{\text{B}};0\rangle$ (Figure 4a and b). This reflects the dynamics of the species outside of the cavity. Placing both the molecules inside the cavity couples the two reactions via the photonic mode enabling the spontaneity of the reaction of molecule A to thermodynamically ‘lift’ B to its unstable configuration $|\text{P}_{\text{B}};0\rangle$ producing mechanical work. However, after molecule A has fully reacted (change in nuclear configuration), inevitably B has to relax again to its stable configuration $|\text{R}_{\text{B}};0\rangle$, completing one cycle of the mechanical motion of B (Figure 3b). The maximum population obtained in $|\text{P}_{\text{B}}\rangle$ before molecule B relaxes back to $|\text{R}_{\text{B}}\rangle$ increases with the light-matter coupling strength ($g_{\text{eff}}$), tunable with the fluence of the driving laser (Figure 4d). For the cycle to be repeated, molecule A needs to be ‘recharged’ or ‘replaced’. To achieve this, we envision a flow setup, as schematically depicted in Figure 4c, that can circulate molecule A inside and out of the cavity. Continued circulation of the A molecules is essential for the molecular machine of B to be oscillating between reactant and product and producing mechanical work (Figure 4a and b). This phenomenon realizes a heat engine producing mechanical work in molecule B, using the (chemical) energy flow from a ‘source’ to a ‘sink’ bath in the form of molecule A [80].