Machine Learning for Polaritonic Chemistry: Accessing Chemical Kinetics

Altering chemical reactivity and material structure in confined optical environments is on the rise, and yet, a conclusive understanding of the microscopic mechanisms remains elusive. This originates mostly from the fact that accurately predicting vibrational and reactive dynamics for soluted ensembles of realistic molecules is no small endeavor, and adding (collective) strong light–matter interaction does not simplify matters. Here, we establish a framework based on a combination of machine learning (ML) models, trained using density-functional theory calculations and molecular dynamics to accelerate such simulations. We then apply this approach to evaluate strong coupling, changes in reaction rate constant, and their influence on enthalpy and entropy for the deprotection reaction of 1-phenyl-2-trimethylsilylacetylene, which has been studied previously both experimentally and using ab initio simulations. While we find qualitative agreement with critical experimental observations, especially with regard to the changes in kinetics, we also find differences in comparison with previous theoretical predictions. The features for which the ML-accelerated and ab initio simulations agree show the experimentally estimated kinetic behavior. Conflicting features indicate that a contribution of dynamic electronic polarization to the reaction process is more relevant than currently believed. Our work demonstrates the practical use of ML for polaritonic chemistry, discusses limitations of common approximations, and paves the way for a more holistic description of polaritonic chemistry.

Delivering a conclusive theoretical understanding for vibrational strong coupling has remained difficult.Especially the experimentally observed resonance dependence in combination with an increase of rate changes for increasing emitter concentration, and clear trends in chemical kinetics are critical features that a theo-retical model should capture.Initial attempts based on standard transition-state theory [49][50][51] failed to reproduce any significant frequency dependence.Active development along the lines of Grote-Hynes [52][53][54] and Pollak-Grabert-Hänggi [55] theory showed some frequency dependence but a connection to experiments has remained unsuccessful.A recent work by Schäfer et al. [13] tackled the problem from first principles by utilizing quantum electrodynamical density-functional theory (QEDFT) [56][57][58][59] in combination with a self-consistent update of the nuclear motion according to Ehrenfest's equation for the experimentally investigated deprotection reaction of 1-phenyl-2-trimethylsilylacetylene (PTA) [14,15].This ab initio theory recovered critical components of the frequency dependence and suggested a microscopic theory based on adjusted vibrational energy redistribution for the cavity induced modification of chemical reactivity.Chen et al. found experimental support for this hypothesis using 2D spectroscopy [10].Nonetheless, a prediction of kinetic quantities remained inaccessible given the computational cost of QEDFT in combination with Ehrenfest dynamics.
In this work, we establish a framework that combines machine learning (ML) models, trained on data from density functional theory (DFT) calculations, with molecular dynamics (MD) simulations to arrive at a more efficient, yet accurate description of the experimentally relevant S N 2 reaction of PTA [14] under strong coupling to a cavity.We find a pronounced frequency dependence for the chemical reaction rate constant along with changes in reaction enthalpy and entropy that are qualitatively consistent with experiment.Interestingly, we discover frequency domains outside the experimentally validated window for which the present ML-accelerated approach predicts a rate constant enhancing character in contrast to earlier fully ab initio simulations [13], which rather suggest inhibition.Here, we tentatively attribute this difference to the simplifications inherent to the present MD approach, suggesting that the latter, despite being widely used, has relevant limitations in its applicability to polaritonic chemistry.Further investigations, beyond the scope of the present work, will be needed to provide a more detailed understanding, which we expect to also generate useful insight into the microscopic mechanisms.

