Ehrenfest+R Dynamics: A Mixed Quantum-Classical Electrodynamics Simulation of Spontaneous Emission

# Ehrenfest+R Dynamics: A Mixed Quantum-Classical Electrodynamics Simulation of Spontaneous Emission

Hsing-Ta Chen Department of Chemistry, University of Pennsylvania, Philadelphia, Pennsylvania 19104, U.S.A.    Tao E. Li Department of Chemistry, University of Pennsylvania, Philadelphia, Pennsylvania 19104, U.S.A.    Maxim Sukharev Department of Physics, Arizona State University, Tempe, Arizona 85287, USA College of Integrative Sciences and Arts, Arizona State University, Mesa, AZ 85212, USA    Abraham Nitzan Department of Chemistry, University of Pennsylvania, Philadelphia, Pennsylvania 19104, U.S.A.    Joseph E. Subotnik Department of Chemistry, University of Pennsylvania, Philadelphia, Pennsylvania 19104, U.S.A.
###### Abstract

The dynamics of an electronic system interacting with an electromagnetic field is investigated within mixed quantum-classical theory. Beyond the classical path approximation (where we ignore all feedback from the electronic system on the photon field), we consider all electron–photon interactions explicitly according to Ehrenfest (i.e. mean–field) dynamics and a set of coupled Maxwell–Liouville equations. Because Ehrenfest dynamics cannot capture certain quantum features of the photon field correctly, we propose a new Ehrenfest+R method that can recover (by construction) spontaneous emission while also distinguishing between electromagnetic fluctuations and coherent emission.

## I Introduction

Light–matter interactions are of pivotal importance to the development of physics and chemistry. The optical response of matter provides a useful tool for probing the structural and dynamical properties of materials, with one possible long term goal being the manipulation of light to control microscopic degrees of freedom. Now, we usually describe light–matter interactions through linear response theory; the electromagnetic (EM) field is considered a perturbation to the matter system and the optical response is predicted by extrapolating the behavior of the system without illumination. Obviously, this scheme does not not account for the feedback of the matter system on the EM field, and many recent experiments cannot be modeled through this lens. For instance, in situations involving strong light–matter coupling, such as molecules in an optical cavity, spectroscopic observations of nonlinearity have been reported as characteristic of quantum effects.(Thompson, 1998; Solano et al., 2003; Fink et al., 2008) As another example, for systems composed of many quantum emitters, collective effects from light–matter interactions lead to very different phenomena from linear response theory, such as coupled exciton–plasma optics(Törmä and Barnes, 2015; Puthumpally-Joseph et al., 2014, 2015; Sukharev and Nitzan, 2017) and superradiance lasers.(Dicke, 1954; Andreev et al., 1980; Oppel et al., 2014)

The phenomena above raise an exciting challenge to existing theories; one needs to treat the matter and EM fields within a consistent framework. Despite great progress heretofore using simplified quantum models,(Dirac, 1927a, b) semiclassical simulations provide an important means for studying subtle light–matter interactions in realistic systems.(Milonni, 1976) Most semiclassical simulations are based on a mixed quantum–classical separation treating the electronic/molecular system with quantum mechanics and the bath degrees of freedom with classical mechanics. While there are many semiclassical approaches for coupled electronic–nuclear systems offering intuitive interpretations and meaningful predictions,(Kapral and Ciccotti, 1999; Tully, 1998, 1990; Wang et al., 1998, 1999) the feasibility of analogous semiclassical techniques for coupled electron–radiation dynamics remains an open question. With that in mind, recent semiclassical advances, including numerical implementations of the Maxwell–Liouville equations,(Ziolkowski et al., 1995; Slavcheva et al., 2002; Fratalocchi et al., 2008; Sukharev and Nitzan, 2011) symmetrical quantum-classical dynamics,(Miller, 1978; Li et al., 2018) and mean-field Ehrenfest dynamics,(Li et al., 2018) have now begun exploring exciting collective effects, even when spontaneous emission is included.

For electron–radiation dynamics, the most natural approach is the Ehrenfest method, combining the quantum Liouville equation with classical electrodynamics in a mean-field manner; this approach should be reliable given the lack of a time-scale separation between electronic and EM dynamics. Nevertheless, Ehrenfest dynamics are known to suffer from several drawbacks. First, it is well-known that, for electronic–nuclear dynamics, Ehrenfest dynamics do not satisfy detailed balance.(Parandekar and Tully, 2006) This drawback will usually lead to incorrect electronic populations at long times. The failure to maintain detailed balance results in anomalous energy flow (that can even sometimes violate the second law of thermodynamics at equilibrium.(Jain and Subotnik, 2018)) For electron–radiation dynamics, this problem may not be fatal since the absorption and emission of a radiation field may be considered relatively fast compared to electronic–nuclear dynamics.

