On-the-fly cavity–molecular dynamics of vibrational polaritons

Abstract

In this work, we combine the density functional tight-binding (DFTB) approach with a light–matter Hamiltonian beyond the long-wavelength approximation to simulate vibrational polaritons formed by coupling molecular vibrations to confined radiation inside a Fabry–Pérot optical cavity. Here, we develop a parallelized propagation scheme with lightweight inter-CPU communication by exploiting the sparse nature of the light–matter interactions in the real space representation. We find that the computationally expensive Born charges required for our propagation can be replaced with the computationally inexpensive Mulliken charges to obtain qualitatively accurate linear spectra especially when the nonlinearity (arising from molecular vibrations) of the light–matter interaction term is not substantial. However, the same approach may not be suitable to be used for studying cavity modification of energy transport or chemical dynamics as this approximation leads to spurious heating of the light–matter hybrid system. We demonstrate the utility of this on-the-fly approach to compute angle resolved polaritonic spectra of water. We implement our approach as an open-source computational package, CavOTF, which is available on GitHub.

---

Arkajit Mandal

Arkajit completed his PhD under the supervision of Prof. Pengfei Huo at the University of Rochester, where he worked on simulating the quantum dynamics of molecules interacting with confined radiation. In 2021, Arkajit joined Columbia University as a postdoctoral research scientist, where he worked with Prof. David R. Reichman. His research at Columbia involved investigating the transport properties of exciton-polaritons formed inside optical cavities. Arkajit began his independent career at Texas A&M in August 2024, where his research group has been investigating cavity modified chemical and physical properties of molecules and materials.

---

1 Introduction

A series of experiments have suggested that ground-state chemical reactions can be modified by coupling molecular vibrations to quantized vacuum radiation inside an optical cavity.^1–9 In this light–matter coupling regime, namely, the vibrational strong coupling (VSC) regime, molecular vibrational modes hybridize with the confined photonic modes to form vibro-polaritons which have been suggested to modify chemical dynamics in the absence of external illumination.¹⁰ Such couplings can lead to either enhancement^3,11,12 or suppression^2,3,13 of bond-breaking processes (or produce no measurable change^14–16). These effects are resonant, emerging when the cavity mode is tuned to a specific molecular vibration, thereby opening up the possibility to achieve chemical selectivity by tuning cavity frequency.^{9,11–13,17–19} Thus, VSC is a fundamentally new and non-invasive approach to controlling chemistry, requiring no external pumping or chemical reagents; instead, it harnesses the ever-present vacuum fluctuations of confined radiation.

However, VSC remains far from a fully reliable or broadly applicable strategy for controlling ground-state chemical reactivity. This is because these remarkable experimental observations currently lack a clear connection to established chemical principles or molecular morphology, or even a semi-empirical framework to guide expectations does not yet exist. Significant challenges also persist in interpretation, particularly in disentangling genuine cavity-induced changes from other experimental factors and in the reproducibility of experimental results.^20–22 These limitations underscore the need for a clear microscopic understanding of cavity-modified ground-state chemical reactivity and the need for accurate modeling of light–matter interactions to unambiguously distinguish the role of cavity radiation in modifying ground-state chemical reactivity.^21–23

Recent theoretical works that study a simple single-molecule-single-cavity-mode model system find resonant modification of chemical kinetics under vibrational strong coupling.^24–28 In contrast, theoretical works that consider an ensemble of molecules coupling to one cavity mode predict little or no change in chemical reactivity under similar conditions.^8,23 This discrepancy raises a fundamental question about what essential physics may be missing in these current theoretical descriptions. Likely limitations include the use of the long-wavelength and single-mode cavity approximations, as well as simplified molecular models that neglect vibrational anharmonicity and bond dissociation. Recently, Li and co-workers have introduced mesoscale simulation approaches that enable atomistic modeling of molecular ensembles within optical cavities.^29–32 However, these methods typically rely on classical force fields, which are inherently nondissociative and thus unable to fully capture the chemical complexity, and reactive character intrinsic to real molecular systems.^28,33–47

