key: cord-0471354-d0m36hjy
authors: Rosenzweig, Bart; Hoffmann, Norah M.; Lacombe, Lionel; Maitra, Neepa T.
title: Analysis of the Classical Trajectory Treatment of Photon Dynamics for Polaritonic Phenomena
date: 2021-11-23
journal: nan
DOI: nan
sha: 3021f14adcb17f0005c3b96c602ade51d8accaaf
doc_id: 471354
cord_uid: d0m36hjy

Simulating photon dynamics in strong light-matter coupling situations via classical trajectories is proving to be powerful and practical. Here we analyze the performance of the approach through the lens of the exact factorization approach. Since the exact factorization enables a rigorous definition of the potentials driving the photonic motion it allows us to identify that the cause of the underestimation of photon number and intensities observed in earlier work is primarily due to an inadequate accounting of light-matter correlation in the classical Ehrenfest force rather than errors from treating the photons quasiclassically per se. The latter becomes problematic when the number of photons per mode begins to exceed a half.

The impressive advances in experimental realizations of strong light-matter coupling induced by confinement continue to open exciting new possibilities for the manipulation of matter [1] [2] [3] [4] . This has led to the establishment of a new field of research in chemistry in the past decade, known as polaritonic chemistry, or molecular polaritonics. The possible applications are everincreasing and just within the last year we have seen experimental demonstrations of new phenomena in different directions, including tunable self-assembled polaritons in microcavities [5, 6] , solar light-harvesting enabled via cavity-coupling [7] , and large enhancement of ferromagnetism from collective strong-coupling [8] . In most cases, to achieve strong light-matter coupling a large number of molecules is required, exploiting the collective nature of the coupling that scales as the square-root of the number of molecules. The coupling strength also scales as the inverse square-root of the volume, so to reach strong-coupling with single molecules, a tighter confinement is required, and here plasmonic cavities have been a successful solution [9] [10] [11] [12] .

A key issue for plasmonic cavities is the large loss entailed from the electron oscillations in the metal, and experimental and theoretical works are on-going to deal with these, and even sometimes take advantage, e.g. [13] [14] [15] [16] [17] [18] [19] . Another challenge in the computational modeling of polaritonic systems in general is the efficient accounting of the many photon modes in the cavity that can influence the dynamics. In this work, we focus on this latter issue; in particular, on the treatment of photons via classical trajectory dynamics.

When more than one cavity mode plays an important role, a full quantum treatment of the photonic degrees of freedom becomes challenging. Different approaches have been explored, including quantum electrodynamical density functional theory [20] [21] [22] , cumulant expansions [23] , tensor networks [24] , and quasinormal modes [25] when accounting also for dissipation. Inspired by mixed quantum-classical methods used for electron-nuclear coupled motion, the multitrajectory Ehrenfest (MTE) approach was proposed [26, 27] and demonstrated to be an efficient and powerful approach. In MTE, the photonic system is represented by the displacement-field coordinate q and its conjugate momentum p, which are evolved using classical equations of motion. The classical photonic trajectories are coupled to the evolution of the matter system via Ehrenfest forces; the matter may be treated quantummechanically which is usually the choice when only electronic degrees of freedom are explicitly considered, or also with mixed quantum-classical approaches where nuclear dynamics are also treated classically. A key aspect is the initial sampling of the photonic degrees of freedom: it has been shown that, starting in a zero-photon state, Wigner-sampling of the initial state is essential to capture the vacuum-fluctuations of the field and approximate photon-matter correlations well enough to describe, for example, the ensuing spontaneous emission of the matter system [26] [27] [28] .