Apart from any concerns about detail balance, Ehrenfest dynamics has a second deficiency related to spontaneous and stimulated emission.(Li et al., 2018) Consider a situation where the electronic system has zero average current initially and exists within a vacuum environment without external fields; if the electronic state is excited, one expects spontaneous emission to occur. However, according to Ehrenfest dynamics, the electron–radiation coupling will remain zero always, so that Ehrenfest dynamics will not predict any spontaneous emission. In this paper, our goal is to investigate the origins of this Ehrenfest failure by analyzing the underlying mixed quantum–classical theory; even more importantly we will propose a new ad hoc algorithm for adding spontaneous emission into an Ehrenfest framework.

This paper is organized as follows. In Sec. II, we review the quantum electrodynamics (QED) theory of spontaneous emission. In Sec. III, we review Ehrenfest dynamics as an ansatz for semiclassical QED and quantify the failure of the Ehrenfest method to recover spontaneous emission. In Sec. IV, we propose a new Ehrenfest+R approach to correct some of the deficiencies of the standard Ehrenfest approach. In Sec. V, we present Ehrenfest+R results for spontaneous emission emanating from a two-level system in 1D and 3D space. In Sec. VI, we discuss extensions of the proposed Ehrenfest+R approach, including applications to energy transfer and Raman spectroscopy.

Regarding notation, we use a bold symbol to denote a space vector in Cartesian coordinate. Vector functions are denoted as and denotes the corresponding quantum operator. We use for integration over 3D space. We work in SI units.

## Ii Review of Quantum theory for spontaneous emission

Spontaneous emission is an irreversible process whereby a quantum system makes a transition from an excited state to the ground state, while simultaneously emitting a photon into the vacuum. The general consensus is that spontaneous emission cannot fully be described by any classical electromagnetic theory; almost by definition, a complete description of spontaneous emission requires quantization of the photon field. In this section, we review the Weisskopf–Wigner theory(Weisskopf and Wigner, 1930; Scully and Zubairy, 1997) of spontaneous emission, evaluating both the expectation value of the electric field and the emission intensity.

### ii.1 Power-Zienau-Woolley Hamiltonian

Before studying spontaneous emission in detail, one must choose a Hamiltonian and a gauge for QED calculations. We will work with the Power-Zienau-Woolley Hamiltonian(Power and Zienau, 1959; Atkins and Woolley, 1970; Cohen-Tannoudji et al., 1997) in the Coulomb gauge (so that and ) because we believe the method naturally offers a semiclassical interpretation:

 ˆHPZW=ˆHs+ˆHR+ˆH% int, (1)

where

 ˆHR=∫dv{12ϵ0ˆD⊥(r)2+12μ0(∇×ˆA(r))2}, (2)
 ˆHint=−1ϵ0∫dvˆD⊥(r)⋅ˆP⊥(r)+12ϵ0∫dv∣∣ˆP⊥(r)∣∣2. (3)

Here is the vector potential of the EM field and is the transverse field displacement, . is the Hamiltonian of the matter system and will be specified below. Note that the Power-Zienau-Woolley Hamiltonian is rigorously equivalent to the more standard description of QED, but the matter field is now conveniently decomposed into a multipolar form. That being said, in Eq. (1) we have ignored the multipolar magnetization coupling; we are also assuming we may ignore any relativistic dynamics of the matter field.

For QED in the Coulomb gauge, the operator associated with the vector potential can be written in the formalism of second quantization:

 ˆA(r)=i∑kEkωkek(ˆakeik⋅r+ˆa†ke−ik⋅r) (4)

where the matrix element is associated with the frequency , and is the volume of the -dimensonal space. and are the destruction and creation operators of the photon field and satisfy the commutation relations: . The electric and magnetic fields can be obtained according to

 ˆB(r) = ∇×ˆA(r) (5) = i∑kEkk×ekc(ˆakeik⋅r−ˆa†ke−ik⋅r),
 ˆE⊥(r) = −1c∂∂tˆA(r) (6) = i∑kEkek(ˆakeik⋅r−ˆa†ke−ik⋅r),

recalling that in the Coulomb gauge. In term of and , the transverse Hamiltonian of the EM field can be represented equivalently as

 ˆHR=∫dv{ϵ02ˆE⊥(r)2+12μ0ˆB(r)2}=∑kℏωk(ˆa†kˆak+12). (7)

As mentioned above, Eq. (1) has several advantages as a starting point for semiclassical QED. First, the electron–photon coupling is expressed in terms of the electric field and polarization, instead of vector potentials. In fact, the electric dipole Hamiltonian emerges easily if one further makes the long-wavelength approximation.(Cohen-Tannoudji et al., 1997) Second, according to Eq. (1), the displacement is the momentum conjugate to the vector potential , satisfying the canonical commutation relation, . This correspondence will allow us below to use a classical interpretation for propagating the classical EM field.

### ii.2 Electric Dipole Hamiltonian

