Analysis of the classical trajectory treatment of photon dynamics for polaritonic phenomena

The impressive advances in experimental realizations of strong light–matter coupling induced by a confinement continue to open exciting new possibilities for the manipulation of matter.^1–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 ever-increasing, 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,⁷ and large enhancement of ferromagnetism from collective strong coupling.⁸ In most cases, to achieve strong light–matter coupling, a large number of molecules are 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–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 ongoing to deal with these and even sometimes take advantage, e.g., Refs. 13–23. 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. Although polaritonic phenomena are largely driven by resonances between the molecule and cavity modes, self-polarization effects involving all modes could also distort the potential energy landscapes,²⁴ and more needs to be understood about when many modes need to be accounted for. In this work, we focus on a way to treat many modes efficiently 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,^25–27 cumulant expansions,²⁸ tensor networks,²⁹ and quasinormal modes,³⁰ when accounting also for dissipation. Inspired by mixed quantum–classical methods used for the electron-nuclear coupled motion, the multi-trajectory Ehrenfest (MTE) approach was proposed^31,32 and demonstrated to be an efficient and powerful approach. In MTE, each photonic mode is represented by its displacement-field coordinate q and 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 quantum-mechanically, which is usually the choice when only the electronic degrees of freedom are explicitly considered, or also with mixed quantum–classical approaches where the 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 to approximate photon–matter correlations well enough to describe, for example, the ensuing spontaneous emission of the matter system.^24,31,32

This idea of Refs. 31 and 32 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 possible importance of self-polarization when accounting for many modes in the cavity,²⁴ and modifications of chemical reaction rates from vibrational strong coupling^33,34 where only the ground Born–Oppenheimer state is considered. Describing photons quasiclassically 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 bilinear form of the light–matter coupling in the usual nonrelativistic dipole-approximation treatments, makes it highly suitable for quasiclassical approximations. Indeed, classical propagation of the initial phase-space points sampling the initial Wigner distribution gives the exact propagation of the Wigner function when the Hamiltonian is quadratic.³⁵ However, even if the Hamiltonian has up to quadratic terms only in the 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 apparent quadratic form of the Hamiltonian.

A closer look at the earlier studies indicates that the MTE approach, although quite accurate, tends to underestimate the photon number and intensities.^24,31,32 To assess this, having an exact 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. 31). A similar underestimation is seen in the photon number in the single-mode study of cavity-induced suppression of proton-coupled electron transfer in Fig. 2(d) of Ref. 24. 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 the 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 bilinear 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. 24 by a qualitative argument based on the exact factorization approach.^36–39 In this approach, exact time-dependent potentials are defined through exact time-dependent Schrödinger equations for coupled quantum subsystems. The original electron-nuclear formalism was extended to include photons in Ref. 38; depending on how the factorization is done, one can then obtain exact time-dependent potential energy surfaces (TDPES) for the nuclear dynamics,^40,41 the electron dynamics,⁴² or the photon dynamics.³⁸ In this work, we take a closer look at the exact potential driving the photon dynamics, which we call the qTDPES. We show that the underlying force in the MTE treatment of photons significantly differs from the exact, missing terms from the 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 Sec. II by presenting the model we will study and its exact solution. In Sec. 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 a one-photon mode coupled resonantly to a two-level system with the Hamiltonian H = H_m + H_p + H_mp, where