This idea of Ref. [26, 27] to use a distribution of classical trajectories for the photonic system has enabled the investigation of new phenomena arising from strong light-matter coupling, including the importance of self-polarization when accounting for many modes in the cavity [28] , and modifications of chemical reaction rates from vibrational strong-coupling [29, 30] where only the ground Born-Oppenheimer state is considered. Describing photons quasi-classically might seem, on the one hand, very surprising given that the photon is inherently a massless quantum particle. On the other hand, resolving the Hamiltonian in the photon displacement-coordinate representation renders the free-photon Hamiltonian a harmonic form, which, together with the bi-linear form of the light-matter coupling in the usual non-relativistic dipole-approximation treatments, makes it highly suitable for quasiclassical approximations. Indeed, classical Wigner propagation of the initial Wigner distribution is exact for quadratic Hamiltonians [31] . However, even if the Hamiltonian has up to quadratic terms only in photonic displacement coordinate q, coupling to the matter-system dynamics implicitly introduces a more complicated dependence on q beyond quadratic, so we do not expect quasiclassical trajectory calculations to be exact despite the apparently quadratic form of the Hamiltonian.

A closer look at the earlier studies indicate that the MTE approach, although quite accurate, tends to underestimate the photon number and intensities [26] [27] [28] . To assess this, having an exactly numerically-solvable system to check against, as in these works, is extremely valuable. For example, in the studies of spontaneous emission from a two-level or a three-level system coupled to 400 cavity modes initially in vacuum, the intensity of the emitted light in the initial emission event is underestimated by about two-thirds (Figs. 3, 4, 7 and 8 of Ref. [26] ). A similar underestimation is seen in the photon number in the single-mode study of cavityinduced suppression of proton-coupled electron transfer in Fig 2d of Ref. [28] . A question arises: is this due to the quasiclassical approximation itself or from an inadequate accounting of light-matter correlation in the equations of motion driving the photons? Quasiclassical approaches lack phases and interference effects, cannot describe classically-forbidden processes such as tunneling, and are liable not to fully hold onto zero-point energy. The Wigner treatment in particular develops errors from strong anharmonicity. As mentioned earlier, although the terms in the light-matter Hamiltonian used in the computational models that involve the photonic displacement coordinate are purely harmonic with bi-linear coupling, the feedback with the matter system through this coupling effectively induces anharmonic effects on the photonic system.

Whether it is the quasiclassical treatment itself that is largely responsible for the underestimation in the photon number and intensities, or whether it is instead largely due to incorrect forces in the classical equations driving the photon dynamics has not yet been determined. The latter reason was suggested in Ref. [28] by a qualitative argument based on the exact factorization approach [32] [33] [34] [35] . In this approach exact timedependent potentials are defined through exact timedependent Schrödinger equations for coupled quantum subsystems. The original electron-nuclear formalism was extended to include photons in Ref. [34] : depending on how the factorization is done, one can then obtain exact time-dependent potential energy surfaces (TDPES) for the nuclear dynamics [36, 37] , for the electron dy-namics [38] , or for the photon dynamics [34] . In the present work, we take a closer look at the exact potential driving the photon dynamics which we call the qT-DPES. We show that the underlying force in the MTE treatment of photons significantly differs from the exact, missing terms from light-matter correlation. In particular, the neglect of the photonic coordinate dependence of the coefficients that weigh the effective cavity-Born-Oppenheimer surfaces in evaluating the Ehrenfest forces has paramount importance. On the other hand, it is the quasiclassical approximation itself that leads to errors when the photon number per mode exceeds about a half.

We begin in section II by presenting the model we will study and its exact solution. In section III we introduce the particular flavor of the exact factorization approach used in this work, different approximate quasiclassical analyses of it and the MTE, before presenting the results and analysis in Sec. IV. The implications for more realistic systems are discussed in Sec V, with prospects for effective practical methods based on the exact factorization approach.

The particular model we consider is simply one photon mode coupled resonantly to a two-level system, with Hamiltonian H = H m + H p + H mp where H m = g |g g| + e |e e|, H p = ω 2 c q 2 − ∂ 2 /∂q 2 /2 H mp = r eg λω c q (|e g| + |g e|)

Here |g(e) are the two electronic states, while q represents the displacement-field coordinate operator of the photonic system. Eq. (1) follows from the nonrelativistic Pauli-Fierz Hamiltonian [20] [21] [22] 39] , under the condition that the electronic states are parity eigenstates, and neglecting any dependence on nuclear variables. (No rotating wave approximation is made). The coupling matrix element r eg = e|r|g wherer is the (possibly many-)electron position operator. The Pauli-Fierz Hamiltonian also contains a self-polarization term, but in this simple two-level case, it can be absorbed into the matter part H m , as a renormalization of the energies g and e . We choose a cavity resonant with the electronic energy difference, ω c = e − g . The exact solution may be expressed as the expansion