In practice, for atomic problems, we often consider an electronic system with a spatial distribution on the order of a Bohr radius interacting with an EM field which has a wavelength much larger than the size of the system. In this case, we can exploit the long-wavelength approximation and recover the standard electric dipole Hamiltonian (i.e. a generalized Göppert-Mayer transformation(Cohen-Tannoudji et al., 1997)):

 ˆHint≈−ˆd⋅ˆE⊥(r=0). (8)

In this representation, the coupling between the atom and the photon field is simple: one multiplies the dipole moment operator, , by the electric field evaluated at the origin (where the atom is positioned). This bi-linear electric dipole Hamiltonian is the usual starting point for studying quantum optical effects, such as spontaneous emission.

### ii.3 Quantum Theory of Spontaneous Emission

For a quantum electrodynamics description of spontaneous emission, we need to consider only a simple two-level system

 ˆHs=ε0|0⟩⟨0|+ε1|1⟩⟨1| (9)

which is coupled to the photon field. We assume and . The electronic dipole moment operator takes the form of

 (10)

where is the transition dipole moment of the two states. Using Eq. (6), with a dipolar approximation, the coupling between the two level system and the photon field can be expressed as

 (11)

where the matrix element is given by . We assume the initial wavefunction is and the density matrix is .

Based on the generalization of Weisskopf–Wigner theory (see Appendix A), we can write down the excited state population as

 ρ11(t)=ρ11(0)e−κt (12)

in a coarse-grained sense (meaning that we average over a time scale satisfyting ). The decay rate for a three-dimensional system is given by the Fermi’s golden rule (FGR) rate(Nitzan, 2006)

 κ3D=|μ01|2Ω33πℏϵ0c3. (13)

Similarly, for an effectively one-dimensional system, we imagine a uniform charge distributions in the plane and a delta function in the direction. The effective dipole moment in 1D is defined as . The decay rate for this effectively 1D case is

 κ1D=μ201Ωℏϵ0c. (14)

Eqs (13) and (14) are proven in Ref Li et al., 2018, as well as in Appendix A. Below, we will use to represent the FGR rate for either or depending on context. Note that, in general, Fermi’s golden rule is valid in the weak coupling limit (), which is also called the FGR regime.

We assume that the initial condition of the photon field is a vacuum, i.e. there are no photons at . For a given initial state, , the expectation value of the observed electric field for an effectively 1D system is given by

 ⟨E⊥(x,t)⟩=|c0(t)|c1(0)×R(x,t)sinΩ(t−|x|/c) (15)

where

 R(x,t)=Ωμ01cϵ0e−κ2(t−|x|c)×θ(ct−|x|) (16)

Note that contains an event horizon () for the emitting radiation. The observed electric field represent the coherent emission at the frequency . In a coarse-grained sense, since , the coherent emission has a magnitude given by

 ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯⟨E⊥(x,t)⟩2=|c0(t)|2|c1(0)|2×R(x,t)22. (17)

We note that the coherent emission depends on the instantaneous population of the ground state .

The expectation value of the intensity distribution can be obtained as

 ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯⟨E2⊥(x,t)⟩=|c1(0)|2×R(x,t)22, (18)

which conserves energy of the total system. Note that the variance of the observed electric field (i.e. the fact that ) reflects a quantum mechanical feature of spontaneous emission. For proofs of Eqs. (1518), see Appendix A.

## Iii Ehrenfest Dynamics as ansatz for quantum electrodynamics

Ehrenfest dynamics provides a semiclassical ansatz for modeling QED based on a mean-field approximation together with a classical EM field and quantum matter field.(Li et al., 2018) In general, a mean-field approximation should be valid when there are no strong correlations among different subsystems. In this section, we review the Ehrenfest approach for treating coupled electron–radiation dynamics, specifically spontaneous emission.

### iii.1 Semiclassical Hamiltonian

According to Ehrenfest dynamics, we define a semiclassical Hamiltonian, that approximates Eq. (1):(Mukamel, 1999)

 ˆHSC=ˆHs−∫dvE⊥(r,t)ˆP(r)+∫dv{ϵ02E⊥(r,t)2+12μ0B(r,t)2}. (19)

Here, we treat the EM field operators as classical electrodynamics parameters and , so that we can commute with in the coupling term. Note that the average polarization is defined as . From the semiclassical Hamiltonian, we identify an electronic Hamiltonian

 (20)

which depends only parametrically on the EM field. As a sidenote, we mention that, in Eqs. (13), we have neglected a formally infinite self-interaction energy. If we include such a term, we can argue that, for a single charge center, one can write a slightly different electronic Hamiltonian namely

 (21)

All numerical results presented below are nearly identical using either Eq. (20) or Eq. (21).

### iii.2 Ehrenfest dynamics

Within Ehrenfest dynamics, the electronic system is described by the electronic density matrix while the EM field, and , are propagated by classical electrodynamics. The electronic density matrix evolves according to the Liouville equation,

 ∂∂tˆρ(t)=−iℏ[ˆHel(E),ˆρ(t)]. (22)

