Effect of Many Modes on Self-Polarization and Photochemical Suppression in Cavities

The standard description of cavity-modified molecular reactions typically involves a single (resonant) mode, while in reality the quantum cavity supports a range of photon modes. Here we demonstrate that as more photon modes are included, physico-chemical phenomena can dramatically change, as illustrated by the cavity-induced suppression of the important and ubiquitous process of proton-coupled electron-transfer. Using a multi-trajectory Ehrenfest treatment for the photon-modes, we find that self-polarization effects become essential, and we introduce the concept of self-polarization-modified Born-Oppenheimer surfaces as a new construct to analyze dynamics. As the number of cavity photon modes increases, the interplay between photon emission and absorption inside the increasingly wide bands of these surfaces, together with their deviations from the cavity-free Born-Oppenheimer surfaces, leads to enhanced suppression. The present findings are general and will have implications for the description and control of cavity-driven physical processes of molecules, nanostructures and solids embedded in cavities.

The interaction between photons and quantum systems is the foundation of a wide spectrum of phenomena, with applications in a range of fields. One rapidlyexpanding domain is cavity-modified chemistry, by which we mean here nuclear dynamics concomitant with electron dynamics when coupled to confined quantized photon modes [1][2][3][4]. The idea is to harness strong light-matter coupling to enhance or quench chemical reactions, manipulate conical intersections, selectively break or form bonds, control energy, charge, spin, or heat transfer, and reduce dissipation to the environment, for example. This forefront has has been strongly driven by experiments [2,[5][6][7][8][9][10][11], with theoretical investigations revealing complementary insights [4,[12][13][14][15][16][17][18][19][20][21][22][23][24][25][26][27][28][29][30][31]. However, apart from a handful of exceptions [32][33][34][35][36][37][38] the simulations of cavity-modified chemistry largely involve coupling to only one (resonant) photon mode, and the vast majority uses simple model systems for the matter part. The modeling of realistic cavity set-ups requires coupling to multiple photon modes that are resonant in the cavity even if they are not resonant with matter degrees of freedom, and further, the description should account for losses at the cavity boundaries. Some strategies have been put forward to treat quantized field modes in the presence of dispersive and absorbing materials [39][40][41][42][43] and theories have been developed to treat many modes and many matter degrees of freedom [14,27,30,32,[34][35][36][37][38]44]. So far unexplored however is a demonstration of how the cavity-modified electronic-nuclear dynamics that were simulated using a single loss-less mode change as one increases the number of photon modes.
Molecules coupled to multiple photon modes represent high-dimensional systems for which accurate and computationally efficient approximations beyond model systems are needed. To this end, the Multi-Trajectory Ehrenfest (MTE) approach for light-matter interaction has been recently introduced [33,34], and benchmarked for two-or three-level electronic systems in a cavity. Wigner-sampling the initial photonic state to properly account for the vacuum-fluctuations of the photonic field while using mean-field trajectories for its propagation, this method is able to capture quantum effects such as spontaneous-emission, bound photon states and second order photon-field correlations [33,34]. In particular, as the trajectories are not coupled during their time-evolution the algorithm is highly parallelizable. Therefore, due to the simplicity, efficiency, and especially scalability the MTE approach for photons emerges as an interesting alternative or extension to other multi-mode treatments [27,30,32,34,36,37,42,44]. 1 In this work, we extend the MTE approach to cavitymodified chemistry, and observe for the first time (to our knowledge) the effect that accounting for many photon modes has on coupled electron-ion dynamics. We focus on the process of polaritonic suppression of an important and ubiquitous process in chemistry and biology, the proton-coupled electron transfer [45], finding the electron-nuclear dynamics significantly depends on the number of modes, as sketched in Fig. 1. We neglect (for now) any effects from cavity losses so we can isolate effects purely from having many modes in the cavity rather than a single mode. To validate the MTE treatment of photons, we first study the single-mode case for which exact results can be computed, finding that MTE performs well but tends to underestimate the photon emission and cavity-induced effects. We explain why using the exact factorization approach [46]. Treating also the nuclei classically gives reasonable averaged dipoles, and photon numbers, but a poor nuclear density, as expected. Turning to multi-mode dynamics computed from MTE, we find that as the number of cavity modes increases (without changing the coupling strength), the suppression of both proton transfer and electron transfer significantly increases, the electronic character becomes more mixed throughout, the photon number begins to increase, and the photon frequency acquires a small but growing Lamb-like shift. This suggests that single-mode simulations tend to underestimate the cavity-induced effects on dynamics in realistic cavities. The self-polarization term [19,47,48] in the Hamiltonian that is often neglected in the literature, has an increasingly crucial impact on the dynamics, and we introduce the concept of self-polarizationmodified Born-Oppenheimer (spBO) surfaces as an instructive tool for analysis of chemical processes mediated by cavity-coupling.

A. Self-Polarization-Modified BO Surfaces
Potential energy surfaces play a paramount role in analyzing coupled electron-dynamics: we have Born-Oppenheimer (BO) surfaces for cavity-free dynamics, Floquet [49,50] or quasistatic [51,52] surfaces for molecules in strong fields, cavity-BO [18] or polaritonic surfaces [13] for molecules in cavities and the exact-factorization based time-dependent potential energy surface [45,53,54] for all cases that yields a complete single-surface picture. The surfaces so far explored for molecules in cavities have largely neglected the selfpolarization term, which is typically indeed negligible for single-mode cavities except at ultra-strong coupling strengths [19]; its importance in obtaining a consistent ground-state and maintaining gauge-invariance has also been emphasized [47,48]. In the multi-mode case however, there is a sum over modes in this term that can become as important as the other terms in the Hamiltonian, and, as we shall see below, it cannot be neglected, especially becoming relevant for large modenumbers, contributing forces on the nuclei as the total dipole evolves in time. Therefore, to analyze the dynamics, we define self-polarization-modified Born-Oppenheimer (spBO) surfaces SP BO (R), as eigenvalues  BO . Further, we define 1-photon-spBO surfaces by simply shifting the spBO surfaces uniformly by the energy of one photon, ω α . These can be viewed as approximate (self-polarization modified) polaritonic surfaces, becoming identical to them in the limit of zero coupling. For small non-zero coupling the polaritonic surfaces, defined as eigenvalues ofĤ −T n , whereT n denotes the nuclear kinetic term, resemble the n-photon-spBO surfaces when they are well-separated from each other, but when they become close, the crossings become avoided crossings.
The top middle panel of Fig. 2 shows the spBO surfaces (pink) for the case of a single photon mode at frequency ω α = 0.1 a.u. coupled to our molecule, along with the 1-photon-spBO surfaces (black). Our model molecule consists of one electronic and one nuclear degree of freedom, with the Hamiltonian given in the Methods section, and we truncate the electronic Hilbert space to the lowest two BO-surfaces throughout this paper. For one mode at the coupling strength of λ = 0.005 a.u. (see Methods) the spBO surfaces coincide with the BO surface, i.e. the self-polarization energy is negligible [45].
As the number of cavity modes grows, the spBO surfaces begin to strongly deviate from the BO surfaces. We consider cavities with resonant modes at frequencies ω α = 0.1 + απc L with α = {− M 2 · · · M 2 } with M the number of modes ranging from 0 (single mode), to 10,40,200,440 and L = 50µm is the cavity-length. The black curves in the top panel of Fig. 3 indicate the corresponding spBO surfaces, and clearly show an in-creasing departure from the BO surfaces as the number of photon modes increases. Given that the landscape of such surfaces provides valuable intuition about the nuclear wavepacket dynamics, with their gradients supplying forces, this suggests an important role of the self-polarization term in the dynamics of the nuclear wavepacket, as we will see shortly.
With more modes, the 1-photon-spBO surfaces begin to form a quasi-continuum: band-like structures indicated by the shaded colors in the top row of Fig. 3. The shading actually represents parallel surfaces separated by the mode-spacing 0.00045 a.u. (We note that, as a function of cavity-length, the mode-spacing decreases, approaching the continuum limit as L approaches infinity, however the coupling strength λ α also decreases, vanishing in the infinite-L limit such that the free BO surfaces are recovered). The 1-photon ground-spBO band and 1-photon excited-spBO band show growing width and increasing overlap as the number of photon modes increases, suggesting a nuclear wavepacket will encounter an increasing number of avoided crossings between ground-and excited-polaritonic states as it evolves. Note that as lower frequencies are included in the band, n-photon-spBO (n ≥ 2) states will overlap with the 1-photon-spBO band. For simplicity however we will still refer to these as simply 1-photon-spBO bands with the understanding that they may include some 2-photon and higher-photon-number states for low frequencies. We return to the implications of the spBO bands later in the discussion of the multi-mode cases.

B. Single-Mode Benchmark
First we consider the single-mode case for which we are able to compare the MTE method (see Methods for details) to numerically exact results 2 . The central photon frequency of 0.1 a.u. is chosen to coincide with the BO energy difference at R = −4 a.u., which is where we launch an initial Gaussian nuclear wavepacket on the excited BO surface. We take the initial state as a simple factorized product of the photonic vacuum state ξ 0 (q) for each mode, the excited BO state, and the nuclear Gaussian wavepacket: where q denotes the vector of photonic displacement-field coordinates.
The top panel of Fig. 2 shows the electronic wavefunctions at R = −4 a.u. (left) and R = 4 a.u. (right) in the cavity-free case, showing that the transition of the initial wavepacket to the lower BO surface through non-adiabatic coupling near the avoided crossing re- (2) (r) (2) (r) sults in an electron transfer. Hence the molecule models proton-coupled electron transfer. Ref. [45] found that this proton-coupled electron transfer is suppressed when the molecule is placed in a single-mode cavity resonant with the initial energy difference between the BO surfaces.
The second row of Fig. 2 shows the dynamics of the nuclear wavepacket (see also supplementary materials, movie 1) for the exact cavity-free case (pink), exact single-mode case (black), MTE for photons (blue) and MTE for both photons and nuclei (light blue). As discussed in Ref. [45], the exact dynamics in the cavity shows suppression of proton-coupled electron transfer (compare pink and black dipoles in third panel), due to photon emission at early times (black line in panel (d)) yielding a partially trapped nuclear wavepacket, leading to less density propagating to the avoided crossing to make the transition to the lower BO surface. The BOpopulations in the lowest panel (e) show the initial partial drop to the ground-state surface associated with the photon emission.
Both MTE approaches are able to approximately capture the cavity-induced suppression of the protoncoupled electron transfer, as indicated by the blue and light-blue dipoles and photon-number in panels (b-d), and approximate the BO occupations in panel (e) reasonably well. However both approaches somewhat underestimate the suppression; the photon emission is underestimated by about a third, as is the suppression of the electronic dipole transfer, for example. To understand why, we compare the potentials the MTE photons experience to the exact potential acting on the photons as defined by the exact factorization approach, which was presented in Ref. [46]. In this approach, the total wavefunction of a system of coupled subsystems is factorized into a single product of a marginal factor and a conditional factor, and the equation for the marginal satisfies a Schrödinger equation with potentials that exactly contain the coupling effects to the other subsystem. When the photonic system is chosen as the marginal, one obtains then the exact potential driving the photons, and this was found for the case of an excited two-level system in a single resonant mode cavity in Ref. [46]. It was shown that the potential develops a barrier for small q-values while bending away from an upper harmonic surface to a lower one at large q, creating a wider and unharmonic well. This leads then to a photonic displacement-field density with a wider profile in q than would be obtained via the uniform average of harmonic potentials that underlie the MTE dynamics, i.e. MTE gives lower probabilities for larger electric-field values, hence a smaller photon-number and less suppression compared to the exact.
An additional treatment of the nuclei within MTE yields a spreading of the nuclear wave packet instead of a real splitting ( Fig. 2(a.3)), a well-known problem of Ehrenfest-nuclei. This error is less evident in averaged quantities such as dipoles and BO coefficients.
Having now understood the limitations of MTE, we now apply the MTE framework for photons to the multimode case.

C. MTE Dynamics for Multi-Mode Cases
The top panel of Fig. 3 shows the ground and excited 1-photon spBO bands. As we observed earlier, including more photon modes has two effects on the spBO surfaces. First, the self-polarization morphs them away from the cavity-free BO surfaces, increasing their separation, and what was a narrow avoided crossing in the cavity-free case shifts leftward in R with increased separation. Second, the 1-photon ground and excited spBO bands both broaden with increasing number of crossings with the 0-photon spBO surfaces and with each other in the regions of overlap. As the gradient of these surfaces and the couplings between them are considerably altered, we expect significant differences in the nuclear dynamics when going from the singlemode case to the many-mode case. Indeed, this is reflected in the middle panel of Fig. 3 (panel (b)) and photon number (panel (c)). The corresponding R-resolved BO-occupations of the ground-BO electronic state divided by the nuclear density, |c 1 (R, t)| 2 (as defined in Methods), shown in Fig.4(a), and the R-averaged occupations |C 1,2 (t)| 2 over time plotted in Fig.4(b) also show significant mode-number dependence (A movie is also provided in supplementary materials, movie 2).
Going from a single-mode (black in Fig. 2) to 10modes (green), the spBO surfaces are only slightly distorted from the BO surfaces, but there is an enhancement of the suppression, since the 1-photon ground-spBO band contains 10 surfaces each with slightly shifted crossings with the 0-photon excited spBO surface on which the wavepacket is initially; these crossings become avoided crossings once the matter-photon coupling is accounted for, i.e. in the polaritonic surfaces. This enhances the probability of photon emission (panel (c)) into the narrow band of cavity-modes. This is reflected also in the narrow frequency band of panel (a) in Fig. 5 which provides a spectral decomposition of the occupied photon modes as a function of time. The increased photon emission corresponds to a larger portion of the nuclear wavepacket (panels (a) of Fig. 3) being trapped in the ground electronic state to the left of the avoided crossing than in the single-mode case, while the right-going part continues on the upper electronic-surface. Still, as there is only little distortion of the spBO surfaces, these two branches of nuclear wavepacket follow closely the two branches of the single-mode dynamics. A larger trapped portion of the wavepacket clearly leads to a smaller nuclear dipole moment at larger times but also a smaller electron transfer: the final electron transfer is largely due to the splitting at the electron-nuclear avoided crossing at around R = 2 a.u. to which less nuclear density has reached. The R-resolved BO-occupation of the ground-BO state in Fig.4(a) show that the electronic character throughout the nuclear wavepacket is similar to the single-mode case, especially after the initial interaction region, which is maybe less obvious to discern from the R-averaged occupations in Fig. 4(b) that gives the overall picture from the electronic side over time.
Turning now to the 40-mode case (orange), the distortion of the spBO states from the BO increases, with the avoided crossing shifting a little leftward and widening slightly. Although the overall dynamics follow the 10-mode case closely, the now broadened one-photon bands lead to more and faster initial photon emission compared to the 10-mode case (Fig. 3(c)) , which is also reflected in the more mixed character of the R-resolved BO ground-state population at early times (orange in panel (a.1) in Fig. 4 ). The combined effects of increased early transitions to the electronic ground spBO state and a slightly less sharp electron-nuclear non-adiabatic region, leads to a little more of the nuclear wavepacket being trapped on the left side of the avoided crossing and a reduced electron-transfer, as shown by the electronic and nuclear dipoles and the BO-occupations. A notable difference between the 10-and 40-mode cases is in the spectral decomposition of the occupied photon modes in Fig. 5(b), where a small Lamb-like shift is evident [55] with the center of the dominant band slightly shifting from 0.1 a.u. to 0.102 a.u.. It is important to note that the calculated photon number reflects both a propagating photon (photon emission) as well as a quasi-bound photon component; the latter arises from dressed photon-matter eigenstates where pure BO states get coupled through the counter rotating-wave terms and molecular dipole terms in the Hamiltonian, and, as perturbation theory suggests, grows as the number of photon modes increases. The photon number α a † α a α that we plot in Fig. 4 and its spectrally-resolved version in Fig. 5 do not distinguish these.
In the 200-mode case (red), the self-polarization term distorts the spBO surfaces significantly and shifts the electron-nuclear non-adiabatic region to be centered near R = 0 a.u.. This shift and widening weakens the non-adiabatic coupling at the electron-nuclear avoided crossing significantly, which suggests that the population transfer from the upper to the lower BO state at this crossing would be much reduced, which is in fact the case (see the very gentle slope in panel Fig. 4(a.3)). However, there is actually very little nuclear density reaching this crossing due to the extended overlap of the excited spBO surface and the 1-photon ground spBO band. As a result the photon number rapidly increases from the beginning and almost immediately there is a mixed electronic character throughout the nuclear density, as reflected in Figs. 3(c) and 4(a). The flatter slope of the excited spBO surface together with the increased population in the lower spBO surface (Fig. 4.a), greatly slows the nuclear density down compared to the fewer-  Fig. 3. The photon number continues to grow slightly throughout the evolution. This can also be seen in the spectral decomposition in Fig. 5(c), where for initial times we find a wide band with dominant occupation around a frequency of 0.12 a.u. (inset), which represents the initial spBO energy-difference (0.1 a.u.) with a small Lamb-like shift. As time evolves we see a continual re-absorption and re-emission (yielding the slight constant increase of the photon number) into a wider band building around 0.07 a.u.. Since the energy-difference between the spBO surfaces near the turning point of the nuclear dipole R ≈ −1 a.u. is about 0.05a.u., this frequency can be interpreted as the central frequency of transition between the excited and ground spBO surfaces in the region where the nuclear wavepacket is moving the slowest with the Lamb-like shift again on the order of 0.02a.u..
The 440-mode case (blue) leads to an even stronger suppression of proton-coupled electron transfer. The key feature causing this is the strong deviation of the excited spBO surface such that its gradient slopes back to the left soon after the initial nuclear wavepacket slides down from its initial position at R = −4 a.u., in contrast to the cavity-free excited BO surface. The overlap of the extensively broadened 1(n)-photon-excitedand 1(n)-ground-bands increases significantly creating a near-continuum of avoided crossings. The 0-photon surfaces are everywhere surrounded by near-lying nphoton surfaces with the upper parts of both bands now reaching up into higher energies and the lowest part of the 1-photon ground-state band reaching the fundamental cavity mode of frequency πc/L. Compared to the 200-mode case, even less density reaches the region of closest approach of the two (0-photon) spBO surfaces, which is now even wider. The slopes of bands results in an even slower nuclear dynamics, with the nuclear and electronic dipole returning to their initial positions after only a small excursion away, as evident in Fig. 3. By including the lowest allowable cavity modes, we find a significant increase of the photon number due to the population of low frequency photons. This can be seen in spectral decomposition in Fig.5(d), as we find bright bands rapidly developing at lower frequencies. As expected, we find a larger Lamb-like shift at 0.1276 a.u. at initial times (inset, a cut through the heat map at t = 1.9 fs). However, due to the densely-spaced 0, 1, 2...photon surfaces we observe a quite fast re-absorption and re-emission of photons into a broad band yielding the larger constant increase of the photon number.
Finally, to emphasize the importance of the selfpolarization term on the dynamics, in Fig. 6 we compare the results of the MTE dynamics on the electronic and nuclear dipoles and photon number when this term is neglected (dashed) or included (solid) for 10, 40, 200 and 440 modes. Here we find only small differences for the 10 mode case, however, as anticipated from the discussion above, including more photon modes leads to larger differences in the dynamics. More precisely, the very initial photon emission remains the same with and without self-polarization. However, as more photon modes are accounted for, there are larger deviations as time evolves, especially for the 440 mode case, which yields quantitative deviations up to a factor of 2. The differences in the dynamics are distinct for the electronic and nuclear dipole, where already for the 10 mode case deviations up to 0.3 a.u. (electronic) and 0.2 a.u. (nuclear) are found at later times. The error in neglecting self-polarization becomes especially significant for the 200-and 440-mode cases, where there is qualitatively different behavior in the nuclear dipole. In the 200-mode case, the differences reach 5.8 a.u. for nuclear dipole and 1.8 a.u. for electronic dipole. Indeed, neglecting the self-polarization term leads to an increase of the proton transfer compared to the singlemode case, in contrast to the increased suppression observed when including the self-polarization. Therefore, neglecting the self-polarization term for many photon modes does not only change the quantitative results dramatically, but can also result in overall different physical effects. The nuclear and electronic wavepackets in the 440-mode case becomes delocalized over the entire region, so plotting simply the dipole, an averaged quantity, appears to give more agreement with the self-polarization-neglected dynamics, when in fact the wavepackets look completely different (see also supplementary material, compare movie 2. and 3.).

I. DISCUSSION AND OUTLOOK
Our results suggest that the effect of multiple cavitymodes on the reaction dynamics can be dramatically different than when only a single mode is accounted for. This is particularly true when there are cavity-modes resonant with the matter system. In particular, for the model of cavity-induced suppression of proton-coupled electron transfer investigated here, we find an overall increase of the suppression the more photon modes are accounted for. Two mechanisms are fundamentally responsible for the difference: First, the self-polarization term grows in significance with more modes with the effect that self-polarization-modified BO surfaces are distorted significantly away from their cavity-free shape. Polaritonic surfaces, eigenvalues of H − T n , should include the explicit matter-photon coupling on top of these spBO surfaces. Second, the n-photon-spBO bands become wider and increasingly overlapping, yielding a very mixed electronic character and continual exchange between surfaces. These new dressed potential energy surfaces provide a useful backdrop to analyze the dynamics, and will form a useful tool in analyzing the different surfaces put forward to study coupled photonmatter systems, for example the polaritonic surfaces, and especially the time-dependent potential energy sur-face arising from the exact factorization as this single surface alone provides a complete picture of the dynamics.
The MTE treatment of the photons appears to be a promising route towards treating realistic light-matter correlated systems. In particular, this method is able to capture quantum effects such as cavity-induced suppression of proton-coupled electron transfer, yet overcomes the exponential scaling problem with the number of quantized cavity modes. However, a practical approach for realistic systems will further need an approximate treatment of the matter part. From the electronic side TDDFT would be a natural choice, while a practical treatment of nuclei calls for a classical treatment such as Ehrenfest or surface-hopping in some basis. However, the multiple-crossings inside the n-photon spBO bands suggest that simple surface-hopping treatments based on spBO surfaces should be used with much caution and that decoherence-corrections should be applied, for example those generalized from the exact factorization approach to the electron-nuclear problem [56,57]. Further, the MTE approach could provide a way to accurately approximate the light-matter interaction part of the QEDFT xc functional [4,27,30,44].
Finally, we note that the present findings are general in that the increasing importance of self-polarization with more photon modes is expected to hold for the description and control of cavity-driven physical processes of molecules, nanostructures and solids embedded in cavities in general. These findings could yield a new way to control and change chemical reactions via the self-polarization without the need to explicitly change the light-matter coupling strength itself.

A. Hamiltonian
Here we consider the non-relativistic photon-matter Hamiltonian in the dipole approximation in the Coulomb gauge as [4,18,30,46,58] with the Hamiltonian for the matter in the cavity aŝ Our model is in one dimension, with one electronic coordinate r and one nuclear coordinate R, where the nuclear and electronic kinetic termsT n = − 1 2M ∂ 2 ∂R 2 ,T e = − 1 2 ∂ 2 ∂r 2 , whileĤ SP BO denotes the spBO Hamiltonian, defined by adding the self-polarization term, to the usual BO Hamiltonian. The self-polarization term depends only on matter-operators but scales with the sum over modes of the squares of the photon-matter coupling parameters λ α ; a thorough discussion of this term can be found in Ref. [19,47,48]. Atomic units, in which = e 2 = m e = 1, are used here and throughout. The photon Hamiltonian and photon-matter coupling read as followŝ where α denotes the photon modes,q α = α 1 2ωα (â † α +â α ) is the photonic displacementfield coordinate, related to the electric displacement operator, whilep α is proportional to the magnetic field. We choose the matter-photon coupling strength through the 1D mode function λ α = 2 L 0 sin(k α X) where L denotes the length of the cavity and k α = απ/L the wave vector, and X the total dipole. Here we take X = L/2, assuming that the molecule is placed at the center of the cavity, and that L = 50µm is much longer than the spatial range of the molecular dynamics. In our particular model the matter potentialV m is given by the 1D Shin-Metiu model [59][60][61], which consists of a single electron and proton (Z = 1 above), which can move between two fixed ions separated by a distance L in one-dimension. This model has been studied extensively for both adiabatic and nonadiabatic effects in cavity-free [60][61][62][63] and in-cavity cases [18,45,64]. The Shin-Metiu potential is: We choose here L = 19.0 a.u., a + = 3.1 a.u., a − = 4.0 a.u., a f = 5.0 a.u., and the proton mass M = 1836 a.u.; with these parameters, the phenomenon of proton-coupled electron transfer occurs after electronic excitation out of the ground-state of a model molecular dimer [45].

B. MTE Treatment of Photonic System
A computationally feasible treatment of coupled electron-ion-photon dynamics in a multi-mode cavity calls for approximations. Here we have one electronic and one nuclear degree of freedom but up to 440 photon modes, so we use MTE for the photons, coupled to the molecule treated quantum mechanically. As mentioned above we take the initial state as a simple factorized product of the photonic vacuum state ξ 0 (q) for each mode, the excited BO state, and the nuclear Gaussian wavepacket. More precisely, for the MTE for photons we sample the initial photonic vacuum state from the Wigner distribution given by: ξ 0 (q, p) = α 1 π e − p 2 α ωα −ωαq 2 α . Furthermore, with two electronic surfaces, the equations of motion are as follows, for the lth trajectory: i∂ t C 1 (R, t) C 2 (R, t) = h 11 h 12 h 21 h 22 with the diagonal matrix elements and for i = j, Here the non-adiabatic coupling terms are d ij ij (R) = Φ BO R,i |∂ 2 R Φ BO R,j , and the transition dipole and quadrupole terms r (n) ij = Φ BO R,i |r n |Φ BO R,j . The coefficients C i (R, t) are the expansion coefficients of the electron-nuclear wavefunction in the BO basis: Ψ(r, R, t) = i=1,2 C i (R, t)Φ BO R,i (r). Subsequently the R-resolved and R-averaged BO-populations are defined as |c 1,2 (R, t)| 2 = |C 1,2 (R, t)| 2 /|χ(R, t)| 2 and |C 1,2 (t)| 2 = dR|C 1,2 (R, t)| 2 respectively. In the singlemode case we also present the results for when the proton is also treated by MTE with the nuclear trajectory satisfying MR l (t) = − ∂ R BO (R l ) − α ω α λ α q l α − α λ 2 α Z(Z R l − r l ) . For the photonic system, 10,000 trajectories were enough for convergence for all cases except the 440-mode case which required 50,000 trajectories (the results shown used 100,000 trajectories in all cases).

Data Availability
The data that support the findings of this work are available from the corresponding authors on reasonable request.

ACKNOWLEDGMENTS
We would like to thank Johannes Feist for insightful discussions. Financial support from the US National Science Foundation CHE-1940333 (NM) and the Department of Energy, Office of Basic Energy Sciences, Division of Chemical Sciences, Geosciences and Biosciences under Award de-sc0020044 (LL) are gratefully acknowledged. NMH gratefully acknowledges an IMPRS fellowship. This work was also supported by the European Research Council (ERC-2015-AdG694097), the Cluster of Excellence (AIM), Grupos Consolidados (IT1249-19) and SFB925 "Light induced dynamics and control of correlated quantum systems. The Flatiron Institute is a division of the Simons Foundation.