where |χ d g(e) | 2 represent q-resolved populations associated with the uncoupled ground and excited electronic states respectively. The | notation refers to states in the electronic Hilbert space only. With analogy to the electron-nuclear problem, we may think of χ d g(e) (q, t) as the photonic state coefficients in a "diabatic" expansion of the electronic system, since the electronic basis chosen in Eq. (2) has no dependence on q. Alternatively, one can express the exact solution as

where

are cavity-Born-Oppenheimer states in the electronic Hilbert space, with

Eq. (3) is analogous to the Born-Huang expansion in adiabatic states in the electron-nuclear case.

The exact time-dependent solution is then given either by substituting Eq. (2) or Eq. (3) into the TDSE H|Ψ = i∂|Ψ ∂t:

and analogously for χ d e (q, t), or

where the non-adiabatic coupling matrix elements are

In the MTE approach, the displacement-coordinate of the photonic system is treated as a classical coordinate, and together with its conjugate variable p, define the phase-space for the classical trajectories. An ensemble of independent classical trajectories, initially sampled according to the initial photonic displacement density in q and p, is propagated guided by a mean-field force determined by the electronic system, while the quantum nature of the electronic system is retained, driven by a potential parametrically dependent on the photon coordinate. Analogously to the Ehrenfest equations for the electron-nuclear problem, one obtains the equations for the Ith photonic trajectory in the ensemble:

where |Φ (I) (t) = |Φ q (I) (t) , and the instantaneous classical position q (I) (t) replaces the photonic operator q in the Hamiltonian H qBO . It is convenient to expand |Φ (I) (t) in either the diabatic or Born-Huang basis, as in Eqs.

(2)-(3) where the projected wavefunctions are now purely time-dependent coefficients associated with each trajectory, e.g.

for the diabatic case, replacing the superscript d notation for this basis choice by a tilde, for convenience. Using the diabatic expansion, the equations becomë (11) (where the t-dependence of q (I) , c g,(e) is omitted on the right of the second equation to avoid notational clutter). Instead, using the qBO basis,

we geẗ

with

where the q I -dependence of terms on the right is omitted for notational simplicity).

The Ehrenfest approach for the electron-nuclear problem was derived from treating nuclear degrees of freedom classically starting from a product ansatz of the molecular wavefunction [40, 41] ; here we have applied it to the photonic system. The potentials yielding the forces on the photonic displacement coordinate on the right-hand-side of Eq. (8) have a mean-field character, seen perhaps most clearly in Eq. (13) where the first two terms explicitly show the weighted BO-forces with an additional term coming from non-adiabatic transitions driven by the coupling d eg (q (I) ).

In the electron-nuclear case, the analogous mean-field force on the nuclear coordinate fails to be accurate after the system has propagated through non-adiabatic interaction regions where the electronic eigenstates are strongly coupled. These regions tend to be somewhat localized in space, while we will see in Sec. IV (Fig. 1 ) that in the present electron-photon case, the non-adiabatic coupling extends throughout the photonic displacement wavepacket. The Ehrenfest force is approximate at best, and instead, the exact factorization defines the exact force that drives the photonic system [32] [33] [34] .

In the exact factorization approach, the exact timedependent wavefunction of a system of coupled quantum subsystems is written as a single correlated product of marginal and conditional factors: in the present case,