The EM field is governed by Maxwell’s equations for classical electrodynamics

 ∂∂tB(r,t) = −∇×E⊥(r,t), (23) ∂∂tE⊥(r,t) = (24)

where the average current is generated by the average polarization of the electronic system

 J(r,t)=∂∂tTrs{ˆρ(t)ˆP(r)}. (25)

The total energy of the electronic system and the classical EM field is111If we substitute Eq. (20) for Eq. (21), the conserved energy is

 Utot(ˆρ,E,B)=Trs(ˆρ(t)ˆHs)+∫dv(ϵ02E⊥(r,t)2+12μ0B(r,t)2). (26)

One of the most important strengths of Ehrenfest dynamics is that the total energy () is conserved (as can be shown easily). Altogether, Ehrenfest dynamics is a self-consistent, computationally inexpensive approach for propagating the electronic states and EM field dynamics simultaneously.

### iii.3 Drawback of Ehrenfest Dynamics: Spontaneous Emission

For the purposes of this paper, it will now be fruitful to discuss spontaneous emission in more detail within the context of Ehrenfest dynamics. In the FGR regime, if we approximate the transition dipole moment of the two level system to be a delta function at the origin, we can show that the electric dipole coupling within Ehrenfest dynamics satisfies the relationship

 Hel01=−ℏκImρ01 (27)

for both 1D and 3D systems. For a 1D system, this relation was derived previously in Ref. Li et al., 2018. For a 3D system, this relation can be derived using Jefimenko’s equation for classical electrodynamics with a current source given by Eq. (25) (see Appendix B).

With Eq. (27), we can convert the Liouville equation (Eq. (22)) for Ehrenfest dynamics into a set of self-consistent, non-linear equations of motion for the electronic subsystem:

 ∂ρ11∂t = −2κ(Imρ01)2, (28) ∂ρ01∂t = iΩρ01−iκImρ01(ρ11−ρ00). (29)

In the FGR regime, because , we can approximate the coherence for a time satisfying so that . We may then define an instantaneous decay rate for , satisfying , where

 kEh(t)=2κρ00(t)sin2Ωt. (30)

Note that does not change much within the time scale . To see the population decay in a coarse-grained sense, let us perform a moving average over and denote the average decay rate as

 ¯¯¯¯¯¯¯¯kEh(t)=1τ∫t+τtdt′kEh(t)=κρ00(t); (31)

here we have used .

This analysis quantifies Ehrenfest’s failure to capture spontaneous emission: Eq. (31) demonstrates that Ehrenfest dynamics yields a non-exponential decay and, when , Ehrenfest dynamics does not predict any spontaneous emission. Interestingly, the Ehrenfest decay rate ends up being the correct spontaneous emission rate multiplied by the lower state population at time .

## Iv Ehrenfest+R Method

Given the failure of Ehrenfest dynamics to capture spontaneous emission fully as described above, we now propose an ad hoc Ehrenfest+R method for ensuring that the dynamics of quantum subsystem in vacuum do agree with FGR decay. Our approach is straightforward: we will enforce an additional relaxation pathway on top of Ehrenfest dynamics such that the total Ehrenfest+R emission should agree with the true spontaneous decay rate. We will benchmark this Ehrenfest+R approach in the context of a two-level system in 1D or 3D space. Note that the classical radiation field is at zero temperature, so we may exclude all thermal transitions from to . We begin by motivating our choice of an ad hoc algorithm. In Sec. IV.3, we provide a step-by-step outline so that the reader can easily reproduce our algorithm and data.

### iv.1 The Quantum Subsystem

#### iv.1.1 Liouville equation

As far as the quantum subsystem is concerned, in order to recover the FGR rate of the population in the excited state, we will include an additional relaxation (“+R”) term on top of the Liouville equation,

 ∂ˆρ∂t=−iℏ[ˆHel,ˆρ]+LRˆρ, (32)

where the super-operator enforces relaxation. For a two-level system, the diagonal elements of the super-operator are chosen to be

 [LRˆρ]11=−[LRˆρ]00=−kRρ11, (33)

and the the off-diagonal elements are determined by conserving the purity of the density matrix (if )222When the subsystem sits exactly on the excited state (i.e. ), the coherence has an undetermined phase and the off-diagonal element evolves as

 [LRˆρ]01=[LRˆρ]∗10=12kR(ρ11ρ00−1)ρ01. (34)

The +R relaxation rate in Eq. (33) is chosen as

 kR≡2κ(1−ρ00)Im[ρ01|ρ01|eiϕ]2, (35)

where is the FGR rate. Eq. (35) is similar to Eq. (30) but with an arbitrary phase . Averaging over a time scale (defined in Eq. (31)), we find

 ¯¯¯¯¯¯kR=κ(1−ρ00). (36)