In this work, we develop an on-the-fly molecular dynamics approach in which part of the Hamiltonian is evolved in real space, while the rest is propagated in reciprocal space. This strategy enables the development of a highly parallelized computational algorithm that requires only light inter-CPU communication. In our approach, the electronic structure of matter is solved using the self-consistent-charge density functional tight binding (DFTB-SCC) approach, with each CPU responsible for propagating a distinct portion of the macroscopic system. We find that dipole derivatives (Born charges), required for propagating the dynamics, fluctuate significantly during propagation and cannot be replaced with their time-independent values (i.e., partial charges), as is commonly assumed when using classical force fields, for achieving quantitative accuracy when propagating the cavity-coupled dynamics. Interestingly, we find that for computing qualitatively accurate linear spectra, the computationally cheaper Mulliken charges can be used in place of Born charges, especially when the nonlinearity in the light–matter interaction term is not substantial. However, such a replacement should not be made for analyzing cavity-modified energy transport properties or chemical reactivity as this replacement leads to a significant heating of the system, as expected.

Using this framework, we demonstrate proof-of-principle simulations of a system containing over 8000 atoms, coupled to an ensemble of cavity modes. We present angle-resolved infrared spectra of liquid water under vibrational strong coupling, targeting both the asymmetric stretch and bending modes. We anticipate that this computational tool will allow for in silico experiments providing new opportunities for harnessing confined vacuum radiations in optical cavities to modulate and catalyze chemical reactivity. Based on our computational framework, we develop a computational package, namely “CavOTF”,⁴⁸ which is publicly available on our GitHub repository.

2 Theory

To demonstrate the feasibility of our on-the-fly approach for simulating vibrational polaritons beyond the long-wavelength approximation, we adapt a modified Holstein–Tavis–Cummings Hamiltonian. However, we assume that the field is approximately uniform at the level of an individual matter cell which can be referred to as a local long-wavelength approximation.^31,49,50 Following prior works,^31,51,52 we consider a simple two-dimensional world extending in x and y, with cavity quantization in the y direction and the molecular system extending in the x direction. The cavity radiation modes follow the dispersion relationship , where is the in-place wavevector, with L_s being the length of the supercell extending in the in-place direction x (in this work, we use L_s = 8 × 10⁵ a.u.), L denotes the mirror separation and c is the speed of light in the cavity medium. This simplified setup is used to demonstrate the core concept, while the model can be readily extended to three-dimensional molecular systems and to include multiple polarization directions. To describe the formation and dynamics of vibrational polaritons, we use the following microscopic light–matter Hamiltonian, formulated in the dipole gauge beyond the long-wavelength approximation^51,53,54 and expressed in mass-weighted coordinates within atomic units (ħ = 1) as

Here, the operators â^†_k (â_k) denote the photon creation (annihilation) operator corresponding to the k-th cavity mode with frequency ω_k. Furthermore, Ĥ_M denotes the bare matter Hamiltonian, written as a sum over N independent cells located at R_n = n·a (a is the lattice constant), each containing the same number of molecules. Following recent work,³¹ we treat the cells as non-interacting and denote the total dipole operator of the cell n by M̂_n. The total matter Hamiltonian can be written as

where T̂_n,j is the nuclear kinetic energy operator for the jth atom in the nth box and V̂({r_n,j}) represents the potential energy surface of the molecular subsystem. The light–matter coupling enters through the projected dipole operator, M̂_n = ([small epsilon, Greek, circumflex]·[small mu, Greek, circumflex]_n), where [small mu, Greek, circumflex]_n is the total dipole operator of the nth box and [small epsilon, Greek, circumflex] = x̂ is the polarization direction of the cavity field which has been approximated to be independent of the in-plane wave-vector k since the polariton spectra are probed near k → 0.

For a Fabry–Pérot cavity, the spatial dependence of the cavity field at the cell positions R_n is characterized by the mode function

where g₀ is the light–matter coupling constant. In atomic units, it can be written as , where V is the quantization volume. The light–matter Hamiltonian in eqn (1) can be more explicitly written as

This Hamiltonian is schematically illustrated in Fig. 1a. According to eqn (4), all photonic modes interact with all molecular boxes. Although the full light–matter system described in eqn (4) can be propagated using a parallelized algorithm on a typical high-performance computing cluster (HPC), we pursue a real-reciprocal-space propagation scheme that allows for a straightforward (lightly communicating) hub-and-spoke parallelization using Python's socketserver.

Fig. 1 Schematic illustration of light–matter interactions and the real reciprocal used here for propagating vibro-polariton dynamics. Schematic illustration of light–matter interactions (a) in the standard multi-mode dipole gauge description and (b) in real space. (c) Schematic illustration of the algorithm used here and implemented in our computational package CavOTF.⁴⁸

