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 accounted for, physicochemical 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 increasing deviation of these surfaces from the cavity-free Born–Oppenheimer surfaces, together with the interplay between photon emission and absorption inside the widening bands of these 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.

## I. INTRODUCTION

The interaction between photons and quantum systems is the foundation of a wide spectrum of phenomena, with applications in a range of fields. One rapidly expanding domain is cavity-modified chemistry, by which we mean here nuclear dynamics concomitant with electron dynamics when coupled to confined quantized photon modes.^{1–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, and heat transfer; and reduce dissipation to the environment, for example. This forefront has been strongly driven by experiments,^{2,5–11} with theoretical investigations revealing complementary insights.^{4,12–30} However, apart from a handful of exceptions,^{31–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 setups requires coupling to multiple photon modes that are supported in the cavity even if they are not resonant with matter degrees of freedom, and furthermore, 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–43} and theories have been developed to treat many modes and many matter degrees of freedom.^{14,27,29,31,33–37,44} So far unexplored, however, is an explicit demonstration of how the cavity-modified electronic-nuclear dynamics change as one increases the number of photon modes in the simulation.

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^{32,33} 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 classical 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.^{32,33} In particular, as the trajectories are not coupled during their time-evolution, the algorithm is highly parallelizable. Therefore, due to the simplicity, efficiency, and scalability, the MTE approach for photons emerges as an interesting alternative or extension to other multi-mode treatments.^{27,29,31,33,35,36,42,44,45}

In this work, we extend the MTE approach to cavity-modified chemistry and point out the effect that accounting for many photon modes has on coupled electron–ion dynamics. We focus on the process of polaritonic suppression of the proton-coupled electron transfer,^{46} finding that 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.^{47} 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, the suppression of both proton transfer and electron transfer significantly increases (without changing coupling strength), the electronic character becomes more mixed throughout, and the photon number begins to increase. The results suggest that even when cavity modes are far from the molecular resonances, the chemical properties of the molecule can be dramatically altered by the presence of the cavity even when the coupling strength is not particularly large. The self-polarization term^{19,48,49} in the Hamiltonian that is often neglected in the literature has an increasing impact on the dynamics, and we analyze the error made when neglecting it, depending on the number of photon modes accounted for. To this end, we introduce the concept of self-polarization-modified Born–Oppenheimer (spBO) surfaces as an instructive tool for analysis of chemical processes mediated by cavity-coupling.

## II. HAMILTONIAN

In this work, we consider the non-relativistic photon-matter Hamiltonian in the dipole approximation in the length gauge, i.e., applying the Power–Zienau–Woolley gauge transformation^{50} on the minimal coupling Hamiltonian in the Coulomb gauge, as^{4,18,29,47,51,52}

with the Hamiltonian for the matter in the cavity as

Our model is in one dimension, with one electronic coordinate *r* and one nuclear coordinate *R*, where the nuclear and electronic kinetic terms $T^n=\u221212M\u22022\u2202R2$, $T^e=\u221212\u22022\u2202r2$, while $\u0124BOSP$ 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 *λ*_{α}; thorough discussions of this term can be found in Refs. 19, 48, and 49. Atomic units, in which *ℏ* = *e*^{2} = *m*_{e} = 1, are used here and throughout. The photon Hamiltonian and photon–matter coupling read as follows:

where *α* labels the photon mode, $q^\alpha =12\omega \alpha (\xe2\alpha \u2020+\xe2\alpha )$ is the photonic coordinate, related to the electric field operator, while $p^\alpha =\u2212i\omega \alpha 2(\xe2\u2212\xe2\u2020)$ is proportional to the magnetic field. It is important to note that in this gauge, the photon number is given by

We choose the matter–photon coupling strength through the 1D mode function $\lambda \alpha =2L\u03f50sin(k\alpha X)$, where $L$ denotes the length of the cavity, $k\alpha =\alpha \pi /L$ denotes the wave vector (*α* = 1, 2, 3, …), and *X* denotes the position measured from the center of the cavity. Unless stated otherwise, we take $X=L/2$, assuming that the molecule is placed at the center of the cavity, and $L=12.5\u2009\mu m$, much longer than the spatial range of the molecular dynamics. This cavity-length yields a coupling strength of $\lambda \alpha =(\u22121)\alpha \u2212120.01$ a.u. for modes with odd *α*, *λ*_{α} = 0 for even *α*, and a fundamental cavity mode of frequency *ω*_{0} = 0.0018 a.u., and these parameters are used throughout this paper except for in benchmarking the single-mode results in Sec. IV.

In our particular model, the matter potential $V^m$ is given by the 1D Shin–Metiu model,^{53–55} 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^{54–57} and cavity cases.^{18,46,58} 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.^{46} Furthermore, for computational convenience in the MTE calculations, we truncate the electronic Hilbert space to the lowest two BO-surfaces.

## III. SELF-POLARIZATION-MODIFIED BO SURFACES

Potential energy surfaces play a paramount role in analyzing coupled dynamics: we have Born–Oppenheimer (BO) surfaces for cavity-free dynamics; Floquet^{59,60} or quasistatic^{61,62} surfaces for molecules in strong fields; cavity-BO^{18} or polaritonic surfaces^{13} for molecules in cavities; and exact-factorization based time-dependent potential energy surfaces^{46,63,64} for all cases that yield a complete single-surface picture. The surfaces so far explored for molecules in cavities have largely neglected the self-polarization term, which is often indeed negligible for typical single-mode calculations. However, its importance in obtaining a consistent ground-state and maintaining gauge-invariance has been emphasized.^{48,49} The self-polarization term involves a sum over the number of photonic modes considered.^{65} In the multi-mode case, this sum 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 mode-numbers, contributing forces on the nuclei while the total dipole evolves in time. To analyze the dynamics, we define self-polarization-modified Born–Oppenheimer (spBO) surfaces $\u03f5BOSP(R)$ as eigenvalues of the spBO Hamiltonian $\u0124BOSP\Phi R,BOSP=\u03f5BOSP(R)\Phi R,BOSP$.

Furthermore, we define “*n*-photon-spBO surfaces” by simply shifting the spBO surfaces uniformly by *nℏω*_{α}. We note that this nomenclature should not be taken too literally since in the length gauge, the photon number includes matter-coupling and self-polarization terms on top of the Coulomb-gauge definition of ⟨*a*^{†}*a*⟩, as shown in Eq. (6). That is, a 1-photon-spBO surface does not actually denote a surface where there is one photon in the system, for example. The spBO surfaces 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 $\u0124\u2212T^n$, 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 (0-photon) spBO (pink) and 1-photon spBO (black) surfaces for the single mode case, where the spBO and BO surfaces essentially coincide. The top panel of Fig. 3 shows in black the spBO surfaces for our system for 10, 30, 50, and 70 modes. For 10 modes, they show only a small deviation from the BO surfaces with a small widening and shift of the avoided crossing region. As the number of cavity modes grows, the spBO surfaces clearly show an increasing departure from the BO surfaces. 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.

The band-like structures indicated by the shaded colors in the top row of Fig. 3 represent the 1-photon-spBO surfaces, forming a quasi-continuum. The shading actually represents parallel surfaces separated by the mode-spacing of 0.0018 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 that a nuclear wavepacket will encounter an increasing number of avoided crossings between ground- and excited-polaritonic states as it evolves. It is worth noting that the 1-photon-spBO band overlaps with the (*n* > 1)-photon-spBO band of the lower frequencies, e.g., the 10-photon-spBO curve for the fundamental mode coincides with the 1-photon-spBO curve for the 10th mode. For simplicity, however, we will still refer to these as simply 1-photon-spBO bands with the understanding that they may include some higher-photon-number states for low frequencies. For clarity of the figure, we show only the 1-photon-spBO band, but we note that the (*n* > 1)-photon bands also play a role in the dynamics, in particular, when there is an overlap between the (*n* > 1)-photon spBO ground state and the spBO excited state. We return to the implications of the spBO bands later in the discussion of the multi-mode cases.

Finally, as mentioned above, the spBO surfaces do not incorporate the bilinear light–matter interactions, which if included would turn crossings of the spBO surfaces into avoided crossings of self-polarization modified polaritonic surfaces. Computing these for a large number of photon modes results in a large diagonalization problem. Instead, the spBO surfaces depend on only the matter operators and light–matter coupling strength and so could in principle be computed with a similar computational expense as for BO surfaces while giving already an indication of how chemistry is modified in the cavity, as we will see shortly.

## IV. 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 will consider one electronic and one nuclear degree of freedom but up to 70 photon modes, so we use MTE for the photons, coupled to the molecule treated quantum mechanically.

We launch an initial Gaussian nuclear wavepacket on the excited BO surface at *R* = −4 a.u. 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, $\Psi (r,R,q\u0332,0)=Ne\u2212[2.85(R+4)2]\Phi R,2BO(r)\xi 0(q\u0332)$, where $q\u0332$ denotes the vector of photonic displacement-field coordinates. More precisely, for the MTE for photons, we sample the initial photonic vacuum state from the Wigner distribution given by $\xi 0(q,p)=\u220f\alpha 1\pi e\u2212p\alpha 2\omega \alpha \u2212\omega \alpha q\alpha 2$. Furthermore, with two electronic surfaces, the equations of motion are as follows, for the *l*th trajectory:

with the diagonal matrix elements

and for *i* ≠ *j*,

Here, the non-adiabatic coupling terms are $dij(R)=\u27e8\Phi R,iBO|\u2202R\Phi R,jBO\u27e9$, $dij(2)(R)=\u27e8\Phi R,iBO|\u2202R2\Phi R,jBO\u27e9$, and the transition dipole and quadrupole terms $rij(n)=\u27e8\Phi R,iBO|r^n|\Phi R,jBO\u27e9$. The coefficients *C*_{i}(*R*, *t*) are the expansion coefficients of the electron–nuclear wavefunction in the BO basis, $\Psi (r,R,t)=\u2211i=1,2Ci(R,t)\Phi R,iBO(r)$. Subsequently, *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 single-mode case, we will also present the results for when the proton is also treated by MTE with the nuclear trajectory satisfying $MR\u0308l(t)=\u2212\u27e8\u2202R\u03f5BO(Rl)\u27e9\u2212\u2211\alpha \omega \alpha \lambda \alpha q\alpha l\u2212\u2211\alpha \lambda \alpha 2Z(Z\u27e8R\u27e9l\u2212\u27e8r\u27e9l)$. For all MTE calculations, 20 000 trajectories were enough for the convergence of the results for all cases.

## V. RESULTS

### A. Single-mode benchmark

First, we consider a single-mode case for which we are able to compare the MTE method to numerically exact results.^{66} The cavity-free dynamics of our system shows “proton-coupled electron transfer” in the following sense: 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 nuclear wavepacket to the lower BO surface through non-adiabatic coupling near the avoided crossing results in an electron transfer. Reference 46 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. This energy difference is 0.1 a.u., which would correspond to about the 56th mode in the cavity with the parameters described in Sec. II. A single mode of frequency equal to the fundamental mode of that cavity is so far off the initial resonance that the dynamics is only slightly altered from the cavity-free case (see Sec. IV). Instead, for the purposes of benchmarking the MTE method for cavity-modified dynamics in this section, we use the same parameters as in Ref. 46: cavity frequency of 0.1 a.u. with light–matter coupling strength of *λ* = 0.005 a.u.

The second row of Fig. 2 shows the dynamics of the nuclear wavepacket (see also supplementary material, 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. 46, 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 BO-populations 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 proton-coupled 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. We note that the photon number is by far dominated by the first term in Eq. (6) [compare black solid and gray dashed line in panel (d)]; there is only a single mode at an initially resonant molecular frequency, and the coupling is small enough that the second and third terms are very small. To understand why MTE underestimates the photon number, 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. 47. 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. 47. 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 anharmonic 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 wavepacket 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. We note that an exact treatment of the photons coupled to MTE for only nuclei will not improve this situation. With more photon modes, the polaritonic landscape has even more avoided crossings, which are likely to make the Ehrenfest description for nuclei worse, calling for the development of more advanced propagation schemes.^{67,68}

Having now understood the limitations of MTE, we now apply the MTE framework for photons to the multi-mode case.

### B. MTE dynamics for multi-mode cases

We return to the cavity-parameters of Sec. II, where the fundamental mode has *ω*_{0} = 0.0018 a.u. and the light–matter coupling strength of *λ*_{α} = ±0.01 a.u. for modes that are non-zero at the center of the cavity (see Sec. II). We consider the effect on the dynamics as an increasing number of harmonics of the fundamental are included in the simulation from 1, 10, 30, 50, to 70 modes. Although in principle all modes should be considered, already these cases demonstrate a dramatic impact of the number of modes on the dynamics.

We note that if we had instead used a cavity whose fundamental mode is 0.1 a.u. as in the single-mode demonstration of Sec. V A, then considering the effect of including higher cavity-modes on the dynamics would be more complicated since one rapidly encounters cavity wavelengths short enough that the long-wavelength approximation is broken. Instead, with the parameters of Sec. II that we use here for the many-mode cases, the maximum frequency included is *ω*_{max} = 0.127 (a.u.) in the 70 mode case, which corresponds to a wavelength of *λ*_{max} = 0.057 (*μ*m), much larger than the spatial range of the molecular dynamics.

Another important aspect when including many photon modes is the well-known zero-point energy leakage problem. However, in all cases considered here, we find none (for the 10–50 mode cases) or extremely small leakage for long times and high frequencies (for the 70 mode case) compared to the overall emission and photon number, with no impact on the dynamics. Still, as more modes are included (beyond 70 modes), we anticipate the zero-point energy leakage could become a problem and would need to be addressed carefully.

The top panel of Fig. 3 shows the ground and excited 1-photon spBO bands, as introduced earlier. As mentioned in Sec. III, we do not show the entire (*n* > 1)-photon spBO bands explicitly for clarity, but it is important to note that the 1-photon band does include some (*n* > 1)-photon states of the lower frequency photon modes that are included in each simulation.

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 the 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 single-mode case to the many-mode case.

Indeed, this is reflected in the middle panel of Fig. 3, which shows the nuclear wavepacket at time snapshots 22 fs (a.1), 30 fs (a.2), and 41 fs (a.3), and in the lower panel, showing the nuclear dipole [panel (b)] and electronic dipole [panel (c)]. The corresponding *R*-resolved BO-occupations of the ground-BO electronic state |*c*_{1}(*R*, *t*)|^{2}, shown in Fig. 4(a), and the *R*-averaged occupations |*C*_{1,2}(*t*)|^{2} over time plotted in Fig. 4(b) also show a significant mode-number dependence (a movie is also provided in supplementary material, movie 2).

Dynamics in the single-mode cavity (black line as computed with MTE and gray dashed line for exact, in Fig. 3) is almost identical to the cavity-free case (pink) since the mode is far off the molecular resonance (*ω* = 0.0018 a.u.) for the duration of the dynamics. Differences are seen when the nuclear wavepacket encounters the avoided crossing region, with the single-mode case slightly lagging behind the cavity-free dynamics and with a smaller transfer to the lower electronic state. First, due to the stronger coupling (*λ* = 0.01 a.u.), compared to the single mode benchmarking, the spBO-surfaces (not shown here) already have a very slight distortion from its original BO-form, with a slightly wider and broader avoided crossing region. As a result, the transfer to the ground electronic state is slightly reduced, as evident in 4(a.3), Figs. 4(b), and 3(c). With more population in the upper state, which slopes to the left after the avoided crossing, the wavepacket slows down compared to the cavity-free case [Figs. 3(a.3) and 3(b)]. The 0-photon spBO surfaces at the closest approach have an energy difference of about 0.006 a.u., so the 1-photon ground-state surface does not interact strongly with the excited spBO surface. The overlap of the four- and higher-photon ground-state surfaces with the excited spBO surface leads to a small photon emission, as shown in Fig. 5. We observe that unlike for the parameters of Sec. V A, the larger self-polarization term results in a significant difference between the true photon number of Eq. (6) and the “pseudo-photon-number,” the first term, even for a single-mode, but due to the low frequency of this mode, there is only a limited impact on the energetics of the matter system.

Going now to the 10-mode case (green in Fig. 3), the spBO surfaces are visibly distorted from the BO-surfaces shown in Fig. 2, and we begin to see suppression of both the proton transfer in panels (a) and (b) and more so the electron transfer in panel (c). The largest cavity-frequency has a value of *ω*_{max} = 0.018 a.u., while at the closest approach, the spBO surfaces differ in energy by 0.01 a.u. with their avoided crossing shifting further left to *R* = 1.2 a.u. The suppression of the molecular dynamics begins to occur a little before the wavepacket approaches the avoided crossing between the spBO-surfaces and is due to two effects: first, the gentler slope and weaker crossing of the spBO surfaces causing a weaker effective electron–nuclear non-adiabaticity and, second, the crossings between the 1-photon ground spBO surfaces with the first excited spBO surface causing photon emission. These crossings become avoided crossings once the matter–photon bilinear coupling is accounted for, i.e., in the polaritonic surfaces. At around *R* = 0.6 a.u., a single photon of *ω*_{max} first becomes resonant to the self-polarization modified molecular excitation, enabling transitions from the 0-photon excited surface to the 1-photon ground spBO-band, which continue also at lower frequencies as the wavepacket proceeds to the right and through the avoided crossing at around *R* = 1.2 a.u. This is also reflected in the mixed character of R-resolved BO-population [see Fig. 4(a.2)], showing an increase in the ground electronic state population before reaching the avoided crossing. The part of the wavepacket already transferred to the ground state before reaching the avoided crossing of the spBO surfaces has to climb a potential hill to pass through, hence less reaching the right side. This effect, together with the weakening of the electron-nuclear nonadiabaticity from the self-polarization term distorting the BO surfaces, yields a suppression of both the electron and proton transfer. The photon number, dominated again by the self-polarization contribution [third term in Eq. (6)], shows a corresponding increase at around 23 fs (Fig. 5).

Turning now to the 30-mode case (orange) with *ω*_{max} = 0.055 a.u., the distortion of the spBO states from the BO increases, with the avoided crossing widening slightly and shifting leftward to *R* = 0.5 a.u. The broadened one-photon bands lead to more and earlier photon emission compared to the 10-mode case (Fig. 5). The highest cavity-mode frequencies included now are resonant with the self-polarization modified molecular resonance already at *R* = −1.3 a.u. However, due to the change in the curvature of the excited spBO surface compared to the BO surface, we observe a clear suppression of the dynamics even before the wavepacket reaches this region, as evident from Fig. 4(a.1). Once the wavepacket approaches *R* = −1.3 a.u., the cavity modes become resonant, enabling transitions from the 0-photon excited spBO surface to the broadened 1-photon ground spBO-band; the narrowing of the spBO energy differences as the wavepacket progresses past this point leads to transitions to 1- and (*n* > 1)-photon ground-spBO surfaces of the lower frequency modes. (Again, many crossings of these surfaces become avoided crossings of the polaritonic surfaces.) In fact, we see even earlier an increase in the *R*-resolved BO ground-state population [orange in panel (a.1) in Fig. 4] for the left part of the wavepacket and an increase in the photon number around 20 fs in Fig. 5. Why this happens to the left of the wavepacket rather than the right (the right is less off-resonant than the left) could be due to a stronger photon–matter coupling there from the larger molecular dipole in the left tail of the wavepacket compared to the leading edge. The right part of the nuclear wavepacket shows a more mixed character of the *R*-resolved BO ground-state population at early times. The combined effects of increased early transitions to the electronic ground spBO state and a slightly less sharp electron–nuclear non-adiabatic region lead to a less amount of the nuclear wavepacket reaching the avoided crossing and a reduced electron–proton transfer dramatically, as shown by the electronic and nuclear dipoles and the BO-occupations.

In the 50-mode case (red), the self-polarization term distorts the spBO surfaces further (Fig. 3), shifting the electron–nuclear non-adiabatic region to be centered near *R* = 0 a.u. The 1-photon bands are wider, with *ω*_{max} = 0.091 a.u., and become resonant with the self-polarization modified molecular resonance already at *R* = −3 a.u. This leads to transitions from the initially 0-photon excited spBO surface to the 1-photon ground spBO surfaces already at very short times. This is evident in the almost immediate mixed character of the *R*-resolved BO populations. 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-mode cases and results in a full suppression of both the proton and electron transfer, as evident from panels (a)–(c) of Fig. 3. The wavepacket reflects before appreciably reaching the avoided crossing, and the change in the spatially averaged BO-population in Fig. 4(b) is caused solely by the nuclear wavepacket dropping into the ground-BO state by emitting photons and is not due to the avoided crossing at *R* = 0 a.u. Considering the photon number shown in Fig. 4, although the free photonic field component in panel (b) has a rapid initial increase and then grows throughout the time evolution as in the few-mode cases, while the linear term has a compensating decrease [panel (c)], the total photon number decreases after some time due to the self-polarization contribution. As expected, this term increases with the number of modes at the initial time, but since it is proportional to the total dipole of the system, whose transfer is suppressed, the resulting photon number tracks this behavior and is ultimately reduced compared with the fewer-mode cases.

The 70-mode case (blue) with *ω*_{max} = 0.127 a.u. can be seen as an enhanced 50-mode case and leads to an even stronger suppression of proton-coupled electron transfer, with the same two key features that have been responsible for the cavity-modified dynamics in the 10- and higher-mode cases now having an even greater impact. First, by including the resonant frequency of the initial position of the nuclear wavepacket, part of the wavepacket almost immediately drops to the lower surface by emitting photons (Fig. 5); see also *R*-averaged BO populations at early times [Fig. 4(b)] and the mixed character developing in the *R*-resolved populations of Fig. 4(a). Second, the deviation of the excited spBO surface is now strong enough 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., sloping back to the left, in contrast to the cavity-free excited BO surface. The overlap of the extensively broadened 1-photon-excited- and 1-ground-bands increases significantly, creating a near-continuum of avoided crossings of polaritonic surfaces. The 0-photon surfaces are everywhere surrounded by near-lying *n*-photon surfaces. Compared to the 50-mode case, even less density reaches the region of the avoided crossing, which is now even wider. The slope of the excited spBO-band results in even slower nuclear dynamics, with the nuclear and electronic dipoles returning to their initial positions after only a small excursion away, as evident in Fig. 3. Analogous to the discussion for the 50-mode case, we find a larger initial total photon number, compared to the 50-mode case, followed by a moderate increase and decrease due to the early reflection of the nuclear wavepacket.

Finally, to emphasize the importance of the self-polarization 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, 30, 50, and 70 modes. For the electronic and nuclear dipoles, already for the 10-mode case, deviations up to 1.2 a.u. (electronic) and 0.4 a.u. (nuclear) are found at later times. The error in neglecting the self-polarization term becomes especially notable for the 50- and 70-mode cases, where including the self-polarization term yields a decrease in the proton-transfer (from 50 modes to 70 modes), while not including the self-polarization term yields an increase in the proton transfer. Therefore, neglecting the self-polarization term for many photon modes not only changes the quantitative results dramatically but can also result in overall different physical effects. In fact, the nuclear and electronic wavepackets in the 70-mode case become 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 the supplementary material; compare movies 2 and 3). Turning to the total photon number, we find that when more photon modes are accounted for, the simulations without the self-polarization term first underestimate (10-mode), then coincide (30-modes), and then overestimate the photon number up to a factor of 2.8 (70-modes), compared to the simulations that include the self-polarization term. This can be explained with the trends of the total dipole moment discussed above since a dominant contribution to the photon number is the self-polarization term in Eq. (6), which depends directly on this.

## VI. DISCUSSION AND OUTLOOK

Our results suggest that the effect of multiple cavity-modes on the reaction dynamics can lead to dramatically different dynamics than in the cavity-free case. This is true even when the cavity-modes are far from the electronic resonances encountered in the dynamics and even more so when cavity-modes are 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 in the suppression when 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 with much 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 photon–matter systems, for example, the polaritonic surfaces, and especially the time-dependent potential energy surface 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 toward 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, time-dependent density functional theory (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. 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.^{67,68} Furthermore, the MTE approach could provide a way to accurately approximate the light–matter interaction part of the Quantum-Electrodynamical Density Functional Theory (QEDFT) exchange-correlation functional.^{4,27,29,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. Extensions of these findings to multi-mode cavities suggest a new possibility of controlling and changing chemical reactions via self-polarization without the need to explicitly change the light–matter coupling strength itself.

## SUPPLEMENTARY MATERIAL

See the supplementary material for the following three movies: (1) **MovieSingleMode.mp4** shows the nuclear dynamics for the cavity free case (pink) and for a coupling to one resonant photon mode with a coupling strength of *λ* = 0.005 and resonant frequency of 0.1 for the full quantum solution (black), quantum nuclei with MTE treatment for the photons (blue), and MTE treatment for nuclei and photons (light blue) in comparison. The surfaces represent the BO surfaces (pink) and spBO surfaces (black) with the coupling to the one resonant photon mode. This movie corresponds to the results given in Fig. 2 in the paper. (2) **MovieMultiModeWithSP.mp4** shows the spBO surfaces (upper panel), the nuclear density (second panel), and the *R*-resolved BO-population (third and last panel) over time. Here, we treat the photons with the MTE-approach and the matter part quantum mechanically, and all simulations include the self-polarization term. Here, we choose the coupling strength to be *λ* = 0.01 for all cases and compare the dynamics of the cavity free case (pink), the coupling to a single resonant mode (black), and the coupling to 10- (green), 30- (orange), 50- (red), and 70- (blue) photon modes. This movie corresponds to Figs. 3 and 4 in the paper. (3) **MovieMultiMode.mp4** shows the same calculations, dynamics, and comparison as (2), however, now without taking the self-polarization into account, and the upper panel now shows the BO-surfaces. This movie gives more details for the discussion of Fig. 6.

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding authors upon reasonable request.

## ACKNOWLEDGMENTS

We would like to thank Johannes Feist for insightful discussions and the anonymous referees for very helpful comments. Financial support from the US National Science Foundation under Grant No. CHE-1940333 (N.T.M.) and the Department of Energy, Office of Basic Energy Sciences, Division of Chemical Sciences, Geosciences, and Biosciences under Award No. DE-SC0020044 (L.L.) is gratefully acknowledged. N.M.H. gratefully acknowledges an IMPRS fellowship. This work was also supported by the European Research Council (Grant No. 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.

## REFERENCES

This includes Quantum-Electrodynamical Density Functional Theory (QEDFT),^{4,27,30,44} which is an exact non-relativistic generalization of time-dependent density functional theory that dresses electronic states with photons and allows the electronic properties of real materials to be retained in a computationally efficient way.

We note that the two-mode case can also be solved exactly numerically, but the single-mode comparison here already illustrates the main points.