Thus, the average total decay rate predicted by Eq. (32) is

 κ=¯¯¯¯¯¯¯¯kEh+¯¯¯¯¯¯kR. (37)

In other words, Eqs. (3235) should recover the true FGR rate and correct Ehrenfest dynamics.

The phase in Eq. (35) can be chosen arbitrarily without affecting the total decay rate in a coarse-grained sense (i.e. if we perform a moving average over ). In what follows, we will run multiple trajectories (indexed by ) with chosen randomly. The choice of a random allows us effectively to introduce decoherence within the EM field, so that we may represent the time/phase uncertainty of the emitted light as an ensemble of classical fields. Each individual trajectory still concurrently carries a pure electronic wavefunction. Note that a random phase does not affect the FGR decay rate of the quantum subsystem.

#### iv.1.2 Practical Implementation

In practice, we will work below with the wavefunction , rather than the density matrix . For each time step , the wavefunction is evolved with a two-step propagation scheme:

 (38)

The operator carries out standard propagation of the Schrödinger equation with the electronic Hamiltonian given by Eq. (20). The quantum transition operator implements the additional +R population relaxation from Eqs. (33) and (35). Explicitly, the transition operator is defined by

 (c′0c′1)=ˆT0←1[kR](c0c1) (39)

where

 c′1=c1e−kRdt/2≈c1|c1|√|c1|2−kR|c1|2dt, (40)

and if ,

 c′0=c0 ⎷1+|c1|2|c0|2(1−e−kRdt)≈c0|c0|√|c0|2+kR|c1|2dt (41)

Note that the phase of the wavefunction is preserved by the transition operator in our implementation. In other words, the transition operator simply transfers of population from to without changing the phase of the wavefunction.

Note that, if the subsystem happens to begin purely on the excited state (i.e. or ), there is an undetermined phase in the wavefunction representation. In other words, we can write say and choose randomly. In this case, the transition operator is defined as

 c′1 = eiθe−κdt/2≈eiθ√1−κdt, (42) c′0 = √1−e−κdt≈√κdt. (43)

As emphasized in Ref. Li et al., 2018 and Sec. III, for these initial conditions, and so that the +R relaxation must account for all of the required spontaneous decay.

#### iv.1.3 Energy Conservation

While Ehrenfest dynamics conserves the total energy of the quantum subsystem together with the EM field, our proposed extra +R relaxation changes the energy of the quantum subsystem by an additional amount:

 dUsdt=Tr{ˆHsLRˆρ}=−ΩkRρ11. (44)

Thus, during a time step , the change in energy for the radiation field is

 δUR=ΩkRρ11dt. (45)

For the Ehrenfest+R approach to enforce the energy conservation, this energy loss must flow into the EM field in the form of light emission. In other words, we must rescale the and fields.

### iv.2 The Classical EM fields

At every time step, with the +R correction of the quantum wavefunction, we will rescale the Ehrenfest EM field ( and ) for each trajectory () as follows:

 EℓEh+R=EℓEh+αℓδER, (46)
 BℓEh+R=BℓEh+βℓδBR, (47)

or, in matrix notation,

 (EℓEh+RBℓEh+R)=R[δUℓR](EℓEhBℓEh). (48)

Here, the coefficients and depend on the random phase from Sec. IV.1. In choosing the rescaling function , there are several requirements:

1. and must be transverse fields.

2. Since the +R correction enforces the FGR rate, it is crucial that the rescaled EM field does not interfere with propagating the quantum subsystem. Therefore, the spatial distribution of and must be located outside of the polarization distribution. In other words, , ensuring the electronic Hamiltonian, Eq. (20), does not change much after we rescale the classical EM field.

3. The magnitude of must be equal to times the magnitude of for all in space so that the emission light propagates only in one direction.

4. The directional energy flow must be outward, i.e. the Poynting vector, must have for all (assuming the light is emanating from the origin).

5. On average, we must have energy conservation, i.e. the energy increase of the classical EM field must be equal to the energy loss of the quantum subsystem described in Eq. (45).

Unfortunately, it is very difficult to satisfy all of these requirements concurrently, especially (c), (d), and (e). Nevertheless, we will make an ansatz below which we believe will be robust.

Given a polarization distribution , the rescaling functions for our ansatz are picked to be of the form

 δER = ∇×∇×P−gP⊥, (49) δBR = −∇×P−h(∇×)3P, (50)

where and are chosen to best accommodate requirements (b)–(d). Note that Eqs. (49) and (50) are both transverse fields. The rescaling direction in Eqs. (49) and (50) are motivated by a comparison of the electrodynamical quantum–classical Liouville equation (QCLE) and Ehrenfest dynamics in the framework of mixed quantum-classical theory (to be published). In 3D space, we simply choose , but the dynamics in 1D are more complicated. (In Appendix C, we show numerically that and are good directions of the emanated and fields in 3D. For a 1D geometry, we choose and to minimize the spatial overlap of both and . See Appendix C.)