where the marginal χ(q, t) represents the photonic wavefunction and |Φ q (t) is the electronic state conditionally dependent on the photonic displacement-field coordinate q. The equality on the right is referred to as the "partial normalization condition", and indicates that the conditional wavefunction is normalized to 1 for every q at each time t; the notation ...|... means an inner-product only over the electronic Hilbert space. Despite being a single product (in contrast to Eq. (3)), Eq. (14) represents the exact wavefunction, and the two factors are unique up to a (q, t)-dependent phase factor: e iθ(q,t) can multiply χ(q, t) while e −iθ(q,t) multiplies |Φ q (t) to yield the same |Ψ(q, t) = χ(q, t)|Φ q (t) . The equations for the two components involve terms that depend on the other factor, which intimately represent the electron-photon correlation; these can be found in Ref. [34] and are straightforward generalizations of the electron-nuclear case to the electron-photon one. In particular, with the present choice of factorization, the equation for χ(q, t) has a time-dependent Schrödinger form with a vector potential and a scalar potential.

The phase-freedom appears as a gauge-like freedom in the electron-nuclear coupling potentials, and, in some cases, a gauge can be selected such that only a scalar potential appears. This is always the case when the marginal has only one degree of freedom. (We note that in the general case of N photon modes in the cavity, where the marginal is N -dimensional, the vector potential is also N -dimensional, and has the form A ν = Φ q (t)|i∂ qν Φ q (t) . Even though each mode may be onedimensional, the vector potential may still have a rotational component. ) With only one mode in our model system, we may choose a gauge where A = 0, and then we find that the potential in the Schrödinger equation for χ(q, t) is purely scalar and time-dependent, and we denote it the qTDPES:

and

This allows us to define an exact force driving the photons in an exact-classical-trajectory approach, viä

A question we will ask, is how well do the classical trajectories for the photons evolving under the exact force, Eq. (19) , do, compared with the Ehrenfest photon trajectories, Eq. (13) or Eq. (11) , in reproducing the true photonic dynamics.

To get more insight into the exact force, it is instructive to expand the conditional electronic wavefunction in terms of the qBO wavefunctions that were defined in Eq. (4):

where, due to the partial normalization condition, |C g (q, t)| 2 + |C e (q, t)| 2 = 1 for every q and t. In terms of these coefficients, we observe that the force from the weighted BO contribution to E qTDPES (q (I) , t) somewhat resembles the force in the MTE:

We will see that the two qBO surfaces are almost parallel, both very close to harmonic with deviations only at larger q while the coefficients show a large qdependence near q = 0 (see Sec. IV): as a result, the wBO force is dominated by an almost linear restoring force (first term) together with a large repulsive force away from q = 0 from the last term in Eq. (22) . Comparing now with Eq. (13) we see that the Ehrenfest force has the similar structure of a weighted force from the first term on the right of Eq. (21), but a crucial difference in the last term: while the local shape of the coefficients plays a large role in the wBO force, this is completely absent in the Ehrenfest force, which instead has a term depending on the non-adiabatic coupling between the two electronic states. We will find that this term is much smaller than that from the q-dependence of the coefficients appearing in the wBO term and, in contrast to the large q-dependence of this term, this term in the MTE gives an almost uniform force in q. The Ehrenfest method involves independent trajectories and the notion of how a coefficient associated with one trajectory varies with that of a neighboring trajectory does not arise with these uncoupled trajectories.

The main observable of interest in this study is the photon number, N , defined by N = a † a = ωq 2 + p 2 /ω /2 − 1/2

(As mentioned earlier, there is no contribution from selfpolarization in this model). The distribution of photonic trajectories in the q, p phase space in the MTE calculation directly yields the expectation values on the right from the its variances in position and momentum. When extracting these from the exact factorization framework, a subtlety arises: the marginal wavefunction of the factorization, χ(q, t), reproduces the photonic density and photonic current-density as would be obtained from the full photon-electron wavefunction, |Ψ(q, t) , but for general non-multiplicative observables of the photonic system, the q-dependence of the conditional wavefunction means that there are additional terms to the usual operators. Specifically [33] , for our choice of gauge,

where p 2 χ = − dqχ * (q, t)∂ 2 q χ(q, t) is twice the kinetic energy of the marginal. The second term in Eq. (24) must be added to this in order to get the true expectation value of p 2 . This applies both for cases where a quantum calculation is made for the marginal χ(q, t), as well as when it is computed through a quasiclassical trajectory distribution in (q, p)-phase space; in the latter case the q-density weighted average of E kin must be added to the variance in p.