Since the light–matter coupling term plays a role only around k → 0, we make the simplifying approximation . The photon creation and annihilation operators are then transformed into real-space field operators using .⁵¹ Note that . This transformation localizes the interaction in real space (illustrated in Fig. 1b) and the total Hamiltonian is then expressed as

where denotes the collective light–matter coupling strength, and is the real space photonic coordinate.

In this work, we propagate the dynamics of the nuclear and photonic degrees of freedom classically, i.e. {r̂_n,j, q̂_n} → {r_n,j, q_n}. The equations of motion for both the nuclear and photonic coordinates are derived from Hamilton's equations of motion and integrated in time using a velocity-Verlet-like scheme. The forces F^m_n,j and F^c_n,j acting on the nuclear and photonic degrees of freedom, respectively, are written as

where and . Here, we compute the term −∇_n,jV̂_n({r_n,j}) on-the-fly using the DFTB-SCC approach. We integrate the equation of motion of the photonic and nuclear degrees of freedoms using the following scheme that uses a server-client (hub-and-spoke) architecture (schematically illustrated in Fig. 1c):

Step 0. Initialization. The matter degrees of freedoms {r_n,j(0), p_n,j(0)} at t = 0 are obtained from a long (150 ps) NVT DFTB trajectory. The photonic degrees of freedom {q̃_k, p̃_k} are then initialized performing a constrained molecular dynamics simulation of the Hamiltonian, H_LM({r_n,j = r_n,j(0), p_n,j = p_n,j(0)}) using a simple velocity-Verlet algorithm where {r_n,j, p_n,j} are frozen.

Step 1. Real-space propagation on clients. We evaluate −∇_n,jV̂_n({r_n,j}), μ_n, dμ_n/dr_n,j (computed numerically and updated every three steps). Using these quantities, the nuclear degrees of freedom are fully evolved over a single time step δt, while the photonic momenta are only partially evolved, with their evolution occurring solely due to the localized light–matter coupling term. Specifically, we update the nuclear and photonic degrees of freedom as follows:

Step 2. Receive photonic coordinates from the client and perform Fourier transformation on the server. In this step, the server receives all photonic variables from all the client CPUs and evaluates the Fourier transformed reciprocal-space photonic variables using the following expression:

where we introduce the classical complex variables {a_k} (and their complex conjugates ) which are obtained as

Step 3. Photonic propagation in the reciprocal space on the server. The photonic variables {q̃_k, p̃_k} are propagated analytically using

Step 4. Perform inverse Fourier transformation on the server and broadcast photonic coordinates to clients. Time propagated k-space photonic variables {q̃_k, p̃_k} are inverse Fourier transformed to obtain real space photonic variables {q_n, p̃_n} using the following expressions:

where {a_n} (and their complex conjugates ) are obtained as

The updated values are then sent to all clients. Each client receives only its corresponding pair (q_n, p̃_n).

Step 5. Repeat. Continue from Step 1 for the next time step.

Computational details

Molecular dynamics simulations were performed using the density functional tight binding (DFTB) method with second-order self-consistent charge (SCC) correction. This semiempirical quantum mechanical approach provides a balance between accuracy and computational efficiency, enabling simulations of systems containing hundreds of atoms over picosecond timescales.^55–59 We note that prior work studying light–matter interactions has also used DFTB to obtain the matter's electronic structure.^60,61

In this work, for each Ĥ_Mn, a periodic cubic simulation cell of size 10 Å was constructed, containing 33 water molecules to match the experimental density of bulk liquid water. Periodic boundary conditions were applied in all three spatial directions. Canonical (NVT) molecular dynamics simulations were carried out at 300 K using a Nosé–Hoover thermostat and a time step of 0.5 fs.⁶² Each trajectory was propagated for 150 ps, and atomic coordinates were recorded every 1 fs. Five independent NVT trajectories were generated for statistical sampling. From the equilibrated portion of each trajectory (after 80 ps), 100 configurations with velocities were randomly selected. These were then used to initiate microcanonical (NVE) simulations, each propagated for ∼2 ps with the same time step of 0.3 fs. Electronic structure calculations during the simulations were performed using SCC-DFTB with Brillouin zone integration approximated via Γ-point sampling of a 3 × 3 × 3 folded supercell, implemented using the ‘SupercellFolding’ scheme in DFTB+.⁶³ A detailed list of parameters is provided in the SI.

3 Results and discussion