For a Ehrenfest+R trajectory (labeled by ), the parameters and are chosen to be

 αℓ= ⎷cdtΛδUℓRϵ0∫dv|δER|2×% sgn(Im[ρ01eiϕℓ]) (51)
 (52)

where is self-interference length determined by Eqs. (123124). Note that is associated with the spatial geometry of and . For in the form of a Gaussian distribution (e.g. in a 1D system), we find that the self-interference length is always . By construction, Eqs. (51) and (52) should conserve energy only on average, i.e. an individual trajectory with a random phase may not conserve energy, but the ensemble energy should satisfy energy conservation (see Appendix D).

### iv.3 Step-by-step Algorithm of Ehrenfest+R method

Here we give a detailed step-by-step outline of the Ehrenfest+R method. For now, we restrict ourselves to the case of two electronic states. Given a polarization between the electronic states, before starting an Ehrenfest+R trajectory, we precompute the FGR rate (Eq. (13) or Eq. (14)) and a self-interference length (see Appendix D). At this point, we can initialize an Ehrenfest+R trajectory with a random phase . For time step ,

1. Propagate the wavefunction by and the EM field by Maxwell equations, Eqs. (23) and (24). Here, we denote the EM field as and and is defined by Eq. (20).

2. Calculate the complementary rate (Eq. (35)) and energy change (Eq. (45)).

3. Apply the transition operator (Eq. (39)).

4. Calculate and according to Eq. (51) and Eq. (52) and then rescale the EM field by according to Eq. (4648).

5. Apply absorbing boundary conditions if the classical EM field reaches the end of the spatial grid.

## V Results: Spontaneous Emission

As a test for our proposed Ehrenfest+R ansatz, we study spontaneous emission of a two-level system in vacuum for 1D and 3D systems. We assume the system lies in the FGR regime and the polarization distribution is relatively small in space so that the long-wavelength approximation is valid. For a two-level system with energy difference , we consider two types of initial conditions with distinct behaviors:

1. A superposition state with a fixed relative phase, i.e. where and , :

• The upper state population should decay according to the FGR rate.

• According to Eqs. (15)–(18), the electric field should exhibit coherent emission at frequency .

• The averaged intensity should not equal the coherent emission , i.e. .

2. A pure state with a random phase, i.e. , which corresponds to where is a random phase:

• The upper state population should still decay according to the FGR rate.

• The electric field of each individual trajectory should oscillate at frequency , but the phases of different trajectories should cancel out—so that the ensemble average of the electric field becomes zero, i.e. .

• The averaged intensity should not vanish, i.e. .

Model problems #1 and #2 capture key features when simulating spontaneous emission and can be considered critical tests for the proposed Ehrenfest+R approach. The parameters for our simulation are as follows. The energy difference of the two levels system is . The transition dipole moment is .

For a 1D geometry, we consider a polarization distribution of the form:

 P1D(x)=μ01√aπe−ax2^z, (53)

with and . According to Eq. (53), the polarization is in the direction varying along the direction. For this polarization, the self-interference length is . We use the rescaling function derived in Appendix. C:

 δE1DR(x) = −μ01√aπ4a2x2e−ax2^z, (54) δB1DR(x) = μ01√aπ43a2x3e−ax2^y. (55)

For a 3D geometry, we again assume the polarization is only in the direction, now of the form

 P3D(r)=^zμ012a3/2π3/2e−ar2, (56)

where we use the same parameters for and as for the 1D geometry. The rescaling field in 3D is chosen to be:

 δE3DR(r) = ∇×∇×P3D(r), (57) δB3DR(r) = −∇×P3D(r). (58)

The self-interference length can be obtained numerically as .333Note that the self-interference length strongly depends on dimensionality and is much smaller in 3D than in 1D.

Our simulation is propagated using Cartesian coordinates with for 1D and for 3D. The time step is . Without loss of generality, the random phase for Ehrenfest+R trajectories is chosen from an evenly space distribution, i.e. for .

### v.1 Spontaneous decay rate

Our first focus is an initially coherent state with above. We plot the upper state population and the decay rate of a 1D system in Fig 1(a). As shown in Ref. Li et al., 2018 and summarized in Sec. III above, standard Ehrenfest dynamics does not agree with the FGR decay and cannot be fit to an exponential decay. With Ehrenfest+R dynamics, however, we can quantitatively correct the errors of Ehrenfest dynamics and recover the full spontaneous decay rate accurately. Furthermore, in Fig. 1(b) we prove that, even though Ehrenfest dynamics fails to predict the correct decay rate as a function of initial condition, the decay rate extracted from Ehrenfest+R dynamics agrees very well with the FGR decay rate for all initial conditions. Note that, for the extreme case , Ehrenfest dynamics does not predict any population decay.

Regarding energy conservation, individual Ehrenfest+R trajectories do not conserve energy by design. While the energy loss of the quantum system is roughly the same for every trajectory, the emitted EM energy fluctuates and is not equal to the corresponding quantum energy loss (see Fig. 1(c)). However, an ensemble of trajectories does converse energy on average.