We will now analyze the results of the MTE approach for photons by comparing against quantum and quasiclassical propagation on the exact qTDPES, as well as from propagation on only the wBO component. Due to the simplicity of the model system of Sec. II, the numerically exact solution for |Ψ(q, t) is easily obtained. The exact qTDPES and the wBO component are found by inversion: writing χ(q, t) = |χ(q, t)|e iS(q,t) , where |χ(q, t)| = Ψ(q, t)|Ψ(q, t) , and ∂ q S(q, t) = Im Ψ(q, t)|∂ q Ψ(q, t) /|χ(q, t)| 2 , we can find |Φ q (t) = |Ψ(q, t) /χ(q, t) to insert into Eqs. (15)- (18) . The expression for the phase S(q, t) follows from the gauge choice A = 0, as in Refs. [34, 42] .

We choose parameters for the model of Eq. (1) as e = 0.2, g = −0.2, ω c = 0.4a.u. and r eg λ = 0.01a.u. The resulting dynamics gives Rabi-like oscillations, with a period of T ≈ 2π/0.01.

Before we examine the results of the quasiclassical propagations, we first compare the underlying potentials and forces in MTE dynamics with those arising from the exact factorization.

The top left panel of Fig. 1 shows the two qBO surfaces, which have a form that deviates from harmonic by less than 0.01% in the range shown: diagonalization yields E qBO (q) = ω 2 c q 2 /2 ∓ ω 2 c /4 + (r eg λω c ) 2 q 2 , which shows that the deviation grows with q 2 . The plot also shows the non-adiabatic coupling d eg (q), indicating that the MTE force in Eq. (13) is dominated by the first term, and also that, unlike typical electron-nuclear case, the non-adiabatic coupling in the electron-photon case is delocalized in q-space. The force acting on the MTE trajectories can be compared with that acting on trajectories evolving on the exact qTDPES and wBO surfaces, and to discuss these we must discuss the initial state since these surfaces are both time-dependent and state-dependent.

We begin with the system initially in the state |Ψ(q, 0) = χ g (q)|Φ q,e where the photon ground state is the zero-photon state χ g (q) = ω π 1/4 exp(−ωq 2 /2). Because the conditional electronic state is initially in the excited eigenstate of H qBO , the exact qTDPES and wBO surfaces coincide with the upper qBO surface at the initial time, but quickly begin to deviate from it. After short times, anharmonicity develops at larger q in the wBO and qTDPES, which begin to move towards the lower qBO surface, while at smaller q, the surfaces remain close to the upper qBO surface. This is consistent with the bi-linear nature of the electron-photon coupling scaling with q, and we see that the q-resolved electronic population of the upper qBO surface, shown in the figures, pulls away from the excited state to the lower at larger q. Over time, a barrier builds up near q = 0, at which point the wBO resolutely sticks to the upper qBO surfaces where the electron-photon coupling is locally zero, while everywhere else collapsing to the lower qBO surface, associated with the photon emission. The inclusion of the kinetic and gauge-dependent terms to form the full qTDPES, lead to an even more pronounced and increasingly localized barrier at q = 0, arising primarily from E kin (q, t) [34] , which becomes singular at T /2. During this half-cycle, the photon emission character is evident in |χ(q, t)| 2 which evolves from the harmonic ground-state form to the first-excited state (onephoton state), and also in the population of the excited state |C e (q, t)| 2 which evolves from initially being 1 at all q to 0 at all q except for q = 0.

The correct forces in a quasiclassical description of the photons are thus significantly different to what is in the MTE approach. Clearly the anharmonicity and developing barrier in the wBO and qTDPES surfaces will yield different forces to that of the MTE, leading to a broader distribution in q, which implies greater q 2 and greater photon number. Relating this back to Eq. (22) , where the wBO force is decomposed in terms of qBO components and coefficients, we had noted that the first term is shared with the MTE, the second term is negligible (also clear from the figure), while the third term gives rise to the anharmonicity and barrier: the q-dependence of the electronic coefficient, which intimately reflects the electron-photon correlation, is intimately is missing in MTE.