II. METHODOLOGY
Nonrelativistic quantum electrodynamics commonly starts at the minimal coupling Hamiltonian in the Coulomb gauge [5,[60][61][62], where all charged particles couple via longitudinal Coulomb interaction to each other and to the transverse vector potential (∇ Polaritonic chemistry comes in different flavors that are distinguished by the frequency of the confined modes and the physical nature of the resonator, e.g., plasmonic, Fabry-Pérot, or whispering gallery cavities, which influences their fundamental coupling strength and quality factor.The accurate description of a realistic optical environment is a challenge in itself [63][64][65][66][67], but simple singlemode models often suffice to obtain a qualitative understanding of the relevant emitter dynamics.The goal of this work is to provide a qualitative investigation of vibrational strong coupling and its influence on chemical reactivity for experimentally relevant molecules. In particular, we focus on kinetic changes, its impact on enthalpy and entropy, and the consequences of simplifications in the molecular dynamics.For this reason, we stay conceptually close to the previous work by Schäfer et al. [13] and rely on the Born-Oppenheimer approximation, projecting on the electronic ground-state and ignoring electronic polarizations induced by the cavity, such that the effective nuclear-photonic Hamiltonian takes the form [58,60,62,68] Ĥ = Tnuclei + VPES + ℏω c (â † â + 1  2 ) The nuclei couple self-consistently to a single electromagnetic mode of the cavity, treated classically as the mode displacement q c (t) = ℏ 2ωc ⟨â † + â⟩ in the following, with frequency ω c , effective cavity volume V c , fixed polarization ε c , and molecular dipole moment μ.Taking the classical limit for the nuclei, the system follows standard Hamiltonian mechanics with forces originating from the Poisson bracket {p j , H}.In order to see appreciable changes for the dynamics of a single molecule, we choose large coupling values g 0 = ea 0 ω c /2ℏε 0 V c that are consistent with Ref. 13.We refer the interested reader to Refs.[13,67,[69][70][71] for an extended discussion on the potential motivation of effective single-molecule coupling values.The chosen coupling is considerably larger than experimentally achievable values in Fabry-Pérot cavities but the qualitative agreement with experiments suggests similar microscopic mechanisms.Our results and discussion can be partially transferred to plasmonic cavities which feature substantial single-molecule couplings [21,23].
The electronic force acting on nucleus j are obtained from the potential energy surface (PES) according to and contributes to the total force F j = F j PES + F j c together with the optical force, which can be computed from the derivative of the dipole moment vector The only interaction between light and matter arises then via the time evolution of the photonic mode q c (t) and the gradient of the projected dipole moment.The cavity mode displacement q c (t) depends on the history of the dipole moment (see SI Sec.II A for details) where vanishing of the initial momentum qc (0) = 0 is explicitly enforced.We leverage the modular character of the forces by implementing a custom ("cavity") propagator using the ase Python package [72].The calculator requires merely the initial molecular configuration as well as estimators for the forces arising from the PES (Eq.( 1)) and the dipole moment vector (Eq.( 2)).Possible estimators may include machine-learned models, empirical force fields, and even ab initio calculations based on, e.g., DFT.We note that this approach can also be easily combined with the embedding radiation-reaction approach [67] in the future, thus providing an elegant path for including collective coupling and realistic optical environments.
Obtaining forces and dipole derivatives thus represents the main obstacle in the MD approach, as in reality a molecule can easily pass through tens of thousands of configurations before a reaction occurs.One possible, but practically often too expensive, approach is to perform ab initio electronic structure calculations at each point, Positional information of a structure is translated into descriptor space and neuroevolution potential models for the potential energy surface (PES) and the dipole moment vector µ are trained using supervised learning (left-hand side).The final models are combined to compute the effective forces acting on the nuclei that are then used to propagate the system in time (right-hand side).
usually referred to as ab intio MD [73].While the cost of a single DFT calculation might be comparably low, the sheer quantity of calculations required for a statistically meaningful evaluation, i.e., obtaining thousands of trajectories with tens of thousands of DFT evaluations each, is highly prohibitive and practically limits simulation time as well as system and ensemble size.Another prominent approach is based on empirical force fields [74,75], which are computationally orders of magnitude more efficient than DFT calculations.They are, however, restricted with respect to accuracy as well as availability, and offer limited transferability.Even if we would have a suitable force field for a system, there is no guarantee that, e.g., a force field fitted in an aqueous solution will perform well for simulations in vacuum and vice versa.
As we will show in the following, employing ML techniques in combination with MD provides a feasible, predictive, and scalable path, striking a balance between computational cost and accuracy.

A. Electronic and optical forces using neural networks
In the present work, we developed two different models using the neuroevolution potential (NEP) framework as implemented in the gpumd package [76][77][78].One for predicting the PES V PES (r), as well as the associated force F PES (r), and one for predicting dipole moments µ(r) (Fig. 1; left).We refer the reader to the SI for an extensive discussion of the training and testing procedure of both models.
The NEP models are then combined to obtain the total force used to propagate the system in time with a custom integrator implemented in the ase package [72] (Fig. 1; right).We refer the interested reader to Refs.[78,79] for a more extensive presentation of the NEP framework and its application to tensorial quantities.
The NEP approach employs a simple forward bias multi-layer perceptron with a single hidden layer in combination with a flexible descriptor to predict the atomic energy U i for each atom in a system, by decomposing the total energy into individual contributions from each atom, U = i U i .The model consists of a fully connected network with a single hidden layer, yielding the following expression for the predicted energy, w (1)  µ tanh The two weight matrices w (0) µv and w µ are the weights for the input and hidden layer, with b (0) µ and b (1) being their respective biases, and tanh is used as the activation function for the input layer.The so-called descriptor vector q i v of length N des indexed by v can be seen as representation of the local chemical environment of atom i, is a function of the pairwise distances r ij = r i − r j , and serves as input to the network.If q i v uniquely describes a molecular configuration is determined by a basis expansion over N-body interactions within a cutoff radius r c .The expansion is truncated at 4-body interactions, which is sufficient to accurately describe local changes.
A key feature of the NEP formalism is that these descriptors also contain trainable parameters, which allows the network to tailor the descriptors more individually to different atomic configurations.Predictions for forces and virials can be obtained by computing the gradient of the predicted site energies, i.e., the PES force acting on atom i is Additionally, the NEP formalism may be extended to predict other tensorial properties such as dipole mo-ments, which are obtained as Having replaced the estimates for forces and dipoles with our NEP models, we are now equipped to address the question how an optical resonator might influence the here discussed S N 2 reaction.

III. RESULTS AND DISCUSSION
In a first step, we perform ensemble calculations with preserved particle number, volume, and temperature (NVT) in absence of the cavity to obtain a reference value for the transition-state enthalpy of ∆H ‡ = 0.345 eV = 33.3kJ mol −1 .The latter is in excellent agreement with experimental estimates of ∆H ‡ = 35 ± 4 kJmol −1 [15].DFT calculations using the nudged elastic band approach in combination with a transition-state optimization provide a higher barrier of 0.43 eV .This illustrates the limitation of estimating the enthalpy from the 0 K energy difference between minimum and transition state and the significance of vibrational contributions.Detailed information as well as various benchmarks can be found in the SI.We define a reaction event when the relevant Si-C bond is stretched beyond 3.5 Å, which exceeds the transition-state Si-C distance of approximately 3 Å, to avoid counting eventual recrossing events.If not further specified, all observables are obtained from ensemble averages involving 1000 trajectories.They are initialized with Boltzmann sampled velocities (kept fixed when changing cavity parameters) at the nonequilibrium F -+ PTA state used in Ref. 13, and propagated preserving particle number, volume, and energy (NVE).The cavity displacement is selected such that no electric field exists at time zero, i.e., the cavity force is zero.The initial state features an energy difference to the minimum PTAF -configuration of 1.34 eV.This shortens the necessary calculation time and avoids spurious interplay between thermostat and cavity.One should note, however, that it also limits the transferability of the obtained rate constants to experimental observations.Nonetheless, we can extract changes in rate constant and thus contribute valuable insight to the current hypothesis behind polaritonic chemistry.

A. Strong Coupling
Strong coupling requires the existence of optically active vibrational modes near the cavity frequency.Certainly, the vibrational spectrum is sensitive to temperature, as illustrated in Fig. 2A.The reaction-defining Si-C bond contributes fractionally to most vibrational excitations but primarily at frequencies below 1300 cm −1 .Fig. 2B illustrates the corresponding power spectrum during the reaction process at 400 K with strong coupling to the cavity.We keep the ratio g 0 /ω c constant in our calculations.Changing the cavity length L c (V c = A c L c ), which is the experimental way to tune the frequency of the Fabry-Pérot cavity, leads in a simplified modes picture to g 0 ∝ √ ω c / √ L c and ω c ∝ 1/L c , i.e., a larger distance between the mirrors reduces both frequency and coupling strength with 1/L c .Following the gray-dashed diagonal ω = ω c , we can clearly identify multiple avoided crossings (hybridizations) with substantial Si-C contribution.Each of the avoided crossings contributes with additional low-energy components (following the vertical gray dotted lines), suggesting comparably slow changes, i.e., on the timescale of the reaction.The reorganization of the methyl groups and the proper F-Si-C alignment (bending modes) are critical steps in the reaction pathway.Their interplay with the cavity manifests in the constant contribution at ℏω ≈ 117 cm −1 , which is further detailed in Sect.III E.

B. Resonance dependence
Intuitively, we expect the low-frequency Si-C contribution to play a major role for the reactivity as the defining Si-C bond breaking step requires a few hundred femto seconds.Given the fact that the contribution changes non-monotonically, it is not surprising that the rate constant also changes non-monotonically (Fig. 2C).We observe pronounced regions of inhibited reactivity, especially around 200 cm −1 and the domain including the 770 and 856 cm −1 vibrations.Furthermore, a near twofold enhancement of the rate constant at around 290 and 460 cm −1 is visible.Experimental investigations are only available near 856 cm −1 and support the general inhibiting trend in this domain, albeit featuring an inhibiting effect only in a narrow frequency window.A non-monotonous rate change, however, is in conflict with previous ab initio simulations in Ref. 13.Interestingly, the most pronounced feature at 200 cm −1 seems not related to any optical excitation and is also not connected to the curvature at the transition state, which we estimate with both our NEP model and DFT calculations to be approximately 73 cm −1 .
We will elaborate the conceptual differences to Ref. 13 and possible explanations in Sect.III E, but it should be noted that a strict comparison is difficult due to the difference in observable, temperature, and a substantial difference in statistical sampling.The slow but steady increase in rate constant for large frequencies could partially originate from numerical deviations in the finitedifferences approximation of the dipole gradient (see conservation of energy in SI).With this in mind, let us disentangle the mechanism behind the catalysing and inhibiting effects by estimating enthalpic and entropic changes.I. Negative enthalpy originates from the selected initial state and should only be interpreted relative to the equilibrium enthalpy (see text).The rate constant is calculated as number of products Np per ps.

C. Kinetics
We repeat calculations for the rate constant for various temperatures at four frequencies extracted from Fig. 2C which are characteristic for their respective frequency domains.The results are collected in an Eyring plot (Fig. 2D) and indicate conceptually different mechanisms for inhibition and rate constant enhancement.The corresponding Arrhenius plot is provided in the SI.We would like to emphasize that those changes originate from the independent dynamics of an ensemble of trajectories and the Eyring equation is used only to enable a comparison with the experimentally extracted enthalpy and entropy.Changes in the transmission coefficient κ thus contribute to an altered entropy.
Without cavity environment (ω c = 0, black), we es-timate a weak but negative enthalpic barrier ∆H.This should be seen as the change induced by the chosen initial state and NVE conditions, i.e., the elevated initial state will provide almost no energetic barrier and yet, the cavity will alter this "barrier".We demonstrate in SI Sec.II B that performing rate constant estimates under NVT conditions and sufficient equilibration time provides accurate rates in agreement with the experiment outside the cavity, i.e., our methodology provides the correct kinetics outside the cavity but the chosen initial state shifts the enthalpy up in energy.We collect the changes in enthalpy ∆∆H and entropy ∆∆S in relation to the cavityfree NVE results in Table I.
The ω c = 198 cm −1 excitation, without optically active vibrational mode support, shows a weak increase in entropy, suggesting that the dynamic effect of the cavity is to "simply" raise the transition-state barrier.Even though ω c = 856 cm −1 (green) shows an overall smaller change, entropic changes are large compared to the other frequencies, which suggests that the cavity assigns a slightly stronger dissociative character to the reaction.Both features are in qualitative agreement with experiment [14,15].Taking into account that the correct enthalpy obtained from our NVT calculations in free space is ∆H ‡ = 0.345 eV, the inhibiting frequencies render the reaction more temperature sensitive, which can lead to enhanced rates for large temperatures.Such a trend in temperature sensitivity has been widely observed in experiments [14,15,18].
At the catalysing frequency ω c = 461 cm −1 on the other hand there is almost no change in enthalpy, but a noticeable change in entropy, suggesting that the mechanism dominating here is not related to the experiments which typically showed a clear change in enthalpy.It is important to emphasize, that utilization of the Eyring equation is especially problematic in this domain as a further increase in reaction speed implies that the trajectory will spend less time around the reactant well.This in turn implies that the kinetic arguments underlying the Eyring equation, i.e., a separation of time scales between reactant equilibration and transmission process, become questionable since the transition-state is then part of the equilibration process.Nonetheless, a similar offset in rate constant without change in enthalpy has been observed when employing Grote-Hynes rate-theory [80].All three domains relate to different kinetic changes and thus suggest slightly different mechanisms.

D. Vibrational Dynamics
Let us shine a bit more light on the mechanistic differences between the chosen frequency domains that catalyze (461 cm −1 ) or inhibit reactions without (198 cm −1 ) or with (856 cm −1 ) vibrational support.Fig. 3 illustrates the accumulated difference in normal mode occupation between a given cavity frequency and free space during the reaction, averaged over the full ensemble.The corresponding occupation differences are presented in the SI, where the overall structure for ω c = 856 cm −1 is comparable to Ref. [13].Optically relevant domains around 300, 550, 770, and 1200 cm −1 are noticeably affected more strongly when selecting the cavity frequencies ω c = 461 cm −1 and ω c = 856 cm −1 .Choosing a specific cavity frequency affects the vibrational modes in energetic proximity more strongly.This effect is especially apparent for ω c = 856 cm −1 which involves the C=C stretching mode around 1200 cm −1 (see green bars in Fig. 3).We quantify the changes in cross-correlation between the IR spectrum and the accumulated (absolute) difference (Fig. 3 top (bottom)) in Table II with the help of the relative difference between the cross-correlation for ω c ∈ 198, 461, 856 cm −1 and the high-frequency value ω c = 1251 cm −1 .The resulting change indicates how strongly differences in normal-mode occupation correlate with infrared activity.Changes in the normal mode occupation for cavity frequencies with optical support, i.e., ω c = 461 cm −1 and ω c = 856 cm −1 , feature larger correlation.

TABLE II. Relative difference of the cross-correlation function δδC
− 1, where we define , where A(A)D represents the accumulated (absolute) difference.Only the frequency domain illustrated in Fig. 3 was utilized in the calculation of the cross-correlation function.Overall stronger correlation for ω c = 461, 856 cm −1 suggests that the microscopic mechanism is more strongly characterized by redistribution of vibrational energy between optically active modes and is thus optically mediated, as suggested in Ref. 13.In other words, the cavity facilitates energy exchange between optically active modes, especially within an energy-window around the cavity frequency.Since the Si-C bond is an essential ingredient in the reaction process, a stronger involvement of Si-C bond stretching (gray circles shown in Fig. 3) in the affected normal modes will result in a larger impact on the chemical reaction.The vibrational analysis supports then also the previous hypothesis that chemical changes without support by optically active modes follow, to a certain degree, a different mechanism [13].Surprisingly, however, even though ω c = 461, 856 cm −1 seem to share a very similar mechanism based on the ML+MD analysis, their effect on the enthalpy is qualitatively different (catalysing vs inhibiting).Albeit previous ab initio calculations did not provide access to kinetic changes, the corresponding analysis provided no indication of a qualitative difference between ω c = 461 cm −1 and ω c = 856 cm −1 [13].Let us reflect in the following section on the underlying approximation of our molecular dynamics simulations in order to understand this contradicting observation.

E. Limitations of simplified cavity molecular dynamics
The experimentally relevant domain around 856 cm −1 [14,15]  culations as well as in experiment and the recent ab initio QEDFT calculations with nuclear motion according to the Ehrenfest equations of motion [13].The existing information to changes in chemical rates is thus consistent at ω c = 856 cm −1 .On the other hand, the ML+MD calculations show enhanced rate constants at 461 cm −1 , an effect that has not been observed in the QEDFT calculations.A direct comparison of Fig. 2 to Ref. [13] is, however, problematic since the latter used an incomplete sampling with a strong bias towards the high-energy tail of the Boltzmann distribution.In other words, the QEDFT results have been obtained at an effectively higher temperature.
To allow for a more reliable comparison and shed light on the apparent discrepancy between the approaches, we recalculated the rate constant changes with our ML+MD approach for the same initial velocities as the QEDFT calculations.Fig. 4 sets the newly obtained rate constants from our ML+MD calculations (black dots) in contrast with the average Si-C distances calculated with QEDFT (blue stars) and taken from Ref. [13].We demonstrate in SI Sec.II D that Si-C distance and rate constant are well correlated in this case.
Ignoring the offset, a slight frequency shift, and the qualitatively different behaviour near ω c = 0, the overall shapes of the reaction rate constant profiles shown in Fig. 4 are consistent.Given identical initial conditions, ML+MD and QEDFT provide thus a similar profile in the intermediate frequency domain but this profile is elevated into the catalysing domain for the ML+MD calculations.Looking back at Fig. 2C, the major catalysing feature of the ML+MD calculations at 461 cm −1 is consistent between proper (Fig. 2C) and incomplete sam-pling (Fig. 4).Especially low-frequency features are considerably shifted and altered in strength, potentially originating from the quicker average reaction time when sampling from the high-energy tail of the Boltzmann distribution, which implies that the cavity has less time to influence the reaction.We recall from Sect.III C, that shorter reaction times emphasize dynamic contributions, calling the concept of a kinetic reaction rate further into question.It is thus plausible that the overwhelming catalysing strength is partially an artifact and we therefore suggest to focus on the qualitative trend only.This leaves but one question: Can we draw any conclusion about the mechanism of vibrational strong coupling from the observed discrepancy?

Relevance of dynamic electronic polarization
As we established at the beginning, effective nuclear forces are comprised of the adiabatic electronic Born-Oppenheimer forces and the dynamic optical forces mediated via dipolar changes induced by nuclear displacement.This implies that the static polarization of the electric system due to the instantaneous cavity field as well as its non-adiabatic corrections, quantum nuclear, and quantum light-matter effects are absent -as is the case in most available theoretical investigations for chemical reactivity affected by strong coupling.Non-adiabatic electron-nuclear effects are expected to play a minor role as the electronic excited space is separated by about 3 eV at minimum and transition state.While there has been recent discussions about the potential need to consider the full quantum light-matter interaction in model sys- (Top) Black dots: Change in rate constant compared to freespace for the unidirectional reaction PTAF − → FtMeSi+PA − using the 30 initial configurations used in Ref. [13] and g0/ωc = 1.132, but propagated with our NEP based molecular dynamics calculator.Blue stars: Change in average Si-C distance obtained from the QEDFT calculations in Ref. [13] amplified by a factor 50. Transmission spectrum obtained from DFT at 0K with harmonic approximation consistent with Ref. [13] (blue dashed), and using our NEP model and GPUMD at 400K NVE conditions (red solid).Vertical lines indicate characteristic features observed in Ref. [13].We show in SI Sec.II D that Si-C distance and rate constant are clearly correlated, i.e., the qualitative trend of QEDFT and ML+MD can be set in relation.(Bottom) ML+MD calculations repeated for different fundamental light-matter coupling strength.
tems to recover resonant features in the cavity modified reactivity [81,82], it remains up to debate if this is true for realistic systems under standard ambient conditions.Such a question requires a nuanced discussion based on the specific system at hand and we expect that the answer will vary strongly between collective and single molecular coupling.
Lastly, and probably most importantly, nuclear motion will induce strong optical fields which in turn polarize the electronic system.This happens similar to a static external potential ∝ g 0 q c ε c • μe , or via the selfpolarization contributions It is important to note that this polarization, although being instantaneous in the sense of the Born-Oppenheimer approximation, is inherently dynamic due to its origin from the nuclear dynamics.We will label this effect in the following dynamic electronic polarization, to emphasize that the effect does not originate from the hybridization of the light-matter ground state nor from electronic transitions (non-adiabatic coupling elements).This form of dynamic electronic polarization can be formally incorporated via the cavity Born-Oppenheimer approximation [83,84], where the photonic displacement q c is considered as an additional parametric variable and self-polarization contributions are added to the electronic structure calculations.When included, the dynamic electronic polarization can lead to notable asymmetries in the vibrational polaritons [13,84] and increase the effective reaction barrier, thus reducing the chemical reactivity [71,85].Our current ML+MD calculations are lacking the possibility to describe this dynamic electronic polarization and will thus tend towards a higher reactivity and more symmetric Rabi-splittings (see Fig. 2 B).
Consulting perturbation theory, cavity induced changes in the electronic ground-state energy scale approximately as where ∆E e ≫ ω c represents the electronic excitation energy for the dominant dipole transitions.Perturbation theory should provide an adequate estimate for electronic changes since the expansion parameter ∝ 1/ √ V c is small and only g 0 /ω c takes appreciable values.Fixing again g 0 /ω c = const.and noticing a larger dipole near the transition state (see Note [87]), an additional increase of the barrier ∝ ω c would be missing in our MD calculations compared to the QEDFT calculations.As a result, our ML+MD calculations provide sensible values near ω c = 0 and for large frequencies, but have the tendency to be overly reactive in the intermediate domain as they lack an additional suppression from dynamic electronic polarization.
Larger cavity frequencies ω c are driven resonantly only by vibrational modes of comparable frequency.If such a vibration contributes to the reaction, i.e., whether it will be noticeably excited during the reaction, depends largely on its Si-C contribution.If the cavity can exert notable effects on a reaction event will therefore depend on three time scales: (i) The strongly occupied reactive modes of frequency ω SiC , (ii) the cavity frequency ω c , and (iii) the frequency with which the PES is modulated ω modPES (dynamic electronic polarization) as consequence of the nuclear dynamics.Since ω modPES depends partially on the optical mode q c , which oscillates with ω c , it seems intuitive that the largest effect of dynamic electronic polarization on the reactive modes ω SiC can be expected when all time scales are comparable.As emphasized by Fig. 2B and Fig. 3, most of the optically active modes with relevant Si-C contribution are located below 840 cm −1 .The dynamic contribution of electronic polarization will then play a decreasing role at higher frequencies.How the precise interplay between nuclear motion, strongly coupled cavity, and cavity-modulated electronic polarization affects the reactivity goes beyond the scope of this work.Our results emphasize the role of dynamic electronic polarization but also indicate that it is unlikely to be the only relevant contribution.Am-ple experimental work, however, demonstrated changes in solute-solvent interaction [18,88] under vibrational strong coupling which indicates the involvement of dispersive interactions mediated by electronic polarization.

Additional considerations
Let us consider the analogy of our MD calculations as the self-consistently driven dynamic of a ball on a highdimensional energy surface.If we intend to cross a specific barrier, but let the cavity periodically remove kinetic energy before inserting it back at a later time, we can imagine that the likelihood to cross the barrier is modulated by the frequency with which the cavity is oscillating.We can indeed observe such an effect at low cavity frequencies where reaction events appear in bursts that are related to the cavity frequency.Fig. 5 shows the Fourier transform of the auto-correlation function of the change in product of the S N 2 reaction.The linear dispersion at low frequencies clearly shows that the bursts in reaction events are correlated with the cavity frequency, a feature that is absent in the QEDFT calculations of Ref. [13].This trend continues up to a bending mode at 117 cm −1 , which contributes to the necessary rearrangement of the methyl groups.Bending modes have recently been identified as a relevant component in cavity-enhanced charge transfer [38], which might suggest that the observed interplay between linear dispersion and bending mode could be of wider relevance.Surely, this simplified picture of modulated reactivity can only hold if vibrational energy redistribution into other modes is unlikely -a condition that is rarely fulfilled and explains why the effect disappears quickly.
It seems intuitive to assume that the here observed discrepancy, supposedly scaling with g 2 0 , should only matter for sizeable coupling strengths, and the change in rate constant should approach the cavity-free reference value monotonously for decreasing coupling.However, this would either imply that the ML+MD calculation would need to qualitatively change from catalysis to inhibition while approaching k g0=0 , which is not observed in Fig. 4 (bottom), or that the QEDFT calculations should become catalysing during this process, for which there was no indication in Ref. 13.One might thus draw the conclusion that the dynamic electronic polarization, missing in ML+MD, is required for the qualitative prediction of chemical changes via single-molecule vibrational strong coupling at most frequencies and for all coupling values.Furthermore, such features might even play a role in collectively coupled systems where individual molecules can exhibit large dynamic electronic polarization in response to collective vibrational dynamics [70,71].SI Sec.II H includes the same investigation performed in this section using a second NEP model based on DFT data calculated with a smaller electronic basis-set.The overall trend up to 700 fs (as in Ref. 13) is consistent with the here presented ML+MD calculations.Since both ML+MD investigations show a consistent qualitative deviation from the QEDFT calculations, they both support the argument of a dynamic electronic polarization.In addition, the smaller basis-set results in a smaller reaction barrier which reduces the significance of the catalysed trajectories at longer times.The most characteristic feature remaining after 2 ps is a strong inhibition around 200 cm −1 , consistent with Fig. 2C.This emphasizes that estimates of kinetic quantities require sufficient sampling.
Sampling, calculating, and learning potentials for every frequency and coupling, in order to obtain the full cavity Born-Oppenheimer surface, is in principle possible.If we intend to widely scan those parameters or would like to include more than a single mode, it becomes, however, practically unfeasible.A (self-consistent) perturbative correction [84,85,89] to the electronic quantities might suffice to suppress some of the problematic features.Such a correction could be constructed via static electronic polarizabilities, which can be predicted with an additional neural network [79], to correct the electronic energies and dipole moments according to the adiabatic cavity field.We envision to combine such a corrected treatment with a full GPU implementation of our approach, which would dramatically reduce the required computational time and pave the way towards explicit ensembles, solvents, and accurate estimates of the chemical kinetics.

IV. CONCLUSION
Recent years have seen a rapid development of theoretical models for the description of strong light-matter coupling, and polaritonic chemistry in particular.Most theoretical work focuses on model systems, which is natural given the young age of the field, and yet poses a major challenge as overly simplified models are unable to truly connect to experiment.Here, we illustrate a combination of ab initio trained machine learning models and modular cavity molecular dynamics that aims to describe realistic molecules in realistic optical environments.This work allows, for the first time, to directly relate theoretical predictions to the experimentally measured changes in chemical kinetics.
We describe theoretically the appearance of singlemolecule strong coupling and its influence on the Si-C bond for the experimentally investigated S N 2 reaction [14].A clear frequency dependence of the rate constant, a critical aspect of polaritonic chemistry, is observed and translates into changes in enthalpy and entropy that are consistent with experimental observations.Interestingly, we observe inhibiting and catalysing effects for the same reaction, the latter of which stand in contrast to previous ab initio calculations [13].In total, three different regimes can be identified that are set apart by differences in the kinetic changes.(i) A strongly inhibiting effect without clear vibrational contribution results in a strong increase in enthalpy but relatively small increase in entropy.(ii) Vibrationally supported catalysis predominantly increases the entropy and only slightly lowers the enthalpy, which results effectively in a simple shift in the Eyring plot, but we emphasize, that the shorter reaction time likely results in an overestimation of this effect.(iii) Vibrationally supported inhibition raises the enthalpy and results in the strongest increase in entropy.The latter observation is qualitatively consistent with experiment and suggests a slight change in the chemical character of the reaction.Vibrationally supported rate changes (ii-iii) are accompanied by a more pronounced change in normal-mode occupation in optically active domains, suggesting that the microscopic mechanism is caused by a stronger interplay of optically active modes via the cavity, i.e., by cavity mediated changes in the redistribution of vibrational energy.
The discrepancy with ab initio calculations, although sharing comparable patterns when scanning the cavity frequency and comparing identical ensembles, suggests that dynamic changes in electronic polarization induced by nuclear motion and mediated by the cavity play a considerable role at the selected coupling strength.This might explain why many simplified molecular dynamics simulations aiming to understand vibrational strong coupling have been able to capture some frequency dependence but often showed strong detuning, or simultaneous catalysing and inhibiting features for the same reaction, which, to best of our knowledge, is in conflict with current experimental work.Our work demonstrates therefore the importance of ab initio QED in the future of polaritonic chemistry and emphasizes the significance of theory that is tailored to describe experimentally relevant reactions.Future development will focus on (self-consistent) perturbative corrections and the full transfer of the established framework to graphical processor units (GPUs), thus providing access to realistic (optical) environments and explicit description of solute-solvent ensembles.The latter is essential for investigating potential modification of solvation dynamics induced by strong coupling [18,88,90].
Polaritonic chemistry remains an equally fascinating and puzzling domain of research.While major questions are yet to be answered [2][3][4]91], especially the connection between local chemistry and collective coupling as well as the interplay with solvation, the continuous growth of theoretical methodology and additional experiments draw an optimistic picture for the future of polaritonics.Our work adds to this a new facet and a clear perspective for possible future development, providing valuable insight that can be experimentally validated.

SUPPORTING INFORMATION
Information on training procedure and performance estimates for the NEP models, Numerical and Simulation details, Derivation and Implementation of optical forces, NVT calculations with corresponding estimates for enthalpy and entropy, Details for kinetic estimates from NVE calculations (including Arrhenius plots), Details and results for rate and Si-C distance estimates for consistency check, Details for normal mode projections and Si-C stretching contributions, Units, A second evaluation of the consistency check based on a different NEP model trained from DFT calculations using the smaller 6-31G* basis, Extensive information and data on the calculation of vibrational frequencies (intermediate and transition state).

ACKNOWLEDGMENTS
We thank Zheyong Fan, Magnus Rahm, Stefano Corni, and Göran Johansson for insightful discussions.C.S. acknowledges support from the Swedish Research Council through Grant No. 2016-06059 and funding from the Horizon Europe research and innovation program of the European Union under the Marie Sk lodowska-Curie grant agreement no.101065117.J.F., E.L., and P.E.acknowledge funding from the Knut and Alice Wallenberg foundation through Grant No. 2019.0140,funding from the Swedish Research Council through Grant No. 2020-04935 as well as the Swedish Foundation for Strategic Research via the SwedNESS graduate school (GSn15-0008).The computations were enabled by resources provided by the National Academic Infrastructure for Supercomputing in Sweden (NAISS) at NSC, PDC, and C3SE partially funded by the Swedish Research Council through grant agreement no.2022-06725.
Partially funded by the European Union.Views and opinions expressed are, however, those of the author(s) only and do not necessarily reflect those of the European Union or REA.Neither the European Union nor the granting authority can be held responsible for them.
We use an active learning approach in this work.An initial training set includes unrelaxed structures from the potential energy surface published in Ref. [1] and additional structures generated by rattling the minimum energy PTAF − structure.We then used the ORCA code version 5.0 [2] (PBE, 6-31G * basis) to obtain energies, forces, and dipoles.We trained a first version of the NEP model, using the NEP-3 potentials in GPUMD [3].With this first generation model, we performed molecular dynamic simulations and selected trajectories with bad performance, adding those structures to the training set and repeating the procedure for a total of 7 generations.It should be noted that GPUMD has been undergoing changes since those initial attempts.
The comparably small 6-31G * basis provided reasonable energies but showed limited reliability for the dipole moments.We prepared a new dataset that included all rattling and molecular dynamics structures as well as a randomly sampled set of structured from the potential energy surface.We performed DFT PBE def2-TZVP calculations (using tight SCF convergence) to generate in total 20170 structures to train the final dipole NEP model with the following parameters.The energy/force model has been trained on the same dataset with the input parameters.Both models are sufficiently converged and their performance is illustrated in the following sections.

B. Initial performance estimates
Let us start our model evaluation with the simple scatter plot Fig. 1 for the expected and predicted energies.The model is well converged within its training set.Additional tests for an independent test-set follow later.
We compared our dipole model against the established symmetry adapted Gaussian Process Regression (SA-GPR) employed by TENSOAP [4] for which we randomly selected a set of 993 structures.Fig. 2 shows dipole predictions for GPUMD and TENSOAP.SA-GPR is considerably slower and practical application for molecular dynamics is limited as it requires at each step to build kernel elements between the test and training set.That said, they require typically less data.For this reason, we decided to use less data for SA-GPR which keeps the training time and memory requirements low and allows for somewhat comparable evaluation times, would one perform molecular dynamics with both models.
Our NEP model performs overall well and is quicker to evaluate than SA-GPR, making it the more convenient choice for our purpose.This short discussion is not suited to make a general claim about the superiority of one of the models but merely serves as sanity check among the existing approaches.

C. Validation of the NEP models compared to ab initio calculations
To validate our NEP model against data outside the training set, we select 99 configurations for various trajectories, at various times, from our MD simulations performed in SI Sec.II D. This provides an estimate for the quality of the model in general and the specific reliability of our reruns from which we draw the conclusion that electronic polarization might play a more relevant role than currently thought.For each configuration we evaluate the potential and kinetic energies of the electronic system, the dipole, and the total force (both electronic and cavity contributions) on the Si-C bond using the NEP model, and using ORCA (Figure 3).NEP model and ORCA are overall in good agreement, and we can expect that our NEP model is able to accurately reproduce the trajectories that one would obtain from ab initio MD.There is one outlier among the data for the kinetic energy featuring a large Si-C bond distance (bond is broken).ML models are excellent at interpolation but configurations that push beyond the tightly sampled do- main will suffer from reduced accuracy.However, those configurations do not play any role in our analysis since the bond is considered as broken once it crosses a Si-C distance of 3.5 Å, i.e., any data generated for larger values do not enter into our analysis.
We perform additional consistency checks between the DFT codes ORCA and NWChem with the def2-TVZP basis set (Figure 4) and for two different basis sets def2-TVZP and 6-31G* within NWChem (Figure 5).ORCA and NWChem are in good agreement, besides a small and irrelevant constant shift in the potential energy.However, the difference between the basis sets is larger than the difference between NEP model and reference calculations, which provides further evidence for the quality of our model in the relevant domain.

D. Numerical Details
All ASE calculations use the Velocity Verlet propagator with a time-step of 0.1 fs.We obtained the Jacobian of the dipole moment with a 2nd-order central-difference approximation using displacements h = 10 −4 Å, implemented in the Python package calorine [5].This value is located at the beginning of a stable plateau illustrated in Fig. 6, i.e., smaller steps are not useful.
The ensemble averaged energy loss during propaga- tion is given in Fig. 7.It increases with increasing frequency since the ratio g/ω is kept constant and a larger cavity frequency leads thus to stronger cavity induced forces which in turn increase the accumulated error due to the finite-difference approximation.We plan to extend GPUMD with analytic derivatives in the future which would entirely mitigate this error.
Our NVE calculations set a temperature by sampling initial velocities from a Boltzmann distribution and removing the center of mass momentum.Ref. [1] showed that considering a solvent resulted in notable changes of the relative infrared activity of vibrational excitations but the good agreement in enthalpy (using NVT conditions) suggests that the effect on the reaction is small enough to draw relevant conclusions from our investigations.

II. SUPPLEMENTARY INFORMATION
A. Obtaining total forces from electronic forces and dipoles Splitting the Hamiltonian in matter Ĥ0 and lightmatter component Supplementary Figure 5. Parity plots for ab initio calculations using the def2-TVZP versus 6-31G* basis sets.Calculations were performed using NWChem.Comparison of (a) potential energies of the electronic system (we have subtracted −22 214 eV from the values), (b) kinetic energies of the electronic system, (c) dipole moments and (d) the (electronic + cavity) force acting on the Si-C bond projected on the bond vector.and introducing the simplified classical limit q, p, μ → q, p, µ, we obtain the classical Hamilton function The careful reader will notice that the classical nuclear limit does not strictly imply μ → µ since the total dipole moment includes electronic and nuclear contributions.The electronic remainder, especially (ε c • μe ) 2 + ε c • μe ε c • µ n , polarizes the electronic system and thus influences the nuclear forces.We discuss the potential consequences of this subtlety and the formally correct treatment in terms of the cavity Born-Oppenheimer [6] in further detail in the main text.Following classical Hamilton mechanics for the canonical displacement mode of the cavity oscillator provides, due to the equivalence of kinetic and canonical momentum in Power-Zienau-Wooley gauge ∂ t q = {q, H LM } = ∂H ∂p = p, the mode-resolved Maxwell equation where the trapezoidal rule is used to approximate the integrals.
We emphasise that µ and F PES can be obtained either directly from DFT, or our NEP models, but only the latter is computationally tractable for MD simulations.

B. NVT reference calculation
Theoretically predicted rates require sufficient time for thermalization to reach a statistically meaningful distribution near an equilibrium state of the system.This requires long propagation times and suggests the use of NVT conditions.Both aspects are problematic in calculations involving the cavity for two major reasons.First, the interplay between thermostat and cavity might give rise to spurious features that misguide our interpretation.Second, observing an appreciable number of reactions under such conditions requires long propagation times.The latter is not an issue for calculations on the GPU, but the current CPU based cavity calculator is certainly limited in this aspect.
Fig. 8 presents an example of the number of reactant molecules over time at 400K.We observe two characteristic domains, an initial burst of reactions, and a second domain for the properly thermalized reactant.We extract the rates as linear fits to the latter and calculate a transition-state enthalpy of ∆H ‡ = 0.345 eV from the Eyring plot in Fig. 9, which is in agreement with experimental estimates ∆H ‡ = 35 ± 4 kJ/mol.[7] It should be noted that we only estimate the Si-C breaking reaction step here, as the correct fluoride attacking and final protonation steps would require an explicit treatment of the solvent.

C. Rate calculations and thermodynamics from NVE
A reaction event took place when the Si-C bond is stretched beyond a value of 3.5 Å.This is located about 0.5 Å behind the transition state, allowing us to account IR transmission (arb.units) Supplementary Figure 11.Top: Rate for the unidirectional reaction P T AF − → F tM eSi + P A − using the 30 initial configurations used in Ref. [1] and g/ω = 1.132, but propagated with our NEP based molecular dynamics calculator.Transmission spectrum obtained from Octopus at 0K, harmonic approximation (red solid), and using our NEP model and GPUMD at 400K NVE conditions (blue dashed).Vertical lines indicate characteristic features observed in Ref. [1].Middle: Trajectory and time-averaged Si-C distance using all 30 trajectories for different time-intervals.Bottom: Trajectory and time-averaged Si-C distance for the subset of reactive trajectories for different time-intervals.This corresponds to the observable shown in Ref. [1].The relatively clear correlation with the rate supports the reliability of the averaged Si-C distance as indicator for rate changes.A sanity check for the difference between middle and bottom picture is, that the average Si-C distance at ωc = 0 (outside cavity) is clearly higher for the subset of reactive trajectories which increases their average value.The overall differences are small, rate and Si-C distance show a strong correlation.
from the transition-state curvature should be correctly accounted for.Future work should investigate how the cavity Born-Oppenheimer approximation modifies our observation, as we would expect then a closer agreement with previous model and ab initio calculations [1,10].
E. Normal mode occupations Fig. 12 illustrates the difference in normal mode occupation for 3 different choices of cavity frequency.The normal mode occupation o is calculated by normalizing the force extracted from the trajectory at a given time f(t) = F(t)/∥F(t)∥ 2 and projecting it onto the orthornormal set of normal mode forces o Notice that the y-axis uses the frequency and we stretch the normal-mode occupation accordingly.The overall structure for ω c = 856 cm −1 is comparable to Ref. [1].For comparison, one should take into account, that the experimentally most relevant normal mode is located at 849 cm −1 in Ref. [1] while it is located at 831 cm −1 when using our NEP potential and the ASE internal vibrational mode calculator (0 K). Noticeable is the rather unaffected region around 750-1000 cm −1 from which only the optically active modes at 770 and 831 cm −1 stand out.The domain between 450-550 cm −1 is quite pronounced (similar to Ref. [1]).

F. Si-C stretching contribution in vibrational modes
Since the reactive step involves breaking the Si-C bond, the contribution of Si-C stretchings in the normal modes is an important indicator for the expected impact on the reaction when energy is redistributed between optically active modes.Figure 13 demonstrates that, in agreement with Fig. 2C from the main manuscript, the Si-C stretching contributions are located foremost in the energy-window between 160 cm −1 and 840 cm −1 , with the exception of the C=C bond stretching around 1200 cm −1 and a high-energy mode beyond above 2000 cm −1 .Nonetheless, this supports our argumentation that the reactive modes that could be effected of dynamic electronic polarization are localized below 840 cm −1 which explains why changes at higher frequencies disappear.

G. Unit conversion between atomic units and ASE units
The dimensionless ratio g 0 /ω c determines our coupling strength.It should be noted that using the implemented cavity-calculator requires the coupling λ = 1 versions between units) which can be directly related to Ref. [1] by taking

H. NEP model using a smaller electronic basis
As detailed in Sec.I A, we started training a second NEP model based on DFT calculations (using ORCA) with the smaller 6-31G* basis set.Dipole moments and stretched configurations are less reliable when using the small 6-31G* basis, such that the following presentation should be consider with caution.As shown in Figure 14, the trained NEP models based on the 6-31G* calculations accurately reproduce the forces obtained from ORCA.
The smaller basis results in a lower reaction barrier and thus a quick saturation of the limited set of trajectories.Nonetheless, we can use those NEP models to investigate if the dynamic electronic polarization remains an important components resulting in comparable deviations.
Figure 15 presents reaction rate constants and Si-C distances obtained from starting ML+MD calculation using the 6-31G* trained NEP models.Two aspects require additional discussion: 1.The overall shape is consistent with Figure 11 for the average change in Si-C distance in the first 700 fs (blue lines in middle and bottom plot).The catalysing character around 450/cm is smaller and seems to be broader, such that the domain around 800/cm obtains additional weight, better emphasizing the resonant dependence.The overall effect remains small at larger frequencies.
2. All frequencies besides 200/cm quickly reach comparable number of products as obtained at the catalysing frequencies explained in the first aspect, i.e., the lower barrier allows the non-catalyzed trajectories to catch up.The only feature that is clearly standing out after 2 ps is the strong inhibition at 200/cm (in Si-C averages and rate).
The basic qualitative behavior of both ML+MD investigations is consistent in the first 700 fs.However, the quick increase at all frequencies (besides 200/cm) implies that only few features remain to play a role after longer propagation time or broader sampling, i.e., only the features reported in Fig. 2C of the paper can be expected to remain relevant after proper sampling.Nonetheless, that the qualitative trend of both ML+MD investigations (up to 700 fs) is consistent and qualitatively contradicts the QEDFT calculations (in the same 700 fs domain) supports our conclusions in the main paper.

I. Vibrational frequencies
Figure 16 presents the normal-mode frequencies and their difference obtained at the intermediate state (blue) and the transition-state (orange) when using our NEP model or ORCA (both based on def2-TZVP basis).The modes are sorted by their frequency.
The ORCA transition-state (TS) geometry was obtained by reoptimizing an initial gues using the "OptTS" flag.The hessian was calculated at the beginning of the optimization to ensure reliable convergence.Using our NEP model, the TS and its frequencies have been calculated by using the final TS structures from ORCA, fixing the positions of the two relevant Si-C atoms, performing a BFGS optimization to reduce forces down to a maximum value of 10 −12 , and calculating the vibrational modes of this relaxed TS structure.The lowest TS frequency of NEP model and ORCA calculation are with 69.5 and 73.38 cm −1 in close agreement.
In addition, we provide below the explicit normal-mode frequencies.

FIG. 1 .
FIG.1.Methodological flow-chart.Density functional theory (DFT) is used to calculates energies, forces, and dipoles.Positional information of a structure is translated into descriptor space and neuroevolution potential models for the potential energy surface (PES) and the dipole moment vector µ are trained using supervised learning (left-hand side).The final models are combined to compute the effective forces acting on the nuclei that are then used to propagate the system in time (right-hand side).

FIG. 2 .
FIG.2.Modification of reaction rate constants by coupling to the cavity.(A) Temperature-dependent infrared spectrum outside the cavity, projected on the cavity polarization.(B) Difference in power spectra of the Si-C bond for PTAF − coupled to the cavity with g0/ωc = 1.132 at 400 K vs the same system outside the cavity.Gray dotted lines (ℏω = 290, 520, 770, 856 cm −1 ) serve as guides to the eye.(C) Rate constant (black dots) and standard error (black bars) for the unidirectional reaction PTAF − → FtMeSi + PA − at 400 K and g0/ωc = 1.132.The rate constant is calculated as number of products after 2 ps.Transmission spectrum at 400 K (red line, defined as negative of absorption spectrum).Vertical magenta-colored dashed lines indicate relevant frequencies used for the kinetic estimates in (D).(D) Eyring plot in free space and for relevant cavity resonances given in inverse cm with g0/ωc = 1.132.The extracted change in enthalpy and entropy is collected in TableI.Negative enthalpy originates from the selected initial state and should only be interpreted relative to the equilibrium enthalpy (see text).The rate constant is calculated as number of products Np per ps.

)FIG. 4 .
FIG.4.Impact of dynamic electronic polarization.(Top) Black dots: Change in rate constant compared to freespace for the unidirectional reaction PTAF − → FtMeSi+PA − using the 30 initial configurations used in Ref.[13] and g0/ωc = 1.132, but propagated with our NEP based molecular dynamics calculator.Blue stars: Change in average Si-C distance obtained from the QEDFT calculations in Ref.[13] amplified by a factor 50. Transmission spectrum obtained from DFT at 0K with harmonic approximation consistent with Ref.[13] (blue dashed), and using our NEP model and GPUMD at 400K NVE conditions (red solid).Vertical lines indicate characteristic features observed in Ref.[13].We show in SI Sec.II D that Si-C distance and rate constant are clearly correlated, i.e., the qualitative trend of QEDFT and ML+MD can be set in relation.(Bottom) ML+MD calculations repeated for different fundamental light-matter coupling strength.

