The exact time-dependent potential energy surface driving the nuclear dynamics was recently shown to be a useful tool to understand and interpret the coupling of nuclei, electrons, and photons in cavity settings. Here, we provide a detailed analysis of its structure for exactly solvable systems that model two phenomena: cavity-induced suppression of proton-coupled electron-transfer and its dependence on the initial state, and cavity-induced electronic excitation. We demonstrate the inadequacy of simply using a weighted average of polaritonic surfaces to determine the dynamics. Such a weighted average misses a crucial term that redistributes energy between the nuclear and the polaritonic systems, and this term can in fact become a predominant term in determining the nuclear dynamics when several polaritonic surfaces are involved. Evolving an ensemble of classical trajectories on the exact potential energy surface reproduces the nuclear wavepacket quite accurately, while evolving on the weighted polaritonic surface fails after a short period of time. The implications and prospects for application of mixed quantum-classical methods based on this surface are discussed.

## I. INTRODUCTION

The interesting fact that simply confining a quantum system increases its coupling to vacuum fluctuations of the radiation field so much as to change the system’s properties has been known since the early days of quantum mechanics. A key realization was that spontaneous emission rates of a two-level system can be increased by placing it in a resonant cavity due to the enhanced light–matter coupling strength.^{1} The recent burst of activity in cavity quantum electrodynamics is largely due to impressive experimental advances in manipulating matter at the meso- and nano-scales, which have driven the investigations toward exploring the possibilities of new cavity-mediated phenomena. The playing field is rich when considering enhanced light–matter interaction on top of the already complex interactions of real matter, as opposed to the simple two-level-like systems that were considered earlier.^{2–10}

Polaritonic phenomena involve a complex interplay of photonic, electronic, and nuclear degrees of freedom, but often we are particularly interested in the dynamics of one particular subsystem. In polaritonic chemistry applications, it is typically the nuclear (or ionic) system we are interested in, for example, in investigating changes in chemical reactivity when molecules are confined. The exact factorization (EF) approach^{11–16} offers a natural context to study this problem since it provides a rigorous and unique definition of potentials acting on a subsystem, which exactly incorporate the effects of coupling to the other subsystems. In particular, one obtains a time-dependent Schrödinger equation (TDSE) for the nuclear wavefunction $\chi (R\u0332\u0332,t)$, in which the potentials contain the complete effect of coupling to the photonic and electronic systems. The full, exact wavefunction of the system has the form of a single correlated product, $\chi (R\u0332\u0332,t)\Phi R\u0332\u0332(r\u0332\u0332,q\u0332\u0332,t)$, where $r\u0332\u0332,q\u0332\u0332$ represent the electronic and photonic coordinates, respectively. This is quite in contrast to another exact representation for the full wavefunction, the Born–Huang expansion generalized to include photon degrees of freedom, which contains an infinite number of terms of correlated products, $\u2211j\chi j(R\u0332\u0332,t)\Phi R\u0332\u0332j(r\u0332\u0332,q)$; here, $\Phi R\u0332\u0332j(r\u0332\u0332,q)$ are the polaritonic eigenstates. Instead of the single TDSE for the nuclear system that appears in EF, one has a set of coupled equations to solve, and it is not possible to identify a potential that drives the nuclear motion. The concept of a single equation for the nuclear wavefunction, as arises in EF, is especially useful when turning to mixed quantum-classical methods where being able to unambiguously define a force acting on a classical nuclear trajectory at each point is desirable.

The EF was originally derived with the static coupled electron-nuclear problem in mind,^{11,12,17} while the time-dependent problem was introduced later.^{13,14} Extensions of the EF idea to a variety of different settings have been made, e.g., the static purely electronic problem^{18} with connections to density functional theory,^{19} the time-dependent electronic problem with connections to the single-active electron approach,^{20} and recently, the electronic Fock space to derive a new quantum embedding method.^{21–23} A reverse factorization for electron-nuclear problems has also been explored where the nuclear wavefunction is conditionally dependent on the electronic coordinate, which has an advantage when one is most interested in the electronic system.^{24}

References 15, 16, and 25 extended the EF approach to polaritonic systems, focusing on phenomena that are either primarily interesting for the photonic system,^{15} the electronic system,^{25} or the nuclear system.^{16} The present work expands on the latter. Reference 16 studied the exact potential driving a model of cavity-induced suppression of proton coupled electron-transfer. The model molecule was photo-excited into the first excited Born–Oppenheimer state, and when outside the cavity, the nuclear wavepacket evolves toward a narrow avoided crossing where there is some transfer of population to the ground-electronic state. Since, in this state, the electronic dipole is localized on the opposite side to where it was localized initially, the process is known as proton-coupled electron-transfer. The enhanced spontaneous emission rate from the excited molecule when placed in an optical cavity led to a significant partial photon emission when the cavity resonance was tuned to the initial excitation frequency. The interplay between the partial photon emission, the partial electronic de-excitation, and coupling to nuclear motion led to less of the nuclear wavepacket evolving to the narrow avoided crossing and, hence, a suppression of the proton-coupled electron-transfer. While the shape of the polaritonic surfaces suggests that such a phenomenon could occur, they alone could not predict the suppression, as a weaker light–matter coupling strength led to practically identical surfaces but with no suppression dynamics. What would be needed to determine the dynamics would be coupling matrix elements between the surfaces. On the other hand, the exact potential arising from the EF approach could distinguish clearly the suppression phenomenon.