In this work, the primary computational bottleneck arises from evaluating the derivative of the dipole moment dμ_n/dr_n,j (Born charges), with respect to nuclear motion. In typical classical force fields, atomic partial charges are fixed parameters; consequently, the Born charges reduce to a time-independent constant value dμ_n/dr_n,j ≈ z_n,j. In our DFTB-based simulations, atomic charges are time-dependent and fluctuate along the trajectory. In Fig. 2b, we show the time-dependent Born charges on the O and H atoms, labeled as ∇_Oμ and ∇_Hμ, respectively. We calculate these Born charges using the simple forward finite difference approach, i.e., .

Fig. 2 (a) Born charge distribution of oxygen and hydrogen atoms in liquid water and (b) their time-dependent fluctuations.

Fig. 2a presents the Born charge distribution for the O and H atoms sampled from a dynamical evolution of 500 fs. We observe a broad distribution of O's Born charges, with the majority lying between −1 and 0 with a mode near −0.6. H atoms show a much sharper distribution, with its peak appearing near +0.3 and most values distributed between −0.2 and +0.8.

In Fig. 2b, we present time-dependent Born charges as well as their statistical averages 〈∇_Oμ〉 and 〈∇_Hμ〉. Clearly, the time-dependent fluctuations are large. These large variations illustrate a possible limitation of the presently employed bi-linear (with a linear dipole moment) light–matter couplings in prior theoretical works, which have had limited success in explaining cavity modified ground-state chemical kinetics in the collective coupling regime.^8,23 Neglecting the fluctuations in the Born charges underestimates the strong anharmonic nature of molecular vibrations, misrepresenting the actual light–matter interaction under VSC. While outside the scope of the present work, our computational tool enables investigation of non-linear and anharmonic effects in cavity modification of ground-state chemical reactivity, which will be pursued in our future work.

In this work, we study the modification of the molecular and photonic spectra in the VSC regime. The cavity modified molecular spectra were computed based on linear response theory, where the spectral intensity I_m(ω) is computed from the dipole moment autocorrelation function:

where 〈⋯〉 denotes statistical averaging and [scr S, script letter S] is a smoothening function (a weighted moving average) which allows us to obtain reasonably converged spectra with a small number of trajectories, as was done in ref. 64. We used 40 individual trajectories to compute the molecular spectra which couple to single cavity mode.

Fig. 3c shows the simulated molecular spectra featuring a peak corresponding to the H–O–H bending mode in bulk water outside of the cavity (schematically illustrated in Fig. 3a). This particular peak appears at ∼0.19 eV with a narrow profile and is in good agreement with past experimental and theoretical works.^65,66

Fig. 3 Schematic figure of water (a) outside and (b) inside an optical cavity. Simulated molecular spectra of H–O–H bending motion (c) outside and (d) inside the cavity under vibrational strong coupling. The effects of using Mulliken vs. Born charges on (e) temperature and (f) cavity modified molecular spectra.

In Fig. 3d, we illustrate the cavity modified molecular spectra when considering a single radiation mode coupled to a single box of water (99 atoms). The light–matter Hamiltonian for this reduced system is written as

Here, we have adjusted the mirror spacing to tune the cavity frequency ħω₀ to 0.19 eV, to be in resonance with the H–O–H bending mode. As expected, light–matter interactions lead to a Rabi-splitting (∼0.024 eV for η = 2η₀ a.u. chosen here) with two peaks corresponding to the upper (+) and lower (−) vibro-polaritons.

In this calculation, we used the computationally expensive Born charges in propagating the dynamics of the light–matter hybrid system. In Fig. 3f, we explore the possibility of using the computationally inexpensive Mulliken charges as an approximate replacement for the Born charges in order to compute the cavity modified spectra. Recent works^67,68 that study laser–matter interactions using DFTB used Mulliken charges for propagating dynamics. Fig. 3f shows that the Mulliken charges do provide reasonably accurate polaritonic spectra illustrating an efficient route to computing vibro-polariton linear spectra.

However, such an approximation (replacing Born charges with Mulliken) does lead to long-time heating of the system, making it potentially incompatible for studying cavity modified chemical dynamics or energy transport. Fig. 3e compares the temperature fluctuations over 1 ps. We observed that the temperature fluctuation of water does not exceed ∼307 K, when the Born charges were used (red curve). In contrast, simulations using Mulliken charges exhibit a continuous heating: it exceeds the 310 K within 250 fs and continues to rise over 320 K (orange curve).