1 FIG. 5 .
FIG. 5. Frequency dependence of reaction bursts.Fourier transform of the autocorrelation function of changes in the number of products |F⟨ Ṅp(t) Ṅp(0)⟩| at 400 K and g0/ωc = 1.132where Np is the number of products.The inset shows the time-dependent accumulation of products at low frequencies, in the strongly catalysing domain, and at high cavity frequencies where the 117 cm −1 mode attains a pronounced role during the reaction (see text).

*
Electronic address: christian.schaefer.physics@gmail.comSupplementary Figure1.Simple scatter-plot for the energy model.We use here all training structures but show additional validation checks in the following sections.

Supplementary Figure 2 .
Quality assessment of the prediction for the x,y, and z component (left to right) of the GPUMD and TENSOAP model.We use all 20170 structures for the GPUMD scatter plot.The SA-GPR model of TENSOAP used 800 structures for training and 193 for the here shown scatter plot.

Supplementary Figure 3 .
Parity plots for our NEP model versus ab initio calculations performed with ORCA and the def2-TVZP basis set.Comparison of (a) potential energies of the electronic system (we have subtracted −22 214 eV from the values), (b) kinetic energies of the electronic system, (c) dipole moments and (d) the (electronic + cavity) force acting on the Si-C bond projected on the bond vector.