The lack of this electron-photon correlation barrier and large-q widening in the potential underlying MTE, is expected to lead to an underestimation in photon number. A separate question is then how well does the quasiclassical treatment itself, in contrast to a quantum treatment, perform for dynamics with such anharmonic structures? Even if the exact factorization provides the correct forces that a classical particle would experience, if the classical approximation for the trajectory ensemble is itself poor, it is not of practical use. In Sec. IV B we turn to the propagation itself. 

We now verify the predicted underestimation by running quasiclassical dynamics on the exact factorization surface qTDPES, and the wBO component, and comparing these with MTE. It is also instructive to run a quantum dynamics on the qTDPES and the wBO component: we will find that as the photon emission approaches about half a photon and proceeds beyond that, the quasiclassical approximation itself becomes inaccurate, due to the growing anharmonicity of the potential that includes the increasingly narrow and sharp barrier near q = 0.

The quantum dynamics is computed with the splitoperator method with a time-step of 0.001a.u. and a spatial grid-spacing of 0.1a.u. The quantum dynamics run on the full qTDPES agreed with the results from the original full molecular propagation, which served as a test that our inversion to find the exact qTDPES was accurate. We refer to this below as the exact quantum dynamics. The quasiclassical propagation was carried out by solving Hamilton's equations with a time-step of 0.02 a.u. via leapfrog integration with 10,000 trajectories (initially sampled using the Gaussian initial state).

The top left panel and inset of Fig. 2 shows the photon number during a full Rabi period, computed from the exact quantum dynamics which shows the expected emission of one photon after half a Rabi period. The quantum dynamics on the wBO component appears to be reasonably accurate until about half a photon is emitted, with a slight underestimate, which would be expected from the smaller barrier at q = 0 observed in the potential in Fig. 1 . However after about a quarter Rabi period the quantum propagation on the wBO begins to greatly overestimate N , largely from the E kin contribution to p 2 : the peak in this term at q = 0 grows increasingly sharp and localized as the one-photon state is reached, and because the photonic density remains significant in this region due to the lack of this term in determining the trajectory evolution, it samples the peak and leads to large values. This is also clear from Fig. 3 , where in the second row, the green dashed curve shows the photon displacement-coordinate density obtained from quantum propagation on the wBO surface, and the fourth row shows the wBO surface itself as the green dashed curve, at four time-snapshots. The small barrier at q = 0 can be tunneled through and so leads to only a small dip in the density there at about a quarter-Rabi period. Whether this problem appears in a selfconsistent calculation, where the terms in the qTDPES are computed by a self-consistent conditional electronic state rather than the exact one, is a question for future investigation.

Turning now to the quasiclassical propagations, we see from Fig. 2 that for the first quarter Rabi period or so, quasiclassical propagation on both the full qTDPES and wBO are reasonably accurate, especially for the photon number where there is some error cancellation between a small overestimate of q 2 (t) and underestimate of p 2 (t) . The MTE underestimates the photon number N (t) and the intensities q 2 (t) , p 2 (t) . The time snapshots in Fig. 3 enable a deeper analysis and movies are also provided in the Supplementary Information.

Let us first discuss the MTE case. The top figures in each of the four panels in Fig. 3 shows the MTE photonic displacement-coordinate density, which shows a "breathing" that leads to the larger width, but remain at all times peaked around q = 0 in stark contrast to the exact photonic density. The anharmonicity at large q in the first term of Eq. 13 along with the non-adiabatic coupling cause the density to spread, somewhat similar to a harmonic oscillator that widens (and narrows in the second half of the cycle). But without any bar- rier term, MTE is unable to even remotely capture the correct structure of the density. The increase in the averaged quantities such as q 2 and N as a photon is emitted give an overly rosy picture of the photonic dynamics, which likely lead to errors in other properties of the radiation field. The MTE underestimation was expected from the discussion of Sec. IV A and here we see it explicitly in the dynamics. For the photon number, it is of similar magnitude seen in the previous works on a model molecule and with many modes, as discussed in Sec. I for the specific examples where exact results were available for comparison.