In this paper, we present two further examples of the EF approach applied to polaritonic chemistry. In the first, we revisit the same model of Ref. 16 but begin now in an excited polaritonic state instead of an electronically excited state. The dependence of cavity-altered phenomena on the initial state of polaritonic system has not been widely explored. In particular, when there is an initial excitation in the system, a question is whether this excitation is done in the presence of the cavity such that the initial excitation is a polaritonic one, or whether it is done effectively outside the cavity and then the molecule is inserted into the cavity, in which case the initial excitation is a purely electronic one. For our model of proton-coupled electron-transfer, we find suppression is observed when beginning in a polaritonic state, although not as much as for the initial purely electronically excited state, and we compare the exact potential from EF in the two cases and with the cavity-free case. The second example we consider is one of cavity-induced electronic excitation, where a photo-excited molecule is placed in a cavity whose fundamental frequency is resonant with a molecular frequency at two different geometries: At one nuclear configuration, the mode frequency coincides with the ground and first electronic states, while at the other nuclear configuration, it coincides with the ground and second. Thus, as an excited nuclear wavepacket on the first electronic state evolves from the first nuclear configuration into the second, an electronic excitation into the second excited Born–Oppenheimer state occurs, which is absent in the cavity-free case. With an eye to developing mixed quantum-classical approximations, we plot the force on the proton at different time snapshots due to different components of this exact surface, especially since somewhat localized large step features appear in one of the components, which may exaggerate its effect on the dynamics. A quasiclassical Wigner propagation of nuclear trajectories shows that dynamics on the exact potential reproduces the exact quantum dynamics accurately. We begin in Sec. II with a brief review of the EF approach.

## II. EXACT FACTORIZATION APPROACH

The essential statement of the exact factorization approach is very general: The complete wavefunction for a system of coupled quantum subsystems can be factorized into a marginal wavefunction depending on one (or some) of the subsystem coordinates and a conditional wavefunction of the other coordinates, parametrically dependent on the marginal degrees of freedom. The equations for the two parts depend on the choice of the marginal and conditional and how the Hamiltonian depends on them. For example, consider a system of interacting electrons only. In Ref. 20, the factorization is done in real-space with a chosen number of electronic coordinates associated with the marginal, and the resulting equation for the marginal has a TDSE form. This is not the case in Ref. 22 where instead the factorization is done in Fock space with a chosen set of orbitals included in the marginal; instead, an effective embedding Hamiltonian matrix is deduced, which is of interest to problems in strong correlation.

Consider now any situation for which the components of the Hamiltonian associated with the chosen marginal degrees of freedom have the form of a kinetic term, $T^=\u2212\u2211I\u2207I2/(2MI)$, plus potential terms that are multiplicative operators in the marginal coordinate. Then, the EF equations take the same form as in the original electron-nuclear case.^{13,14} In particular, the equation for the marginal factor is a Schrödinger equation, containing scalar and vector potentials that incorporate exactly the effect of couplings to the other subsystems. Indeed, this is the case for a molecule in a cavity, when we take the non-relativistic photon–matter Hamiltonian in the dipole approximation in the length gauge.^{5,26–28} The full Hamiltonian has the form

where $Tn^=\u2212\u2211I\u2207I2/(2MI)$ is the nuclear kinetic energy, $\u0124BO$ is the Born–Oppenheimer Hamiltonian of the molecule, and the other terms arise from the photonic degrees of freedom. Note that we use atomic units (*e*^{2} = *ℏ* = *m*_{e} = 1) throughout this paper. $\u0124p$ is the “free” photonic Hamiltonian, which has the harmonic form

Here, the sum goes over all the modes *α* in the cavity, but in our examples, we consider only a single mode of frequency *ω*_{α} = *ω*_{c}, resonant with an electronic excitation frequency at a particular molecular geometry; $q^\alpha $ is the displacement field coordinate, and $p^\alpha =\u2212i\u2202\alpha $. Recent work has shown that the dynamics can change significantly as more cavity modes are accounted for,^{29} but here we assume that the harmonics of this fundamental mode are high enough that they do not couple strongly with the molecule. The matter–photon coupling is represented by the term

where *N*_{p} is the number of photon modes, *N*_{n,e} is the number of nuclei and electrons respectively, and $\lambda \alpha \u2192$ is the matter–photon coupling vector. Finally, the self-polarization term is

This has a negligible effect on our examples due to the small value of *λ*_{α} that we use, so we will later safely neglect it. We note that the importance of this term grows as more modes are considered,^{29} and its necessity in obtaining a consistent ground-state and maintaining gauge-invariance has been pointed out.^{30,31} It should be borne in mind that the Hamiltonian above completely neglects any losses at the cavity boundaries, so the phenomena that we will observe will in reality be dampened somewhat.

In the following, we let $R\u0332\u0332$ represent all nuclear coordinates (three degrees of freedom for each nucleus), $r\u0332\u0332$ represent all electronic coordinates, and $q\u0332\u0332$ represent all the displacement field coordinates for the photons (one degree of freedom for each mode). Taking the marginal as the nuclear degrees of freedom,

