Time-ordering in Heisenberg’s equation of motion as related to spontaneous radiation

Despite many years of research into Raman phenomena, the problem of how to include both spontaneous and stimulated Raman scattering into a unified set of partial differential equations persists. The issue is solved by formulating the quantum dynamics in the Heisenberg picture with a rigorous accounting for both time- and normal-ordering of the operators. It is shown how this can be done in a simple, straightforward way. Firstly, the technique is applied to a two-level Raman system, and comparison of analytical and numerical results verifies the approach. A connection to a fully time-dependent Langevin operator method is made for the spontaneous initiation of stimulated Raman scattering. Secondly, the technique is demonstrated for the much-studied two-level atom both in vacuum and in a lossy dielectric medium. It is shown to be fully consistent with accepted theories: using the rotating wave approximation, the Einstein A coefficient for the rate of spontaneous emission from a two-level atom can be derived in a manner parallel to the Weisskopf–Wigner approximation. The Lamb frequency shift is also calculated. It is shown throughout that field operators corresponding to spontaneous radiative terms do not commute with atomic/molecular operators. The approach may prove useful in many areas, including modeling the propagation of next-generation high-energy, high-intensity ultrafast laser pulses as well as spontaneous radiative processes in lossy media.

The question of how to account for spontaneous radiation is closely connected with the foundations of quantum optics. It was this question that first led Einstein in 1917 away from a classical, isotropic model of radiation emission and toward a highly-directional, quantum-mechanical model characterized by the well-known A and B coefficients of spontaneous and stimulated emission 1 . Ten years later Dirac derived the forms of Einstein's A and B coefficients using the new quantum theory developed by Schrödinger and Heisenberg, simultaneously postulating the quantization of the electromagnetic field as a collection of harmonic oscillators for the first time 2,3 . It was found that the harmonic oscillator model was sufficient to explain the light-matter interaction, and it subsequently became a key feature of quantum optics theory.
In the following years, opinions wavered over whether spontaneous radiation resulted primarily from the fluctuations of the quantized electromagnetic field or from the quantum analogue of the classical radiation reaction. In 1963 Jaynes and Cummings showed that spontaneous emission could result from a model in which the fields remain classical and the radiating electron experiences radiation reaction 4 . Ultimately, this model was found to give clearly spurious results for anti-Stokes Raman scattering 5 , since in a semi-classical model the spontaneous transition rate depends on both the initial and final states, while for a fully quantum model the transition rate depends only on the initial state [5][6][7] . In addition, it had been known since 1958 that semi-classical radiation theory can only account for the coherent aspect of spontaneous emission and neglects the incoherent contribution [7][8][9] . In 1973, Senitzky brought the controversy between vacuum fluctuations and radiation reaction to a satisfactory conclusion when he showed that their roles in spontaneous emission are "two sides of the same quantum-mechanical coin". Each phenomenological picture may be brought into clearer focus through different orderings of the relevant quantum operators 10 . Vacuum fluctuations and radiation reaction are therefore related through a fluctuation-dissipation theorem 11 . This revelation, however, does not mean that further insight cannot be gained from semi-classical models of spontaneous radiation 12 .
Even as the theory of Raman scattering was playing a minor role in the elucidation of the nature of spontaneous radiation processes in general, research into Raman phenomena in their own right was developing apace.
Placzek had detailed the quantum theory of spontaneous Raman scattering in 1934 13 , but, with the invention of the laser, focus soon shifted to stimulated Raman processes. Beginning with Shen and Bloembergen in 1965, the theory of stimulated Raman scattering in both the forward and backward directions was formulated [14][15][16][17]  www.nature.com/scientificreports/ works used both classical and semi-classical models to account for pulse propagation effects. As such, the effects of spontaneous Raman scattering were neglected. The first attempts to unify both spontaneous and stimulated Raman scattering within a single theoretical chassis accounting for pulse propagation were undertaken by Raymer and Mostowski in a pair of papers published in 1981 18,19 . Their work assumed a heavily populated ground state and introduced a quantum statistical Langevin operator corresponding to "collision-induced fluctuations" that was shown to reproduce predicted behavior in both the spontaneous and stimulated Raman regimes. This kind of approach is at best incomplete, since it is known both experimentally and quantum-theoretically that, while collisions do induce Raman coherences 20 , an isolated molecule may spontaneously scatter Raman radiation independently of collisions. In subsequent work investigating macroscopic fluctuations in the transient Raman regime [21][22][23][24][25] , they make a distinction between spontaneous Raman radiation induced by vacuum fluctuations and by collisions, and continue to use the Langevin operator as a phenomenological tool. The limitations of this kind of model immediately become apparent as soon as the assumption of a heavily populated ground state is relaxed, as in the case of laser pulses of sufficient intensity to induce a significant fraction of molecules into a vibrationally excited state. Such considerations are especially important in the development of stimulated Raman backscattering amplification (SRBA) of high-energy, high-intensity ultrafast laser pulses, which is regarded as a next-generation technology to surpass the intensity limitations of chirped pulse amplification (CPA). SRBA is envisioned to possibly enable exawatt and even zettawatt laser powers 26 , but two decades of effort in this area has yielded little success 27 .
For laser pulses with this level of intensity, experiment indicates that the interplay between spontaneous and stimulated Raman radiation becomes vitally important. What kind of theoretical model can accurately describe this system? A semi-classical approach readily accommodates the dynamics of pulse propagation but fails to describe the all-important contribution from spontaneous Raman processes in a satisfactory manner. Only a fully quantum-mechanical model, in which the fields are also quantized, will be able to account for the diversity of behavior. Yet, using a conventional Raman Hamiltonian in Heisenberg's equation of motion yields only coherent, stimulated terms.
Spontaneous Raman processes are typically described in the interaction picture through calculation of transition probabilities using the Dyson U-matrix operator 28 . Specifically, the spontaneous terms result from normalordering of field operators within a perturbative time-ordered theoretical framework. The dynamics are usually solved to a few orders of the perturbation, meaning that the solution is valid only for short times. Complete knowledge of the behavior's system for all time requires calculation of an infinite number of perturbative orders.
Such prohibitively arduous calculations can be circumvented through incorporating both time-and normalordering of operators into Heisenberg's equation of motion, from which the approximate, time-ordered perturbative methods were originally derived. In order to do so, it is necessary to point out an insight originally discovered by Feynman 29 that commutators arise at least in part from the chronological order of the operators involved. In other words, for two canonically conjugate operators p and q, or, in the case of creation and annihilation operators â † and â, where in Eqs. (1) and (2) the indices correspond to time-ordered infinitesimal time intervals. It should be noted that the indices need not be sequential; they need only to be in order, with greater values occurring later (to the left) of earlier values (placed to the right). In addition, time indices play a role only for non-commuting operators 30 .
Heisenberg's equation of motion for an operator Â may also be written as an explicit time-ordered expression 29 : In a system of differential equations for non-commuting operators, the correct time-ordering of the operators cannot be neglected when substitutions are made in the course of evaluation. This is shown first for a two-level Raman system. Then, to demonstrate the robust validity of the approach, it is demonstrated for the much-studied two-level atomic system, both in vacuum and in a lossy dielectric.