Here, $g(e)$ are the two electronic states, while q represents the displacement-field coordinate operator of the photonic system. Equation (1) follows from the nonrelativistic Pauli–Fierz Hamiltonian^25–27,43 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 $reg=er̂g$⁠, where $r̂$ 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 $|χg(e)d|2$ represent the q-resolved populations associated with the uncoupled ground and excited electronic states, respectively. The $$ notation refers to the states only in the electronic Hilbert space. With an analogy to the electron-nuclear problem, we may think of $χg(e)d(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 the cavity-Born–Oppenheimer states in the electronic Hilbert space with

Equation (3) is analogous to the Born–Huang expansion in the 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$ as follows:

and analogously for $χed(q,t)$ or

where the nonadiabatic coupling matrix elements are $Dij(q)=Φq,i∂q2Φq,j〉/2$ and $dge(q)=Φq,g∂qΦq,e〉$⁠.

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 as follows:

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 bases, as in Eqs. (2) and (3), where the projected wave functions 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 become

where the t-dependence of q^(I), c_g,(e) is omitted on the right-hand side of the second equation to avoid notational clutter. Instead, using the qBO basis

we get

with $ΔEqBO=EeqBO−EgqBO$⁠, where the q^I-dependence of terms on the right-hand side is omitted for notational simplicity.

The Ehrenfest approach for the electron-nuclear problem was derived from treating the nuclear degrees of freedom classically starting from a product ansatz of the molecular wave function;^44,45 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, shown perhaps most clearly in Eq. (13), where the first two terms explicitly show the weighted BO forces with an additional term coming from nonadiabatic 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 the nonadiabatic 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 nonadiabatic 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.^36–38 

In the exact factorization approach, the exact time-dependent wave function 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 wave function and $Φq(t)$ is the electronic state conditionally dependent on the photonic displacement-field coordinate q. The equality on the right-hand side is referred to as the “partial normalization condition” and indicates that the conditional wave function 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 wave function, 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. 38 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 one-dimensional, 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 as follows:

where

and

This allows us to define an exact force driving the photons in an exact classical trajectory approach via

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 wave function in terms of qBO wave functions that were defined in Eq. (4), which, now again for two states, are

where, due to the partial normalization condition, |C_g(q, t)|² + |C_e(q, t)|² = 1 for every q and t. In terms of these coefficients, we observe that the force from the weighted BO contribution to $EqTDPES(q(I),t)$ somewhat resembles the force in the MTE as follows:

We will see that the two qBO surfaces are almost parallel, both very close to the harmonic with deviations only at larger q, while the coefficients show a large q-dependence 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 a similar structure of a weighted force from the first term on the right-hand side 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 nonadiabatic 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

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

where $⟨p2⟩χ=−∫dqχ*(q,t)∂q2χ(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². 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 the (q and  p) phase space; in the latter case, the q-density weighted average of $Ekin$ must be added to the variance in p. We note that due to the Ehrenfest theorem for a fixed potential, the average q and p are expected to be well-approximated by an ensemble of classical trajectories, but the average values of the squares of q and p may not be well-reproduced.

We will now analyze the results of the MTE approach for photons by comparing against the quantum and quasiclassical propagations on the exact qTDPES and from the 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 $∂qS(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. 38 and 46.

We choose parameters for the model of Eq. (1) as ϵ_e = 0.2, ϵ_g = −0.2, ω_c = 0.4, and r_egλ = 0.01 a.u. The resulting dynamics gives Rabi-like oscillations with a period of T ≈ 2π/0.01 a.u.

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 the harmonic by less than 0.01% in the range shown; diagonalization yields $EqBO(q)=ωc2q2/2∓ωc2/4+(regλωc)2q2$⁠, which shows that the deviation grows with q². The plot also shows the nonadiabatic coupling d_eg(q), indicating that the MTE force in Eq. (13) is dominated by the first term and also that, unlike the typical electron-nuclear case, the nonadiabatic 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 the 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(−ωq2/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 toward the lower qBO surface, while at smaller q, the surfaces remain close to the upper qBO surface. This is consistent with the bilinear 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 leads to an even more pronounced and increasingly localized barrier at q = 0, arising primarily from $Ekin(q,t)$⁠,³⁸ which becomes singular at T/2. During this half-cycle, the photon emission character is evident in |χ(q, t)|², which evolves from the harmonic ground state form to the first-excited state (one-photon state), and also in the population of the excited state |C_e(q, t)|², which evolves from initially being 1 at all q to 0 at all q, except for q = 0.

Thus, the correct forces in a quasiclassical description of the photons are 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²⟩ 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 missing in the 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 the 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 the 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 split-operator method with a time step of 0.001 a.u. and a spatial grid spacing of 0.1 a.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 the leapfrog integration with 10 000 trajectories (initially sampled using the Gaussian initial state), and convergence was checked with respect to both time step and number of trajectories. For the MTE calculation, 5000 trajectories were used, and a time step of 0.1 a.u. was enough for the convergence.

The top left panel and the 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 $Ekin$ contribution to ⟨p²⟩; 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 (multimedia view), where in the second row, the green dashed curve shows the photon displacement-coordinate density obtained from the 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 self-consistent 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, the quasiclassical propagation on both the full qTDPES and wBO are reasonably accurate, especially for the photon number where there is some error cancellations between a small overestimate of ⟨q²(t)⟩ and underestimate of ⟨p²(t)⟩. The MTE underestimates the photon number ⟨N(t)⟩ and the intensities ⟨q²(t)⟩ and ⟨p²(t)⟩. The time snapshots in Fig. 3 (Multimedia View) enable a deeper analysis and movies are also provided in the supplementary material.

Let us first discuss the MTE case. The top figures in each of the four panels in Fig. 3 show 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 nonadiabatic coupling causes the density to spread, somewhat similar to a harmonic oscillator that widens (and narrows in the second half of the cycle). However, without any barrier term, the MTE is unable to even remotely capture the correct structure of the density. The increase in the averaged quantities such as ⟨q²⟩ and N as a photon is emitted to 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 (multimedia view) in the third row, we notice that the density splits too quickly in the early dynamics (t = 100 and 160 a.u. 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 slightly underestimates the outer tail regions of the density 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 over-localize 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 perform 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 the propagation on the wBO than for the propagation on the full qTDPES, as evident from the second row in Fig. 3 (multimedia view). At larger times, there is an underestimate in ⟨q²⟩ that is compensated by an overestimate of ⟨p²⟩, which comes from the $Ekin$ peak in Eq. (24) that is sampled more because the photonic density does not split as much.

Although this 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 the 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.^24,31,32

In contrast, a quasiclassical evolution on the weighted qBO surfaces, wBO, or full qTDPES performs much better for the 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 accurate results at longer times. The sharp barrier forming at q = 0 arises in the gauge-independent term $Ekin$⁠, while in this case, the gauge-dependent term is relatively small.³⁸ Whether there exists a gauge in which the gauge-dependent term has an opposing sharp structure that mitigates the one in $Ekin$ remains to be explored, possibly quasiclassical formulations in such a gauge would be more accurate.

Still, the computational efficiency and reasonable accuracy of MTE treatment of photonic dynamics in cavity-QED are encouraging,^24,31–33 and these 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 the 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 the 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.^39,47–51

Two movies are given in the supplementary material showing the dynamics of the photonic displacement-coordinate density; one movie shows this computed via the MTE against the exact quantum solution, and the other shows this computed via the quasiclassical evolution on the qTDPES and on the wBO component also compared with the exact. The movies correspond to the results shown in Fig. 3.

Financial support from the U.S. National Science Foundation under Grant No. CHE-1940333 (N.T.M. and B.R.) and the Department of Energy, Office of Basic Energy Sciences, Division of Chemical Sciences, Geosciences and Biosciences, under Award No. DE-SC0020044 (L.L.) is gratefully acknowledged. The work of N.M.H. was supported by the German National Academy of Sciences Leopoldina.

The authors have no conflicts to disclose.

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