yields^{13–16,32,33}

and

with the following coupling terms:

the coupling operator between the electron–photon system and the nuclei,

- (ii)
the time-dependent potential energy surface (TDPES),

which has the three components: the (self-polarization-included) weighted-polaritonic component,

where

the kinetic-like term

and the gauge-dependent term

and

- (iii)
the time-dependent vector potential

Since the marginal factor $\chi (R\u0332\u0332,t)$ in Eq. (5) can be multiplied by a phase $ei\theta (R\u0332\u0332,t)$, while the conditional factor $\Phi R\u0332\u0332(r\u0332\u0332,q\u0332\u0332,t)$ is multiplied by the inverse phase $e\u2212i\theta (R\u0332\u0332,t)$, without changing the product $\chi \Phi R\u0332\u0332$, there is a gauge-like dependence to these equations.^{13,14} One can choose different gauges in which $AI(R\u0332\u0332,t)$ is changed by the addition of $\u2207I\theta (R\u0332\u0332,t)$, while $\u03f5GD(R\u0332\u0332,t)$ is changed by the addition of $\u2202t\theta (R\u0332\u0332,t)$. The relation of the coupling terms that appear in the EF to the couplings that appear in a Born–Huang expansion may be found in Ref. 13 for the electron-nuclear case.

We now note three simplifications to these general equations that are relevant for the case studies we will consider. First, our models here have only one nuclear degree of freedom, which means that we can always choose a gauge in which the vector potential $A(R\u0332\u0332,t)=A(R,t)=0$. With this gauge choice, then the exact nuclear wavefunction satisfies a Schrödinger equation (7) with a purely scalar potential of Eq. (9), composed of the three components, Eqs. (10)–(13). Second, the kinetic component *ϵ*_{kin}(*R*, *t*) is negligible compared to the other two terms due to the factor of 1/*M* where *M* is the nuclear mass. Our models consider a proton of mass 1836 times the mass of an electron. This is in contrast to the cases where the electronic^{24,34} or photonic^{15} systems are chosen to be the marginal, where it certainly cannot be neglected. Third, we note that the self-polarization contribution to *ϵ*_{wpol}(*R*, *t*) is negligible for the single-mode cases and coupling strengths we are considering here; this term becomes important on the other hand when many cavity modes are considered.^{29} Thus, the nuclear dynamics is driven by the Schrödinger equation (7), taking *A*(*R*, *t*) = 0 and *ϵ*(*R*, *t*) = *ϵ*_{wpol}(*R*, *t*) + *ϵ*_{GD}(*R*, *t*).

### A. Interpreting the TDPES

The TDPES is composed then of two potentials *ϵ*_{wpol}(*R*, *t*) and *ϵ*_{GD}(*R*, *t*) that completely govern the nuclear motion. Their dependence on the conditional electronic-photonic wavefunction $\Phi R(r\u0332\u0332,q\u0332\u0332,t)$ embodies the correlation between the nuclear and electronic-photonic systems. It is instructive to expand $\Phi R(r\u0332\u0332,q\u0332\u0332,t)$ in terms of the polaritonic eigenstates,

where

defines the polaritonic surfaces *ϵ*_{pol,k}(**R**) and polaritonic eigenstates $\Phi Rk(r\u0332\u0332,q\u0332\u0332)$.^{35,36} Then, due to orthonormality of the polaritonic eigenstates, one can decompose *ϵ*_{wpol} in terms of the polaritonic projections,

The physical meaning of *ϵ*_{wpol} is then clear as a weighted average of the polaritonic surfaces with the weights given by the projections of the conditional electronic-photonic wavefunction on the polaritonic surfaces. With mixed quantum-classical methods in mind, this is somewhat reminiscent of the form of the effective potential the classical nuclei experience in Ehrenfest dynamics, however, with the crucial difference in that the coefficients are here *spatially dependent*. The terms in the electronic equation yield gradients of these terms (with respect to the nuclear coordinate), which allow different trajectories making up the nuclear wavepacket in a multi-trajectory approach to go in different directions. This results in the possibility of wavepacket splitting when there is partial occupation of different electronic surfaces, which does not happen in a multi-trajectory Ehrenfest approach.^{37–40} As in Ref. 38, we will see in Secs. III A and III B that although wavepacket splitting occurs to some extent when *ϵ*_{wpol} is used to propagate, *ϵ*_{GD} greatly increases this effect and is needed to get accurate dynamics.

Similarly, expanding the gauge-dependent term in terms of polaritonic components,

The meaning and role of *ϵ*_{GD} are less directly clear from its definition. Instead, some insight is gained from comparing the total energy of the system $\u27e8\u0124\u27e9$ with the “marginal energy” defined via $\u27e8\u0124nuc\u27e9$ where *H*_{nuc} is the Hamiltonian in the nuclear Schrödinger equation, governing the motion of the nuclear wavepacket. We have

where $\u27e8T^nmarg\u27e9=\u2212\u222b\chi *(R,t)\u220722M\chi (R,t)dR$ is the kinetic energy of the marginal system, which differs from the nuclear kinetic energy $\u27e8T^n\u27e9=\u2212\u222bdr\u0332\u0332dq\u0332\u0332\Psi *(r\u0332\u0332,R,q\u0332\u0332,t)\u2207R22M\Psi (r\u0332\u0332,R,q\u0332\u0332,t)$ by the last term in Eq. (19).^{14} The expectation value of $\xdbep\u2212n$, however, is very small compared with the other terms when the marginal is chosen as the nuclear system due to the 1/*M* factor in Eq. (8)^{37} [note that only the first term of Eq. (8) contributes to its expectation value]. This means that the total energy is essentially

On the other hand, the “marginal energy,” $\u27e8\u0124nuc\u27e9$, is defined as that associated with the Hamiltonian driving the marginal motion,

That is, *ϵ*_{GD}(*R*, *t*) does not contribute to the total energy but instead redistributes energy between the electron–photon system and the nuclear system.

To get a qualitative sense of its structure, we allow ourselves to extend the ideas above to the energy density rather than the total integrated energy. We distinguish two cases. First, when the polaritonic character is distinct in a piecewise localized way, i.e., one polaritonic character dominates the wavefunction in one region of space, while another dominates it in a neighboring region. For definiteness, let us say the polaritonic energy is locally larger on the left than on the right. Then, we expect $\u27e8T^nmarg\u27e9$ would be smaller on the left than on the right, and this is enabled by a step up in *ϵ*_{GD} from left to right in the transition region. This has been observed in the electron-nuclear case,^{38,41} and we will observe it again here in our polaritonic examples in Sec. III. In the second case, where there is a mixed polaritonic character, the expected structure of *ϵ*_{GD}(*R*, *t*) is generally less clear. Here, different parts of the nuclear wavepacket are associated with different polaritonic surfaces in the same region of space. In the cases where the mean-field nature of *ϵ*_{wpol}(*R*, *t*) “washes out” distinct structures, we have found *ϵ*_{GD}(*R*, *t*) can be of paramount importance. For example, when part of the nuclear wavepacket associated with a polaritonic surface begins to reflect from this surface, while the other parts do not, a dynamical step appears in *ϵ*_{GD}(*R*, *t*) to enable the reflection of the part. Because of the average nature of *ϵ*_{wpol}(*R*, *t*), the different curvatures of different polaritonic surfaces “under” the nuclear wavepacket are not always evident.

The equations of the exact factorization approach are no easier to solve than the original full Schrödinger equation; on the contrary, practical numerical schemes to solve them are currently elusive.^{42,43} Rather, the practical power of the approach is due to providing rigorous definitions of coupling potentials, which then offer a clear and well-founded starting point for approximation schemes, as in Refs. 37 and 44–49. Studies of the exact potentials in cases where they could be found^{38,41} were very instructive for these developments. Like in those works, we here obtain the exact potentials for the polaritonic case by inversion. That is, we first compute the full matter–photon wavefunction, which is possible for the systems of one nuclear, one electronic, and one photonic degree of freedom that we study. We then obtain the nuclear wavefunction *χ*(*R*, *t*) = |*χ*(*R*, *t*)|*e*^{iS(R,t)}, which corresponds to the square root of the nuclear density times a phase. The nuclear density is obtained by integration, |*χ*(*R*, *t*)|^{2} = ∫ *drdq*|Ψ(*r*, *R*, *q*, *t*)|^{2}. The phase is fixed by the gauge choice, which for *A*(*R*, *t*) = 0 gives $S(R,t)=\u222bRdR\u2032\u2009I(\u27e8\Psi |\u2202R|\Psi \u27e9)/|\chi (R,t)|2$.^{38} Once *χ*(*R*, *t*) is obtained, the rest is straightforward as Φ_{R}(*r*, *q*, *t*) = Ψ(*r*, *R*, *q*, *t*)/*χ*(*R*, *t*) and we simply use Eqs. (9)–(13) with a simple two-point stencil finite difference in time to compute the gauge-dependent part of the potential.

## III. MODEL HAMILTONIAN

The Hamiltonian for the matter part of our model systems has the Shin–Metiu^{50–52} form

where $T^n=\u221212M\u22022\u2202R2$, $T^e=\u221212\u22022\u2202r2$, and

which describes one electron and one proton moving along the axis between two fixed heavy ions separated by a distance *L*.

We consider two sets of parameters, which will be used to illustrate two effects. The first effect (Sec. III A) is cavity-induced suppression of proton-coupled electron-transfer (PCET), and the parameters of the model are *L* = 19.0 a.u., *a*_{+} = 3.1 a.u., *a*_{−} = 4.0 a.u., *a*_{f} = 5.0 a.u., and *M* = 1836 a.u., the proton mass, which was the set of parameters chosen in Ref. 16. The cavity frequency is chosen as *ω*_{c} = 0.1 a.u., and the photon–matter coupling strength is chosen as *λ* = 0.005 a.u. The second effect (Sec. III B) is cavity-induced electronic excitation (ELEX), and the parameters chosen here are *L* = 19.0 a.u., *a*_{+} = 4.0 a.u., *a*_{−} = 4.0 a.u., *a*_{f} = 5.0 a.u., and *M* = 1836 a.u. with the cavity frequency *ω*_{c} = 0.049 a.u. and again *λ* = 0.005 a.u. These parameters are in a reasonable ballpark for plasmonic nano-gap cavities. The first few BO and polaritonic surfaces are shown in Fig. 1. In all our examples, the initial nuclear wavepacket is a Gaussian centered at R = −4 a.u.

### A. Cavity-induced suppression of PCET

Reference 16 showed that when the Shin–Metiu molecule is placed in a cavity resonant with the initial BO energy difference, the PCET process is suppressed. Part of the nuclear wavepacket never makes it to the electron-nuclear avoided crossing region as it is associated with the emission of a photon. This behavior of the nuclear wavepacket was shown to be directly correlated with the curvature of the exact TDPES as it evolves in time, and impossible to predict from the shape of the static polaritonic surfaces alone. In Ref. 16, the initial state was a Gaussian nuclear wavepacket placed on the second BO surface with no photon in the cavity, which would be the case if the excited molecule was prepared outside the cavity and then placed inside it. Instead, if the molecule is placed inside the cavity before the excitation is done, beginning in a polaritonic state would be more appropriate than the purely electronically excited one. This is the case we consider in the present paper. That is, the initial state has the form $\Psi (r,R,q,0)=Ne\u22122.85(R+4)2\Phi Rk=1(r,q)$, where $\Phi Rk=1(r,q)$ is the first excited polaritonic state [see Eq. (16)].

Figure 2 shows the electronic and nuclear dipoles (i.e., the expectation values of the electronic and nuclear coordinates, respectively) when beginning in the first polaritonic state, compared with beginning in the 0-photon, first-excited electronic state, which we will henceforth call the factorized state, and the cavity-free dipoles. We observe that PCET is suppressed also in the case of beginning in the polaritonic state, although not as much as in the factorized state.

Snapshots of the TDPES are plotted in Fig. 3 along with the nuclear density, in comparison with the cavity-free case and the factorized initial state case. Initially, the exact TDPES coincides with the polaritonic surface in the case of the polaritonic initial state and along the BO excited surface in the case of the factorized state. Here, we note that a global shift to the TDPES has no effect on the dynamics. As we do not set a zero of the potential at any particular *R*, when we say the TDPES coincides with a polaritonic surface, we really mean it is parallel to it. For *R* near −4 a.u., the polaritonic surface has the character of the 1-photon ground electronic BO state to the left and the character of the 0-photon electronically excited BO state to the right. The gentler slope of the polaritonic surface on the left at early times slows the nuclear wavepacket relative to both the dynamics of the factorized state and the cavity-free dynamics. The soft kink in the exact TDPES in the factorized initial-state case, which separates a region where the surface slopes down to the left on the left and down to the right on the right (*R* ≈ −2 a.u. at t = 17.42 fs), develops at somewhat later times in the case of the polaritonic initial state (visible in this case at 21.29 fs). This results in the density splitting into two parts later than for the initially factorized case, and ultimately, less of the wavepacket is trapped in the well on the left, compared to the factorized case. As in Ref. 16, we point out that this behavior cannot be predicted from the polaritonic surfaces alone; the polaritonic surfaces are independent of the choice of the initial state and alone could not predict the distinct behavior resulting from different initial states. On the other hand, the distinct shapes of the TDPES for these two states do correlate with the different nuclear wavepacket dynamics, and their differences with the cavity-free TDPES also correlate with the different dynamics seen for the cavity-free nuclear wavepacket. Similar to the factorized case, the part of the nuclear wavepacket to the left of the soft kink is unable to overcome the energy barrier in the TDPES to reach the region of electron-nuclear avoided crossing (R ≈ 2 a.u.) that leads to the electron-transfer, thus resulting in suppression of PCET. Since less of the nuclear wavepacket is trapped on the left in the case of the initially polaritonic state because of the shape of the TDPES at the earlier times (gentler slope and later development of the soft kink), the suppression of PCET is less than in the factorized initial state case.

Underlying the shapes of the TDPES are the absorption and emission of a photon and electronic transitions between BO states, which, through the conditional wavefunction dependence of the TDPES, control its structure. It can be instructive to consider these conditional quantities when interpreting the dynamics. Movies of the TDPES and its components along with the density are given in the supplementary material, as well as the *R*-resolved BO electronic populations |*C*_{i}(*R*, *t*)|^{2} defined through

Movies of the *n*-photon resolved nuclear densities, $|\chi n\u2212phi|2=|\u27e8\xi n|\Psi (t)\u27e9r,q|2$ where *ξ*_{n}(*q*) are the harmonic oscillator eigenstates of the photonic Hamiltonian, are also provided. For the factorized initial state at early times, partial photon emission occurs throughout the wavepacket, while for the polaritonic initial state, only the right part of the wavepacket becomes correlated with photon emission and electronic de-excitation, while the left part becomes correlated with photon absorption and electronic excitation. As the wavepacket moves toward the right, leaving the resonant region, there has been less time for the part of wavepacket originally on the left to re-emit the photon it absorbed, so more remains on the excited electronic surface compared to the factorized case and less ends up being trapped.

We have found that the essential features of the surface are similar to what was found for the factorized case, but details and timings are different. The gradients of the two dominant contributions to the TDPES, *ϵ*_{wpol} [Eq. (10)] and *ϵ*_{GD} [Eq. (13)], are shown in Fig. 4. At the first time shown, 3.39 fs, the force from the TDPES coincides with the slope of the second polaritonic surface as expected, with only the *ϵ*_{wpol} component contributing, but after very short times, both components begin to come into play. Where different parts of the nuclear wavepacket show a single polaritonic character that differs in different spatial regions, the gradient of *ϵ*_{wpol} follows the gradient of the relevant polaritonic surfaces in a piecewise manner with a peak in the region where they change character, e.g., at times 17.42 fs, 21.29 fs, and 25.56 fs. The gradient of *ϵ*_{wpol} is not simply a weighted average of the gradients of the polaritonic surfaces because the coefficients $Ci(R,t)$ in Eq. (15) vary with *R*, and in particular, these can vary significantly when the nuclear wavepacket correlated with one surface in one region switches to being correlated with another surface. The step-like transition in the coefficient leads to a peak-like structure in its gradient that affects the force the nuclei feel in this region. This structure tends to be partially compensated by a countering step in *ϵ*_{GD}, whose gradient is zero in regions where the nuclear wavepacket is dominated by one surface; the total force from both terms becomes more a step-like feature in the transition region. Typically though, both components need to be considered through most of the nuclear wavepacket. For example, at 21.29 fs, the transition region is quite large, and the *ϵ*_{GD} correction to *ϵ*_{wpol} results in quite a different force throughout the nuclear density except for the right-most part, even reversing the sign of the force. As the nuclear wavepacket separates into two distinct parts evolving on different surfaces, the opposing peaks in *ϵ*_{GD} and *ϵ*_{wpol}, albeit large, have little effect when located in regions of low density in between the parts of the wavepacket (e.g., between *R* ≈ 0 a.u. and 2 a.u. at time 35.80 fs). At this time, we see the gradient of the total surface coincides with the fourth polaritonic surface (parallel to the lowest polaritonic surface with character 1-photon and ground electronic states) in the extreme left part of the wavepacket and is coming entirely from *ϵ*_{wpol} there. On the other hand, the right-part of the wavepacket has a more mixed characters, and really only on its extreme right, it is dominated by a force from one polaritonic surface (the lowest surface). At later times, *ϵ*_{GD} plays a more complex role. For example, consider a situation where there is mixed character of the nuclear wavepacket such that part is correlated with one polaritonic surface, while another part in the same region is correlated with a different polaritonic surface. These two surfaces could have quite different slopes, and *ϵ*_{wpol} typically resembles neither of them. This is particularly noticeable at the later times shown, where for *R* > 0 a.u. reflections from the right-hand-side of the polaritonic surfaces yield complex overlaps of the component densities, and *ϵ*_{GD} can dominate the force in some regions. It sometimes ensures a reflection, which would not happen from purely *ϵ*_{wpol}.

To further demonstrate the importance of *ϵ*_{GD}, and with a view to the development of mixed quantum-classical methods based on the TDPES, we run quasiclassical dynamics on the exact TDPES and compare with such dynamics on just the weighted-polaritonic component. That is, we first sample the initial Wigner phase-space distribution corresponding to the initial Gaussian wavepacket to find 3000 initial positions and momenta^{53} (in fact, the results were well-converged with just 2000 trajectories). We then propagate each of these using classical Hamilton’s equations, $R\u0307=P,P\u0307=\u2212\u2207\u03f5(R,t)$, and make a histogram of the resulting positions. The snapshots of the results are plotted in Fig. 5, and a movie is provided in the supplementary material. On the left are the histograms when the trajectories are classically propagated using the exact TDPES, and we observe the results are remarkably accurate. There appears to be two small differences: The “trapped” part of the quasiclassical evolution appears to travel just a little further to the right, and at later times, the reflection is somewhat too enthusiastic. Still, given that such an evolution lacks adiabatic interference effects, quantum reflection and tunneling, and incurs errors when potentials are far from harmonic,^{53} the results here are very good. (We note that the non-adiabatic interference has been shown to be captured by such an evolution,^{54} but this is not a feature of the dynamics considered here.) In contrast, on the right are the histograms when instead only the gauge-invariant *ϵ*_{wpol} is used for the force on the classical trajectories. The agreement is much worse, not surprising from the discussion of Fig. 4; without the structures of *ϵ*_{GD}, the wavepacket splitting is hardly captured, nor is the detailed reflection at later times. The wavepacket propagates more or less in one piece to the right until it reaches the sharp rise of all the surfaces when the trajectories slow down before turning around.

### B. Cavity-induced ELEX

We now turn to the second set of parameters for our model.

Our initial state is again Gaussian in the nuclear variables, centered at *R* = −4 a.u., and placed on the upper electronic surface, with zero photons in the cavity: $\Psi (r,R,q,t)=Ne\u22122.85(R+4)2\Phi RBO(r)\xi (0)(q)$ where $\xi (0)=(\omega c/\pi )1/4e\u2212\omega cq2/2$ is the zero-photon state in the cavity. We choose our cavity-resonance at *ω*_{c} = 0.049 a.u., resonant with the BO frequency between the ground and first excited electronic states at *R* = ±2.2 a.u. and also resonant between the ground and second excited BO states at *R* = 0 a.u. (see the right-hand-side panel of Fig. 1).

The dynamics outside the cavity is relatively simple: the nuclear wavepacket slides down the first excited BO surface, meeting a weak avoided crossing at the origin, where a little population is transferred to the ground-state. A second transfer occurs later upon reflection back to this region. Figure 6 shows the BO populations, |*C*_{i}(*t*)|^{2}, where |*C*_{i}(*t*)|^{2} = ∫ *dR*|*C*_{i}(*R*, *t*)|^{2} with *C*_{i}(*R*, *t*) defined in Eq. (25). There is a very small transfer from the first to the second excited state, only just discernible in Fig. 6. Now, when the molecule is placed in the cavity, we see a more rapid and larger transfer to the ground-state surfaces along with a much greater rise in the population of the second excited BO state. Thus, this is an example of cavity-induced electronic excitation.

The phenomena can be only partially understood from the polaritonic surfaces shown in Fig. 1 because these surfaces look practically identical to the eye for much smaller coupling strengths where the dynamics is almost the same as outside the cavity. Still, the fact that both the ground-state population and the second excited population are enhanced at some times in the wavepacket evolution can be anticipated from the avoided crossings in the polaritonic surfaces and the nature of these surfaces: The wavepacket starts on the second Born–Oppenheimer surface coinciding with the third polaritonic surface at *R* = −4 a.u. and approaches the narrow avoided crossing with the second polaritonic surface around *R* = −2.2 a.u., suggesting a partial transfer to this surface. The part remaining on the third surface is associated with a 1-photon ground-electronic character to the right of the crossing, as is evident from the shape of that surface, aligning with the ground-state surface shifted up by the energy of one photon in that region. This transfer is consistent with the earlier transfer of population from the electronically excited to the ground electronic state compared with the cavity-free case. The enhancement of the third electronic state population then arises when the part of the wavepacket on the third polaritonic surface (that has the 1-photon ground electronic state character) approaches the narrow avoided crossing with the fourth polaritonic state at *R* = 0. Moving away from the crossing, this fourth surface has the character of the third electronic state.

We now turn to the exact TDPES for this problem, shown in time-snapshots of Fig. 7 along with the cavity-free surface. The effect of the cavity on the nuclear wavepacket dynamics is much smaller than in the PCET example except at later times, but to highlight the electronic character, we plot here also the *R*-resolved electronic populations, defined in Eq. (25). These are plotted as the thin lines on the split scale shown at the bottom of each panel. We see the growth of the coefficient of the third electronic state (red) beginning as the wavepacket goes through *R* = 0 a.u. as expected, and this coefficient is carried along in the tail of the wavepacket as it evolves to the right and then reflects. We see a very complex character develop by the end, especially in the left-hand-side tail of the wavepacket, which has a complex and oscillatory electronic structure. Again, the TDPES shows a piecewise polaritonic character when there are distinct single-surface dominated regions, but again in regions of mixed character, the polaritonic surface structure does not help in predicting the correct force, as given by the TDPES.

Similar observations as for the PCET case can be made for the force plots in Fig. 8, with the salient point being that the force from *ϵ*_{wpol} is, more often than not, countered by an opposing force from *ϵ*_{GD} in the mixed regions and does not resemble a weighted average of forces from polaritonic surfaces. Once again, the quasiclassical propagation on the exact TDPES in Fig. 9 performs remarkably well, while the importance of including the force from *ϵ*_{GD} is evident in its poorer prediction of the dynamics in the right-hand-side panel of Fig. 9.

## IV. CONCLUSIONS AND OUTLOOK

The exact potentials of the EF studied here provide the complete information to propagate the nuclei in cavity-QED problems. Our examples show how their structure correlates directly with the dynamics of the nuclear wavepacket for two different cavity-altered phenomena. Unlike polaritonic surfaces, these exact potentials distinguish initial states of different polaritonic characters and reflect, for example, the different degrees of suppression achieved when an initially polaritonic state is used compared with an initially pure electronic excitation. We note that other possible factorizations may lead to interesting insights, for example, taking the marginal factor to include both the photonic and nuclear degrees of freedom, with the conditional electronic wavefunction parametrically dependent on both degrees of freedom, *χ*(*R*, *q*, *t*)Φ_{R,q}(*r*, *t*). The resulting surfaces would provide corrections to the sometimes used “cavity-BO” surfaces.^{55} For real molecules of interest, finding the exact potentials of the EF method is as hard as solving the full TDSE exactly, so instead approximations need to be made. By studying the exact features in models where they can be found, such as the ones presented here, insight is gained into the effect that the various terms have on the nuclear dynamics and the results provide guidance for building approximate propagation methods. Similar features have been found in the exact surfaces driving the nuclei when coupled to electrons in the absence of cavities or external fields,^{38,41} which have then motivated and justified the mixed quantum-classical propagation methods of Refs. 37, 44, 45, and 56 and Refs. 46–49. This suggests that mixed quantum-classical methods along these lines can be generalized to include photons. The polaritonic surfaces play the role that the underlying Born–Oppenheimer surfaces would in that situation.

In such mixed quantum-classical methods, a key issue is to determine the force on the classical nuclei. While the Ehrenfest method takes a weighted-average over the Born–Oppenheimer surfaces, trajectory surface-hopping uses the gradient of a single surface at any time, switching between surfaces according to a stochastic algorithm. If applied now to problems in polaritonic chemistry, our results here stress that within an Ehrenfest-like approach, a weighted-average of polaritonic surfaces would lead to significantly incorrect dynamics, while if the gauge-dependent term was included to form the total TDPES, the results would improve significantly. However, it is important to realize that the results shown here for the quasiclassical propagation on the weighted polaritonic surfaces are at a higher level than a traditional Ehrenfest calculation since the latter involves an expansion of the conditional electronic wavefunction with nuclear-coordinate-independent coefficients, while the weights used in our calculations include coordinate dependence. The coordinate dependence of the weights affects the forces obtained from the gradients as we explicitly showed. That is, the forces in a traditional Ehrenfest calculation differ from those from the exact weighted-polaritonic surface even if the populations of the polaritonic states averaged over the system are the same.

A traditional surface-hopping scheme, on the other hand, would be even more dogged by the problems of coherence than in the usual cavity-free case due to the increased number of avoided and trivial crossings between the polaritonic surfaces compared to the Born–Oppenheimer ones.^{57} Both decoherence and recoherence effects would need to be captured for an accurate determination of the dynamics. The EF-based surface-hopping approach of Refs. 46–49 as well as the coupled-trajectory scheme of Refs. 37, 44, and 45 can be generalized to the polaritonic case. Here, an additional term, derived from the EF, appears in the electronic equation of motion and how well it captures these effects will be explored in future work, as will other methods of developing EF-based mixed quantum-classical methods for polaritonic chemistry.

## SUPPLEMENTARY MATERIAL

Movies of the dynamics are provided in the supplementary material as follows:

**Movie Fig. 3.mp4**—upper row from left to right:*ϵ*(*R*,*t*),*ϵ*_{wpol}(*R*,*t*), and*ϵ*_{GD}(*R*,*t*) for cavity free (pink), in-cavity with the 0-photon excited electronic initial state (dark blue), and in-cavity with the polaritonic initial state (black). The backdrop of static polaritonic surfaces is shown in light blue for reference (as in Fig. 3). Lower row from left to right: 0-photon resolved density and 1-photon resolved density in the same color code as the upper row, and the*R*-resolved Born–Oppenheimer population for (a) in-cavity with the polaritonic initial state, (b) in-cavity with the 0-photon excited electronic initial state, and (c) cavity-free with |*C*_{1}(*R*,*t*)|^{2}(green), |*C*_{2}(*R*,*t*)|^{2}(orange), and |*C*_{3}(*R*,*t*)|^{2}(red) (the same color code as Fig. 7).**Movie Fig. 5.mp4**—left panel: Quasiclassical propagation of 3000 trajectories (orange) on the total TDPES*ϵ*(*R*,*t*) for the cavity-induced suppression of PCET (case IIIA). These are plotted as a function of*R*(*a*.*u*.) against exact solution (blue), Born–Oppenheimer ground-state population (red), and Born–Oppenheimer first-excited state population (green). Right panel: the same as above, propagated using the weighted-polaritonic PES*ϵ*_{wpol}(*R*,*t*).**Movie Fig7. mp4**—upper row from left to right:*ϵ*(*R*,*t*),*ϵ*_{wpol}(*R*,*t*), and*ϵ*_{GD}(*R*,*t*) for cavity free (pink) and in-cavity with the 0-photon excited electronic initial state (black). The backdrop of static polaritonic surfaces is shown in light blue for reference (as in Fig. 7). Lower row from left to right: 0-photon resolved density and 1-photon resolved density in the same color code as the upper row, and the*R*-resolved Born–Oppenheimer population for (a) with the 0-photon excited electronic initial state and (b) cavity-free with |*C*_{1}(*R*,*t*)|^{2}(green), |*C*_{2}(*R*,*t*)|^{2}(orange), and |*C*_{3}(*R*,*t*)|^{2}(red) (the same color code as Fig. 7).**Movie Fig9. mp4**—left panel: Quasiclassical propagation of 3000 trajectories (orange) on the total TDPES*ϵ*(*R*,*t*) for the cavity-induced electronic excitation (case IIIB). Plotted as a function of*R*(*a*.*u*.) against the exact solution (blue), Born–Oppenheimer ground-state population (red), and Born–Oppenheimer first-excited state population (green). Right panel: the same as above, propagated using the weighted-polaritonic PES*ϵ*_{wpol}(*R*,*t*).

## AUTHORS’ CONTRIBUTIONS

P.M. and B.R. contributed equally to this work.

## ACKNOWLEDGMENTS

Financial support from the US National Science Foundation (Grant No. CHE-1940333) (N.T.M., N.M.H., and B.R.), the Department of Energy, Office of Basic Energy Sciences, Division of Chemical Sciences, Geosciences and Biosciences under Award No. DE-SC0020044(L.L.), and the RISE program at Hunter College (Grant No. 5R25GM060665-20) is gratefully acknowledged.

## DATA AVAILABILITY

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