Effect of Many Modes on Self-Polarization and Photochemical Suppression in Cavities
Norah M. Hoffmann, Lionel Lacombe, Angel Rubio, Neepa T. Maitra
EEffect of Many Modes on Self-Polarization and Photochemical Suppression in Cavities
Norah M. Hoffmann,
1, 2, ∗ Lionel Lacombe, † Angel Rubio,
2, 3, 4, ‡ and Neepa T. Maitra § Department of Physics, Rutgers University at Newark, Newark, NJ 07102, USA Max Planck Institute for the Structure and Dynamics of Matter and Center for Free-ElectronLaser Science and Department of Physics, Luruper Chaussee 149, 22761 Hamburg, Germany Center for Computational Quantum Physics, Flatiron Institute, 162 5th Avenue, New York, NY 10010, USA Nano-Bio Spectroscopy Group and ETSF, Universidad del Pas Vasco, 20018 San Sebastian, Spain (Dated: July 15, 2020)The standard description of cavity-modified molecular reactions typically involves a single (reso-nant) mode, while in reality the quantum cavity supports a range of photon modes. Here we demon-strate that as more photon modes are accounted for, physico-chemical phenomena can dramaticallychange, as illustrated by the cavity-induced suppression of the important and ubiquitous processof 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 thenumber 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 absorp-tion inside the widening bands of these surfaces, leads to enhanced suppression. The present findingsare general and will have implications for the description and control of cavity-driven physical pro-cesses of molecules, nanostructures and solids embedded in cavities.
The interaction between photons and quantum sys-tems is the foundation of a wide spectrum of phenom-ena, with applications in a range of fields. One rapidly-expanding domain is cavity-modified chemistry, bywhich we mean here nuclear dynamics concomitantwith electron dynamics when coupled to confined quan-tized photon modes [1–4]. The idea is to harness stronglight-matter coupling to enhance or quench chemicalreactions, manipulate conical intersections, selectivelybreak or form bonds, control energy, charge, spin, heattransfer, and reduce dissipation to the environment, forexample. This forefront has has been strongly drivenby experiments [2, 5–11], with theoretical investigationsrevealing complementary insights [4, 12–31]. However,apart from a handful of exceptions [32–39] the sim-ulations of cavity-modified chemistry largely involvecoupling to only one (resonant) photon mode, and thevast majority uses simple model systems for the mat-ter part. The modeling of realistic cavity set-ups re-quires coupling to multiple photon modes that are sup-ported in the cavity even if they are not resonant withmatter degrees of freedom, and further, the descrip-tion should account for losses at the cavity boundaries.Some strategies have been put forward to treat quan-tized field modes in the presence of dispersive and ab-sorbing materials [40–44] and theories have been devel-oped to treat many modes and many matter degrees offreedom [14, 27, 30, 32, 34–38, 45]. So far unexploredhowever is an explicit demonstration of how the cavity-modified electronic-nuclear dynamics change as one in- ∗ [email protected] † [email protected] ‡ [email protected] § [email protected] creases the number of photon modes in the simulation.Molecules coupled to multiple photon modes rep-resent high-dimensional systems for which accurateand computationally efficient approximations beyondmodel systems are needed. To this end, the Multi-Trajectory Ehrenfest (MTE) approach for light-matterinteraction has been recently introduced [33, 34], andbenchmarked for two- or three-level electronic systemsin a cavity. Wigner-sampling the initial photonic stateto properly account for the vacuum-fluctuations of thephotonic field while using classical trajectories for itspropagation, this method is able to capture quantumeffects such as spontaneous-emission, bound photonstates and second order photon-field correlations [33,34]. In particular, as the trajectories are not coupled dur-ing their time-evolution the algorithm is highly paral-lelizable. Therefore, due to the simplicity, efficiency, andscalability, the MTE approach for photons emerges as aninteresting alternative or extension to other multi-modetreatments [27, 30, 32, 34, 36, 37, 43, 45]. [46]In this work, we extend the MTE approach to cavity-modified chemistry, and point out the effect that ac-counting for many photon modes has on coupledelectron-ion dynamics. We focus on the process ofpolaritonic suppression of the proton-coupled electrontransfer [47], finding the electron-nuclear dynamics sig-nificantly depends on the number of modes, as sketchedin Fig. 1. We neglect (for now) any effects from cav-ity losses so we can isolate effects purely from havingmany modes in the cavity rather than a single mode. Tovalidate the MTE treatment of photons, we first studythe single-mode case for which exact results can be com-puted, finding that MTE performs well but tends to un-derestimate the photon emission and cavity-induced ef-fects. We explain why using the exact factorization ap-proach [48]. Treating also the nuclei classically gives a r X i v : . [ c ond - m a t . m e s - h a ll ] J u l easonable averaged dipoles, and photon numbers, buta poor nuclear density, as expected. Turning to multi-mode dynamics computed from MTE, we find that asthe number of cavity modes increases, the suppressionof both proton transfer and electron transfer signifi-cantly increases (without changing coupling strength),the electronic character becomes more mixed through-out, and the photon number begins to increase. The re-sults suggest that even when cavity modes are far fromthe molecular resonances, the chemical properties of themolecule can be dramatically altered by the presence ofthe cavity even when the coupling strength is not par-ticularly large. The self-polarization term [19, 49, 50]in the Hamiltonian that is often neglected in the liter-ature, has an increasing impact on the dynamics, andwe analyze the error made when neglecting it, depend-ing on the number of photon modes accounted for. Tothis end, we introduce the concept of self-polarization-modified Born-Oppenheimer (spBO) surfaces as an in-structive tool for analysis of chemical processes medi-ated by cavity-coupling. I. HAMILTONIAN
In this work we consider the non-relativistic photon-matter Hamiltonian in the dipole approximation in thelength gauge, i.e. applying the Power-Zienau-WoolleyGauge transformation [51] on the minimal couplingHamiltonian in Coulomb gauge, as [4, 18, 30, 48, 52, 53] ˆ H = ˆ H SP m + ˆ H p + ˆ V pm , (1)with the Hamiltonian for the matter in the cavity as ˆ H SP m = ˆ T n + ˆ H SPBO where ˆ H SPBO = ˆ T e + ˆ V m + ˆ V SP . (2)Our model is in one dimension, with one electronic co-ordinate r and one nuclear coordinate R , where the nu-clear and electronic kinetic terms ˆ T n = − M ∂ ∂R , ˆ T e = − ∂ ∂r , while ˆ H SPBO denotes the spBO Hamiltonian, de-fined by adding the self-polarization term, ˆ V SP = 12 M (cid:88) α λ α ( Z ˆ R − ˆ r ) , (3)to the usual BO Hamiltonian. The self-polarization termdepends only on matter-operators but scales with thesum over modes of the squares of the photon-mattercoupling parameters λ α ; thorough discussions of thisterm can be found in Ref. [19, 49, 50]. Atomic units, inwhich (cid:126) = e = m e = 1 , are used here and through-out. The photon Hamiltonian and photon-matter cou-pling read as follows ˆ H p ( q ) = 12 M (cid:88) α (cid:0) ˆ p α + ω α ˆ q α (cid:1) (4) ˆ V pm = M (cid:88) α ω α λ α ˆ q α (cid:16) Z ˆ R − ˆ r (cid:17) , (5) Single Photon Mode(a) spBONuclearDensity
Effect of Many Photon Modes(b) spBOPhotonNuclearDensity
FIG. 1. An exemplary sketch of a molecule coupled to manyphoton modes. (a) Sketches the spBO surfaces and the cor-responding nuclear dynamics for a coupling to a single pho-ton mode. (b) Depicts the effect of many photon modes onthe spBO surfaces and the corresponding complete photo-chemical suppression of the proton-coupled electron transfer. where α labels the photon mode, ˆ q α = (cid:113) ω α (ˆ a † α + ˆ a α ) is the photonic coordinate, related to the electric fieldoperator, while ˆ p α = − i (cid:112) ω α (ˆ a − ˆ a † ) is proportional tothe magnetic field. It is important to note that in thisgauge, the photon number is given by ˆ N p = (cid:88) α (cid:0) ˆ a † α ˆ a α + λ α ˆ q α ( Z ˆ R − ˆ r ) + λ α ( Z ˆ R − ˆ r ) / (2 ω α ) (cid:1) (6)We choose the matter-photon coupling strengththrough the 1D mode function λ α = (cid:113) L (cid:15) sin( k α X ) where L denotes the length of the cavity and k α = απ/ L the wave vector ( α = 1 , , ... ) , and X the position mea-sured from the center of the cavity. Unless stated oth-erwise, we take X = L / , assuming that the moleculeis placed at the center of the cavity, and L = 12 . µm ,much longer than the spatial range of the molecular dy-namics. This cavity-length yields a coupling strengthof λ α = ( − α − . a.u. for modes with odd α , λ α =0 for even α , and a fundamental cavity mode of fre-quency ω = 0 . a.u., and these parameters are usedthroughout the paper except for in benchmarking thesingle-mode results in Sec. III.In our particular model the matter potential ˆ V m isgiven by the 1D Shin-Metiu model [54–56], which con-sists of a single electron and proton ( Z = 1 above),which can move between two fixed ions separated bya distance L in one-dimension. This model has beenstudied extensively for both adiabatic and nonadiabaticeffects in cavity-free [55–58] and in-cavity cases [18, 47,29]. The Shin-Metiu potential is: ˆ V m = (cid:88) σ = ± | R + σL | − erf (cid:16) | r + σL a σ (cid:17) | r + σL | − erf (cid:16) | R − r | a f (cid:17) | R − r | (7)We choose here L = 19 . a.u., a + = 3 . a.u., a − =4 . a.u., a f = 5 . a.u., and the proton mass M =1836 a.u.; with these parameters, the phenomenon ofproton-coupled electron transfer occurs after electronicexcitation out of the ground-state of a model molecu-lar dimer [47]. Furthermore, for computational conve-nience in the MTE calculations we truncate the elec-tronic Hilbert space to the lowest two BO-surfaces. II. SELF-POLARIZATION-MODIFIED BO SURFACES
Potential energy surfaces play a paramount role in an-alyzing coupled dynamics: we have Born-Oppenheimer(BO) surfaces for cavity-free dynamics, Floquet [60, 61]or quasistatic [62, 63] surfaces for molecules in strongfields, cavity-BO [18] or polaritonic surfaces [13] formolecules in cavities and the exact-factorization basedtime-dependent potential energy surface [47, 64, 65]for all cases that yields a complete single-surface pic-ture. The surfaces so far explored for molecules in cav-ities have largely neglected the self-polarization term,which is often indeed negligible for typical single-mode calculations. However, its importance in obtain-ing a consistent ground-state and maintaining gauge-invariance has been emphasized [49, 50]. The self-polarization term involves a sum over the number ofphotonic modes considered [66]. In the multi-modecase, this sum can become as important as the otherterms in the Hamiltonian, and, as we shall see below,it cannot be neglected, especially becoming relevant forlarge mode-numbers, contributing forces on the nucleiwhile the total dipole evolves in time. To analyze thedynamics, we define self-polarization-modified Born-Oppenheimer (spBO) surfaces (cid:15)
SPBO ( R ) , as eigenvalues ofthe spBO Hamiltonian: ˆ H SPBO Φ SP R, BO = (cid:15) SPBO ( R )Φ SP R, BO .Further, we define ” n -photon-spBO surfaces” by sim-ply shifting the spBO surfaces uniformly by n (cid:126) ω α . Wenote that this nomenclature should not be taken too lit-erally since in the length gauge, the photon number in-cludes matter-coupling and self-polarization terms ontop of the Coulomb-gauge definition of (cid:104) a † a (cid:105) as shownin Eq. (6). That is, a -photon-spBO surface does not ac-tually denote a surface where there is one photon in thesystem, for example. The spBO surfaces can be viewedas approximate (self-polarization modified) polaritonicsurfaces, becoming identical to them in the limit of zerocoupling. For small non-zero coupling the polaritonicsurfaces, defined as eigenvalues of ˆ H − ˆ T n , resemble the n -photon-spBO surfaces when they are well-separatedfrom each other, but when they become close, the cross-ings become avoided crossings. The top middle panel of Fig. 2 shows the ( -photon)spBO (pink) and -photon spBO surfaces (black) for thesingle mode case, where the spBO and BO surfaces es-sentially coincide. The top panel of Fig. 3 shows in blackthe spBO surfaces for our system for 10, 30, 50, and 70modes. For 10 modes, they show only a small deviationfrom the BO surfaces with a small widening and shiftof the avoided crossing region. As the number of cav-ity modes grows, the spBO surfaces clearly show an in-creasing departure from the BO surfaces. Given that thelandscape of such surfaces provides valuable intuitionabout the nuclear wavepacket dynamics, with their gra-dients supplying forces, this suggests an important roleof the self-polarization term in the dynamics of the nu-clear wavepacket, as we will see shortly.The band-like structures indicated by the shaded col-ors in the top row of Fig. 3 represent the -photon-spBOsurfaces, forming a quasi-continuum. The shading actu-ally represents parallel surfaces separated by the mode-spacing . a.u. (We note that, as a function of cavity-length, the mode-spacing decreases, approaching thecontinuum limit as L approaches infinity, however thecoupling strength λ α also decreases, vanishing in theinfinite- L limit such that the free BO surfaces are recov-ered). The -photon ground-spBO band and -photonexcited-spBO band show growing width and increas-ing overlap as the number of photon modes increases,suggesting a nuclear wavepacket will encounter an in-creasing number of avoided crossings between ground-and excited- polaritonic states as it evolves. It is worthnoting that the -photon-spBO band overlaps with the ( n > -photon-spBO band of the lower frequencies, e.g.the 10-photon-spBO curve for the fundamental modecoincides with the 1-photon-spBO curve for the 10thmode. For simplicity however we will still refer to theseas simply -photon-spBO bands with the understand-ing that they may include some higher-photon-numberstates for low frequencies. For clarity of the figure, weshow only the -photon-spBO band, but we note that the ( n > -photon bands also play a role in the dynamics,in particular when there is overlap between the ( n > -photon spBO ground state and the spBO excited state.We return to the implications of the spBO bands later inthe discussion of the multi-mode cases.Finally, as mentioned above, the spBO surfaces do notincorporate the bilinear light-matter interactions, whichif included would turn crossings of the spBO surfacesinto avoided crossings of self-polarization modified po-laritonic surfaces. Computing these for a large num-ber of photon modes results in a large diagonalizationproblem. Instead, the spBO surfaces depend on only thematter operators and light-matter coupling strength socould in principle be computed with a similar computa-tional expense as for BO surfaces while giving alreadyan indication of how chemistry is modified in the cavity,as we will see shortly.3 II. MTE TREATMENT OF PHOTONIC SYSTEM
A computationally feasible treatment of coupledelectron-ion-photon dynamics in a multi-mode cavitycalls for approximations. Here we will consider oneelectronic and one nuclear degree of freedom but up to70 photon modes, so we use MTE for the photons, cou-pled to the molecule treated quantum mechanically.We launch an initial Gaussian nuclear wavepacket onthe excited BO surface at R = − a.u. We take the ini-tial state as a simple factorized product of the photonicvacuum state ξ ( q ) for each mode, the excited BO state,and the nuclear Gaussian wavepacket: Ψ( r, R, q,
0) = N e − [2 . R +4) ] Φ BO R, ( r ) ξ ( q ) , where q denotes the vec-tor of photonic displacement-field coordinates. Moreprecisely, for the MTE for photons we sample the ini-tial photonic vacuum state from the Wigner distributiongiven by: ξ ( q, p ) = (cid:81) α π e (cid:20) − p αωα − ω α q α (cid:21) . Furthermore,with two electronic surfaces, the equations of motion areas follows, for the l th trajectory: ¨ q lα ( t ) = − ω α q lα − ω α λ α ( Z (cid:104) R (cid:105) l − (cid:104) r (cid:105) l ) , (8) i∂ t (cid:18) C ( R, t ) C ( R, t ) (cid:19) = (cid:18) h h h h (cid:19) (cid:18) C ( R, t ) C ( R, t ) (cid:19) , (9)with the diagonal matrix elements h ii = (cid:15) i BO ( R ) − M ∂ R + (cid:88) α (cid:16) λ α ω α q lα ( ZR − r ii ( R ))+ λ α · (( ZR ) − ZRr ii ( R ) + r (2) ii ) (cid:17) (10)and for i (cid:54) = j , h ij = − M d ij ( R ) ∂ R − d (2) ij ( R )2 M − (cid:88) α λ α ω α q lα r ij ( R ) + (cid:88) α (cid:18) λ α · (cid:16) − ZRr ij ( R ) + r (2) ij ( R ) (cid:17)(cid:19) (11)Here the non-adiabatic coupling terms are d ij ( R ) = (cid:104) Φ BO R,i | ∂ R Φ BO R,j (cid:105) , d (2) ij ( R ) = (cid:104) Φ BO R,i | ∂ R Φ BO R,j (cid:105) , and the transi-tion dipole and quadrupole terms r ( n ) ij = (cid:104) Φ BO R,i | ˆ r n | Φ BO R,j (cid:105) .The coefficients C i ( R, t ) are the expansion coefficientsof the electron-nuclear wavefunction in the BO basis: Ψ( r, R, t ) = (cid:80) i =1 , C i ( R, t )Φ BO R,i ( r ) . Subsequently, R -resolved and R -averaged BO-populations are definedas | c , ( R, t ) | = | C , ( R, t ) | / | χ ( R, t ) | and | C , ( t ) | = (cid:82) dR | C , ( R, t ) | respectively. In the single-mode casewe will also present the results for when the pro-ton is also treated by MTE with the nuclear trajectorysatisfying M ¨ R l ( t ) = −(cid:104) ∂ R (cid:15) BO ( R l ) (cid:105) − (cid:80) α ω α λ α q lα − (cid:80) α (cid:0) λ α Z ( Z (cid:104) R (cid:105) l − (cid:104) r (cid:105) l ) (cid:1) . For all MTE calculations20,000 trajectories were enough for convergence the re-sults for all cases. RESULTSA. Single-Mode Benchmark
First we consider a single-mode case for which weare able to compare the MTE method to numerically ex-act results [67]. The cavity-free dynamics of our systemshows ”proton-coupled electron transfer” in the follow-ing sense: The top panel of Fig. 2 shows the electronicwavefunctions at R = − a.u. (left) and R = 4 a.u.(right) in the cavity-free case, showing that the transitionof the initial nuclear wavepacket to the lower BO sur-face through non-adiabatic coupling near the avoidedcrossing results in an electron transfer. Ref. [47] foundthat this proton-coupled electron transfer is suppressedwhen the molecule is placed in a single-mode cavity res-onant with the initial energy difference between the BOsurfaces. This energy difference is . a.u. which wouldcorrespond to about the 56th mode in the cavity withthe parameters described in Sec. I. A single mode offrequency equal to the fundamental mode of that cav-ity is so far off the initial resonance that the dynam-ics is only slightly altered from the cavity-free case (seeSec. III). Instead, for the purposes of benchmarking theMTE method for cavity-modified dynamics in this sec-tion, we use the same parameters as in Ref [47]: cavityfrequency of . a.u. with light-matter coupling strength λ = 0 . a.u.The second row of Fig. 2 shows the dynamics ofthe nuclear wavepacket (see also supplementary mate-rials, movie 1) for the exact cavity-free case (pink), ex-act single-mode case (black), MTE for photons (blue)and MTE for both photons and nuclei (light blue). Asdiscussed in Ref. [47], the exact dynamics in the cavityshows suppression of proton-coupled electron transfer(compare pink and black dipoles in third panel), due tophoton emission at early times (black line in panel (d))yielding a partially trapped nuclear wavepacket, lead-ing to less density propagating to the avoided crossingto make the transition to the lower BO surface. The BO-populations in the lowest panel (e) show the initial par-tial drop to the ground-state surface associated with thephoton emission.Both MTE approaches are able to approximatelycapture the cavity-induced suppression of the proton-coupled electron transfer, as indicated by the blue andlight-blue dipoles and photon-number in panels (b–d),and approximate the BO occupations in panel (e) rea-sonably well. However both approaches somewhat un-derestimate the suppression; the photon emission is un-derestimated by about a third, as is the suppression ofthe electronic dipole transfer, for example. We note thatthe photon number is by far dominated by the first termin Eq. (6) (compare black and dashed gray line in panel(d)); there is only a single mode at an initially resonantmolecular frequency, and the coupling is small enoughthat the second and third terms are very small. To un-derstand why MTE underestimates the photon number,4 -4(2) (r)-12 0 12r[a.u.] Φ -4(1) (r) 0.10.20.30.4 -6 -4 -2 0 2 4 6 E n e r g y [ a . u . ] R[a.u.] Φ (r)-12 0 12r[a.u.] Φ (r)00.20.40.60.81 -4 -2 0 2 4 N u c l e a r D e n s i t y [ a . u . ] R[a.u.](a.1) -4 -2 0 2 4R[a.u.](a.2) -4 -2 0 2 4R[a.u.](a.3)-6-4-20246 0 10 20 30 40 E l e c t r o n i c D i p o l e [ a . u . ] time[fs](b) -4-20246 0 10 20 30 40 N u c l e a r D i p o l e [ a . u . ] time[fs](c) 00.10.20.30.40.50.6 0 10 20 30 40 P h o t o n N u m b e r [ a . u . ] time[fs](d)00.20.40.60.81 0 5 10 15 20 25 30 35 40 | C , | [ a . u . ] time[fs](e) FIG. 2. Single-Mode Case: The top panel shows the ground(lower) and excited (upper) BO wavefunctions at R = − a.u.(left) and at R = 4 a.u. (right) and the spBO surfaces (pink)and one-photon spBO surfaces (black). The spBO surfacesare essentially identical with the BO surfaces however for thesingle-mode case. The second panel depicts the nuclear den-sity for cavity-free (pink), full quantum treatment (black), MTEtreatment of the photons only (blue) and MTE treatment ofboth photons and nuclei (light blue) at time snapshots t = 22 fs(a.1) , t = 30 fs (a.2) and t = 38 fs (a.3). The third panel showsthe electronic (b) and nuclear (c) dipole and the photon num-ber (d). In panel (d) the dashed-grey line shows the photonnumber Eq. (6) while the black shows the first term only. Thelowest panel depicts the BO occupations, | C , ( t ) | . we compare the potentials the MTE photons experienceto the exact potential acting on the photons as defined bythe exact factorization approach, which was presentedin Ref. [48]. In this approach, the total wavefunction ofa system of coupled subsystems is factorized into a sin-gle product of a marginal factor and a conditional factor,and the equation for the marginal satisfies a Schr ¨odingerequation with potentials that exactly contain the cou-pling effects to the other subsystem. When the photonicsystem is chosen as the marginal, one obtains then theexact potential driving the photons, and this was found for the case of an excited two-level system in a single res-onant mode cavity in Ref. [48]. It was shown that the po-tential develops a barrier for small q -values while bend-ing away from an upper harmonic surface to a lowerone at large q , creating a wider and anharmonic well.This leads then to a photonic displacement-field densitywith a wider profile in q than would be obtained via theuniform average of harmonic potentials that underliethe MTE dynamics, i.e. MTE gives lower probabilitiesfor larger electric-field values, hence a smaller photon-number and less suppression compared to the exact.An additional treatment of the nuclei within MTEyields a spreading of the nuclear wave packet insteadof a real splitting (Fig. 2(a.3)), a well-known problemof Ehrenfest-nuclei. This error is less evident in av-eraged quantities such as dipoles and BO coefficients.We note that an exact treatment of the photons cou-pled to MTE for only nuclei will not improve this sit-uation. With more photon modes, the polaritonic land-scape has even more avoided crossings, which are likelyto make the Ehrenfest description for nuclei worse, call-ing for the development of more advanced propagationschemes [68, 69].Having now understood the limitations of MTE, wenow 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.I, where thefundamental mode has ω = 0 . a.u. and the light-matter coupling strength λ α = ± . a.u. for modeswhich are non-zero at the center of the cavity (see Sec. I).We consider the effect on the dynamics as an increasingnumber of harmonics of the fundamental are includedin the simulation: from 1, 10, 30, 50, to 70 modes. Al-though in principle all modes should be considered, al-ready these cases demonstrate a dramatic impact of thenumber of modes on the dynamics.We note if we had instead used a cavity whose funda-mental mode is . a.u. as in the single-mode demon-stration of the previous section, then considering theeffect of including higher cavity-modes on the dynam-ics would be more complicated since one rapidly en-counters cavity wavelengths short enough that the long-wavelength approximation is broken. Instead, with theparameters of Sec. I that we use here for the many-modecases, the maximum frequency included is ω max =0 . a.u. ] in the mode case, which corresponds to awavelength of λ max = 0 . µm ] , much larger than thespatial range of the molecular dynamics.Another important aspect when including many pho-ton modes is the well-known zero-point energy leakageproblem. However, in all cases considered here we findnone (for the 10 - 50 mode cases) or extremely smallleakage for long times and high frequencies (for the 70mode case) compared to the overall emission and pho-5on number, with no impact on the dynamics. Still, asmore modes are included (beyond 70 modes) we antici-pate the zero-point energy leakage could become a prob-lem and would need to be addressed carefully.The top panel of Fig. 3 shows the ground and excited -photon spBO bands, as introduced earlier. As men-tioned in Sec. II, we do not show the entire ( n > -photon spBO bands explicitly for clarity, but it is impor-tant to note that the -photon band does include some ( n > -photon states of the lower frequency photonmodes that are included in each simulation.As we observed earlier, including more photon modeshas two effects on the spBO surfaces. First, the self-polarization morphs them away from the cavity-free BOsurfaces, increasing their separation, and what was anarrow avoided crossing in the cavity-free case shiftsleftward in R with increased separation. Second, the -photon ground and excited spBO bands both broadenwith increasing number of crossings with the -photonspBO surfaces and with each other in the regions ofoverlap. As the gradient of these surfaces and the cou-plings between them are considerably altered, we expectsignificant differences in the nuclear dynamics when go-ing from the single-mode case to the many-mode case.Indeed, this is reflected in the middle panel of Fig. 3which shows the nuclear wavepacket at time snapshots fs (a.1), fs (a.2), fs (a.3) and in the lower panel,showing the nuclear dipole (panel b) and electronicdipole (panel (c)). The corresponding R -resolved BO-occupations of the ground-BO electronic state | c ( R, t ) | ,shown in Fig.4(a), and the R -averaged occupations | C , ( t ) | over time plotted in Fig.4(b) also show signif-icant mode-number dependence (A movie is also pro-vided in supplementary materials, movie 2).Dynamics in the single-mode cavity (black as com-puted with MTE and gray dashed for exact, in Fig. 3)is almost identical to the cavity-free case (pink), sincethe mode is far off the molecular resonance ( ω =0 . a.u.) for the duration of the dynamics. Differ-ences are seen when the nuclear wavepacket encoun-ters the avoided crossing region, with the single-modecase slightly lagging behind the cavity-free dynamics,and with a smaller transfer to the lower electronic state.First, due to the stronger coupling ( λ = 0 . a.u.),compared to the single mode benchmarking, the spBO-surfaces (not shown here) already have a very slight dis-tortion from its original BO-form, with a slightly widerand broader avoided crossing region. As a result, thetransfer to the ground electronic state is slightly re-duced, as evident in Fig. 4(a.3, b) and Fig. 3c). Withmore population in the upper state which slopes to theleft after the avoided crossing, the wavepacket slowsdown compared to the cavity-free case (Fig. 3a.3,b). The -photon spBO surfaces at closest approach have anenergy difference of about . a.u. so the -photonground-state surface does not interact strongly with theexcited spBO surface. The overlap of the - and higher-photon ground-state surfaces with the excited spBO sur- E n e r g y [ a . u . ] R[a.u.]10 -6-4-2 0 2 4 6R[a.u.]30 -6-4-2 0 2 4 6R[a.u.]50 -6-4-2 0 2 4 6R[a.u.]700.20.61 N u c l e a r D e n s i t y [ a . u . ] (a.1)0.20.61 (a.2)00.20.61 -6 -4 -2 0 2 4 6R[a.u.] (a.3)-4-2 0 2 4 6 0 10 20 30 40 ( n ) - D i p o l e [ a . u . ] time[fs](b) -8-6-4-2 0 2 4 6 0 10 20 30 40 ( e ) - D i p o l e [ a . u . ] time[fs](c) FIG. 3. The ground- and excited -photon spBO bands, rep-resenting surfaces separated by 0.0018 a.u. (see text) for 10modes (green), 30 modes (orange), 50 modes (red) and 70modes (blue). The middle panel depicts the nuclear density attime snap shots t = 22 fs (a.1), t = 30 fs (a.2) and t = 41 fs (a.3)in the same color code along with the single mode case com-puted within MTE-for-photons (black), exact (gray dashed)and cavity-free (pink) . The lowest panel shows the nucleardipole (b) and electric dipole (c). face leads to a small photon emission as seen in Fig. 5.We observe that unlike for the parameters of the previ-ous section, the larger self-polarization term results ina significant difference between the true photon num-ber of Eq. (6) and the “pseudo-photon-number”, the firstterm, even for a single-mode, but due to the low fre-quency of this mode there is only a limited impact onthe energetics of the matter system.Going now to the 10-mode case (green in Fig. 3), thespBO surfaces are visibly distorted from the BO-surfacesshown in Fig.2, and we begin to see suppression of boththe proton transfer in panels (a) and (b), and more sothe electron transfer in panel (c). The largest cavity-frequency has a frequency ω max = 0 . a.u., while atclosest approach the spBO surfaces differ in energy by . a.u. with their avoided crossing shifting further leftto R = 1 . a.u. The suppression of the molecular dy-6 | c ( R , t ) | (a.1)00.51 | c ( R , t ) | (a.2)00.51 -4 -2 0 2 4 | c ( R , t ) | R[a.u.](a.3) 0 0.2 0.4 0.6 0.8 1 0 5 10 15 20 25 30 35 40 | C , ( t ) | time[fs](b) FIG. 4. Groundstate BO-surface population (a) at time snapshots t = 22 fs (1), t = 30 fs (2) and t = 41 fs (3) over R andthe averaged population over time (c) in the same color codeas Fig.4. namics begins to occur a little before the wavepacketapproaches the avoided crossing between the spBO-surfaces, and is due to two effects: first, the gentlerslope and weaker crossing of the spBO surfaces caus-ing a weaker effective electron-nuclear non-adiabaticity,and second, the crossings between the -photon groundspBO surfaces with the first excited spBO surface caus-ing photon emission. These crossings become avoidedcrossings once the matter-photon bilinear coupling is ac-counted for, i.e. in the polaritonic surfaces. At around R = 0 . a.u., a single photon of ω max first becomes reso-nant to the self-polarization modified molecular excita-tion, enabling transitions from the 0-photon excited sur-face to the -photon ground spBO-band, which continuealso at lower frequencies as the wavepacket proceeds tothe right and through the avoided crossing at around R = 1 . a.u.. This is also reflected in the mixed charac-ter of R-resolved BO-population (see Fig.4(a.2)), show-ing an increase of the ground electronic state popula-tion before reaching the avoided crossing. The part ofthe wavepacket already transferred to the ground statebefore reaching the avoided crossing of the spBO sur-faces has to climb a potential hill to pass through, henceless reaches the right side. This effect, together withthe weakening of the electron-nuclear nonadiabaticityfrom the self-polarization term distorting the BO sur-faces, yields a suppression of both the electron and pro-ton transfer. The photon number, dominated again bythe self-polarization contribution (third term in Eq. 6),shows a corresponding increase at around f s , Fig.5. Turning now to the 30-mode case (orange) with ω max = 0 . a.u., the distortion of the spBO states fromthe BO increases, with the avoided crossing wideningslightly and shifting leftward to R = 0 . a.u.. The broad-ened one-photon bands lead to more and and earlierphoton emission compared to the 10-mode case (Fig. 5).The highest cavity-mode frequencies included now areresonant with the self-polarization modified molecularresonance already at R = − . a.u.. However, dueto the change in curvature of the excited spBO surfacecompared to the BO surface, we observe a clear sup-pression of the dynamics even before the wave packetreaches this region, as evident from Fig. 4.(a.1). Oncethe wave packet approaches R = − . a.u. the cavitymodes become resonant enabling transitions from the -photon excited spBO surface to the broadened -photonground spBO-band; the narrowing of the spBO energydifferences as the wavepacket progresses past this pointleads to transitions to - and ( n > -photon ground-spBO surfaces of the lower frequency modes. (Again,the many crossings of these surfaces become avoidedcrossings of the polaritonic surfaces). In fact we seeeven earlier an increase in the R -resolved BO ground-state population (orange in panel (a.1) in Fig. 4) for theleft part of the wave packet, and an increase in the pho-ton number around fs in Fig 5. Why this happensto the left of the wavepacket rather than the right (theright is less off-resonance than the left), could be due toa stronger photon-matter coupling there from the largermolecular dipole in the left tail of the wavepacket com-pared to the leading edge. The right part of the nuclearwave packet shows a more mixed character of the R -resolved BO ground-state population at early times. Thecombined effects of increased early transitions to theelectronic ground spBO state and a slightly less sharpelectron-nuclear non-adiabatic region, leads to a less ofthe nuclear wave packet reaching the avoided crossingand a reduced electron-proton-transfer dramatically, asshown by the electronic and nuclear dipoles and the BO-occupations.In the 50-mode case (red), the self-polarization termdistorts the spBO surfaces further (Fig. 3), shifting theelectron-nuclear non-adiabatic region to be centerednear R = 0 a.u.. The -photon bands are wider, with ω max = 0 . a.u., and become resonant with the self-polarization modified molecular resonance already at R = − a.u. This leads to transitions from the initially0-photon excited spBO surface to the -photon groundspBO surfaces already at very short times. This is evi-dent in the almost immediate mixed character of the R -resolved BO populations. The flatter slope of the excitedspBO surface together with the increased population inthe lower spBO surface (Fig. 4.a), greatly slows the nu-clear density down compared to the fewer-mode cases,and results in a full suppression of both the proton andelectron transfer, as evident from panels (a-c) in Fig. 3.The wave packet reflects before appreciably reachingthe avoided crossing, and the change in the spatially-7 T o t a l P h o t o n N u m b e r Time[fs](a) 0 1 2 3 4 5 6 0 10 20 30 40 P s e u d o P h o t o n N u m b e r Time[fs](b)-5-4-3-2-1 0 1 0 10 20 30 40 B ili n e a r C o u p li n g T e r m Time[fs](c) 01234 0 10 20 30 40 S e l f - P o l a r i z a t i o n T e r m Time[fs](d)
FIG. 5. The different components of the photon number: (a)Total photon number of Eq. (6) (b) Pseudo-photon number(first term in Eq. (6)) (c) Bilinear coupling term (second termin Eq. (6)), (d) Self-polarization term (third term in Eq. (6)), us-ing the same color code as Fig.4. averaged BO-population in Fig.4.(b) is caused solely bythe nuclear wave packet dropping into the ground-BOstate by emitting photons, and is not due to the avoidedcrossing at R = 0 a.u. Considering the photon numbershown in Fig. 4, although the free photonic field com-ponent in panel (b) has a rapid initial increase and thengrows throughout the time evolution as in the few-modecases, while the linear term has a compensating decrease(panel c), the total photon number decreases after sometime, due to the self-polarization contribution. As ex-pected, this term increases with the number of modes atthe initial time, but since it is proportional to the the totaldipole of the system, whose transfer is suppressed, theresulting photon number tracks this behavior and is ul-timately reduced compared with the fewer-mode cases.The 70-mode case (blue) with ω max = 0 . a.u. canbe seen as an enhanced 50-mode case and leads to aneven stronger suppression of proton-coupled electrontransfer, with the same two key features that have beenresponsible for the cavity-modified dynamics in the 10-and higher-mode cases now having an even greaterimpact. First, by including the resonant frequency ofthe initial position of the nuclear wave packet, partof the wave packet almost immediately drops to thelower surface by emitting photons Fig. 5; see also R -averaged BO populations at early times, Fig. 4.(b) andthe mixed character developing in the R -resolved pop-ulations of Fig. 4.(a). Second, the deviation of the ex-cited spBO surface is now strong enough that its gra-dient slopes back to the left soon after the initial nu-clear wavepacket slides down from its initial positionat R = − a.u., sloping back to the left, in contrast tothe cavity-free excited BO surface. The overlap of theextensively broadened -photon-excited- and -ground- bands increases significantly creating a near-continuumof avoided crossings of polaritonic surfaces. The -photon surfaces are everywhere surrounded by near-lying n -photon surfaces. Compared to the 50-modecase, even less density reaches the region of the avoidedcrossing, which is now even wider. The slope of the ex-cited spBO-band results in even slower nuclear dynam-ics, with the nuclear and electronic dipole returning totheir initial positions after only a small excursion away,as evident in Fig. 3. Analogously to the discussion forthe 50-mode case, we find a larger initial total photonnumber, compared to the 50-mode case, followed by amoderate increase and decrease due to the early reflec-tion of the nuclear wave packet.Finally, to emphasize the importance of the self-polarization term on the dynamics, in Fig. 6 we comparethe results of the MTE dynamics on the electronic andnuclear dipoles and photon number when this term isneglected (dashed) or included (solid) for 10, 30, 50 and70 modes. For the electronic and nuclear dipole alreadyfor the 10 mode case deviations up to . a.u. (electronic)and . a.u. (nuclear) are found at later times. The errorin neglecting self-polarization becomes especially no-table for the 50- and 70-mode cases, where including theself-polarization yields a decrease of the proton-transfer(from 50 modes to 70 modes), while not including theself-polarization yields an increase of the proton trans- P h o t o n N u m b e r [ a . u . ] time[fs](c.1) (b.4)(b.3)(b.2)-4-2024 N u c l e i D i p o l e [ a . u . ] (b.1) (a.4)Without SP (dashed) With SP (solid)70(a.3)50(a.2)30-8-6-4-2024 E l e c t r o n i c D i p o l e [ a . u . ] (a.1)10 FIG. 6. Difference of the photon number (upper panel),nuclear dipole (middle panel) and electronic dipole (lowerpanel) without self-polarization term (dashed) and with self-polarization term (solid) in the same color code as Fig.4. . (70-modes), compared tothe simulations that include self-polarization. This canbe explained with the trends of the total dipole momentdiscussed above, since a dominant contribution to thephoton number is the self-polarization term in Eq. (6)which depends directly on this. IV. DISCUSSION AND OUTLOOK
Our results suggest that the effect of multiple cavity-modes on the reaction dynamics can lead to dramati-cally different dynamics than the cavity-free case. Thisis true even when the cavity-modes are far from theelectronic resonances encountered in the dynamics, andeven more so when cavity-modes are resonant with thematter system. In particular, for the model of cavity-induced suppression of proton-coupled electron trans-fer investigated here, we find an overall increase ofthe suppression the more photon modes are accountedfor. Two mechanisms are fundamentally responsible forthe difference: First, the self-polarization term grows insignificance with more modes with the effect that self-polarization-modified BO surfaces are distorted signif-icantly away from their cavity-free shape. Polaritonicsurfaces, eigenvalues of H − T n , should include the ex-plicit matter-photon coupling on top of these spBO sur-faces. Second, the n -photon-spBO bands become widerand increasingly overlapping, yielding a very mixedelectronic character with much exchange between sur-faces. These new dressed potential energy surfaces pro-vide a useful backdrop to analyze the dynamics, andwill form a useful tool in analyzing the different surfacesput forward to study coupled photon-matter systems,for example the polaritonic surfaces, and especially thetime-dependent potential energy surface arising fromthe exact factorization as this single surface alone pro-vides a complete picture of the dynamics.The MTE treatment of the photons appears to be apromising route towards treating realistic light-mattercorrelated systems. In particular, this method is able to capture quantum effects such as cavity-induced sup-pression of proton-coupled electron transfer, yet over-comes the exponential scaling problem with the num-ber of quantized cavity modes. However, a practical ap-proach for realistic systems will further need an approx-imate treatment of the matter part. From the electronicside TDDFT would be a natural choice, while a prac-tical treatment of nuclei calls for a classical treatmentsuch as Ehrenfest or surface-hopping in some basis. Themultiple-crossings inside the n -photon spBO bands sug-gest that simple surface-hopping treatments based onspBO surfaces should be used with much caution andthat decoherence-corrections should be applied, for ex-ample those generalized from the exact factorization ap-proach to the electron-nuclear problem [68, 69]. Fur-ther, the MTE approach could provide a way to accu-rately approximate the light-matter interaction part ofthe QEDFT xc functional [4, 27, 30, 45].Finally, we note that the present findings are generalin that the increasing importance of self-polarizationwith more photon modes is expected to hold for thedescription and control of cavity-driven physical pro-cesses of molecules, nanostructures and solids embed-ded in cavities in general. Extensions of these findingsto multi-mode cavitites suggest a new possibility of con-trolling and changing chemical reactions via the self-polarization without the need to explicitly change thelight-matter coupling strength itself. ACKNOWLEDGMENTS
We would like to thank Johannes Feist for insightfuldiscussions and the anonymous referees for very help-ful comments. Financial support from the US NationalScience Foundation CHE-1940333 (NM) and the Depart-ment of Energy, Office of Basic Energy Sciences, Divi-sion of Chemical Sciences, Geosciences and Biosciencesunder Award de-sc0020044 (LL) are gratefully acknowl-edged. NMH gratefully acknowledges an IMPRS fel-lowship. This work was also supported by the EuropeanResearch Council (ERC-2015-AdG694097), the Cluster ofExcellence (AIM), Grupos Consolidados (IT1249-19) andSFB925 ”Light induced dynamics and control of corre-lated quantum systems. The Flatiron Institute is a divi-sion of the Simons Foundation.
DATA AVAILABILITY STATEMENT
The data that support the findings of this study areavailable from the corresponding authors upon reason-able request.9
1] T. W. Ebbesen, Accounts of chemical research , 2403(2016).[2] J. George, T. Chervy, A. Shalabney, E. Devaux, H. Hiura,C. Genet, and T. W. Ebbesen, Phys. Rev. Lett. , 153601(2016).[3] H. Hiura, A. Shalabney, and J. George, ChemRxiv (2018).[4] M. Ruggenthaler, N. Tancogne-Dejean, J. Flick, H. Appel,and A. Rubio, Nature Reviews Chemistry , 0118 (2018).[5] C. Riek, D. V. Seletskiy, A. S. Moskalenko, J. F. Schmidt,P. Krauspe, S. Eckart, S. Eggert, G. Burkard, and A. Leit-enstorfer, Science , 420 (2015).[6] A. S. Moskalenko, C. Riek, D. V. Seletskiy, G. Burkard,and A. Leitenstorfer, Phys. Rev. Lett. , 263601 (2015).[7] T. Byrnes, N. Y. Kim, and Y. Yamamoto, Nature Physics , 803 (2014).[8] J. Kasprzak, M. Richard, S. Kundermann, A. Baas,P. Jeambrun, J. Keeling, F. Marchetti, M. Szyma ´nska,R. Andre, J. Staehli, et al. , Nature , 409 (2006).[9] S. Schmidt, Physica Scripta , 073006 (2016).[10] A. Thomas, E. Devaux, K. Nagarajan, T. Chervy,M. Seidel, D. Hagenm ¨uller, S. Sch ¨utz, J. Schachen-mayer, C. Genet, G. Pupillo, et al. , arXiv preprintarXiv:1911.01459 (2019).[11] A. Thomas, L. Lethuillier-Karl, K. Nagarajan, R. M. Ver-gauwe, J. George, T. Chervy, A. Shalabney, E. Devaux,C. Genet, J. Moran, et al. , Science , 615 (2019).[12] J. Feist and F. J. Garcia-Vidal, Phys. Rev. Lett. , 196402(2015).[13] J. Galego, F. J. Garcia-Vidal, and J. Feist, Phys. Rev. X ,041022 (2015).[14] J. Flick, M. Ruggenthaler, H. Appel, and A. Rubio, Pro-ceedings of the National Academy of Sciences , 15285(2015).[15] J. Galego, F. J. Garcia-Vidal, and J. Feist, Nature Commu-nications , 13841 EP (2016).[16] J. Schachenmayer, C. Genes, E. Tignone, and G. Pupillo,Phys. Rev. Lett. , 196403 (2015).[17] M. Cirio, S. De Liberato, N. Lambert, and F. Nori, Phys.Rev. Lett. , 113601 (2016).[18] J. Flick, H. Appel, M. Ruggenthaler, and A. Rubio, Jour-nal of Chemical Theory and Computation , 1616 (2017),pMID: 28277664.[19] C. Sch¨afer, M. Ruggenthaler, H. Appel, and A. Rubio,Proceedings of the National Academy of Sciences ,4883 (2019).[20] R. F. Ribeiro, L. A. Mart´ınez-Mart´ınez, M. Du, J. Campos-Gonzalez-Angulo, and J. Yuen-Zhou, Chem. Sci. , 6325(2018).[21] M. Kowalewski, K. Bennett, and S. Mukamel, The Journalof Physical Chemistry Letters , 2050 (2016).[22] J. F. Triana, D. Pel ˜A¡ez, and J. L. Sanz-Vicario, The Jour-nal of Physical Chemistry A , 2266 (2018), pMID:29338227.[23] T. Szidarovszky, G. J. Hal´asz, A. G. Cs´asz´ar, L. S. Ceder-baum, and A. Vib´ok, The Journal of Physical ChemistryLetters , 6215 (2018).[24] G. Groenhof and J. J. Toppari, The Journal of PhysicalChemistry Letters , 4848 (2018), pMID: 30085671.[25] H. L. Luk, J. Feist, J. J. Toppari, and G. Groenhof, Jour-nal of Chemical Theory and Computation , 4324 (2017),pMID: 28749690. [26] O. Vendrell, Chemical Physics , 55 (2018).[27] J. Flick, C. Sch¨afer, M. Ruggenthaler, H. Appel, andA. Rubio, ACS Photonics , 992 (2018).[28] J. F. Triana, D. Pel´aez, and J. L. Sanz-Vicario, The Journalof Physical Chemistry A , 2266 (2018).[29] J. Flick, Exact nonadiabatic many-body dynamics: Electron-phonon coupling in photoelectron spectroscopy and light-matter interactions in quantum electrodynamical density-functional theory , Ph.D. thesis, Humboldt-Universit¨at zuBerlin Berlin (2016).[30] M. Ruggenthaler, J. Flick, C. Pellegrini, H. Appel, I. V.Tokatly, and A. Rubio, Phys. Rev. A , 012508 (2014).[31] C. Sch¨afer, M. Ruggenthaler, and A. Rubio, Physical Re-view A , 043801 (2018).[32] M. S´anchez-Barquilla, R. Silva, and J. Feist, (2019),arXiv:1911.07037.[33] N. M. Hoffmann, C. Sch¨afer, A. Rubio, A. Kelly, andH. Appel, Phys. Rev. A , 063819 (2019).[34] N. M. Hoffmann, C. Sch¨afer, N. S¨akkinen, A. Rubio,H. Appel, and A. Kelly, The Journal of Chemical Physics , 244113 (2019).[35] J. Flick, M. Ruggenthaler, H. Appel, and A. Rubio, Pro-ceedings of the National Academy of Sciences , 3026(2017).[36] J. del Pino, F. A. Y. N. Schr¨oder, A. W. Chin, J. Feist, andF. J. Garcia-Vidal, Phys. Rev. Lett. , 227401 (2018).[37] S. Franke, S. Hughes, M. K. Dezfouli, P. T. Kristensen,K. Busch, A. Knorr, and M. Richter, Phys. Rev. Lett. ,213901 (2019).[38] J. Flick, D. M. Welakuh, M. Ruggenthaler, H. Appel, andA. Rubio, ACS photonics , 2757 (2019).[39] K. B. Arnardottir, A. J. Moilanen, A. Strashko, P. T¨orm¨a,and J. Keeling, arXiv preprint arXiv:2004.06679 (2020).[40] B. Huttner and S. M. Barnett, Phys. Rev. A , 4306 (1992).[41] S. Y. Buhmann, Dispersion forces I: Macroscopic quantumelectrodynamics and ground-state Casimir, Casimir–Polderand van der Waals Forces , Vol. 247 (Springer, 2013).[42] S. Buhmann,
Dispersion Forces II: Many-Body Effects, Ex-cited Atoms, Finite Temperature and Quantum Friction , Vol.248 (Springer, 2013).[43] S. Y. Buhmann and D.-G. Welsch, Physical Review A ,012110 (2008).[44] S. Scheel and S. Buhmann, Acta Physica Slovaca. Reviewsand Tutorials , 675 (2008).[45] C. Pellegrini, J. Flick, I. V. Tokatly, H. Appel, and A. Ru-bio, Phys. Rev. Lett. , 093001 (2015).[46] This includes Quantum-Electrodynamical Density Func-tional Theory (QEDFT) [4, 27, 30, 45], which is an exactnon-relativistic generalization of time-dependent densityfunctional theory that dresses electronic states with pho-tons and allows to retain the electronic properties of realmaterials in a computationally efficient way.[47] L. Lacombe, N. M. Hoffmann, and N. T. Maitra, Phys.Rev. Lett. , 083201 (2019).[48] N. M. Hoffmann, H. Appel, A. Rubio, and N. T. Maitra,The European Physical Journal B , 180 (2018).[49] C. Sch¨afer, M. Ruggenthaler, V. Rokaj, and A. Rubio, ACSphotonics , 975 (2020).[50] V. Rokaj, D. M. Welakuh, M. Ruggenthaler, and A. Ru-bio, Journal of Physics B: Atomic, Molecular and OpticalPhysics , 034005 (2018).
51] E. A. Power and S. Zienau, Philosophical Transactions ofthe Royal Society of London. Series A, Mathematical andPhysical Sciences , 427 (1959).[52] I. V. Tokatly, Phys. Rev. Lett. , 233001 (2013).[53] A. Mandal, S. M. Vega, and P. Huo, arXiv preprintarXiv:2005.00201 (2020).[54] S. Shin and H. Metiu, The Journal of Chemical Physics , 9285 (1995).[55] J.-Y. Fang and S. Hammes-Schiffer, The Journal of Chem-ical Physics , 8442 (1997).[56] J.-Y. Fang and S. Hammes-Schiffer, The Journal of Chem-ical Physics , 5727 (1997).[57] A. Abedi, F. Agostini, Y. Suzuki, and E. K. U. Gross, Phys.Rev. Lett. , 263001 (2013).[58] F. Agostini, A. Abedi, Y. Suzuki, S. K. Min, N. T. Maitra,and E. K. U. Gross, J. Chem. Phys. , 084303 (2015).[59] J. Galego, C. Climent, F. J. Garcia-Vidal, and J. Feist, Phys.Rev. X , 021057 (2019).[60] H. Sambe, Phys. Rev. A , 2203 (1973).[61] T. Fiedlschuster, J. Handt, and R. Schmidt, Phys. Rev. A , 053409 (2016). [62] M. Thachuk, M. Y. Ivanov, and D. M. Wardlaw, The Jour-nal of Chemical Physics , 4094 (1996).[63] Y. Sato, H. Kono, S. Koseki, and Y. Fujimura, Journal ofthe American Chemical Society , 8019 (2003).[64] A. Abedi, N. T. Maitra, and E. K. Gross, Physical ReviewLetters , 123002 (2010).[65] A. Abedi, N. T. Maitra, and E. Gross, The Journal ofChemical Physics , 22A530 (2012).[66] D. De Bernardis, T. Jaako, and P. Rabl, Physical Review A , 043820 (2018).[67] We note that the two-mode case can also be solved ex-actly numerically, but the single-mode comparison herealready illustrates the main points.[68] J.-K. Ha, I. S. Lee, and S. K. Min, The Journal of Phys-ical Chemistry Letters , 1097 (2018), pMID: 29439572,https://doi.org/10.1021/acs.jpclett.8b00060.[69] S. K. Min, F. Agostini, and E. K. U. Gross, Phys. Rev. Lett. , 073001 (2015)., 073001 (2015).