### v.2 Emission Fields in 1D

We now turn our attention to the coherent emission and the intensity of the EM field. We start by considering a 1D geometry. According to Eq. (15), for a given time , the electric field of spontaneous emission can be expressed as a function of and shows oscillatory behavior proportional to for short times. Also, an event horizon is observed at , i.e. no electric field should be observed for because of causality.

We find that the electric field obtained by an individual Ehrenfest+R trajectory shows the correct oscillations at frequency with an additional phase shift. For an initially coherent state, the ensemble average of Ehrenfest+R trajectories agrees with Eq. (15) very well (see Figs. 2(a) and 2(b) for two cases with different initial conditions.) When the initial state is exclusively the excited state, the ensemble average of Ehrenfest+R trajectories vanishes by phase cancellation and we recover (see Fig. 2(c)).

Now we compare the emission intensity and the magnitude of the coherent emission . On the right panels of Fig. 2, we plot the coarse-grained behavior of Ehrenfest+R trajectories. We show that Ehrenfest+R can accurately recover the spatial distribution of both and , as well as the event horizon. Note that in Fig. 2, the electric field and the intensity at large corresponds to emission at earlier times. If we start with a coherent initial state, the relative proportion of coherent emission is given by , see Eqs. (15) and (17). For , the coherent emission is responsible for 50% of the total energy emission at early times (), and the coherent emission dominates later (). Obviously, if we begin with a wavefunction prepared exclusively on the excited state, there is no coherent emission due to phase cancellation among Ehrenfest+R trajectories. In the end, using an ensemble of trajectories with random phases , Ehrenfest+R is effectively able to introduce some quantum decoherence among the classical trajectories and can recover both and .

This behavior of Ehrenfest+R dynamics should be contrasted with the behavior of standard Ehrenfest dynamics, where we run only one trajectory and we observe only coherent emission with . Although the coherent emission obtained by standard Ehrenfest dynamics is close to the quantum result when is small (see Fig. 2(a)), the magnitude of the coherent emission is incorrect in general. The electric field does oscillate at the correct frequency.

### v.3 Emission Fields in 3D

For a 3D geometry, for reasons of computational cost, we propagate the dynamics of spontaneous emission for short-times only (). Our results are similar to the 1D case and are plotted in Fig. 3. For a coherent initial state (Fig. 3 (a), (c)), each Ehrenfest+R trajectory yields an electric field and EM intensity oscillating at frequency , and these features are retained by the ensemble average. For the case of dynamics initiated from the excited state only (Fig. 3 (b), (d)), each trajectory still oscillates at frequency , but the average electric field is zero () so that the peak of the spectrum vanishes.

In Fig. 3, we also compare our result versus the well-known classical Poynting flux of electric dipole radiation. In Fig. 3 (a), our reference is

 I(r,t)=μ0c2Ω4μ20116π2sin2θr2sin2Ω(t−rc), (59)

and, in Fig. 3 (b), our reference is the mean electromagnetic energy flux

 ¯¯¯I(r)=μ0c2Ω4μ20132π2sin2θr2. (60)

In general, Ehrenfest+R dynamics yields a similar distribution as the classical dipole radiation. When initiated from a coherent state, both methods behave as ; when initiated from the excited state, Ehrenfest+R method shows dependence while Ehrenfest dynamics does not yield any emission (not shown in the plot.) However, we note that the intensity of the Ehrenfest+R results is slightly larger than that of classical dipole radiation. This difference is attributed to the fact that the classical dipole radiation includes only coherent emission, which is captured by standard Ehrenfest dynamics. By contrast, Ehrenfest+R dynamics can also yield so-called incoherent emission (), which is effectively a quantum mechanical feature with no classical analogue.

## Vi Conclusions and Future work

In this work, we have proposed a heuristic, new semiclassical approach to quantum electrodynamics, based on Ehrenfest dynamics and designed to capture spontaneous emission correctly. Our ansatz is to enforce extra electronic relaxation while also rescaling the EM field in the direction and . Our results suggest that this Ehrenfest+R approach can indeed recover the correct FGR decay rate for a two-level system. More importantly, both intensity and coherent emission can be accurately captured by Ehrenfest+R dynamics, where an ensemble of classical trajectories effectively simulates the statistical variations of a quantum electrodynamics field. Obviously, our approach here is not unique; a more standard approach would be to explicitly model the EM vacuum fluctuations with a set of harmonic oscillators. Nevertheless, by avoiding the inclusion of high frequency oscillator modes, our ansatz eliminates any possibility of artificial zero point energy loss or other anomalies from quasi-classical dynamics.(Peslherbe and Hase, 1994; Brieuc et al., 2016)

