{'fraction_non_alphanumeric': 0.053880116740252934, 'fraction_numerical': 0.019207263837960536, 'mean_word_length': 4.367519078859285, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Description of spontaneous photon emission and local density of states in presence of a lossy polaritonic inhomogenous medium

We provide a description of spontaneous emission in a dispersive and dissipative linear inhomogeneous medium based on the generalized Huttner-Barnett model [Phys. Rev. A 46, 4306 (1992)]. Our discussion considers on an equal footing both the photonic and material fluctuations which are necessary to preserve unitarity of the quantum evolution. Within this approach we justify the results obtained in the past using the Langevin noise method that neglects the removal of photonic fluctuations. We finally discuss the concept of local density of states (LDOS) in a lossy and dispersive inhomogeneous environment that provides a basis for theoretical studies of fluorescent emitters near plasmonic and polaritonic antennas.


I. INTRODUCTION
In the recent years the theoretical problem of describing the coupling of single fluorescent quantum emitters with a metallic nano-particle supporting surface plasmon (SP) modes has become very urgent due to many applications envisioned with photonic and quantum information processing technologies at the nanoscale. In particular the concept of local density of states (LDOS) [1,2] is central since it provides a figure a merit for quantifying the coupling of quantum emitters to plasmonic systems and it plays a central role in recent studies using nearfield optical microscopes [1,[3][4][5][6][7]. However, one of the main issue with plasmonic systems is that they are intrinsically dissipative and that a self-consistent quantum electrodynamics (QED) description of plasmons, i.e. respecting rigorously the unitarity of time evolution, involves necessarily the inclusion of additional degrees of freedom associated with fluctuating currents and dipoles in the metal. This problem has been generally studied using the Langevin noise method advocated originally by Gruner and Welsch [8][9][10] and it has been intensively used meanwhile in QED in macroscopic media [11][12][13][14][15][16][17][18][19][20], e.g., for calculating the coupling of quantum fluorescent emitters to fluctuating fields near plasmonic antennas [21][22][23][24][25][26][27][28][29]. The justification and consistency of the approach (which intuitively generalizes earlier fundamental results obtained by Rytov and Lifschitz in the context of Casimir force calculations [30][31][32][33][34][35][36][37][38][39][40][41]) has been however the subject of some controversies in the past and a rigorous mathematical derivation based on a QED hamiltonian formalism valid for the most general inhomogeneous media has been looked for during years. The main issue is that the Langevin noise approach was originally motivated by the Huttner-Barnett Hamiltonian formalism [42][43][44][45][46][47][48] which for a homogeneous dissipative and dispersive dielectric material, respecting Kramers-Krönig causality constraints, generalizes the historical Hopfield-Fano model for polaritons in bulk media [49,50]. The Huttner-Barnett formalism leads directly to the results postulated by Gruner and Welsch in their seminal work and in particular it was used in the calculation of the emission rate of fluorescent dipoles imbedded in a lossy dielectric matrix [42][43][44][45][46][47][48]. However, the rigorous equivalence between the two formalisms in the more general inhomogeneous situation was still lacking for years.
Very recently, based on important calculation by Suttorp and coworkers, several works [52][53][54][55][56][57][58] proposed a mathematical justification of the Langevin noise method based on the direct generalization of the Huttner-Barnett model for inhomogeneous media. However, we showed in two publications [59], and [60] that these earlier derivations overlooked the role of boundary conditions concerning the value taken by dielectric susceptibilities at spatial infinity. We showed that to preserve the fundamental unitarity of the quantum evolution one must necessarily add a contribution associated with vacuum photonic fluctuation to the fluctuating source current terms calculated by Gruner and Welsch [8][9][10]. Actually, it implies that, rigorously speaking, the formalism introduced by Gruner and Welsch [8][9][10] relies on the assumption that the surrounding dielectric environment is necessarily lossy and this, even at spatial infinity. Contrarily to a widespread belief this situation is not generally applicable to nanoparticle antennas which are by definition spatially localized and very often surrounded by vacuum in the calculations.
The main objective of the present work is to leverage on the analysis started by us in [59,60] in order to describe the coupling of quantum emitters to fluctuating currents and fields. In the present paper we will consider specifically the regime of spontaneous emission and the coupling of transition dipoles to coherent laser sources (i.e., the derivation of the optical Bloch equations). We will demonstrate that the most general formalism derived in [60] actually allows us to generalize the results obtained for spontaneous emission that were only rigorously derived for a dipole in vacuum. We will show that for this specific problem the dynamical equations and emission rates actually agree with previous results obtained using the Gruner and Welsch method without the assumptions made in this phenomenological approach.
The layout of this paper is a follows: In section II we provide a general description of Huttner-Barnett formalism [43,59,60] using the dual Lagrangian method given in [59,60]. We include in the model the interaction between a fluorescent molecule and the photonic and dielectric environment. In section III we analyze the spatio-temporal evolution of the electric field operator and separate the contributions associated with free photons from fluctuating currents (associated with the dielectric medium and the molecule). In sect IV we study within the Wigner-Weisskopf method the spontaneous emission of a two-level atom in a general lossy and dissipative inhomogeneous environment. We conclude with a discussion in section V about LDOS and photonic wave functions associated with spontaneously emitted photons and we derive the dynamics associated with the optical Bloch equations.

II. THE DYNAMICAL EQUATIONS OF A TRANSITION DIPOLE IN A FLUCTUATING PHOTONIC AND DIELECTRIC ENVIRONMENT
We start with the dual-Lagrangian density for the coupled system and where we use the Coulomb gauge constraint ∇ · F(x, t) = 0 for the electric potential. The electric potential F(x, t) is the dual of the usual magnetic potential A(x, t) and we showed that it is a specifically adapted choice for quantization of electromagnetic problems involving dipole densities P(x, t) [59,60]. In this description the material part associated with the inhomogeneous dielectric medium reads and corresponds to the Huttner-Barnett model [43,59,60]. In this model the field X ω (x, t) completely characterizes the fully causal dielectric environment satisfying Kramers-Krönig relations [43,59,60]. The contribution [61] is associated with the Schr'odinger matter fields Ψ a (x, t) and Ψ * a (x, t) which describe the atomic dipole of mass M in the external potential V (x). In the second quantization formalism we expand the atomic wave function operator as Ψ a (x, t) = m b m (t)ψ m (x) where ψ m (x) (labeled by m) are energy eigenstates of the time independent Schrodinger equation E m ψ m = (− 2 ∇ 2 2m + V )ψ m defined for the energy E m and forming a complete orthogonal set such that d 3 xψ * m (x)ψ n (x) = δ n,m . The usual second quantization procedure leads to the (fermionic) anti-commutators . This is clearly consistent with the usual canonical quantization procedure since the Lagrangian function which together with the canonical anti-commutation rules {b m , Π bn } = i δ n,m , etc... imply Eq. 2.4. Moreover, in this model the dipole density P(x, t) includes contributions from the surrounding dielectric medium P diel. (x, t) and of the Schrodinger field P Ψ (x, t).
For the dielectric part we use in agreement with Huttner and Barnett (see [43,59,60]): For the Schrodinger field we here model the dipole density fluid using the formula: where µ m,n = e d 3 ξψ * m (ξ)ξψ n (ξ) = µ * n,m denote some transition dipoles and e < 0 is the electron charge. The model relies on the Openheimer approximation which separates the dynamical evolution into a center of mass motion x(t) (supposed here irrelevant since the atom is not moving) and a relative motion ξ(t) corresponding approximately to the electron motion with respect to the nuclei. Here we introduced ∆(x − x 0 ) which is a narrow peaked distribution centered on the coordinate x 0 associated with the center of mass of the atomic system and normalized as d 3 x∆(x − x 0 ) = 1. ∆(x − x 0 ) reduces to the Dirac distribution δ 3 (x − x 0 ) in the point-like dipole limit. Still, in order to avoid singular divergences we will here keep ∆(x − x 0 ) finite.
Inserting the total dipole density P = P Ψ + P diel. into the Lagrangian Eq. 2.1 leads to the Euler-Lagrange equations [59,60] corresponding to Maxwell's equations in the medium with full polarization P.
For the Schrodinger fields we must be more careful. Indeed by using the Euler-Lagrange classical method we should deduce the evolution equationṡ is a constant of motion. However, at that stage Eq. 2.10 is not fully quantum since in the present formalismĒ(t) and b n (t) do not commute and the order of operators must be specified. To derive the equation of motion for the matter field we should therefore introduce the Hamiltonian of the system. Using the previous Lagrangian we derive the full Hamiltonian H(t) reads: is clearly an integral of motion as it can be for instance proven by using the Poynting theorem discussed in [59] and which reads: Integrating the dissipated power d 3 xJ · E leads, after some simple calculations which will not be repeated here, to d dt H(t) = − Σ∞ dΣ · cE × B → 0. This shows that the energy E = H(t) is naturally conserved in the limit of an infinite integration volume and supposing that the field decays fast enough at infinity (an hypothesis which makes sense only with some additional physical assumptions). Moreover, to build a quantized version of this field theory we insert in Eq. 2.11 the usual normal order product convention : [...] : [59,60] which removes some infinite spurious quantities in the energy (here the normal order means that we let bosonic and fermionic operators to commute or respectively anticommute, additionally the fermionic and bosonic operators mutually commute). Using the Heisenberg evolution law i d dt O(t) = [O(t), H(t)] for any quantum operator O(t) we thus easily deduce once again Eqs. 2.8 and 2.9. Furthermore, a rigorous application of commutation and anticommutation rules shows that Eq. 2.10 is actually correct if the electric field operator is indeed positioned before the annihilation operator b n (taking into account this constraint it is also easy to prove the constancy of N (t)).
We will now in the following consider only the case of the idealized two-level atom. For this we introduce the ground state m = 1 and the excited state m = 2 such as and we can therefore write the Schrodinger part of the Hamiltonian H Ψa (t) as We get for the evolution equation of σ z (t): ). Importantly the field operatorĒ(t)+P Ψ (t) commutes with atomic operators such as σ(t), σ z (t) and this makes the calculation easier.
We similarly introduce the evolution equation for σ(t) as :σ which again involves the fieldĒ(t) +P Ψ (t).

III. DESCRIPTION OF THE FULL ELECTRIC FIELD OPERATOR IN PRESENCE OF A POLARIZABLE MEDIUM AND A FLUORESCENT DIPOLE
At that stage and before to solve the previous dynamical problem we should give a reminding concerning the formal structure of the fully quantized electric field operator discussed in [59,60].
We showed in [60] that the total electric field evaluated at time t and at point x is given by In this equation t 0 is an initial time (which can be sent to the remote past denotes the field associated with free-space propagating photons existing in absence of any dielectric surrounding medium and atoms [59,60]. The quantization of this free field expanded into plane wave modes (labeled by the wave-vector k α and the transverse polarizationǫ α,j , with j = 1, 2) leads to the following operator expressions [60]: V (V is the infinite box volume of the Born von Karman modal expansion method) and where the modal expansion coefficients c . Furthermore, as shown in [59,60] we have c where t 0 is the initial time mentioned before and where c α,j (t 0 ) is an operator associated with the total field acting at time t 0 . In the rest of this work we will write Eq. 3.2 in the more compact form Here it is explicitly obtained as an inverse Fourier integral [60]: which vanishes for τ < 0 and involves the knowledge of the usual Green dyadic function [2,60,62] G v (x, x ′ , ω) in vacuum, i.e., solution of the equation Using the definition for the Fourier transform t)e +iωτ (and an equivalent formula for the dipole density P(x, ω)) it is actually easier to write Eq. 3.1 as . Now in the present problem involving a polarizable medium with complex permittiv- where P (0) (x, ω) is the fluctuating dipole density introduced by Gruner and Welsch [8][9][10] and given by [60]: ω (x, 0) a fluctuating dipolar term associated with the dielectric medium (we have f Therefore, it is specially convenient to introduce the total Green dyadic function [2,60,62] G(x, x ′ , ω) in the polarizable medium solution of the equation We thus rewrite Eq. 3.6 as Maxwell's equations in the dielectric medium in absence of P eff. (x ′ , ω) and we have the integral relation [60] where in the last equality we used the important Lippman-Schwinger integral relations Therefore, we see that the electric field E (0) (x, ω) is completely determined by the knowledge of the vacuum electric field E (v) (x, ω) through a linear equation. Clearly E (0) (x, ω) can be also rewritten as where the functions E (0) α,j (x) are the classical electric fields which are solutions of the scattering problem of a plane wave E (v) α,j (x) with pulsation ω α by the polarizable medium. Therefore, we have We point out that contrarily to E do not constitute in general an orthogonal mode basis but this does not prevent us to use it for expanding the field E (0) (x, ω) as in Eq. 3.13. Before leaving this subsection it is important to give the field expression in the time domain which reads [60]: where like for the vacuum case the causal propagator ∆(τ, x, x ′ ) is given by the inverse Fourier transform which vanishes for τ < 0 [60].

IV. THE SIPE APPROACH AND THE WIGNER-WEISSKOPF APPROXIMATION
In order to solve the system of dynamical equations we will first use the Wigner-Weisskopf approximation [63] as analyzed by J. Sipe [64]. For this purpose we first consider the total Hamiltonian H(t) written as We now expand the quantum states |Ψ(t) into the mode basis associated with the non interacting part H (0) (t). We will consider the problem of spontaneous emission of light by a two-level system and following Sipe [64] we postulate that the quantum state at time t 0 = 0 reads where 0 1 means that no electron is in the lower energy state E 1 and 1 2 implies that there is an electron in the upper energy state E 2 (in the following the notation |0 1 , 1 2 , 0 m , 0 P ; t 0 and other similar ones mean that the vector is actually an eigenstate of a complete set of operator at time t 0 = 0). Since N (t) = N (0) we see that for t ≥ 0 the system is in an eigenstate of N (t) corresponding to the eigenvalue N = +1. Similarly 0 m and 0 P mean that there is no polarization and photon excitation in the system at the initial time. At time t > t 0 the atomic system evolves to its ground state and we approximately have: with S(0) = 1 and b α,j (0) = 0, b ω (x, 0) = 0. We also define a second possible state which we call the ground state of the complete system and that we write neglecting dressing and following Sipe [64] as: which characterizes a system with fundamental energy E 1 (i.e., neglecting dressing). To solve Eq. 2.14 coupled to Maxwell's equations we here consider the matrix elements G(0)|A(t)|I(0) = G(t)|A(0)|I(t) associated with the operator A(t) is expressed in the Heisenberg picture. We also remind that if U (t, 0) denotes the unitary evolution operator associated with the full Hamiltonian H(t), we have |I(t) = U (t, 0)|I(0) , |G(t) = U (t, 0)|G(0) and, therefore we have A(t) = U −1 (t, 0)A(0)U (t, 0). Using Eq. 2.14 we get In order to solve this equation we should evaluate the different matrix elements involved. First, we have where we have used Eq. 4.4. Second, the matrix element can be rewritten as −e i E 1 t G(0)|(Ē(0) +P Ψ (0))|I(t) ≃ G(0)|(Ē(t)+P Ψ (t))|I(0) since the quantum state |G(0) corresponds to the eigenvalue σ z = −1. Regrouping the terms we finally get the dynamical equation Now, to evaluate G(0)|Ē(t)|I(0) we need to insert the electric field solution of Maxwell's equations evaluated at the dipole position. As discussed in the previous subsection the total field is given by Eq. 3.15. By using the form for G(0)| and |I(t) it is not difficult [65] to show that we have Therefore, the only contribution to the electric field matrix element comes out from the field generated by the atomic dipole itself with: where we have used Regrouping all these expressions, Eq. 4.8 finally reads: where ω 0 = ω 21 − |µ 1,2 | 2 ∆(0) is a modified pulsation, due to the dipole G(0)|P Ψ (t)|I(0) . This result is central for the present analysis since it shows that the vacuum photon field and the material fluctuating currents are not playing an effective role in the dynamical equation. Therefore, the transition dynamics is driven by the self-interaction of the source electromagnetic field. We emphasize that in Eq. 4.11 the Green dyadic tensor ∆(τ, x, x 0 ) is actually a highly singular function both in the spatial (near x 0 ) and in the time domain (near τ = 0). Actually, we showed in [60], using the Laplace transform method, that ∆(τ, x, is a distribution that is regular in the time domain at τ = 0. Introducing this definition in Eq. 4.11 shows that the integral t 0 dτ Q(τ, x, x 0 )S(t − τ ) actually vanishes for t → 0 + and that the dipole correction to ω 21 seen in ω 0 compensates exactly for the additional dipole term coming from the equality Eq. 4.11 can be solved more easily using the Laplace transform formalism applied to the field S(t). Using Eq. 4.11 one obtains: Therefore the following solution holds: (4.13) The calculation of this integral is given in the appendix A using the Wigner-Weisskopf approach. In short, the idea is to assume for long time t ≫ 0 the exponential decay law S(t) = S(0)e −iω0t , whereω 0 is a complex frequency defined asω 0 = ω 0 − iΓ/2 + δ, where Γ ≥ 0 and δ are real numbers (we have also S(0) = 1 by definition of the operator and quantum state). Now if we suppose that to a good approximation (called the polar approximation) we havẽ which indeed justifies the decay law (see Appendix A). This allows us to write the decay rate and the lamb shift where we used the Kramers-Kronig relation at the end. We point out that from the symmetry Im The first integral term is the correct Lamb shift obtained in the rotating wave approximation and the Wigner-Weisskopf theory. The second integral has actually the wrong sign (see the analysis of the problem in [34]) and is clearly non resonant. Only a more precise theory going beyond the Wigner-Weisskopf polar approximation would justify this erroneous (small) value and this will not be considered here.
Moreover, the total field at point x would be in principle calculated using: However, the polar approximation allows us to define in a simpler form this matrix element associated with the spontaneously emitted photon electric field: The explicit calculation of E(x, t) can only be done with the same conditions as used for evaluating S(t) in Eq. 4.13. For Eq. 4.13 it was necessary to suppose t ≫ t 0 . This is however not sufficient here and we will consider the far-field far away from the source region where the propagator reads asymptotically : is a smoothly varying form factor characterizing the emission profile (for the permittivity we here suppose a background but this could be vacuum. Using this robust far-field approximation we have where the Heaviside function is reminiscent of the causal nature of the single photon emission (since the photon emission starts at t 0 = 0 no light exists outside the future-oriented light cone with apex located at x 0 , t 0 ) and is here justified by the nature of the Bromwich integral.

A. Local density of states and polaritonic wave functions
Some important remarks should be done here concerning the above derivation and its meaning. First, as observed by Sipe, E(x, t) is defining together with B(x, t) := G(0)|B(x, t)|I(0) a wave function for the single emitted photon. More precisely, starting from Maxwell's quantum equations for operators E(x, t) and B(x, t) we can define some Maxwell's equations for the complex fields E(x, t) and B(x, t) which reads: is defined as G(0)|D(x, t)|I(0) and involves the complex polarization field P(x, t) := G(0)|P(x, t)|I(0) . From [59,60] and Eq. 4.9 we get where χ(x, τ ) is the local linear susceptibility of the inhomogeneous medium defined in [59]. From Eq. 4.9 we have and S(t) = e −iω0t is the complex valued dipole amplitude given in Eq. 4.13. In other words, if we insert the source term given by Eq. 5.3 in the Maxwell equations Eq. 5.1 we can solve the problem directly using the propagator ∆(τ, x, x ′ ) defined previously for the inhomogeneous dielectric problem. This solution is essentially classical and will automatically lead to Eq. 4.20 in the far-field of the quasi point-like dipole µ 1,2 S(t) associated with the polarization density given by Eq. 5.3. The methods is associated with the first quantization approach of photon proposed by Sipe [64] and Bialinicky-Birula [66] in which D and B define a wave function for the emitted photon. A second, remark connected to the first one deals with the energy conservation and the meaning of Γ in Eq. 4.15. Indeed, the structure of this mathematical expression for Γ is reminiscent of a classical calculation for the power radiated by an oscillating point-like dipole [2]. This is clear since we can write Γ as with µ 1,2 = |µ 1,2 |n. This is rigorously equivalent to the classical formula obtained for the power P 0 of a radiating dipole at the pulsation ω 0 which reads The last expression is identical to Eq. 5.4 if we identify the radiative power P 0 and the rate ω 0 Γ. In order to give a justification to this identification we start from Eq. 5.1 and we obtain a complexified version of the Poynting theorem which reads: with J = ∂ t P the complex dipolar current associated with Eq. 5.2: Moreover, by integration over the volume we can define the dissipated power inside the particle as W Ψ (t) : t)]. Using Eq. 5.3 we get: This is with a minus sign the total radiated power P 0 discussed previously but weighted by the exponential decay factor e −Γt . Integrating −W Ψ (t) over time we get the total energy emitted by the dipole from the initial time which approaches ω 0 if t → +∞. In evaluating W Ψ (t) we used the fact that while G(x 0 , x 0 , ω 0 + i0 + ) is a badly mathematically defined quantity this is not so for Im[G(x 0 , x 0 , ω 0 + i0 + )] which can be easily obtained by contour integration in the complex plane (see appendix B) and leads in the homogeneous surrounding medium case to: Re[n 0 (ω 0 )]I. (5.10) where n 0 (ω 0 ) is the surrounding medium optical index. This allows us to define the LDOS in the general case and to justify directly Eq. 5.5. Furthermore, we have also where we have used a standard separation [62] of the Green tensor as G = G ref + G 0 , where G 0 is a contribution of the bulk medium of permittivityε(ω 0 ) = n 2 0 (ω 0 ) and G ref is an additional contribution originating from the inhomogeneities and various interfaces present in the system.
It is also important for the present study to make a comment concerning the theory of intensity measurement proposed by Glauber [34,68,69]. We remind that following the theory of Glauber the photon detection rate I(x, t) at point x and time t should generally be expressed as a convolution between the temporal response of the detector M (τ ) and the first-order correlation function of the electric field γ(x, t, τ ) = E (−) (x, t)E (+) (x, t − τ ) , i.e., (5.12) In the formula for γ(x, t, τ ), E (+) (x, t) and E (−) (x, t) are respectively the positive and negative frequency operator parts of the electric field containing as usual only annihilation and creation operator for the photon field. In the broadband detector limit usually considered the formula simplifies and we get I( which is the standard formula of Glauber [68]. We stress that in order to derive Eq. 5.12 a dipolar coupling with the detector was taken into account using a Hamiltonian interaction of the usual form H int = −p · E where p is a dipole operator for the detector. However, from the point of view of the present dual formalism the principal field to be coupled to the detector is not the electric field E but the displacement D = E + P with P the total dipole density of the medium (which in the interaction picture does not include the detector dipole contribution). As we showed in [59,60] this displacement field is properly quantized by introducing a plane wave expansion with the general form (compare with Eq. 3.2): where c α,j (t) and c † α,j (t) are respectively the annihilation and creation operators associated with the photons in this dual formalism obeying usual commutation relations for bosons [59,60]. Comparing with Eq. 3.2 for E (v) (x, t) we see that the time dependency of c α,j (t) is not in general harmonic due to the coupling with the dipolar sources present [59,60]. Using this description the positive frequency part D (+) (x, t) of the displacement field operator D(x, t) is clearly defined as and D (−) (x, t) = D (+) † (x, t) as usual. In the dual formalism the interaction Hamiltonian for the detection process actually reads H new int = −p · D and therefore the single photon rate correlation function is still given by Eq. 5.12 with the correlation function now replaced by γ(x, t, τ ) = D (−) (x, t)D (+) (x, t − τ ) . For all practical needs in the laboratory the use of D instead of E will not change anything since most single photon detectors are located in the far-field region, i.e., generally speaking in the air with D ≃ E. Still in the near-field regime the new formalism is in principle more powerful since it includes from the ground the lossy and dispersive dielectric environment.
Furthermore, for the single photon process considered before we can write where we used the fact that the ground state is supposed here to be approximately the same at time t = 0 and time t. Therefore, the recorded single photon intensity in the far-field by an idealized broadband detector requires only the knowledge of I(x, t) ∝ D (−) (x, t)D (+) (x, t) = | G(0)|D (+) (x, t)|I(0) | 2 which from the previous analysis and Eq. 4.20 is given by |E(x, t)| 2 . In agreement with Sipe's analysis in vacuum [64] we thus finally obtain a description of photon detection in terms of a single photon wave function in presence of a dielectric environment.

B. The continuous regime and the optical Bloch equations
Before to conclude it is here important to study the continuous excitation regime when a laser mode interacts with the two-level atom considered previously. In order to find the optical Bloch equations in this regime we go back to Eqs. 2.13, 2.14 and study the dynamics of σ(t) , and σ z (t) where the average is taken on an arbitrary initial state for the two-level atom. More specifically, we are interested in the evolution of σ(t) , and σ z (t) with time under the influence of a quasi-classical electromagnetic wave characterized by an harmonic electric field E (0) (ω L is the pulsation of the quasi-classical laser field). This electric field is a solution of the homogeneous classical Maxwell equations in presence of the dielectric medium. Therefore, from the point of view of the QED approach considered here it will be necessary to include a contribution of the vacuum electric operator E (0) (x, t) in Eqs. 2.13, 2.14. It is also clear that without the introduction of the operator E (0) (x, t) to preserve unitarity it would be also impossible to describe the excitation by an incident laser mode considered as a pure photonic state. With our formalism it is thus possible to describe the interaction process in complete analogy with what is done in the literature for an atom excited by a laser beam in vacuum (i.e. without a dielectric lossy and dispersive surrounding). Now, lets start with Eq. 2.14 and consider the average Here comes a difficulty because to solve this equation one must specify the operator ordering in (Ē(t) + P Ψ (t))σ z (t) . This is a central issue which is well documented in the case of an atom in vacuum [34]. The usual trick, that we should apply here as well, is to take a normal ordering in which positive frequency part of the electric field are positioned to the right of σ z while the negative frequency part of the electric field operator is positioned to the left of σ z . This is allowed because atomic and field operators defined at the same time commute. If we can do that we will remove the contributions from E (0) (x, t) associated with vacuum fluctuations and only study the effect of the (classical) external interacting field and of the radiation-reaction. We mention, that there are mathematical subtleties in the definition of positive and negative frequency parts and it exist actually two ways to define it which are not rigorously equivalent.
On the one side we could be tempted to consider a Fourier transformÃ(ω) of any operator A(t) and thus define the positive frequency part as . This way of defining A (±) (t) is actually correct if we have no interaction. However, the canonical approach [69,70] is to use the separation between annihilation and creation operators for the fields and to define the positive frequency part by using only annihilation operators (respectively the negative frequency part is defined using only the creation operators for the fields). This is clearly the definition used for D (+) (x, t) in Eq. 5.14 but now this should be generalized for taking into account the relation E(t)+P Ψ (t) =D(t)−P diel. (t). We give in Appendix C a detailed discussion of this important point in the present dual formalism. We now write:  Ψ (t)) † . We have also for the electric field operators [59]). As shown in Appendix D a rigorous application of the operator L (±) t is in general difficult and an exact result is only obtained in the vacuum. Fortunately, here we are interested in the dynamics at long time and we can use the approximation (see Eqs. 3.11, and 3.12): together with the relation E (0,−) (x, t) = (E (0,+) (x, t)) † , ∆ (−) (τ, x, x ′ ) = (∆ (+) (τ, x, x ′ )) * . The next step is to remove the excitation from the initial state at time t 0 → −∞. For this we use [71] for every operator the unitary transformation A new (t) = T A(t)T −1 where T is the displacement operator defined α,j (t) = f α,j (t 0 )e −iωα(t−t0) the modal coefficients in the expansion of the free laser field. More precisely, in analogy with Eq. 3.13 the laser field is written E (5.20) and where E (0) α,j (x) is defined by Eq. 3.14 from the plane wave modes E (v) α,j (x) (for the particular example used in this section the laser is monochromatic so that we have necessarily ω α = ω L ). Moreover, within this formalism the application of T on the initial coherent state |L, t 0 leads to the photon vacuum: T |L, t 0 = |0, t 0 and since T acts only on the operator field If we suppose that the initial quantum state is |L, t 0 ⊗ |atom, t 0 (where |atom, t 0 is any coherent superposition of the atomic states |1, t 0 and |2, t 0 ) then the coherent laser field can be removed from the quantum state (which now reads |0, t 0 ⊗ |atom, t 0 ) and we should in turn add a classical laser field in the dynamical Eq. 5.16. The previous analysis therefore generalizes the usual method for removing coherent states. However, here the trick is now valid in presence of lossy and dispersive media.
The rest of the derivation is more conventional. We write σ(t) = S(t)e −iωLt and using the rotating wave approximation we neglect contributions from E (0,−) L and σ † . We get after some calculations: −∞ e i(ωL−ω)τ S(t) . (5.22) In the last equality we used the Markovian approximation [34,70] σ z (t)S(t − τ ) ≃ σ z (t)S(t) = − S(t) (since σ z (t)σ(t) = −σ(t) by definition). In the long time limit with t 0 → −∞ we have is the quantum rate at the laser frequency (compare with Eq. 4.15) and is the new Lamb shift. We emphasize that for most applications the difference between Γ and Γ ′ , δ and δ ′ can be neglected. Therefore, within the general approach considered we obtained the first Bloch equation which is now valid in presence of lossy and dispersive linear media. It is possible using the same procedure to deduce the second optical Bloch equation which reads within tthe same approximations: (5.26) Therefore, we can by using the generalized Huttner Barnett model justify the use of optical Bloch equations which were often introduced with the more phenomenological Langevin noise approach [21][22][23][24][25][26][27][28][29].

VI. FINAL REMARKS AND CONCLUSION
To conclude, we provided a description of spontaneous emission for a fluorescent two level atom using the generalized Hutner-Barnett approach given in [59,60]. We showed that within this Hamiltonian description it is clearly possible to analyze rigorously spontaneous emission in a lossy and dispersive inhomogeneous dielectric environment. Importantly, we showed that while our description used the complete electromagnetic field including photon vacuum fluctuations and Langevin's noise current associated with the dielectric environment the spontaneous emission process can be understood as resulting from a self-coupling of the fluorescent dipole. In this perspective the LDOS appears thus as a consequence of the classical radiation reaction (in agreement with semi-classical approaches neglecting the quantization of losses [2,3]). On the one side, spontaneous emission is thus interpreted rather classically as a radiation reaction due to the Lorentz force [34]. However, on the other side the full unitarity of quantum mechanics is respected in our formalism in order to preserve the canonical equaltime commutation relations. Therefore, as already discussed by Milonni in the context of photon/atoms coupling in vacuum [34] the description is not univocal and depends on the order we introduce operators in the dynamical equations. Here, the choice was done in order to favor the classical radiation force interpretation but other choices are clearly possible and all of them are equivalent. Furthermore, this fact can be seen as a direct consequence of the preservation of unitarity in our description. Without the inclusion of both photon vacuum and fluctuating currents the alternative representations would not exist and the full unitarity would be broken. The present work justifies semi-classical results [2,3,72] and alternative quantum approaches based on the Langevin's noise method [21][22][23][24][25][26][27][28][29] which neglected the role of photon vacuum. We think that the present work will motivate further studies in order to analyze other regimes of coupling between emitters and dielectric media and will impact our description of quantum polaritonic an plasmonic physics in the quantum regime (e.g., with near-field optical microscopes involving singlephoton emitters [6,7,[73][74][75].

Appendix A: The Wigner-Weisskopf approximation
The evaluation of the integral Eq. 4.13 is in general difficult and we will here use the method proposed by Wigner and Weisskopf [63]. For this we introduce the is regular at time τ = 0. The ansatz made by Wigner and Weisskopf is to suppose for time t → +∞ an exponential decay S(t) = S(0)e −iω0t whereω 0 is a complex frequency which is defined asω 0 = ω 0 − iΓ/2 + δ with Γ ≥ 0 and δ real numbers. This leads to the relation in the limit γ → 0 + . This integral is calculated by contour integration in the complex plane and taking into account that the poles of B(ω) are all located in the lower frequency half-plane (the derivation is identical to the one forε(ω) as shown in [60]). Using the residue theorem we get where the sum is taken over the residues or poles Ω m of B(ω). This equality cannot be valid at every time since the left hand side is independent of t while the right hand side depends explicitly on t. Actually, the exponential decay law is only valid for long time, i.e., t ≫ 0. However, the equality between the right hand side and To conclude this evaluation we observe that in the limit where there is a continuum of poles Ω m we can with a good approximation [76] write Furthermore, in the limit where the poles are near the real axis and where Γ/2 ≫ −Ω ′′ m we have: Inserting Eqs. A6 and A7 in Eq. A5 leads to: However from the definition of B(z), with z a complex number, we have B(−z * ) = B * (z). Therefore, for ω 0 ± i0 + we have after separating the real from the imaginary part B ′ (ω 0 + i0 + ) = B ′ (ω 0 − i0 + ) and B ′′ (ω 0 + i0 + ) = −B ′′ (ω 0 − i0 + ). This implies (not the difference of sign with Eq. A7) and we get after comparison wih Eq. A8 which is the final result.
act on the material oscillator fields X ω (x, t) in order to separate the contribution containing annihilation operators f ω (x, t) from the contribution containing only creation operators f † ω (x, t). Moreover, as shown in [59,60] we have f ω (x, t) = i∂tXω (x,t)+ωXω(x,t) √ 2 ω and Therefore if we define and use the canonical commutation relations [59] we get: which is actually equivalent to Therefore, after integration in the complex plane we get and x, x ′ )) * . We emphasize that the longitudinal part G v,|| (x, x ′ , ω) doesn't contribute to ∆ (±) v (t − t ′ , x, x ′ ) as it can be shown directly. Furthermore, if we take the imaginary part of Eqs. D2 and D3 and use once again the separation we can directly demonstrate the rigorous equivalence in which the contribution of G v,|| (x, x ′ , ω) vanishes once again. We emphasize that, despite some similarities, this result is different from Eq. 5.19 (in particular due to the presence of the imaginary part and the Heaviside function θ(t − t ′ ) in Eq. D9). In order to evaluate asymptotically Eq. D9 for long time t − t ′ we use the definition ω 2 c 2 G v (x, x ′ , ω) = ∇ × ∇ × [G v (x, x ′ , ω)I] where G v (x, x ′ , ω) = e iω/c|x−x ′ | 4π|x−x ′ | is the standard scalar Green function of the Helmholtz equation in vacuum [60]. Now, we consider the integral We remind that following Feynman [78] we have also the integral δ + (x) := +∞ 0 dω π e −iωx = 1 iπ 1 x−i0 + = 1 iπ P.V.[1/x]) + δ(x) for x real. Therefore, we deduce We are interested in the regime τ → +∞ and from its definition δ + (τ + |x − x ′ |/c) can be neglected. We thus obtain Finally, we deduce the asymptotic result which is valid in the limit τ → +∞. In this formula the absence of a contribution like 0 −∞ dω[...] clearly means that such a term is non resonant. Moreover, from Eq. 3.11 and 3.12 we know that generally speaking E (0) (x, ω) and G(x, x ′ , ω), i.e., the relevant fields for the inhomogeneous problem, can be calculated by using Lippman-Schwinger integrals which depend on the knowledge of the dyadic Green function in vacuum G v (x, x ′ , ω). Using the previous asymptotic of L it is thus not difficult to deduce the generality of Eq. 5.19 in the long time limit τ → +∞.