In Fig. 4a–c, we illustrate the effect of total light matter coupling η given in eqn (5) on molecular spectra. Here, we have tuned the cavity frequency ħω to be in resonance with the H–O–H bending mode. It shows that as η increases from η₀ to 3η₀ a.u. (where η₀ = 0.0001 a.u.), the upper and lower polariton peaks separate, reflecting an increase in the Rabi-splitting from 0.012 eV to 0.036 eV, proportional to the total light matter coupling (∝η), as is expected from a simple Jaynes–Cummings model^5,6,51 (note that η is the total coupling between light and matter and it does not depend on the number of water boxes or number of molecules). This indicates the lack of anharmonic or nonlinear effects in this particular molecular vibration which also explains the success of using the Mulliken charges in this specific instance.

Fig. 4 Cavity modified molecular spectra and angle resolved photonic spectra under vibrational strong coupling. Panels (a)–(d) show cavity modified molecular spectra at various light–matter coupling strengths for η₀ = 0.0001 with ω₀ = 0.19 eV, and panels (f)–(i) show the corresponding spectra for η₀ = 0.00034 with ω₀ = 0.43 eV. Angle resolved photonic spectra are shown in panel (e) for ω₀ = 0.19 eV and in panel (j) for ω₀ = 0.43 eV.

In Fig. 4f–h, we show the effect of η (chosen as multiples of η₀ = 0.00034 a.u., a constant) on molecular spectra when tuning the cavity frequency to 0.43 eV. This frequency is nearly in resonance with the H–O–H symmetric stretching mode that appears at ∼0.43 eV, as shown in the bare molecular spectrum in Fig. 4i. In the bare molecular spectrum, the symmetric and asymmetric H–O–H stretching vibrations can be observed at ∼0.43 eV and ∼0.47 eV, respectively. Unlike in H–O–H bending mode, the vibro-polariton modes formed here combine two molecular vibrations, with one cavity radiation mode leading to more complex polariton modified molecular spectra. With an increase in η, from η₀ to 3η₀ (here η₀ = 0.00034 a.u.), we observe these two peaks moving farther apart.

In addition to these single cavity mode calculations, we now present the angle resolved vibro-polariton spectra which are computed from the photonic autocorrelation function I_k(ω) which is expressed as

where q̃_k is the photonic position in the reciprocal space. We computed this dispersion using 5 independent trajectories (which were found to converge the polariton spectra for this specific setup). Each trajectory includes 81 localized cavity radiation modes, with each mode coupling to a box of water containing 33 water molecules as schematically shown in Fig. 1b. In our simulation, we kept the total length constant, L_s = 80 × 10⁴ a.u. Each water box was placed at a uniform distance, resulting in a lattice constant of 0.987 × 10⁴ a.u., which is appropriate for computing spectra. However, when studying energy transport or chemical reactivity, the value of a must be carefully chosen and should be used as a convergence parameter.

Fig. 4e shows the angle resolved vibro-polariton spectra as a function of the in-plane cavity wave-vector k, with the k = 0 cavity mode tuned to 0.19 eV. Thus, k = 0 cavity mode is in resonance with the H–O–H bending vibration, forming the upper and lower polariton bands. When k increases, the cavity photon energy rises, while the molecular vibrational energy remains flat. This causes the upper polariton branch to become increasingly photonic, while the lower branch becomes increasingly matter-like at higher k.

In Fig. 4j, we tuned the cavity mirror spacing to adjust the cavity frequency of the k = 0 cavity mode to be 0.43 eV. Here, we can clearly see two flat bands, corresponding to symmetric and asymmetric stretching modes simultaneously hybridized with the cavity band, leading to the formation of three polaritonic bands.

Finally, in Fig. 5, we compute the angle resolved spectra and benchmark the usage of Mulliken charges as an approximate replacement of Born charges. Here, we tune the photon frequency to resonate with the H–O–H bending (0.19 eV) or H–O–H stretching frequencies (0.43 eV) as well as two off-resonance frequencies, 0.16 eV and 0.22 eV. Fig. 5a–c show that the spectra computed using Mulliken charges closely match those obtained with Born charges, capturing all qualitative features of the spectra. In contrast, in Fig. 5(d), we observe noticeable deviations in the angle resolved spectra when Born charges are replaced with Mulliken charges. This is because the cavity modes are simultaneously coupling to the H–O–H symmetric and asymmetric stretching which are known to exhibit substantially more nonlinearities in comparison to the bending motion.⁶⁹ This illustrates that in the presence of substantial molecular nonlinearity, the use of Mulliken charges is insufficient and should be approached with caution.