Results
The two-level Raman system. Here, the technique is demonstrated for spontaneous Raman scattering from a two-level molecular system defined by a ground state �b⟩ and an excited vibrational state �a⟩ , with molecular operators ̂ that obey the Heisenberg equation of motion and are defined by The system is governed by an unperturbed molecular Hamiltonian Ĥ M , a Raman interaction Hamiltonian Ĥ I , and a field Hamiltonian Ĥ F such that (4) nm = �n⟩⟨m� www.nature.com/scientificreports/ where and the coupling factor G is given by where ij = e ij is the dipole moment between levels �i⟩ and �j⟩, is the polarization vector of field mode , and p and s designate the pump and Stokes fields, respectively. Taking into account the time-ordering aspects of the Heisenberg equation of motion as outlined above (and using overset indices to indicate time-order), the equations of motion for the molecule and field are then where ̂i j =̂i j e −i ij t , â =â e i t , and Δ ≡ ab − ( p − s ) is the two-photon detuning. Eqs. (11)-(16) introduce the convention of initially placing the molecular operators to the right of the field operators.
The time-and normal-ordered solution of these equations can be illustrated in a straightforward manner. For example, the operator ̂a a can be written as the sum of two functions, each operating in single time bin only such that Allowing the other operators to be expressed in the same manner yields, for ̇̂ ( 2) aa , www.nature.com/scientificreports/ In going from Eqs. (21) and (22), it is assumed that the molecular operators are normal-ordered functions of the field operators. In Eq. (21), the operators are time-ordered but not normal-ordered, while in Eq. (22) the commutation relations have been used. Since in Eq. (22) the operators are both time-and normal-ordered, the time indices may be omitted. The partial derivative terms may be evaluated through simple integration of Eqs. (11)- (16), and, although implicit dependence on the field operators may be accounted for, it is shown below that differentiation with respect to explicit field operator dependence is sufficient to achieve accurate results.
For example, an analytical perturbative solution in powers of a parameter can be found if, for any operator Â and the coupling constant G, the substitutions are made. In this case, using the procedure outlined above and the initial conditions ̂a a (0) = 0 , ̂b b (0) = 1 , it is straightforward to show that up to second order in where n p and n s are the number of pump and Stokes photons, respectively. Equation (27) is the probability of finding the molecule in the excited vibrational state �a⟩ as a function of time. It is identically equal to the probability of a spontaneous Raman transition from state �b⟩ to �a⟩ as calculated with the conventional U-matrix method up to second order using the conventional interaction Hamiltonian V = −̂ ⋅̂ .
Again, such a perturbative solution is accurate only for small times. As the time increases beyond the domain of the perturbation, a numerical calculation using time-and normal-ordered Eqs. (11)-(16) yields greater accuracy, since in fact it takes account of all perturbative orders. This kind of numerical solution (conventional 4th order Runge-Kutta, and where the two-photon detuning Δ = 0 ) is shown in Figs. 1 and 2 for the case of a single pump photon scattering from a molecule and spontaneously generating a Stokes photon. Figure 1 shows that, for The end result is that the pump photon is annihilated, a Stokes photon is created, and the molecule has transitioned to the excited vibrational state. These results clearly show the advantage of the approach outlined above when the value of t increases beyond perturbative domains. Of interest in Fig. 2 is the fact while the Stokes intensity increases by one photon, the Stokes field remains at zero. Physically, this corresponds to the unpredictability of the phase of a spontaneously generated photon. Spontaneous Raman scattering is an incoherent process. If, alternatively, the incoming Stokes field were nonzero, the scattering process would contain both coherent and incoherent components, resulting in a corresponding change in the Stokes field amplitude. As they stand, Eqs. (11)-(16), when both time-and normal-ordered, model stimulated (coherent) and spontaneous (incoherent) scattering processes separately. In order to model spontaneously initiated stimulated Raman scattering, an additional equation is required to model the stochastic manner in which incoherent, spontaneously generated Stokes photons contribute to coherent processes. This is generally done using a Langevin operator formalism 19,31 . The quantum-mechanical formulation of the Langevin approach is described in detail in the work of Melvin Lax (see 32 and references therein). The equation of motion for a quantum operator Â (t) subjected to a random Markoffian noise force F (t) is given as 32 : where For a system with Gaussian noise statistics, the moments in Eqs. (29) and (30) are sufficient to describe the dynamics. Of key importance is the definition of the generalized time-dependent Einstein relation for the diffusion matrix 32,33 The quantity 2⟨D (t)⟩ in Eq. (31) measures the extent to which the chain rule for differentiating a product has been violated 32 . Since the time-ordered method described above assumes the molecular operators are functions of the fields, it is appropriate to calculate the relationship between the intensity and field operators: Figure 2. Full evolution of operators in a one-photon spontaneous Raman scattering process. Blue curves correspond to pump operators, while red curves correspond to Stokes operators. Solid curves correspond to intensities â † a , while dotted curves correspond to fields â . Curves corresponding to ̂b b and ̂a a follows those of â † p a p and â † s a s , respectively. The dotted black curve corresponds to the perturbative analytical solution for ̂a a in Eq. (27). Here, the two-photon detuning Δ = 0. www.nature.com/scientificreports/ The modified equations of motion for the fields are then Although the above has been calculated for the Stokes field, the same procedure may be applied to the pump field, as well. It may also be applied to the two-level atomic system outlined in the next section in order to characterize the onset of superfluorescence [34][35][36] .
Of significant importance is the fact that the diffusion matrix relating the intensity to the fields is fully timedependent and derived from the total system Hamiltonian, unlike in similar approaches 19,31,35,36 . This has great consequences when higher-order interaction terms are included in the Hamiltonian. These higher-order interaction terms may very well become important for ultra-intense laser pulses, specifically of the kind used for stimulated Raman backscattering amplification (SRBA) in plasma 27 , and may affect estimates of parasitic instabilities 37 .
The two-level atomic system in vacuum. Traditionally, stimulated and spontaneous radiation have been investigated assuming a model two-level atom in which a single electron transitions between upper and lower energy states, and there is a large body of significant, historical literature detailing the development of scientific understanding of this system (see the references listed in the "Introduction").
The aim of this section is to show in a heuristic manner how the present technique fits into this historical narrative. In particular, it is shown how the Einstein A coefficient that characterizes the rate of spontaneous emission may be derived in a straightforward manner that parallels the Weisskopf-Wigner approximation 38 . The Lamb frequency shift is also calculated.
A two-level atom is assumed with electronic ground state �b⟩ and excited state �a⟩ , with atomic operators ̂i j that obey the Heisenberg equation of motion. The dynamics of the system are governed by an unperturbed Hamiltonian Ĥ A , an interaction Hamiltonian Ĥ I , and a field Hamiltonian Ĥ F such that where where P ij ≡ ⟨i� ij ⋅ �j⟩ , ij ≡ e ij , E ≡ ℏ ∕2 0 V 1∕2 , and V is the volume of integration over the field mode. Letting â =â exp −i t and ̂i j =̂i j exp i ij t and taking the rotating wave approximation, the technique outlined above for the Raman system yields the equations of motion www.nature.com/scientificreports/ where n =â † a is the photon number operator. In order to make spontaneous radiative terms explicit in Eqs.  www.nature.com/scientificreports/ are entirely spontaneous, the field operators are identically zero and do not change. Physically, it corresponds to the fact that the phase of the spontaneously emitted radiation cannot be predicted. Taking this into account, it is then straightforward to show that the equation of motion for the population of the excited state �a⟩ is given by Assuming closely-spaced field modes, where the factor of 2 out front accounts for two polarization states. Since |P ab | 2 = P 2 ab cos 2 , where cos = ij ⋅ , Eq. (59) can be written as Parallel to the Weisskopf-Wigner approximation 38,39 , it can be seen that the time integral is non-negligible only in the region = ab . 3 can therefore be replaced with 3 ab , and the lower limit on the frequency integral is extended to −∞ . Furthermore, since Equation (61) can be written as where is the well-known Einstein coefficient that characterizes the rate of spontaneous emission. Following a similar argument, it is straightforward to show that for the population of the ground state and therefore Equations (63), (65), and (66) show that, while the dynamics of the ground state population depend on the excited state population, the dynamics of the excited state population are independent of the ground state population. Consequently, an atom initially in the excited state can decay spontaneously to the ground state, but an atom initially in the ground state can never spontaneously transition to the excited state. This is in good agreement with energy conservation and all accepted theories of spontaneously emitted radiation 11 .
Finally, a total photon number operator N can be defined such that Using Eq. (58), it can easily be shown that which again is in agreement with accepted theories 39 . The Lamb shift can be calculated through analysis of the time-dependence of the atomic de-excitation operator ̂b a and takes the form of a slight frequency shift 40 . The equation of motion for the atomic de-excitation operator in the absence of an externally applied field is Since ̂b a is assumed to be a slowly varying function of time, it may be taken outside the integral. As time advances 40 ,  The two-level atomic system embedded in a lossy dielectric. Spontaneous radiation in lossy media is an ongoing area of research [43][44][45][46] . Here, the simple case of an infinite, homogeneous, and isotropic medium is treated, with the understanding that the interested reader may generalize the approach to any Hamiltonian that sufficiently describes the system under investigation. The canonical quantization of the electromagnetic field within a dispersive and lossy dielectric medium obeying the Kramers-Kronig relations is performed following the approach of Huttner and Barnett 47 , resulting in boson creation and annihilation operators Ĉ † and Ĉ which correspond to polaritons instead of free-space photons. Significantly, in this approach, the one-to-one relation between the wave vector and frequency of the field in the dielectric is lost: in a lossy dielectric medium, the fields are expressed as continuous double integrals over both the wave vector and the frequency, which are treated as independent variables 47 .
The polariton creation and annihilation operators obey the commutation relation and the transverse electric field has the form where E ≡ ℏ 2 c ∕2 0 V 1∕2 , V is the volume of integration over the field mode, 2 c = 2 ∕ 0 , is the charge density, and is the effective mass density associated with the harmonic polarization (see 47,48 for details). The function f ( , ) is defined as where the dimensionless function ( ) characterizes the frequency dependence of the coupling between the electromagnetic and dressed matter fields. A local field correction factor ( + 2)∕3 has also been included (consistent with the approach in 48,49 ). The local field correction factor can in general be derived from first principles through a more detailed consideration of Umklapp processes 50 . Also, the complex dielectric function = � + i �� obeys the Kramers-Kronig relations and is defined as A two-level impurity atom is assumed to be embedded within the medium at the origin and has electronic ground state �b⟩ and excited state �a⟩ , with atomic operators ̂i j that obey the Heisenberg equation of motion. The dynamics of the impurity atomic system are governed by an unperturbed Hamiltonian Ĥ A , an interaction Hamiltonian Ĥ I , and a field Hamiltonian Ĥ F such that where where P ij ≡ ⟨i� ij ⋅ �j⟩ , ij ≡ e ij , Ĉ ( , ) =Ĉ( , ) exp (−i t) , and ̂i j =̂i j exp i ij t . Taking the rotating wave approximation, the equations of motion for the impurity atomic system are where the modified spontaneous decay rate A L in a lossy, dispersive dielectric medium is which is identically equal to accepted values 48 . www.nature.com/scientificreports/ The Lamb shift can also be calculated in an analogous manner to the case of an atom in vacuum. The equation of motion for the atomic de-excitation operator in the absence of an externally applied polariton field is where Λ L is the modified Lamb frequency shift given by