As far as computational cost is concerned, one Ehrenfest+R trajectory costs roughly the same amount as one standard Ehrenfest trajectory, and all dynamics are numerically stable. Implementation of Ehrenfest+R dynamics is easy to parallelize and incorporate within sophisticated numerical packages for classical electromagnetics (e.g. FDTD(Taflove and Taflove, 1998)).

Given the promising results presented above for Ehrenfest+R, we can foresee many interesting applications. First, we would like to include nuclear degrees of freedom within the quantum subsystem to explicitly address the role of dephasing in spontaneous and stimulated emission. Second, we would like to study more than two states. For instance, a three-level system with an incoming EM field can be employed for studying inelastic light scattering processes, such as Raman spectroscopy. Third, we would also like to model multiple spatial separated quantum emitters, such as resonance energy transfer.

At the same time, many questions remain and need to be addressed:

1. The current prescription for Ehrenfest+R approach is fundamentally based on enforcing the FGR rate. However, in many physical situations, such as molecules in a resonant cavity or near a metal surface, the decay rate of the quantum subsystem can be modified by interactions with environmental degrees of freedom. How should we modify the Ehrenfest+R approach to account for each environment?

2. For a quantum subsystem interacting with a strong incoming field, including the well-known Mollow triplet phenomenon(Mollow, 1969) and other multi-photon processes, EM field quantization can lead to complicated emission spectra involving frequencies best described with dressed states. Can these effectively quantum features be captured by Ehrenfest+R dynamics?

3. Finally, and most importantly, it remains to test how the approach presented here behaves when there are many quantum subsystems interacting, leading to coherent effects (i.e. plasmonic excitations). Can our approach simulate these fascinating experiments?

These questions will be investigated in the future.

## Appendix A Generalized Weisskopf–Wigner Theory of Spontaneous Emission

Consider the electric dipole Hamiltonian given by Eq. (11). For comparison with semiclassical dynamics in Sec. V we will now derive the exact population dynamics and the emission EM field of a two level system in vacuum based on Weisskopf–Wigner theory and a retarded Green’s function approach.

### a.1 Dressed state representation

Let be a state of the EM field with one photon of mode , as expressed in a Fock space representation. Let us denote the vacuum state as . For a system composed of an atom interacting with the EM field, the dressed state representation has the following basis (including up to a single photon per mode)(Scully and Zubairy, 1997; Nitzan, 2006)

 |j;k⟩ = |j⟩|0,⋯,1k,⋯,0⟩ (61) |j;0⟩ = |j⟩|{0}⟩ (62)

Here are the wavefunctions for the two level system. For such a setup, the total wavefunction in the dressed state representation must be of the form:

 |ψ(t)⟩=C0,0(t)|0;0⟩+C1,0(t)|1;0⟩+∑kC0,k(t)|0;k⟩+∑kC1,k(t)|1;k⟩. (63)

For spontaneous emission, let the initial wavefunction of the two-level system in vacuum be written as

 |ψ(0)⟩=C0|0;0⟩+C1|1;0⟩ (64)

with . We would like to propagate and calculate as a function of time. We emphasize that, in Eqs. (63) and (64), the Hilbert space is restricted to one photon states.

For visualization purpose, it is helpful to write down the electric dipole Hamiltonian explicitly in matrix form in the dressed state representation,

 H = H0+V =

Here the set is an infinite set of matrices with exclusively diagonal elements for . is an infinite row with corresponding elements

 Vk=iμ01⋅ek√ℏωk2ϵ0Ln (66)

between the vacuum state and a one-photon state with mode . Let us denote the diagonal part of the matrix as the unperturbed Hamiltonian and the off-diagonal part as the coupling Hamilton . Note that the two quantum states in vacuum ( and ) are coupled to two different continuous manifolds and , respectively.

Given that , the manifold will always include a quantum state that is energetically resonant with the state. However, the the manifold will always be off-resonant with for all . Therefore, as the lowest order approximation, we can assume and

 C0,0(t)≈|C0,0(t)|e−iε0t/ℏ (67)

for short time behavior.

### a.2 Retarded Green’s function formulation

We employ a retarded Green’s function formulation(Nitzan, 2006) to obtain the time evolution of and . The retarded Green’s operators are for the full Hamiltonian and for the unperturbed Hamiltonian where is a positive small quantity (). Using Dyson’s identity , we can obtain the retarded Green’s function in a self-consistent expression

 G11(ε) = (68) Gk1(ε) = Vkε−ε0−ℏωk+iη1ε−ε1+iη+i2Γ(ε), (69)

where the self energy is . The self energy can be evaluated by a Cauchy integral identity (ignoring the principle value part). For 1D, we can consider a dipole moment and use the density of states of a 1D system to obtain the self energy as

 Γ1D(ε) = 2iL2π∫dkμ201E2kε−ε0−ℏωk+iη = iμ2012πϵ0ℏc[−iπ(ε−ε0)] = μ201ϵ0ℏc(ε−ε0)

Here, . For 3D, we consider a dipole moment so that