Fig. 5 Angle-resolved photonic spectra under vibrational strong coupling when using Mulliken charges vs. Born charges at (a) ω₀ = 0.16 eV, (b) ω₀ = 0.19 eV, (c) ω₀ = 0.22 eV, and (d) ω₀ = 0.43 eV. Here, we use the light–matter coupling η = η_c·ω₀^3/2 with η_c = 0.00105.

4 Conclusions

In this work, we propagate the dynamics of strongly coupled vibro-polaritons, formed by coupling molecular vibrations with confined radiation modes inside an optical cavity, where molecular potentials, dipoles and their gradients are computed on-the-fly. In our approach, we adapt a real-space description of the light–matter interaction and exploit its sparsity to devise a low-communication, hub-and-spoke server–client parallelization scheme. Using this approach, we develop an open-source code CavOTF⁴⁸ which is available on GitHub.

We illustrate the utility of our approach by computing the cavity modified molecular and angle resolved photonic spectra. Obtaining the dipole gradients, namely, the Born charges, is the computational bottleneck of our approach as they are evaluated numerically. We find that Mulliken charges can often be used as a computationally inexpensive replacement for computing qualitatively accurate linear spectra. However, when substantial molecular nonlinearity persists, this approach produces less accurate linear spectra. As expected, using Mulliken charges leads to spurious heating of the system and therefore should not be used when studying cavity modified chemical dynamics or cavity modified energy transport.

Conflicts of interest

There are no conflicts to declare.

Data availability

The code is available on Github⁴⁸ with the data corresponding to the figures (as well as the codes to generate the figures) also available in a separate repository.⁷¹

Supplementary information (SI): details of the parameters used in this work; figure illustrating rapid convergence with respect to the number of trajectories. See DOI: https://doi.org/10.1039/d5cp04879f.

Acknowledgements

This work was supported by the Texas A&M startup funds and the Strategic Transformative Research Program (STRP) grant provided by the Texas A&M University. This work used TAMU FASTER and ACES at Texas A&M University through allocation PHY230021, TAMU LAUNCH at Texas A&M University, and SDSC Expanse at the San Diego Supercomputer Center through allocations CHE250156 and CHEM250162 from the Advanced Cyberinfrastructure Coordination Ecosystem: Services⁷⁰ & Support (ACCESS) program, which is supported by U.S. National Science Foundation grants #2138259, #2138286, #2138307, #2137603, and #2138296. A. M. appreciates discussion with Xin Yan, Gerrit Groenhof, Pengfei Huo and Elious Mondal. The authors thank Michael Fowler for his assistance in preparing Fig. 3.

Notes and references

