{'fraction_non_alphanumeric': 0.04750878475819065, 'fraction_numerical': 0.018077305872077672, 'mean_word_length': 4.312830473723967, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Mixed Quantum-Classical Electrodynamics: Understanding Spontaneous Decay and Zero Point Energy

The dynamics of an electronic two-level system coupled to an electromagnetic field are simulated explicitly for one and three dimensional systems through semiclassical propagation of the Maxwell-Liouville equations. We consider three flavors of mixed quantum-classical dynamics: the classical path approximation (CPA), Ehrenfest dynamics, and symmetrical quantum-classical (SQC) dynamics. The CPA fails to recover a consistent description of spontaneous emission. A consistent"spontaneous"emission can be obtained from Ehrenfest dynamics--provided that one starts in an electronic superposition state. Spontaneous emission is always obtained using SQC dynamics. Using the SQC and Ehrenfest frameworks, we further calculate the dynamics following an incoming pulse, but here we find very different responses: SQC and Ehrenfest dynamics deviate sometimes strongly in the calculated rate of decay of the transient excited state. Nevertheless, our work confirms the earlier observations by W. Miller [J. Chem. Phys. 69, 2188-2195, 1978] that Ehrenfest dynamics can effectively describe some aspects of spontaneous emission and highlights new possibilities for studying light-matter interactions with semiclassical mechanics.


I. INTRODUCTION
Understanding the dynamics of light-matter interactions is essential for just about any flavor of physical chemistry; after all, with a few exceptions, photons are the most common means nowadays to interrogate molecules and materials in the laboratory. Today, it is standard to study molecules and materials with light scattering experiments, absorption spectroscopy, pump-probe spectroscopy, etc. For a chemist, the focus is usually on the matter side, rather than the electromagnetic (EM) field side: one usually pictures an incoming EM field as a time-dependent perturbation for the molecule. Thereafter, one calculates how the molecule responds to the perturbation and, using physical arguments and/or semiclassical insight, one extrapolates how the molecular process will affect the EM field. For instance, in an absorption experiment, we usually assume linear response theory [1] when calculating how much energy the molecule absorbs. More precisely, one calculates a dipole-dipole correlation function and then, after Fourier transform, one can make an excellent prediction for the absorption pattern. For weak electric fields, this approach often results in reliable data.
However, in many situations involving strong light/matter interactions (e.g. laser physics), the states of the radiation field and the material sub-systems have to be considered on equal footing. An example of strong recent interest is the host of observed phenomena that manifest strong exciton-photon coupling. [2][3][4] Closely related, and also in recent focus, are observations and models pertaining to strong interactions between molecules and electromagnetic modes confined in optical cavities. [5][6][7] As another example, recent studies by Mukamel [8], Bucksbaum [9] and coworkers who have explored the proper interpretation of x-ray pump-probe scattering experiments and, in particular, the entanglement between electrons, nuclei and photons. Beyond the analysis of simplified quantum models, the important tools in analyzing many of these phenomena are variants of coupled Maxwell and Schrödinger (or, when needed, quantum-Liouville) equations, where the radiation field is described by classical Maxwell equations while the molecular system is modeled with a handful of states and described quantum mechanically. [3,[10][11][12][13][14][15][16][17]. A classical description of the radiation field is obviously an important element of simplification in this approach, which makes it possible to simulate the optical response of realistic model systems. However, open questions remain in this area, in particular: • How does spontaneous emission emerge, if at all, in semiclassical calculations?
• How do we best describe computationally the possibly simultaneous occurrence of absorption, scattering, fluorescence and non-linear optical response following a pulse or CW excitation of a given molecular system that may interact with its environment?
• How do we treat both quantum-mechanical electron-electron interactions (e.g. spinorbit coupling) and classical electronic processes (e.g. electronic energy transfer) in a consistent fashion?
In the future, our intention is to address each and every one of these questions. For the present article, however, our goal is to address the first question. We note that spontaneous emission rates can be evaluated from the rate of energy emission by a classical dipolar antenna [18]. An important quantification of this observation has been provided by Miller [19] who has shown that apart from semiclassical corrections, spontaneous decay rates can be ascertained from classical dynamics. Indeed, for a dipolar harmonic oscillator Miller has shown that a semiclassical decay rate can be obtained from classical dynamics exactly. His treatment [19], however, raises several questions. First, in Ref. 19, the molecular system is represented by a classical harmonic oscillator rather than a 2-level system. How will the observations made by Miller be affected with a proper quantum-mechanical treatment?
What will be the performance of mixed semiclassical treatments for spontaneous emission, and which semiclassical treatment will perform best? Second, in Ref. 19, no explicit light pulses are applied to the electronic system, but one can ask: If a pulse of light is applied to the system, and we use mixed quantum-classical dynamics, is the propagated photon field consistent with the ensuing molecular dynamics? With an external temperature, do we recover detailed balance? In this article, we will address most of these questions, paying special attention to the recent symmetrical quantum-classical (SQC) dynamics protocol of Cotton and Miller [20,21]. This article is arranged as follows. In Sec. II, we briefly review the theory of spontaneous decay. In Sec. III, we introduce the semiclassical Hamiltonian in our model. In Sec. IV, we implement Ehrenfest dynamics, CPA and SQC. In Sec. V, simulation details are given. In Sec. VI, we compare results for spontaneous decay. In Sec. VII, we simulate and analyze two cases: (i) the arrival of an incoming pulse and (ii) dephasing effects. We conclude in Sec. VIII.
For notation, we use the following conventions: ω 0 is used to represent the energy difference between the excited state |e and the ground state |g ; ω k (or ck ) is used to represent the energy of the photon with wave vector k ; µ 12 is the electric transition dipole moment of the molecule; 1 a represents the molecular size so that the transition dipole moment with a characteristic charge q is approximately µ 12 ≈ 1 a q; ς is used to represent a dephasing rate; U 0 denotes the total energy of an incident pulse; k 0 denotes the peak position, in Fourier space, of an incident pulse; b is a parameter fixing the width of an incident pulse in space; and c is the speed of light. We work below in SI units.

II. THEORY OF SPONTANEOUS EMISSION
For completeness, and because we will work in both one and three dimensions, it will be convenient to briefly review the theory of spontaneous emission and dipole radiation.
Consider a molecular species in an excited state |e which can decay to the ground state |g by emitting a photon spontaneously.

A. The Fermi's Golden Rule (FGR) Rate
Let the vacuum state for the radiation field be |0 . Suppose that initially the system is in state |e ⊗ |0 . At long times, we expect to observe spontaneous emission, so that the final state will be |g ⊗ a † q,s |0 . Here, a † q,s creates a photon with wave vector q and polarization s. We now apply Fermi's Golden Rule (FGR) for the emission rate. We further make the dipole approximation, so that the interaction Hamiltonian for a molecule sitting at the origin is H int = −qˆ r ·ˆ E(0), where q is the electronic charge,ˆ r is the position operator for the quantum system, andˆ E(0) is the electric field at the origin. In such a case, the decay rate k in 3D can be calculated as follows [22]: Here, µ 12 = | e|qˆ r|g | is the three-dimensional transition dipole moment of the molecule, k , s is the a unit vector in direction of the electric field indexed by the wave vector k and the polarization vector s, and ω 0 is the energy difference between |e and |g . Eqn. (1a) is the usual FGR expression. In Eqn. (1b), if we replace the discrete k with the continuous dϕdθdk sin θk 2 ρ(k ), where ρ(k ) = V /(2π) 3 is the three-dimensional density of states (DOS) for the photons, we recover Eqn. (1c).
In what follows below, it is useful to study EM radiation in 1D as well as in 3D. To that end, we will imagine charge distributions that are function of x only, i.e. they are uniform in y and z directions. In 1D, the density of states (DOS) for the photon field is ρ( k ) = L x /2π.
Therefore, the decay rate in 1D is: Using V = L x L y L z and defining the one-dimensional dipole moment |µ 1D 12 | 2 = |µ 3D 12 | 2 /L y L z , we can rewrite the final 1D rate as Below, we will use µ 12 to represent either µ 1D 12 and µ 3D 12 depending on context. Note that, in 1D, the spontaneous decay rate k 1D depends linearly on the frequency ω 0 and quadratically on the transition dipole moment µ 12 . In 3D, however, k 3D depends cubically on ω 0 instead of linearly, but still quadratically on µ 12 . Note that, for Eqns. (1c) and (3) to apply, two conditions are required: (i) The dipole approximation must be valid, i.e. the wavelength of the spontaneous light must be much larger than width of molecule.
(ii) The coupling between molecule and radiation field must be weak to ignore any feedback of the EM field, i.e. ω 0 must be much larger than the inverse lifetime.

B. The Abraham-Lorentz Rate
While FGR is the standard protocol for modeling spontaneous emission with quantum mechanics, we can also recover a similar decay rate with classical mechanics by using the Abraham-Lorentz equation [23] . For a classical charged harmonic oscillator moving in the x direction with mass m, the Abraham-Lorentz equation reads where τ = q 2 /6π 0 c 3 m has the dimension of time. The last term in Eqn. (4) represents the recoil force on a particle as it feels its own self-emitted EM field. Since τ 1/ω 0 , we can assume the damping effect is small and so we replace mτ ...
Eqn. (5) represents a damped harmonic oscillator, which has a well-know solution since ω 0 τ 1. In Eqn. (6), the amplitude x 0 and the phase φ will depend on the initial conditions, and the decay rate k AL is At this point, we can write down the total energy of the harmonic oscillator: To relate the Abraham-Lorentz rate k AL to the FGR rate in 3D, we require a means to connect a classical system with mass m to a pair of quantum mechanical states. To do so, we imagine the oscillator is quantized and that the motion is occurring in the ground state, where x 2 0 = /2mω 0 . This is equivalent to asserting that the initial energy of the dipole is 1 2 ω 0 , which we set equal to the total dipole energy, mω 2 0 x 2 0 . If we further assert that the dipole operator is off-diagonal (as in Eqn. (16)), we may substitute qx 0 ≈ µ 12 , which leads to the following Abraham-Lorentz rate (k AL ) With this ansatz, the Abraham-Lorentz decay rate k AL is equal to the FGR rate in 3D. Note that several ad hoc semiclassical assignments must be made for this comparison, and it is not clear how to generalize the Abraham-Lorentz approach to treat more than two electronic states in a consistent fashion.

C. The Asymptotic Electromagnetic Field
Below, we will analyze different schemes for solving Maxwell's equations coupled together with the Liouville equation, and it will be helpful to compare our results with the standard theory of dipole radiation. According to classical electrodynamics, if a dipole is located at the origin and is driven by an oscillating field, the electromagnetic (EM) field is generated with the energy density (at time t and position r) given in the far-field by [24] u( r, t) = µ 0 c 2 Here, without loss of generality, we assume that the dipole is pointing in the z−direction, so that θ is the polar angle from the z-axis. r is the distance from the observer to the dipole (sitting at the origin). Eqn. (10) predicts that, for the energy density, there is sin 2 θ dependence on the polar angle θ and 1/r 2 dependence on the distance r. Note that Eqn.
(10) is valid in the far-field when r λ d, where λ is the wavelength of EM field and d is the size of the dipole.

III. THE SEMI-CLASSICAL HAMILTONIAN
We consider the problem of a two-level system coupled to a radiation field. After a Power-Zienau-Woolley transformation [25,26] is applied, the Hamiltonian reads as follows: Here,B = ∇ ×Â,D ⊥ = 0Ê +P ⊥ .Â is the vector potential for the EM field andP ⊥ is the polarization operator for the matter. For the EM field, the relevant commutators Hamiltonian of the electronic system, which will be defined below. We ignore all magnetic moments in Eqn. (11).
Eqn. (11) is a large Hamiltonian, written in the context of a quantum field. For semiclassical dynamics, it is convenient to extract the so-called "electronic Hamiltonian" that depends only parametrically on the EM field. Following Mukamel [25], one route to achieve such a semiclassical Hamiltonian is to consider the equation of motion for an observable of the matterQ: If we approximate that the E-field is classical, so that we may commuteÊ ⊥ with all matter operators, we find the following semiclassical electronic Hamiltonian: With only one charge center, however, we will not need to distinguish between the longitudinal and perpendicular components, and so we will drop the ⊥ notation below.
For this paper, we consider the simplest case of two electronic states: the ground state |g and the electronic excited state |e . Thus, we represent H s as follows: Furthermore, we assume that (a) these states carry no permanent dipole and (b) the transition between them is characterized by two single electron orbitals ψ g and ψ e and an effective charge q such that the transition dipole density is given by with a corresponding polarization operator: For example, in 3D, in the common case that ψ e ( r) is a p z orbital ( 2 1/2 a 5/4 π 3/4 ze −ar 2 /2 ) and ψ g ( r) is an s orbital (( a π ) 3/4 e −ar 2 /2 ), ξ( r) would be ξ 3D ( r) = q r 2 1/2 a 2 π 3/2 ze −ar 2 If we consider a charge distribution that is effectively 1D, changing along in the x direction but polarized in the z direction, the reduced form of ξ( r) would be The magnitude of ξ( r) is related to the magnitude of the total transition dipole moment, µ 12 : Eqn. (19) guarantees that, when the width ofP( r) approaches 0, Eqn. (13) becomes the standard dipole Hamiltonian,Ĥ el =Ĥ s − µ 12 · E(0). This definition allows us to rewrite Eqns. (17)(18) above, as follows: Note that ξ 3D and ξ 1D have different units.
In Appendix A we will show that under the point dipole limit -where the width of ξ( r) is much smaller than the wavelength of EM field, so that ξ( r) can be treated as a delta function -some analytic results can be derived for the coupled electronic-photons dynamics.

IV. METHODS
Many mixed quantum-classical semiclassical dynamics tools have been proposed over the years to address coupled nuclear-electronic dynamics, including wave packet dynamics [27,28], Ehrenfest dynamics [29], surface-hopping dynamics [30,31], multiple spawning dynamics [32], and partially linearized density matrix dynamics (PLDM) [33]. Except for the Ehrenfest (mean-field) dynamics, other methods are usually based on the Born-Oppenheimer approximation, which relies on the timescale separation between (slow) classical and (fast) quantum motions. Such methods cannot be applied in the present context because the molecular timescales and the relevant photon periods are comparable. [34] The Ehrenfest approximation relies on the absence of strong correlations between interacting subsystems, and may be valid under more lenient conditions. We therefore limit the following discussion to the application of the Ehrenfest approximation and its variants [35].

A. Ehrenfest Dynamics
According to Ehrenfest dynamics for a classical radiation field and a quantum molecule, the molecular density operatorρ(t) is propagated according to while the time evolution of the radiation field is given by the Maxwell's equations Here, the current density operator,ˆ J = dP /dt, is replaced by its expectation value: If we substitute Eqns. (16) and (21) into Eqn. (23), the current density J( r) can be simplified to where ρ 12 is the coherence of the density matrix ρ.
Two points are noteworthy: First, because Eqn. (21) does not include any dephasing or decoherence, there is also an equivalent equation of motion for the electronic wavefunction Here H el ij is a matrix element of the operatorĤ el =Ĥ s − d r E( r) ·P( r). Second, under the dynamics governed by Eqns. (21) and (22), the total energy of the system U tot is conserved, where Altogether, Eqns. (21), (22), and (23) capture the correct physics such that, when an electron decays from the excited state |e to the ground state |g , an EM field is generated while the total energy is conserved.

Advantages and disadvantages of Ehrenfest dynamics
The main advantage for Ehrenfest dynamics is a consistent, simple approach for simulating electronic and EM dynamics concurrently.
Several drawbacks, however, are also apparent for Ehrenfest dynamics. First, consider Eqn. (24). Certainly, if the initial electronic state is an eigenstate of H s , i.e. (C 1 , C 2 ) = (0, 1), then ρ 12 (t = 0) = C 1 C * 2 = 0 and there will be no current density J( r) if there is no EM field initially in space. Thus, in disagreement with the exact quantum result, there is no spontaneous emission: the initial state (0, 1) will never decay. According to Ehrenfest dynamics, spontaneous emission can be observed only if C 1 = 0 and C 2 = 0, i.e., if the initial state is a linear combination of the ground and excited states.
Second, it is well known that, for finite temperature, Ehrenfest dynamics predicts incorrect electronic populations at long time: the electronic populations will not satisfy detailed balance [36]. Here, finite temperature would correspond to a thermal distribution of photon modes at time t = 0, representing the black-body radiation. However, for the purposes of fast absorption and/or scattering experiments, where there is no equilibration, this failure may not be fatal.

B. The Classical Path Approximation (CPA)
If Ehrenfest dynamics provides enough accuracy for a given simulation, the relevant dynamics can actually be further simplified and reduced to the standard "classical path approximation (CPA)" [37]. To make this reduction, note that the EM field can be considered the sum of 2 parts: (i) the external EM field E ext ( r) that represents a pulse of light approaching the electronic system and (ii) the scattered EM field E scatt ( r) generated from spontaneous or stimulated emission from the molecule itself. Thus, at any time, where we impose free propagation for the external EM field, i.e., E ext ( r, t) = E ext ( r − ctr ext , 0). Herer ext represents the unit vector in the propagation direction of the external EM field.
According to the CPA, we ignore any feedback from electronic evolution upon the EM field, i.e., we neglect the d r E scatt ( r) ·P( r) term of Eqn. (21). Thus, the electronic dynamics now obey while photon dynamics still obeys Eqn. (22). This so called classical path approximation underlines all usual descriptions of linear spectroscopy, and should be valid when | E scatt | | E ext |. In such a case, the coherence ρ 12 and current density J are almost unchanged if we neglect the d r E scatt ( r) ·P( r) term.

Advantages and disadvantages of the CPA
Obviously, the advantage of Eqn. (27) over Eqn. (21) is that we can write down an analytical form for the light-matter coupling ( d rE( r)P ( r)), since E ext propagates freely.
That being said, the disadvantage of the CPA is that one cannot obtain a consistent description of spontaneous emission for the electronic degrees of freedom, because the total energy is not conserved; see Eqns. (22) and (27). As such, the classical path approximation would appear reasonably only for studying the electronic dynamics; EM dynamics are reliable only for short times.

C. Symmetrical Quasi-classical (SQC) Windowing Method
As discussed above, the Ehrenfest approach cannot predict exponential decay (i.e. spontaneous emission) when the initial electronic state is (0, 1). Now, if we want to model spontaneous emission, the usual approach would be to include the vacuum fluctuations of the electric field, in the spirit of stochastic electrodynamics [38]. That being said, however, there are other flavors of mean-field dynamics which can improve upon Ehrenfest dynamics and fix up some failures. [33,39] (i.e., the inability to achieve branching, the inability to recover detailed balance, etc.) Miller's symmetrical quasi-classical (SQC) windowing [21] is one such approach.
The basic idea of the SQC method is to propagate Ehrenfest-like trajectories with quantum electrons and classical photons (EM field), assuming two modifications: (a) one converts each electronic state to a harmonic oscillator and includes the zero point energy (ZPE) for each electronic degree of freedom (so that one samples many initial electronic configurations and achieves branching); and (b) one bins the initial and final electronic states symmetrically (so as to achieve detailed balance). We note that SQC dynamics is based upon the original Meyer-Miller transformation [40], which was formalized by Stock and Thoss [41], and that there are quite a few similar algorithms that propagate Ehrenfest dynamics with zero-point electronic energy [39]. While Cotton and Miller have usually propagated dynamics either in action-angle variables or Cartesian variables, for our purposes we will propagate the complex amplitude variable C 1 , C 2 so as to make easier contact with Ehrenfest dynamics [42]. For- For completeness, we will now briefly review the nuts and bolts of the SQC method for a two-level system coupled to a bath of bosons.
1. Standard SQC procedure for a two-level system coupled to a EM field 1. At time t = 0, the initial complex amplitudes C 1 (0) and C 2 (0) are generated by Eqn.
Here, RN is a random number distributed uniformly between [0, 1] and n j = 0, 1 is the action variable for electronic state j. n j = 0 implies that state j is unoccupied while n j = 1 implies state j is occupied. θ j = 2πRN is the angle variable for electronic state j. Note that |C 1 | 2 + |C 2 | 2 = 1, but rather, on average |C 1 | 2 + |C 2 | 2 = 1 + 2γ, such that γ is a parameter that reflects the amount of zero point energy (ZPE) included. Originally, γ was derived to be 1/2 [40], but Stock et al. [43] and Cotton and Miller [21] have found empirically that 0 < γ < 1/2 often gives better results.
4. At each time t, one may calculate raw populations (before normalization) as follows: Here, N is the number of trajectories and W 1 is the window function for the ground state |g , centered at (n 1 , n 2 ) = (1, 0); W 2 is the window function for the excited state |e , centered at (n 1 , n 2 ) = (0, 1). (l) means the lth trajectory.
5. The true density matrix at time t is calculated by normalizing Eqn. (30) in the following manner: Miller and Cotton have also proposed a protocol to calculate coherences and not just populations [44], but we have so far been unable to extract meaningful values from this approach. Future work exploring such coherences would be very interesting.

Choice of window function and initial distribution
Below, we will study a two-level system weakly coupled to the EM field, i.e. the polarization energy will be several orders less than ω 0 . For such a case, one must be very careful about binning. Cotton and Miller [45] have suggested that triangular window functions with γ = 1/3 perform better than square window functions in this regime. Therefore, we have invoked the triangular window function in Eqn. (32) with γ = 1/3 below.
Here, h(x) is Heaviside function. Fig. 1 gives a visual representation of the triangular window function in Eqn. (32). The bottom and upper pink triangles represent areas where W 1 = 0 and W 2 = 0 respectively. To be consistent with the choice of triangular window functions, one must modify the standard protocol in Eqn. (28). Instead of the standard square protocol, assuming we start in excited state |e , one generates a distribution of initial action variables (n 1 (0), n 2 (0)) within the area where W 2 = 0 (see Eqn. 32) uniformly. Visually, this initialization implies a distribution of (n 1 (0), n 2 (0)) inside a triangle centered at (0, 1) in the (n 1 , n 2 ) configuration space, as demonstrated in Fig. 1. The protocol for initializing angle variables is not altered: one sets θ j = 2πRN, j = 1, 2.

Advantages and disadvantages of SQC dynamics
Compared with Ehrenfest dynamics, one obvious advantage of SQC dynamics is that the latter can model spontaneous emission when the initial electronic state is (0, 1). Moreover, the SQC approach must recover detailed balance in the presence of a photonic bath at a given temperature [46] -provided that the parameter γ is chosen to be small enough for the binning [42].
At the same time, the disadvantage of the SQC method is that all results are sensitive to the binning width γ. γ should be big enough to give enough branching, but also should be small enough to enforce detailed balance [42]. As a result, one must be careful when choosing γ. Although not relevant here, it is also true that SQC can be unstable for anharmonic potentials. [42] Lastly, as a practical matter, we have found SQC requires about 1000 times more trajectories than Ehrenfest dynamics.

D. Classical Dynamics with Abraham-Lorentz Forces
Although (as shown above) classical electrodynamics with Abraham-Lorentz forces can be useful to model self-interaction, we will not analyze Abraham-Lorentz dynamics further in this paper. Because the correspondence between Ehrenfest dynamics and Abraham-Lorentz dynamics is not unique or generalizable, we feel any further explanation of Abraham-Lorentz equation would be premature. While a Meyer-Miller transformation [40] can reduce a quantum mechanical Hamiltonian into a classical Hamiltonian, the inverse is not possible.
Thus, it is not clear how to run classical dynamics with Abraham-Lorentz forces starting from an arbitrary initial superposition state (C 1 , C 2 ) in the {|g , |e } basis. For instance, following the approach above in Section II B, we might set mω 2 0 x 2 = |C 2 (0)| 2 ω 0 /2. However, doing so leads to a rate of decay equal to k FGR /|C 2 (0)| 2 . This result goes to infinity in the limit C 2 → 0; see Fig. 11. Future work may succeed at finding the best correspondence between semiclassical dynamics and the Abraham-Lorentz framework, but such questions will not be the focus of the present paper.  f Eqns. (34)(35) located far away at time zero).

B. Propagation procedure
Equations of motion (Eqns. (21), (22)) are propagated with a Runge-Kutta 4th order solver, and all spatial gradients are evaluated on a real space grid with a two-stencil in 1D and a six-stencil in 3D. Thus, for example, if we consider Eqn. (22) in 1D, in practice we approximate: etc. Here (i) is a grid index. This numerical method to propagate the EM field (Eqn. (22)) is effectively a finite-difference time-domain (FDTD) method [47,48].

C. Absorbing boundary condition (ABC)
To run calculations in 3D, absorbing boundary condition (ABC) are required to alleviate the large computational cost. For such a purpose, we invoke a standard, one-dimensional smoothing function [49,50] S(x): In 1D, by multiplying the E and B field with S(x) after each time step, we force the E and B fields to vanish for |x| > R 1 .
In 3D, we choose the corresponding smoothing function to be of the form of Eqn. (35), where S(x), S(y) or S(z) is exactly the same as Eqn. (34). Note that this smoothing function has cubic (rather than spherical) symmetry. for the 3D dynamics in Fig. 7.

D. Extracting Rates
Our focus below will be on calculating rates of emission; these rates will be subsequently compared with FGR rates. To extract a numerical rate (k) from Ehrenfest or SQC dynamics, we simply calculate the probability to be on the excited state as a function of time (P 2 (t)) and fit that probability to an exponential decay: P 2 (t) ≡ P 2 (0)e −kt . For Ehrenfest dynamics, all results are converged using the default parameters in Table I. For SQC dynamics, longer simulation times are needed (to ensure P 2 (t end ) < 0.02); in practice, we set t end = 150 fs.
Note that, for SQC dynamics, P 2 (t) in SQC is calculated by Eqn. (31b) and we sample 2000 trajectories.

VI. RESULTS
We now present the results of our simulations and analyze how Ehrenfest and SQC dynamics treat spontaneous emission. The initial state is chosen to be (C 1 , C 2 ) = ( 1/2, 1/2) for Ehrenfest dynamics. We begin in one-dimension.
A. Ehrenfest Dynamics: 1D In Fig. 2, we plot P 2 (t) for the default parameters in Table I. Clearly, including ABC's has no effect on our results. For this set of parameters, Ehrenfest dynamics predicts a decay rate that is ∼ 1/3 slower than Fermi's Golden Rule (FGR) in Eqn. (3).
In Fig. 3, we now examine the behavior of Ehrenfest dynamics across a broader parameter regime. In Fig. 3a and 3b, we plot the dependence of the decay rate on the energy difference of electronic states, ω 0 , and the dipole moment, µ 12 . Ehrenfest dynamics correctly predicts linear and quadratic dependence, respectively, in agreement with FGR in 1D (see Eqn. (3)).
Generally, the fitted decay rate from Ehrenfest dynamics is ∼ 1/3 slower than FGR. As far as the size of the molecule is concerned, in Fig. 3c, we plot the decay rate k as a function of the parameter a (in Eqn. 20b). Note that our results are independent of molecular size when a > 0.05 nm −2 . This independence underlies the dipole approximation: when the width of the molecule is much smaller than wavelength of light, 1/a c/ω 0 , the decay rate should not be dependent on the width of molecule. Note that ω 0 = 16.46 eV for these simulations, which dictates that results will be dependent on a for a < 0.05 nm −2 . Finally,

Initial Conditions
The results above were gathered by setting C 1 = 1/2. Let us now address how the initial conditions affect the Ehrenfest rate of spontaneous decay. In Fig. 4 we plot k vs.
Here, we differentiate how k is extracted, either from a (a) a fit of the long time decay (t end = 99 fs) or (b) a fit of the short time decay (t end = 5 fs). Clearly, the decay rates in Fig. 4a and 4b are different, suggesting that the decay of P 2 is not purely exponential (see detailed discussion in Appendix); the decay constant is itself a function of time. Moreover, according to Fig. 4b, the short time decay rate appears to be linearly dependent on |C 1 (0)| 2 and, in the limit that |C 1 (0)| 2 → 1, both fitted decay rates k approach the FGR result.  Table I. The initial electronic state is (C 1 , C 2 ) = ( 1/2, 1/2). Note that Ehrenfest dynamics captures most of the correct FGR physics.
These results suggest that the fitted decay rate k satisfies where k FGR is the FGR decay rate. In fact, in the Appendix, we will show that Eqn. (36) can be derived for early time scales (2π/ω 0 t 1/k FGR ) under certain approximations.
We also mention that the same failure was observed previously by Tully when investigating the erroneous long time populations predicted by Ehrenfest dynamics. [46,51,52]

Distribution of EM field
Beyond the electronic subsystem, Ehrenfest dynamics allows us to follow the behavior of the EM field directly. In Fig. 5, we plot the distribution of the EM field at times 3.00 fs (a-b) , 30.00 fs (c-d ), and 99.00 fs (e-f ) with two methods: Ehrenfest (red lines) and the CPA (light blue lines). On the left hand side, we plot the electric field in real space (E z (x)); on the right hand side, we plot the EM field in Fourier space (E z (k x )). Here, the Fourier transform is performed over the region x > 0, which corresponds to light traveling exclusively to the right. In the insets on the right, we zoom in on the spectra in a small neighborhood of ω 0 (here, 16.46 eV).
From Fig. 5, we find that Ehrenfest dynamics and the CPA agree for short times. However, for larger times, only Ehrenfest dynamics predicts a decrease in the EM field (corresponding to the spontaneous decay of the signal). This decrease is guaranteed by Ehrenfest dynamics because this method conserves energy. By contrast, because it ignores feedback and violates energy conservation, the CPA does not predict a decrease in the emitted EM field as a function of time (or any spontaneous decay). Thus, overall, as shown in Fig. 5f, the long time EM signal will be a Lorentzian according to Ehrenfest dynamics or a deltafunction according to the CPA. These conclusions are unchanged for all values of the initial  Table I have been used here. Note that Ehrenfest dynamics and the CPA agree for short times but only Ehrenfest dynamics predicts a decrease in the EM field for larger times, which is a requirement of energy conservation.

B. SQC: 1D
The simulations above have been repeated with SQC dynamics. In Fig. 6a, we plot P 2 (t) for a single trajectory that begins on the excited state (C 2 = 1) for the default parameters (see Table I). The remaining three sub-figures in Fig. 6   All unreported parameters are set to their default values in Table I. The initial electronic state is (C 1 , C 2 ) = (0, 1). Note that the SQC decay rates are very close to the FGR rates (less than 10 % difference), whereas Ehrenfest dynamics completely fail and predicts k = 0 for this case (when C 2 = 1 initially). For these simulations, we apply ABC's. C 2 = 1, the overall SQC decay rate k is almost the same as FGR (less than 10 % difference), whereas Ehrenfest dynamics completely fails and predicts k = 0. [53] C. Ehrenfest Dynamics: 3D Finally, all of the Ehrenfest simulations above have been repeated in 3D. Overall, as shown in Fig. 7, the results are qualitatively the same as in 1D. However, as was emphasized in Sec. II, the decay rate now depends cubically (and not linearly) on ω 0 .
Concerning the radiation of EM field in 3D, in Fig. 8 All unreported parameters are set to their default as in Table I. The initial electronic state is (C 1 , C 2 ) = ( 1/2, 1/2). The Ehrenfest decay rates in 3D depend correctly only a, ω 0 and µ 12 and match FGR. For these simulations, we apply ABC's.
θ through sin 2 θ. These results are in very good agreement with theoretical dipole radiation (black line, Eqn. 10). Lastly, in Fig. 9, we plot the energy density as a function of the radial distance r from the molecule, while keeping the polar angle fixed at θ = π/2 (a) and θ = π/4 (b   Table I. ABCs are not applied here. The initial electronic state is (C 1 , C 2 ) = ( 1/2, 1/2). Note the strong and perhaps surprising agreement between Ehrenfest/CPA dynamics and the classical dipole radiation; this agreement depends on the choice of initial electronic states, as is proven in the Appendix.

VII. DISCUSSION
The results above suggest that, for their respective domains of applicability, both Ehrenfest dynamics and SQC can recover spontaneous emission. We will now test this assertion by investigating the response to (i) photo-induced dynamics and (ii) dephasing. are the same as in Fig. 8. The radial distribution of EM energy density is the same for Ehrenfest and the CPA at short times and, just as in Fig. 8, these radial distributions agree with the classical dipole radiation result (provided the initial electronic state is (C 1 , C 2 ) = ( 1/2, 1/2)).

A. An incoming pulse in one dimension
To address photo-induced dynamics, we imagine there is an incident pulse at t = 0 of the form: Here, A(b, k 0 , x 0 ) is an normalization coefficient with value The total energy of incident pulse is U 0 . The parameter b determines the width of the pulse in real space. k 0 defines the peak of the pulse in reciprocal space. x 0 represents the center of pulse at t = 0.
At time zero, the Fourier transform of E z (x) is: resembles a single large packet at k x = 0. For resonance with the molecule, |E z (k x )| should be large at ck x = ω 0 (16.46 eV by default).

Electronic dynamics
In Fig. 10, we plot the electronic population of the excited state as a function of time after exposure to incident pulses of different intensity (U 0 ) and wavevector (k 0 ); see Eqn. (37).
We plot short and long times, on the left and right hand sides, respectively. For strong, resonant pulses, (U 0 = 19.7 keV, k 0 = 0.013 nm −1 ), there is obviously a strong response (see a-b). For strong, off-resonant pulses (U 0 = 19.7 keV, k 0 = 0.334 nm −1 ), obviously the response is weaker. In both situations, SQC (green line) and Ehrenfest dynamics (red line) agree almost exactly for short times. At longer times, however, the SQC P 2 (t) value decays ∼ 2 times faster than the Ehrenfest dynamics result.
Let us consider now weak pulses. In Fig. 10e-h, we plot the excited state population when the incident pulse is weak (U 0 = 3.29 keV), keeping all other parameters unchanged. Now, there is much less agreement between SQC and Ehrenfest dynamics, especially for long times. Generally, SQC predicts a faster decay rate for P 2 (t) than Ehrenfest dynamics for small |E z (ω 0 /c)|.
The statement above is quantified in Fig. 11. Here, we vary U 0 , which results in a change in the initial absorption (which is quantified by 1 − P 2 (t = 0.5 fs) on the x-axis). This graph quantifies how the population decay on the excited state depends on the initial condition: the decay of P 2 decreases when the initial excited state population decreases. Obviously, this Ehrenfest data is in complete agreement with Fig. 4. Now, the new piece of data in Fig. 11 is the SQC data. Here, we see that SQC behaves in  Table I. The initial electronic state (C 1 , C 2 ) = (1, 0). Two methods are compared: Ehrenfest dynamics (red line) and SQC (green line). Note that SQC and Ehrenfest dynamics disagree for long times, especially for weak pulses. See Fig. 11. For these simulations, we apply ABC's. Numerical results for Ehrenfest dynamics show that enforcing ABC's does not make any difference at all.    Fig. 10. Note that Ehrenfest and and CPA dynamics agree at short times but disagree at long times when energy conservation becomes important.

Distribution of the EM field
At this point, we should also comment on the EM field that is produced following incident radiation for the two-level system. Effectively, our results are consistent with Fig. 5 above.
In Fig. 12, on the left, we plot E z (x) versus x in space at times 14.00 fs (a), 55.99 fs (c) and 149.00 fs (e). On the right hand side, we plot the Fourier transform E z (k x ) versus photon energy ck x . As above, we find that, for short times, Ehrenfest dynamics (red lines) and the CPA (light blue lines) are in good agreement. Thereafter, however, the agreement Note that, for Fig. 13, we are operating in the linear response regime: the incoming pulse energy U 0 is relatively weak. In Appendix C, we plot the absorption spectra for a few different incoming fields and demonstrate that the results are linear with U 0 . We also show that standard linear response theory yields a good estimate of the overall lineshape.

B. Dephasing effects
In the present article, we have now shown that semiclassical theories -Ehrenfest and SQC -can both recover some elements of spontaneous emission, which is mostly thought to be a quantum effect [14,19]. With this claim in mind, however, there is now one final subject that must be addressed, namely the role of dephasing. After all, in a large simulation with an environment, dephasing can and will occur; therefore one must wonder whether or not such dephasing will affect the rate of spontaneous emission.
To answer this question, we have run several simple calculations that replace Eqn. (21) by Eqn. (39), Thus, we have propagated electron-photon dynamics by altering the electronic equation of motion but keeping the classical EM equations the same. ς in Eqn. (39) is an empirical dephasing rate: when ς = 0, there is no dephasing and when ς > 0 there is a finite rate of coherence loss between the two electronic states.
In Fig. 14a, we plot the rate of spontaneous emission k as a function of the dephasing rate ς. When dephasing increases, the coherence between the electronic states is expected to decrease, and so the current should decrease, and thus the rate of spontaneous emission is expected to decrease as well. However, perhaps surprisingly, the fitted rate for establishing equilibrium also increases.
Most importantly, in Fig. 14b, we plot the final population of the excited state P 2 (t end ).
As should be expected, the long term excited state population increases (does not reduce to zero) when dephasing increases with either SQC or Ehrenfest dynamics. This graph highlights the limitations of semiclassical methods: as currently implemented, one cannot include both spontaneous emission and dephasing.

VIII. CONCLUSION
In this article, we have simulated the semiclassical dynamics of light coupled to a two-level electronic system with three different methods: Ehrenfest, the CPA and SQC. Most results have been reported in one dimensional, but we have also considered Ehrenfest dynamics in 3D with absorbing boundary conditions. As far as spontaneous emission is concerned, the CPA cannot consistently recover the effect and violates the energy conservation. That being said, Ehrenfest dynamics do predict spontaneous decay consistently, but only provided that we start in a non-trivial superposition state (with C 1 , C 2 = 0). Using electronic ZPE, SQC dynamics predicts spontaneous decay even with C 1 = 0. Both latter methods yield results fairly close to the correct FGR rate. In all cases, unfortunately, spontaneous emission 14. An analysis of the effects of dephasing on spontaneous emission for Ehrenfest and SQC methods in 1D. A plot of (a)the fitted decay rate k as a function of the dephasing rate ς; (b) the normalized long time population of P 2 , P 2 (t end )/P 2 (t = 0) as a function of ς. All simulation parameters are set to their default values in Table I; t end = 400 fs. The initial electric population for the excited state is set to P 2 (0) = 1/2 for Ehrenfest dynamics and P 2 (0) = 1 for SQC dynamics.
Note that both methods fail to recover spontaneous emission in the presence of strong dephasing.
For these simulations, we apply ABC's. Numerical results for Ehrenfest dynamics show that enforcing ABC's does not make any difference at all.
is destroyed when dephasing is introduced, which represents a fundamental limitation of semiclassical dynamics.
Perhaps most interestingly, we have also studied photo-initiated excited dynamics and, in this case, we find very different dynamics as predicted by the different semiclassical methods.
First, as far the EM field is concerned, we have demonstrated that Ehrenfest dynamics can recover the correct absorption spectra, at least qualitatively; at the same time, however, CPA dynamics gives qualitatively incorrect spectra because the method ignores feedback and does not conserve energy. Second, and equally interesting, Ehrenfest dynamics predicts that the overall stimulated decay rate will depend smoothly on initial state (C 1 , C 2 ) but will approach the FGR rate in the weak resonance regime. Vice versa, SQC recovers FGR when (C 1 , C 2 ) = (0, 1) but overestimates the stimulated decay rate, sometimes by as much as a factor of 10 in the weak coupling limit. These SQC anomalies should be very important for designing improved binning protocols in the future [45]. At present, because the cost of SQC dynamics is roughly 1000 times greater than Ehrenfest dynamics and because the method appears to fail for low intensity applied fields, further modification will likely be required before the method can be practical for large-scale simulations.
Looking forward, many questions remain. (i) There are many other semiclassical methods for studying coupled nuclear electronic dynamics [21,[54][55][56][57]; will these methods give us new insight into electrodynamics? (ii) Might we learn more about spontaneous emission by considering ZPE effects through RPMD-like algorithms [58]? (iii) Will different semiclassical methods behave similarly or differently with more than two electronic states? (iv) Can we converge multiple-spawning [59][60][61][62][63] and/or MC-TDH [64][65][66] calculations and generate exact quantum electrodynamical trajectories so that, in the future, we may benchmark other, less exact, semi-classical approximations? And lastly, (iv), are there other, new and non-intuitive features that will emerge when we study multiple pulses incoming upon a molecule? These questions will be answered in the future.
using the well known solution for a wave equation with a source: Here,ρ 12 is the time derivative of ρ 12 . If we average over many different initial electronic populations with different phases, Imρ 12 (0) = 0, the average coupling is simpler: Here, we have denoted ξ(x) = | ξ(x)|. Now, for simplicity, suppose the width of the molecule is infinitely small (i.e., a point-dipole approximation), ξ(x) ≈ µ 12 δ(x). In such a case, Eqn.
(41) can be simplified as: and therefore, from Eqn. (21), At this point, we make the weak coupling approximation, and assume that the off-diagonal terms in H el are infinitely small, so that ρ 12 (t) ≈ √ P 1 P 2 e −iω 0 t is a meaningful first order approximation. Eqn. (43) then reads: where k FGR = ω 0 c 0 |µ 12 | 2 is the FGR spontaneous decay rate in 1D (see Eqn. 3). From Eqn. (44), we can derive the instantaneous transfer rate plus an analytical solution for all times as follows.
Also, τ is much larger than the phase oscillating period (2π/ω 0 τ ), therefore sin 2 (ω 0 t ) can be viewed as a rapid oscillation and we approximate the integral by Then we have ln P 2 (t + τ ) P 2 (t) ≈ −k FGR P 1 (t)τ As a result, we write Ehrenfest dynamics for P 2 in the form of an exponential decay where the instantaneous decay rate is time-dependent On the one hand, for short times, the decay rate is proportional to the initial population k FGR P 1 (0) as shown in Fig. 4b, and we can conclude that Ehrenfest dynamics recovers the FGR rate when P 1 (0) → 1.
On the other hand, we may recast Eqn. (44) in terms of the population difference, ∆P = P 2 − P 1 , Just as above, the instantaneous behavior within the time scale τ can be obtained by Hence, we find an analytical form for P 2 according to Ehrenfest dynamics: For short times, we take t → 0 and find that the instantaneous decay rate is also proportional to P 1 (0). For the initial population P 1 (0) = P 2 (0) = 1/2, as was considered in Fig. 2, the analytical solution becomes P 2 (t) = e −k FGR t /(1 + e −k FGR t ). This formula agrees with the numerical result in Fig. 2.

B. Connecting Ehrenfest Dynamics with classical dipole radiation in 3D
Here, we show that Ehrenfest dynamics agrees with classical dipole radiation at short times assuming that the initial conditions satisfy (|C 1 |, |C 2 |) = ( 1/2, 1/2). First, consider classical dipole radiation, and let the oscillating dipole (in the z-direction) be situated at the origin. The current takes the form I = −qω sin(ωt + φ)ê z and if the dipole width d is small enough, the current density is J( r) = lim d→0 d · I δ( r) = −µ 12 ω sin(ωt + φ)δ( r)ê z This is the source that acts as input for Maxwell's equations and yield classical dipole radiation.
Thus, this initial electronic state guarantees that Eqns. (53) and (54) will be identical at short times: the EM field from Ehrenfest dynamics will agree with classical dipole radiation exactly. This exact agreement will fail for other initial states or at long times. Even though both methods have the same geometric form, in general, Ehrenfest dynamics would need to be rescaled to match classical dipole radiation in absolute value.

C. Absorption spectra with different incoming field intensities
In this subsection, we plot the absorption lineshape for a variety of different incoming fields and prove that the data in Fig. 13 is occurring in the linear regime. Indeed, according to Fig. 15, the overall absorption signal is linearly proportional to the incoming energy U 0 .
The absorption lineshape can be recovered approximately by simply assuming a Lorentzian signal with width k FGR and a uniform fitting for the total norm. Note that there is a small shift in the maximal signal location: according to Ehrenfest dynamics, the peak is centered at ω 2 0 + ∆ 2 (rather than ω 0 ) where ∆ is the time-averaged off-diagonal coupling in the