Supplementary Figure 4 .
Parity plots for ab initio calculations performed with ORCA versus NWChem.The def2-TVZP basis was used in both cases.Comparison of (a) potential energies of the electronic system (we have subtracted −22 214 eV from the values), (b) kinetic energies of the electronic system, (c) dipole moments and (d) the (electronic + cavity) force acting on the Si-C bond projected on the bond vector.

Supplementary Figure 14 .
Parity plots for the alternative NEP model trained on the 6-31G* basis set versus ab initio calculations performed with ORCA and the 6-31G* basis set.Comparison of (a) potential energies of the electronic system (we have subtracted −22 206.5 eV from the values), (b) kinetic energies of the electronic system, (c) dipole moments and (d) the (electronic + cavity) force acting on the Si-C bond projected on the bond vector.

Supplementary Figure 15 .
Rate and Si-C averages corresponding to Figure11but obtained using the 6-31G* trained NEP models.Top: Rate for the unidirectional reaction P T AF − → F tM eSi + P A − using the 30 initial configurations used in Ref.[1]  and g/ω = 1.132, but propagated with our NEP based molecular dynamics calculator.Transmission spectrum obtained from Octopus at 0K, harmonic approximation (red solid), and using our NEP model and GPUMD at 400K NVE conditions (blue dashed).Vertical lines indicate characteristic features observed in Ref.[1].Middle: Trajectory and time-averaged Si-C distance using all 30 trajectories for different time-intervals.Bottom: Trajectory and time-averaged Si-C distance for the subset of trajectories used in Figure11.

TABLE I .
Change in enthalpy and entropy compared to freespace Eyring result (black data in Fig.2D).
is inhibited in our ML+MD cal-FIG.3.Change in mode occupation due to coupling to cavity.Accumulated (top) and accumulated absolute (bottom) difference in normal mode occupation for different cavity frequencies (given in cm −1 ) with g0/ωc = 1.132 vs free space at 400 K.The corresponding differences in occupation and numerical details are presented in the SI.We quantify the relative changes in cross-correlation between the IR spectrum and accumulated (absolute) difference in TableII.Gray circles show the contribution of Si-C stretching to the vibrational normal mode (see text and SI Sec.II F).Notice that we average here over a time-domain of 2 ps, and not 0.7 ps as in Ref.13.