1. A. Thomas, L. Lethuillier-Karl, K. Nagarajan, R. M. Vergauwe, J. George, T. Chervy, A. Shalabney, E. Devaux, C. Genet and J. Moran, et al., Science, 2019, 363, 615–619.
2. W. Ahn, J. F. Triana, F. Recabal, F. Herrera and B. S. Simpkins, Science, 2023, 380, 1165–1168.
3. K. Nagarajan, A. Thomas and T. W. Ebbesen, J. Am. Chem. Soc., 2021, 143, 16877–16889.
4. G. Yin, T. Liu, L. Zhang, T. Sheng, H. Mao and W. Xiong, Science, 2025, 389, 845–848.
5. A. Mandal, M. A. Taylor, B. M. Weight, E. R. Koessler, X. Li and P. Huo, Chem. Rev., 2023, 123, 9786–9879.
6. T. E. Li, B. Cui, J. E. Subotnik and A. Nitzan, Annu. Rev. Phys. Chem., 2022, 73, 43–71.
7. M. Ruggenthaler, D. Sidler and A. Rubio, Chem. Rev., 2023, 123, 11191–11229.
8. J. A. Campos-Gonzalez-Angulo, Y. R. Poh, M. Du and J. Yuen-Zhou, J. Chem. Phys., 2023, 158, 230901.
9. J. Lather and J. George, J. Phys. Chem. Lett., 2020, 12, 379–384.
10. S. Kéna-Cohen and J. Yuen-Zhou, ACS Cent. Sci., 2019, 5, 386–388.
11. J. Lather, P. Bhatt, A. Thomas, T. W. Ebbesen and J. George, Angew. Chem., Int. Ed., 2019, 58, 10635–10638.
12. A. Thomas, J. George, A. Shalabney, M. Dryzhakov, S. J. Varma, J. Moran, T. Chervy, X. Zhong, E. Devaux and C. Genet, et al., Angew. Chem., 2016, 128, 11634–11638.
13. K. Hirai, R. Takeda, J. A. Hutchison and H. Uji-i, Angew. Chem., 2020, 132, 5370–5373.
14. J. C. Nelson and M. L. Weichman, J. Chem. Phys., 2024, 161, 074304.
15. L. Chen, A. P. Fidler, A. M. McKillop and M. L. Weichman, Nanophotonics, 2024, 13, 2591–2599.
16. A. P. Fidler, L. Chen, A. M. McKillop and M. L. Weichman, J. Chem. Phys., 2023, 159, 164302.
17. R. M. Vergauwe, A. Thomas, K. Nagarajan, A. Shalabney, J. George, T. Chervy, M. Seidel, E. Devaux, V. Torbeev and T. W. Ebbesen, Angew. Chem., Int. Ed., 2019, 58, 15324–15328.
18. J. Lather, A. N. Thabassum, J. Singh and J. George, Chem. Sci., 2022, 13, 195–202.
19. F. J. Garcia-Vidal, C. Ciuti and T. W. Ebbesen, Science, 2021, 373, eabd0336.
20. M. A. Michon and B. S. Simpkins, J. Am. Chem. Soc., 2024, 146, 30596–30606.
21. G. D. Wiesehan and W. Xiong, J. Chem. Phys., 2021, 155, 241103.
22. M. V. Imperatore, J. B. Asbury and N. C. Giebink, J. Chem. Phys., 2021, 154, 191103.
23. L. P. Lindoy, A. Mandal and D. R. Reichman, Nanophotonics, 2024, 13, 2617–2633.
24. L. P. Lindoy, A. Mandal and D. R. Reichman, Nat. Commun., 2023, 14, 2733.
25. S. M. Vega, W. Ying and P. Huo, J. Am. Chem. Soc., 2025, 19727.
26. W. Ying, M. A. Taylor and P. Huo, Nanophotonics, 2024, 13, 2601–2615.
27. W. Ying and P. Huo, Commun. Mater., 2024, 5, 110.
28. J. Sun and O. Vendrell, J. Phys. Chem. Lett., 2023, 14, 8397–8404.
29. T. E. Li, A. Nitzan and J. E. Subotnik, Angew. Chem., Int. Ed., 2021, 60, 15533–15540.
30. X. Ji and T. E. Li, J. Phys. Chem. Lett., 2025, 16, 5034–5042.
31. T. E. Li, J. Chem. Theory Comput., 2024, 20, 7016–7031.
32. T. E. Li, A. Nitzan and J. E. Subotnik, Nat. Commun., 2022, 13, 4203.
33. D. S. Wang, J. Flick and S. F. Yelin, J. Chem. Phys., 2022, 157, 224304.
34. D. S. Wang and S. F. Yelin, ACS Photonics, 2021, 8, 2818–2826.
35. C. Schäfer, J. Flick, E. Ronca, P. Narang and A. Rubio, Nat. Commun., 2022, 13, 7817.
36. S. Mondal, D. S. Wang and S. Keshavamurthy, J. Chem. Phys., 2022, 157, 244109.
37. W. Liu, J. Chen and W. Dou, J. Phys. Chem. C, 2024, 128, 12544–12550.
38. X. Li, A. Mandal and P. Huo, Nat. Commun., 2021, 12, 1315.
39. A. Mandal, X. Li and P. Huo, J. Chem. Phys., 2022, 156, 014101.
40. J. A. Campos-Gonzalez-Angulo, R. F. Ribeiro and J. Yuen-Zhou, Nat. Commun., 2019, 10, 4685.
41. B. Xiang, R. F. Ribeiro, M. Du, L. Chen, Z. Yang, J. Wang, J. Yuen-Zhou and W. Xiong, Science, 2020, 368, 665–667.
42. Z. T. Brawley, S. Pannir-Sivajothi, J. E. Yim, Y. R. Poh, J. Yuen-Zhou and M. Sheldon, Nat. Chem., 2025, 17, 439–447.
43. Y. R. Poh, S. Pannir-Sivajothi and J. Yuen-Zhou, J. Phys. Chem. C, 2023, 127, 5491–5501.
44. P.-Y. Yang and J. Cao, J. Phys. Chem. Lett., 2021, 12, 9531–9538.
45. B. S. Simpkins, A. D. Dunkelberger and J. C. Owrutsky, J. Phys. Chem. C, 2021, 125, 19081–19087.
46. T. Liu, G. Yin and W. Xiong, Chem. Sci., 2025, 16, 4676–4683.
47. Q. Yu and J. M. Bowman, J. Chem. Theory Comput., 2024, 20, 4278–4287.
48. https://github.com/mandalgrouptamu/cavOTF, Accessed: 2025-12-01.
49. O. Dmytruk and M. Schiró, Phys. Rev. B, 2021, 103, 075131.
50. J. Li, D. Golez, G. Mazza, A. J. Millis, A. Georges and M. Eckstein, Phys. Rev. B, 2020, 101, 205140.
51. D. Jasrasaria, A. Mandal, D. R. Reichman and T. C. Berkelbach, J. Chem. Phys., 2025, 162, 014109.
52. R. H. Tichauer, J. Feist and G. Groenhof, J. Chem. Phys., 2021, 154, 104112.
53. A. Mandal, D. Xu, A. Mahajan, J. Lee, M. Delor and D. R. Reichman, Nano Lett., 2023, 23, 4082–4089.
54. N. C. Benson and V. Daggett, J. Phys. Chem. B, 2012, 116, 8722–8731.
55. M. Elstner and G. Seifert, Philos. Trans. R. Soc., A, 2014, 372, 20120483.
56. M. Gaus, Q. Cui and M. Elstner, J. Chem. Theory Comput., 2011, 7, 931–948.
57. G. Zheng, M. Lundberg, J. Jakowski, T. Vreven, M. J. Frisch and K. Morokuma, Int. J. Quantum Chem., 2009, 109, 1841–1854.
58. S. Wickramasinghe, A. Hoehn, S. T. Wetthasinghe, H. Lin, Q. Wang, J. Jakowski, V. Rassolov, C. Tang and S. Garashchuk, J. Phys. Chem. B, 2023, 127, 10129–10141.
59. J. Jakowski, B. Hadri, S. J. Stuart, P. Krstic, S. Irle, D. Nugawela and S. Garashchuk, Proceedings of the 1st Conference of the Extreme Science and Engineering Discovery Environment: Bridging from the eXtreme to the campus and beyond, 2012, pp. 1–7.
60. D. Sidler, C. M. Bustamante, F. P. Bonafé, M. Ruggenthaler, M. Sukharev and A. Rubio, Nanophotonics, 2025, 14, 4941–4955.
61. C. M. Bustamante, F. P. Bonafé, M. Sukharev, M. Ruggenthaler, A. Nitzan and A. Rubio, J. Chem. Theory Comput., 2025, 21, 9823–9831.
62. G. J. Martyna, M. E. Tuckerman, D. J. Tobias and M. L. Klein, Mol. Phys., 1996, 87, 1117–1157.
63. B. Hourahine, B. Aradi, V. Blum, F. Bonafe, A. Buccheri, C. Camacho, C. Cevallos, M. Deshaye, T. Dumitrica and A. Dominguez, et al., J. Chem. Phys., 2020, 152, 124101.
64. https://github.com/TaoELi/cavmd_examples_h2o/blob/main/water_VUSC/collect_all_data_N.py, Accessed: 2025-05-26.
65. C.-C. Yu, K.-Y. Chiang, M. Okuno, T. Seki, T. Ohto, X. Yu, V. Korepanov, H.-O. Hamaguchi, M. Bonn and J. Hunger, et al., Nat. Commun., 2020, 11, 5977.
66. W. B. Carpenter, J. A. Fournier, R. Biswas, G. A. Voth and A. Tokmakoff, J. Chem. Phys., 2017, 147, 084503.
67. Q. Xu, D. Weinberg, M. S. Okyay, M. Choi, M. Del Ben and B. M. Wong, J. Am. Chem. Soc., 2024, 146, 35313–35320.
68. Q. Xu, M. Del Ben, M. Sait Okyay, M. Choi, K. Z. Ibrahim and B. M. Wong, J. Chem. Theory Comput., 2023, 19, 7989–7997.
69. M. Vinaykin and A. V. Benderskii, J. Phys. Chem. Lett., 2012, 3, 3348–3352.
70. T. J. Boerner, S. Deems, T. R. Furlani, S. L. Knuth and J. Towns, Practice and Experience in Advanced Research Computing (PEARC '23), 2023, pp. 1–4.
71. https://github.com/amiraminitamu/Data-cavOTF, Accessed: 2025-12-09.

---

Footnote

† These authors contributed equally to this work.

---