Discussion and conclusion
It may be instructive to discuss the above method in view of the widely used techniques outlined by Shaul Mukamel that calculate contributions to spontaneous radiative processes, as described in 51 . In the latter approach, the calculation of the spontaneous light emission signal mirrors the calculation for the nonlinear response function S (3) (t 3 , t 2 , t 1 ) , which itself results from a perturbative expansion of the density matrix operator ̂(t) according to the U-matrix method. The difference in the case of spontaneous radiation is that the emitted field has no photons initially (it is a vacuum state) and therefore must be treated quantum-mechanically instead of classically, as is commonly done in evaluating S (3) (t 3 , t 2 , t 1 ) . In practice this involves evaluating the action of field operators on field modes, for example â † �n⟩ = √ n + 1�n + 1⟩ and â�n⟩ = √ n�n − 1⟩ . Such operations are the origin, mathematically, of the spontaneous radiation in the U-matrix method, as well. Specifically, it is the +1 term in √ n + 1 that corresponds to spontaneous radiation 52 . In fact, this +1 term results directly from the evaluation of the commutator â,â † , such that it is rigorously correct to state that â † �n⟩ = � n + �â ,â † � �n + 1⟩ (see Sect. 2.3 of 53 ).
Mathematically, then, spontaneous radiation results directly from the noncommutativity of the field creation and annihilation operators. In order to introduce spontaneous radiation into a physical theory, evaluation of field commutators must be correctly implemented. Here, in this work, instead of calculating the action of operators on field modes as in 51 , evaluation of commutators occurs in normal-ordering of time-ordered operator products. The advantage of working directly with the time-ordered Heisenberg equation of motion is that its validity for a given Hamiltonian remains valid for all times. Current widely-used approaches that use a perturbative expansion "may be valid at short times but will always break down at longer times" 51 , and this is clearly shown in Fig. 2. This is not to say that, in many circumstances, the power and versatility of such perturbative approaches (and their variants 54,55 ) may not prove to be more convenient and of sufficient accuracy for small t.
It should be noted that after the initial expression of the equations of motion in Eqs. (11)-(16), (45)- (51), and (81)-(84), it is assumed that the chemical (either molecule or atom) operators are functions of the field operators and therefore do not commute with them. Conventionally, it is assumed without nuance that the chemical and field operators commute with each other 56 . This is not the case, however, for the chemical's interaction with spontaneously emitted or scattered radiation. On the contrary, it has been shown that as time progresses, the chemical source operators that give rise to the spontaneous radiation evolve into the joint chemical/field Hilbert space 56,57 , and thus fail to commute. In this sense, it is not surprising that in Eqs. (11)-(16), (45)-(51), and (81)-(84), only the spontaneous terms are affected by an assumption of non-commutativity between the chemical and field operators, while the coherent terms remain unaffected.
In addition, it has long been known that, in the case where chemical operators initially commute, the complementary contributions of radiation reaction and vacuum fluctuations to spontaneous radiative processes can be apportioned in differing amounts by changing the order of the operators. Consequently 7 , it is possible, for example, to express the total Lamb frequency shift Λ as a sum of contributions from both radiation reaction (RR) and vacuum fluctuations (VF) such that Λ = Λ RR + Λ VF . In the case where the operators are normal-ordered, Λ RR = Λ and Λ VF = 0 . When the operators are anti-normal-ordered, Λ RR = −Λ while Λ VF = 2Λ . In the case of symmetric-ordering, Λ RR = 0 and Λ VF = Λ . Again, assuming the chemical and field operators initially commute, it is possible to shuffle between these different orderings to emphasize various physical interpretations 7 .
In the method described above, however, the chemical and field operators do not commute, because the chemical operators are assumed to be functions of the fields. This assumption is necessary in order to make use of the initial conditions ̂a a (0) , ̂b b (0) . Doing so effectively restricts the ordering of operators to normal-order, since in this context anti-normal ordering yields an excited state population ̂a a that is dependent on the ground state population ̂b b , which is unphysical. In normal-ordering, only radiation reaction contributes to spontaneous radiation. In anti-normal-ordering, both radiation reaction and vacuum fluctuations contribute to spontaneous radiation. In symmetric-ordering, spontaneous radiation is entirely due to vacuum fluctuations, and therefore it is expected that, like anti-normal-ordering, symmetric-ordering will also yield unphysical results. As such, the method outlined in this work attributes spontaneous radiation entirely to radiation reaction.
In the absence of freedom to choose the ordering of chemical and field operators, however, the method described above allows a freedom in the temporal extension of the operators. For example, the convention of initially placing the chemical operators to the right of all the field operators has been used, and as a result the temporal extension of the chemical operators is forward in time. However, equivalent results may also be obtained if the chemical operators are initially placed to the left of all the field operators, and the temporal extension of the chemical operators is backward in time: www.nature.com/scientificreports/ Either of the options in Eq. (97) may be used to yield physically meaningful results, as long as the chemical operators are assumed to be normal-ordered functions of the field operators.
In conclusion, it is evident from the above that the Heisenberg equation of motion can account for spontaneous radiative processes in a straightforward way when both time-and normal-ordering of the operators are taken into account. Moreover, this is done within the context of a unified set of partial differential equations, which are amenable to pulse propagation calculations and simulations. While the technique is immediately useful for investigating spontaneous radiative processes in lossy media as well as modeling the interaction between spontaneous and stimulated Raman scattering from high-intensity laser pulses, its simplicity and generality suggest that it may yield fruitful results in many other areas of physics.