We now consider the quasiclassical dynamics on the exact qTDPES. Looking at the time-snapshots in Fig. 3 in the third row, we notice that the density splits too quickly in the early dynamics (t = 100, 160au in the figure). This may be explained by the fact that the trajectories are pushed away from the potential barrier developing in the middle, while the quantum wavepacket can still live in the classically forbidden region. The quasiclassical underestimates the outer tail regions of the density slightly but overestimates the peaks in order to compensate the density that has been shunned away from the barrier. When the potential barrier becomes larger, we observe that the classical trajectories tend to overlocalize in each (partial) well. The two errors somewhat compensate in the averaged quantities -less density in the classically-forbidden region near q = 0, but also less density in the tails -leading to the reasonably accurate averages observed earlier in Fig. 2 .

The average N (t) from the quasiclassical propagation on the wBO surface appears to performs about the same for much of the quarter Rabi period, as that on the full qTDPES, due to an error cancellation effect: the underlying potential displays a smaller barrier at the origin so the trajectories are not repelled away as much. In fact the quasiclassical q-resolved displacement density is more accurate for propagation on wBO than for propagation on the full qTDPES, as evident from the second row in Fig. 3 . At larger times, there is an underestimate in q 2 which is compensated by an overestimate of p 2 which comes from the E kin peak in Eq. (24) which is sampled more because the photonic density does not split as much.

Although the present work considered only the very simplest case of a two-level system coupled to a single photon-mode, we expect that the essential results generalize to realistic systems with many photon-modes, many electronic levels, and nuclear degrees of freedom: an MTE treatment of photons tends to underestimate the photon number and intensities, due to an inadequate accounting of matter-photon coupling in the treatment of the independent-trajectory Ehrenfest coefficients. The underlying photonic displacement-coordinate density does not show the character expected from the partial emission of a photon. The underestimation observed in the simple model has about the same magnitude as that observed in more complex systems where exact results were available for comparison [26] [27] [28] .

In contrast, a quasiclassical evolution on the weighted qBO surfaces, wBO, or full qTDPES performs much better for photon numbers of less than a half per mode, while for larger photon numbers per mode, large errors arise from the quasiclassical approximation itself that performs poorly when the underlying potential changes rapidly in space and becomes far from locally quadratic. The quasiclassical calculations were performed here within a particular gauge choice for the exact factorization (zero vector potential), and an open question is whether a different gauge choice would result in more FIG. 3 . Time snapshots of the photonic displacement coordinate density, computed via MTE (top), quasiclassical (orange) and quantum (green dashed) on the wBO surface along with the exact quantum (black) for reference (second row), and the quasiclassical (orange) and quantum (black) on the full qTDPES (third row), shown at a) t = 100au, b) t = 160au, c) t = 240au, d) t = 350au. The fourth row shows the exact qTDPES (black) along with the wBO component (green dashed).

accurate results at longer times. The sharp barrier forming at q = 0 arises in the gauge-independent term E kin however while in this case the gauge-dependent term is relatively small [34] .

Still, the computational efficiency and reasonable accuracy of MTE treatment of photonic dynamics in cavity-QED is encouraging [26] [27] [28] [29] , and the present results should be taken more as an explanation of its errors and as a note of caution that in such calculations the photon numbers are likely underestimated. Global quantities obtained from averaging over the Ehrenfest trajectories may appear more accurate than the underlying photonic displacement-coordinate density, which likely has a qualitatively incorrect structure. (We also note that with MTE methods, zero point energy leakage from higher modes to lower modes may occur, and should be checked; this was not a problem with the single-mode calculation here).

Again, it is important to note that the quasiclassical calculations on the exact-factorization surfaces are not self-consistent ones, i.e. we are using the surfaces extracted from inversion of an exact calculation and simply running quasiclassical dynamics on them. We are using them as an analysis tool rather than as a method in itself. A self-consistent mixed quantum-classical calculation based on the exact factorization will be investigated in the future, for example extending those developed for the electron-nuclear problem [35, [43] [44] [45] [46] [47] . 

Quantum electrodynamics in the coulomb gauge