Aurélien Drezet

Description of spontaneous photon emission and local density of states in presence of a lossy polaritonic inhomogenous medium

Keywords: numbers: 42.50.Ct, 41.20.Jb, 73.20.Mf

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,[START_REF] Novotny | Principles of Nano-Optics[END_REF] 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,[START_REF] Colas Des Francs | [END_REF][4][START_REF] Colas Des Francs | Relationship between scanning near-field optical images and local density of photonic states[END_REF][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][START_REF] Knöll | QED in dispersing and absorbing media[END_REF][START_REF] Matloob | [END_REF][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][START_REF] Ginzburg | Theoretical Physics and Astrophysics[END_REF][START_REF] Rytov | Theory of Electrical Fluctuations and Thermal Radiation[END_REF][START_REF] Callen | [END_REF][START_REF] Milonni | The Quantum Vacuum: An Introduction to Quantum Electrodynamics[END_REF][START_REF] Rosa | [END_REF][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][START_REF] Barnett | The quantum optics of dielectrics[END_REF] 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][START_REF] Barnett | The quantum optics of dielectrics[END_REF]. 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

L = B 2 -D 2 2 + F • ∇ × P - P 2 2 + L M + L Ψa (2.1)
where by definition B(x, t) = 1 c ∂ t F(x, t), and D(x, t) = ∇ × F(x, t) 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

L M = +∞ 0 dω (∂ t X ω ) 2 -ω 2 X 2 ω 2 . (2.2)
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]

L Ψa = i Ψ * a ∂ t Ψ a -V Ψ * a Ψ a - 2 2M ∇Ψ * a • ∇Ψ a (2.3)
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

{b n (t), b † m (t)} = δ n,m , (2.4)
and

{b n (t), b † m (t)} = 0 = {b n (t) † , b † m (t)}
. This is clearly consistent with the usual canonical quantization procedure since the Lagrangian function

L Ψa = d 3 xL Ψa (x, t) = i m b † m ḃm -m E m b † n , b † m implies the canonical momenta: Π bm = ∂L Ψa ∂ ḃm = i b n (t) † , Π b † m = ∂L Ψa ∂ ḃ † m = 0 (2.5)
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]):

P diel. (x, t) = +∞ 0 dω 2σ ω (x) π X ω (x, t). (2.6)
For the Schrodinger field we here model the dipole density fluid using the formula:

P Ψ (x, t) ≃ e d 3 ξΨ * a (ξ, t)ξΨ a (ξ, t)∆(x -x 0 ) = m,n b † m (t)b n (t)µ m,n ∆(x -x 0 ) (2.7)
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 ∆(xx 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∆(xx 0 ) = 1. ∆(xx 0 ) reduces to the Dirac distribution δ 3 (xx 0 ) in the point-like dipole limit. Still, in order to avoid singular divergences we will here keep ∆(xx 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.

∇ × B(x, t) = 1 c ∂ t D(x, t), ∇ • D(x, t) = 0 ∇ × E(x, t) = - 1 c ∂ t B(x, t), ∇ • B(x, t) = 0 (2.8)
where E = D -P Ψ -P diel. . Similarly, we deduce using the Euler-Lagrange method a dynamical equation for the dielectric field X ω (x, t):

∂ 2 t X ω (x, t) + ω 2 X ω (x, t) = 2σ ω (x) π E(x, t) (2.9)
with here E = D -P Ψ -P diel. . For the Schrodinger fields we must be more careful. Indeed by using the Euler-Lagrange classical method we should deduce the evolution equations

ḃm (t) = -i E m b m (t) + i Ē(t) • n µ * n,m b n (t) (2.10)
with Ē(t) = d 3 E(x, t)x∆(xx 0 ) ≃ E(x 0 , t) the averaged field applied on the atomic dipole. A similar equation is obtained for the complex conjugate (Hermitian) field b † m (t). From this equation we deduce

d( m b † m (t)bm(t)) dt
= 0 which means that the observable

N (t) = m b † m (t)b m (t)
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:

H(t) = d 3 x : B(x, t) 2 + E(x, t) 2 2 : +H M (t) + H Ψa (t) (2.11)
where

H M (t) = d 3 x +∞ 0 dω ωf † ω (x, t)f ω (x, t) and H Ψa (t) = m E m b † m (t)b m (t). H(t)
is clearly an integral of motion as it can be for instance proven by using the Poynting theorem discussed in [59] and which reads:

-∂ t ( B 2 +E 2 2 ) = ∇ • (cE × B) + J • E where J = ∂ t P.
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

E 2 -E 1 = ω 21 > 0. We have N (t) = b † 1 (t)b 1 (t) + b † 2 (t)b 2 (t) = N (t 0
) and we can therefore write the Schrodinger part of the Hamiltonian H Ψa (t) as

H Ψa (t) = E 1 b † 1 (t)b 1 (t) + E 2 b † 2 (t)b 2 (t) = E 1 + E 2 2 N (t) + ω 21 2 σ z (t) (2.12)
with by definition

σ z (t) = b † 2 (t)b 2 (t) -b † 1 (t)b 1 (t).
We get for the evolution equation of σ z (t):

σz (t) = -2i( Ē(t) + PΨ (t)) • µ 1,2 σ(t) + 2i( Ē(t) + PΨ (t)) • µ * 1,2 σ † (t) (2.13)
with by definition

σ(t) = b † 1 (t)b 2 (t) and σ † (t) = b † 2 (t)b 1 (t).
We emphasize that in this equation it is the field Ē(t) + PΨ (t) = D(t) -Pdiel. (t) which plays a central role, with PΨ (t) = (µ 1,2 σ(t)

+ µ * 1,2 σ † (t))∆(0) (since P Ψ (x, t) = (µ 1,2 σ(t)+µ * 1,2 σ † (t))∆(x-x 0 )
). 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 :

σ(t) = -iω 21 σ(t) - i( Ē(t) + PΨ (t)) • µ * 1,2 σ z (t) (2.14)
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

E(x, t) = E (v) (x, t) + t-t0 0 dτ d 3 x ′ ∆ v (τ, x, x ′ ) • P(x ′ , t -τ ). (3.1)
In this equation t 0 is an initial time (which can be sent to the remote past t 0 → -∞ if needed), E (v) (x, t) = D (v) (x, t) 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]:

E (v) (x, t) = α,j - ω α 2 c (v) α,j (t) kα × ǫα,j Φ α (x) + hcc. (3.2)
with Φ α (x) = e ikα•x / √ V (V is the infinite box volume of the Born von Karman modal expansion method) and where the modal expansion coefficients c

(v) α,j (t) = c (v)
α,j (0)e -iωαt (with ω α = c|k α |) satisfy the usual commutation relations for bosons (i.e. [c

(v) α,j (t), c (v) † α ′ ,j ′ (t)] = δ α,α ′ δ j,j ′ etc...)
. Furthermore, as shown in [59,60] we have c

(v) α,j (0) = c α,j (t 0 )e -iωαt0
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

E (v) (x, t) = α,j E (v) α,j (x)c (v) α,j (0)e -iωαt + hcc. (3.3)
which involves the transverse electric mode profiles E (v) α,j (x). We emphasize that the set of functions E (v) α,j (x) constitute an orthogonal mode basis satisfying the condition

d 3 xE (v) α,j (x) • E (v) * α ′ ,j ′ (x) = ωα 2 δ α,α ′ δ j,j ′ Moreover, the tensor ∆ v (τ, x, x ′ ) seen in Eq. 3.
1 is a retarded dyadic Green propagator which was written -

∂ 2 τ Uv (τ,x,x ′ ) c 2 in [60].
Here it is explicitly obtained as an inverse Fourier integral [60]:

∆ v (τ, x, x ′ ) = +∞ -∞ dω 2π e -iωτ ω 2 c 2 G v (x, x ′ , ω) (3.4)
which vanishes for τ < 0 and involves the knowledge of the usual Green dyadic function [START_REF] Novotny | Principles of Nano-Optics[END_REF]60,62] G v (x, x ′ , ω) in vacuum, i.e., solution of the equation

∇ × ∇ × G v (x, x ′ , ω) - ω 2 c 2 G v (x, x ′ , ω) = Iδ(x -x ′ ). (3.5)
Using the definition for the Fourier transform

E(x, ω) = +∞ -∞ dt 2π E(x, t)
e +iωτ (and an equivalent formula for the dipole density P(x, ω)) it is actually easier to write Eq. 3.1 as

E(x, ω) = E (v) (x, ω) + d 3 x ′ ω 2 c 2 G v (x, x ′ , ω • P(x ′ , ω), (3.6)
where

Ẽ(v) (x, ω) = α,j [E (v) α,j (x)c (v) α,j (0)δ(ω -ω α ) + E (v) * α,j (x)c (v) † α,j (0)δ(ω + ω α ). Now in the present problem involving a polarizable medium with complex permittiv- ity ε(x, ω) = ε ′ (x, ω) + i ε ′′ (x, ω) we have P(x, ω) = P Ψ (x, ω) + P (0) (x, ω) +( ε(x, ω) -1) E(x, ω) (3.7)
where P (0) (x, ω) is the fluctuating dipole density introduced by Gruner and Welsch [8][9][10] and given by [60]:

P (0) (x, ω) = ε′′ (x, ω) π f (0) ω (x, 0)Θ(ω) + ε′′ (x, -ω) π f (0) † -ω (x, 0)Θ(-ω). (3.8) with f (0)
ω (x, 0) a fluctuating dipolar term associated with the dielectric medium (we have f

(0) ω (x, t) = f (0) ω (x, 0)e -iωt and f (0) ω (x, 0) = f (0) ω (x, t 0 )e iωt0
). Therefore, it is specially convenient to introduce the total Green dyadic function [START_REF] Novotny | Principles of Nano-Optics[END_REF]60,62] G(x, x ′ , ω) in the polarizable medium solution of the equation

∇ × ∇ × G(x, x ′ , ω) - ω 2 c 2 ε(x, ω)G(x, x ′ , ω) = Iδ(x -x ′ ). (3.9)
We thus rewrite Eq. 3.6 as

E(x, ω) = E (0) (x, ω) + d 3 x ′ ω 2 c 2 G(x, x ′ , ω) • P eff. (x ′ , ω), (3.10)
where

P eff. (x, ω) = P Ψ (x, ω) + P (0) (x, ω).
Here the new electric field operator E (0) (x, ω) is a solution of Maxwell's equations in the dielectric medium in absence of P eff. (x ′ , ω) and we have the integral relation [60]

E (0) (x, ω) = E (v) (x, ω) + d 3 x ω 2 c 2 G v (x, x ′ , ω) •( ε(x, ω) -1) E (0) (x ′ , ω) = E (v) (x, ω) + d 3 x ω 2 c 2 G(x, x ′ , ω) •( ε(x, ω) -1) E (v) (x ′ , ω)(3.11)
where in the last equality we used the important Lippman-Schwinger integral relations

G(x, x ′ , ω) = G v (x, x ′ , ω) + d 3 x ω 2 c 2 G v (x, x ′ , ω) •( ε(x, ω) -1)G(x, x ′ , ω) = G v (x, x ′ , ω) + d 3 x ω 2 c 2 G(x, x ′ , ω) •( ε(x, ω) -1)G v (x, x ′ , ω).
(3.12)

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

E (0) (x, ω) = α,j [E (0) α,j (x)c (v) α,j (0)δ(ω -ω α ) +E (0) * α,j (x)c (v) † α,j (0)δ(ω + ω α ) (3.13)
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

E (0) α,j (x) = E (v) α,j (x) + d 3 x ω 2 α c 2 G v (x, x ′ , ω α ) •( ε(x, ω α ) -1) E (0) α,j (x ′ , ω) = E (v) α,j (x) + d 3 x ω 2 α c 2 G(x, x ′ , ω α ) •( ε(x, ω α ) -1) E (v) α,j (x ′ ). (3.14)
We point out that contrarily to

E (v) α,j (x) the fields E (0) α,j (x)
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]:

E(x, t) = E (0) (x, t) + t-t0 0 dτ d 3 x ′ ∆(τ, x, x ′ ) • P eff. (x ′ , t -τ ) (3.15)
where like for the vacuum case the causal propagator ∆(τ, x, x ′ ) is given by the inverse Fourier transform

∆(τ, x, x ′ ) = +∞ -∞ dω 2π e -iωτ ω 2 c 2 G(x, x ′ , ω), (3.16)
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

H(t) = d 3 x : B(x, t) 2 + D(x, t) 2 2 : +H M (t) + H Ψa (t) -d 3 xD(x, t) • P(x, t) (4.1)
H(t) can be formally separated into a non-interacting part : +H M (t) + H Ψa (t) and a coupling term

H (0) = d 3 x : B(x,t) 2 +D(x,t)
H (I) (t) = -d 3 xD(x, t) • P(x, t).
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

|I(0) = |0 1 , 1 2 , 0 m , 0 P ; t 0 (4.2)
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:

|I(t) ≃ S(t)e -i E 1 t |0 1 , 1 2 , 0 m , 0 P ; t 0 + d 3 x +∞ 0 dωb ω (x, t)f ω (x, 0) † |1 1 , 0 2 , 0 m , 0 P ; t 0 + α,j b α,j (t)c α,j (0) † |1 1 , 0 2 , 0 m , 0 P ; t 0 + ... (4.3)
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:

|G(t) ≃ e -i E 1 t |1 1 , 0 2 , 0 m , 0 P ; t 0 = e -i E 1 t |G(0) (4.4)
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

d dt G(0)|σ(t)|I(0) = -iω 21 G(0)|σ(t)|I(0) - iµ * 1,2 • G(0)|σ z (t)( Ē(t) + PΨ (t))|I(0) (4.5)
In order to solve this equation we should evaluate the different matrix elements involved. First, we have

G(0)|σ(t)|I(0) = G(t)|σ(0)|I(t) ≃ e i E 1 t G(0)|σ(0)|I(t) = S(t) (4.6)
where we have used Eq. 4.4. Second, the matrix element

G(0)|σ z (t)( Ē(t) + PΨ (t))|I(0) = G(t)|σ z (0)( Ē(0) + PΨ (0))|I(t) ≃ e i E 1 t G(0)|σ z (0)( Ē(0) + PΨ (0))|I(t) (4.7)
can be rewritten as

-e i E 1 t G(0)|( Ē(0) + PΨ (0))|I(t) ≃ G(0)|( Ē(
+ iµ * 1,2 • G(0)|( Ē(t) + PΨ (t))|I(0) . (4.8)
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 [START_REF]The operators Ē(0) (t) and P(0) (t) can be directly expressed as function of operators cα,j (0)[END_REF] to show that we have

G(0)| Ē(0) (t)|I(0) = 0, G(0)| P(0) (t)|I(0) = 0.
Therefore, the only contribution to the electric field matrix element comes out from the field generated by the atomic dipole itself with:

G(0)| PΨ (t)|I(0) = (µ 1,2 G(0)|σ(t)|I(0) +µ * 1,2 G(0)|σ † (t)|I(0) )∆(0) ≃ µ 1,2 S(t)∆(0). (4.9)
where we have used

G(0)|σ † (t)|I(0) ≃ e i E 1 t G(0)|σ † (0)|I(t) = 0. (4.10)
Regrouping all these expressions, Eq. 4.8 finally reads:

d dt S(t) = -iω 0 S(t) + t 0 dτ iµ * 1,2 • ∆(τ, x 0 , x 0 ) • µ 1,2 •S(t -τ ) (4.11)
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,

x 0 ) = [Q(τ, x, x 0 ) -Iδ 3 (x -x 0 )δ(τ )]Θ(τ )
where Q(τ, x, x 0 ) 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

t 0 dτ ∆(τ, x 0 , x 0 ) • µ 1,2 S(t -τ ) = t 0 dτ Q(τ, x 0 , x 0 ) • µ 1,2 S(t -τ ) -G(0)| PΨ (t)|I(0) .
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:

pS(p) -S(0) = -iω 0 S(p) - iµ * 1,2 • p 2 c 2 G(x 0 , x 0 , ip) • µ 1,2 S(p) (4.12)
Therefore the following solution holds:

S(t) = γ+i∞ γ-i∞ idp 2π
e pt S(0)

p + iω 0 + iµ * 1,2 • p 2 c 2 G(x0,x0,ip)•µ 1,2 . (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 have

ω0 = ω 0 - µ * 1,2 • ω 2 0 c 2 G(x 0 , x 0 , ω 0 + i0 + ) • µ 1,2 (4.14)
which indeed justifies the decay law (see Appendix A). This allows us to write the decay rate

Γ = 2Im[ µ * 1,2 • ω 2 0 c 2 G(x 0 , x 0 , ω 0 + i0 + ) • µ 1,2 ] (4.15)
and the lamb shift

-δ = Re[ µ * 1,2 • ω 2 0 c 2 G(x 0 , x 0 , ω 0 + i0 + ) • µ 1,2 ] = +∞ -∞ dω π P.V.[ Im[ µ * 1,2 • ω 2 c 2 G(x0,x0,ω)•µ 1,2 ] ω -ω 0 ] (4.16)
where we used the Kramers-Kronig relation at the end. We point out that from the symmetry Im[G(x, x ′ , -ω)] = -Im[G(x, x ′ , ω)] we have:

-δ = +∞ 0 dω π P.V.[ Im[ µ * 1,2 • ω 2 c 2 G(x0,x0,ω)•µ 1,2 ] ω -ω 0 ] + +∞ 0 dω π Im[ µ * 1,2 • ω 2 c 2 G(x0,x0,ω)•µ 1,2 ] ω + ω 0 ]. (4.17)
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 [START_REF] Milonni | The Quantum Vacuum: An Introduction to Quantum Electrodynamics[END_REF]) 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:

E(x, t) := G(0)|E(x, t)|I(0) = +∞ -∞ dω 2πi ω 2 c 2 G(x, x 0 , ω) • µ 1,2 e -iωt S(0) ω 0 -ω -i0 + - ω 2 0 c 2 µ * 1,2 • G(x, x 0 , ω) • µ 1,2 (4.18)
However, the polar approximation allows us to define in a simpler form this matrix element associated with the spontaneously emitted photon electric field:

E(x, t) := G(0)|E(x, t)|I(0) = +∞ -∞ dω 2πi ω 2 c 2 G(x, x 0 , ω) ω0 -ω -i0 + • µ 1,2 e -iωt S(0) (4.19)
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 :

G(x, x 0 , ω 0 + i0 + ) ≃ F(x, x 0 , ω 0 + i0 + )e iω0 √ ε(ω0)R/c
where R = |xx 0 | ≫ c/ω 0 and where F(x, x 0 , ω 0 + i0 + ) 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

E(x, t) ≃ ω 2 0 c 2 F(x, x 0 , ω 0 + i0 + ) • µ 1,2 e -iω0(t- √ ε(ω0)R/c) S(0)Θ(t -ε(ω 0 )R/c), (4.20)
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.

V. DISCUSSIONS

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:

∇ × B(x, t) = 1 c ∂ t D(x, t), ∇ • D(x, t) = 0 ∇ × E(x, t) = - 1 c ∂ t B(x, t), ∇ • B(x, t) = 0 (5.1)
In these equations D(x, t) 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

P(x, t) = P eff. (x, t) + t 0 χ(x, τ )dτ E(x, t -τ ), (5.2)
where χ(x, τ ) is the local linear susceptibility of the inhomogeneous medium defined in [59]. From Eq. 4.9 we have

P eff. (x, t) ≃ µ 1,2 S(t)∆(x -x 0 ).
(5.3) 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 [START_REF] Bialynicki-Birula | [END_REF] 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 [START_REF] Novotny | Principles of Nano-Optics[END_REF]. This is clear since we can write Γ as

Γ = π 3 ω 0 |µ 1,2 | 2 ρ LDOS (x 0 ) (5.4)
where

ρ LDOS (x 0 ) = 6ω 0 πc 2 Im[n * • G(x 0 , x 0 , ω 0 + i0 + ) • n (5.5)
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

P 0 = π 3 ω 2 0 |µ 1,2 | 2 ρ LDOS (x 0 ). (5.6)
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:

-∂ t (|B| 2 + |E| 2 ) = 2c∇ • (Re[E × B * ]) + 2Re[J • E * ] , (5.7)
with J = ∂ t P the complex dipolar current associated with Eq. 5.2:

J (x, t) = ∂ t P eff. (x, t) + t 0 χ(x, τ )dτ ∂ t E(x, t -τ ) +χ(x, t)dτ E(x, 0) (5.8)
Moreover, by integration over the volume we can define the dissipated power inside the particle as W Ψ (t) :=

d 3 x2Re[∂ t P * eff. (x, t) • E(x, t)].
Using Eq. 5.3 we get:

W Ψ (t) ≃ -2Im[µ * 1,2 • ω 3 0 c 2 G(x 0 , x 0 , ω 0 + i0 + ) •µ 1,2 ]e -Γt = -ω 0 Γe -Γt . (5.9)
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

t 0 = 0 to time t: δE = - t 0 dt ′ W Ψ (t ′ ) = ω 0 (1 -e -Γt )
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:

Im[G(x 0 , x 0 , ω 0 + i0 + )] = ω 0 6πc
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

ρ LDOS (x 0 ) = ω 2 0 π 2 c 3 Re[n 0 (ω 0 )] + 6ω 0 πc 2 Im[n * • G ref (x 0 , x 0 , ω 0 + i0 + ) • n] (5.11)
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 [START_REF] Milonni | The Quantum Vacuum: An Introduction to Quantum Electrodynamics[END_REF][START_REF] Glauber | [END_REF]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.,

I(x, t) = 2Re[ +∞ 0 M (τ )γ(x, t, τ )].
(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(x, t) ∝ E (-) (x, t)E (+) (x, t) which is the standard formula of Glauber [START_REF] Glauber | [END_REF]. 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):

D(x, t) = α,j - ω α 2 c α,j (t) kα × ǫα,j Φ α (x) + hcc.
(5.13)

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

D (+) (x, t) = α,j - ω α 2 c α,j (t) kα × ǫα,j Φ α (x) (5.14)
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

G(0)|D(x, t)|I(0) ≃ G(0)|D (+) (x, t)|I(0) (5.15)
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 (x, t) = E 0 L (x)e -ωLt + cc. (ω 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

d dt σ(t) = -iω 21 σ(t) - i ( Ē(t) + PΨ (t))σ z (t) • µ * 1,2
(5.16)

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 [START_REF] Milonni | The Quantum Vacuum: An Introduction to Quantum Electrodynamics[END_REF]. 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

A (+) (t) = +∞ 0 dω Ã(ω)e -iωt (similarly A (-) (t) = 0 -∞ dω Ã(ω)e -iωt
). 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 Ē(t)+ PΨ (t) = D(t)-Pdiel. (t). We give in Appendix C a detailed discussion of this important point in the present dual formalism. We now write:

( Ē(t) + PΨ (t))σ z (t) = [( Ē(+) (t) + P(+) Ψ (t))σ z (t) + σ z (t)( Ē(-) (t) + P(-) Ψ (t))] (5.17)
with for the dipole field operator P(+) Ψ (t) = µ 1,2 σ(t) and P(-) Ψ (t) = ( P(+) Ψ (t)) † . We have also for the electric field operators (5.18) where the field [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):

E (±) (x, t) = E (0,±) (x, t) + t-t0 0 dτ d 3 x ′ ∆ (±) (τ, x, x ′ ) • P eff. (x ′ , t -τ ).
E (0,±) (x, t) = L (±) t [E (0,±) (x, t)] and the dyadic ∆ (±) (τ, x, x ′ ) = L (±) τ [∆(τ, x, x ′ )] are defined by applying the operator L (±) t = 1 2 [1 ± i∂t c √ -∇ 2 ] (see Ap- pendix C and
E (0,+) (x, t) ≃ +∞ 0 dω E (0) (x, ω)e -iωt ∆ (+) (τ, x, x ′ ) ≃ +∞ 0 dω 2π e -iωτ ω 2 c 2 G(x, x ′ , ω) (5.19)
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 [START_REF] Cohen-Tannoudji | Photons and Atoms: an Introduction to Quantum electrodynamics[END_REF] for every operator the unitary transformation A new (t) = T A(t)T -1 where T is the displacement operator defined

as T = α,j e [f (v) * α,j (t)c (v) α,j (t)-f (v) α,j (t)c (v) † α,j (t)] with f (v)
α,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 (0) L (x, t) = +∞ -∞ dω E (0) L (x, ω)e -iωt with E (0) L (x, ω) = α,j [E (0) α,j (x)f (v) α,j (t 0 )e iωαt0 δ(ω -ω α ) + cc.
(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

E (0) (x, t) we deduce T E (0,+) (x, t)T -1 = E (0,+) (x, t) + E (0,+) L (x, t), T E (0,-) (x, t)T -1 = E (0,-) (x, t) + E (0,-) L (x, t).
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:

d dt S(t) = -i(ω 0 -ω L ) S(t) - i Ē(0) L • µ * 1,2 σ z (t) + N (t) (5.21)
where

N (t) = -i +∞ 0 dω 2π ω 2 c 2 µ * 1,2 • G(x 0 , x 0 , ω) • µ 1,2 × t-t0 -∞ e i(ωL-ω)τ σ z (t)S(t -τ ) ≃ +i +∞ 0 dω 2π ω 2 c 2 µ * 1,2 • G(x 0 , x 0 , ω) • µ 1,2 × t-t0
-∞ e i(ωL-ω)τ S(t) . (5.22) In the last equality we used the Markovian approximation [START_REF] Milonni | The Quantum Vacuum: An Introduction to Quantum Electrodynamics[END_REF]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

N (t) = i µ * 1,2 •G(x0,x0,ω)•µ 1,2 S(t) = (-Γ ′ /2 -iδ ′ ) where Γ ′ = 2Im[ µ * 1,2 • ω 2 0 c 2 G(x 0 , x 0 , ω L + i0 + ) • µ 1,2 ] (5.23)
is the quantum rate at the laser frequency (compare with Eq. 4.15) and

-δ ′ = +∞ -∞ dω π P.V.[ Im[ µ * 1,2 • ω 2 c 2 G(x0,x0,ω)•µ 1,2 ] ω -ω L ] (5.24)
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

d dt S(t) = -i(ω 0 -δ -ω L ) S(t) - i Ē(0) L • µ * 1,2 σ z (t) -Γ/2 S(t) (5.25)
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:

d dt σ z (t) = -2i Ē(0) * L • µ 1,2 S(t) + 2i Ē(0) L • µ * 1,2 S † (t) -Γ(1 + σ z (t) ).
(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 [START_REF] Novotny | Principles of Nano-Optics[END_REF][START_REF] Colas Des Francs | [END_REF]). On the one side, spontaneous emission is thus interpreted rather classically as a radiation reaction due to the Lorentz force [START_REF] Milonni | The Quantum Vacuum: An Introduction to Quantum Electrodynamics[END_REF]. 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 [START_REF] Milonni | The Quantum Vacuum: An Introduction to Quantum Electrodynamics[END_REF] 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 [START_REF] Novotny | Principles of Nano-Optics[END_REF][START_REF] Colas Des Francs | [END_REF][START_REF] Drezet | [END_REF] 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].

VII. ACKNOWLEDGMENTS

This work was supported by Agence Nationale de la Recherche (ANR), France, through the PLACORE (ANR-13-BS10-0007) grant. The author gratefully acknowledges discussions with S. Huant.

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 notation β(p) = -

µ * 1,2 • p 2 c 2 G(x,x ′ ,ip)•µ 1,2
and we rewrite Eq. 4.12 as Ṡ(t) = -iω 0 S(t) + i t 0 dτ β(τ )S(tτ ) = -iω 21 S(t)+i t 0 dτ q(τ )S(t-τ ) where q(τ ) = β(τ )+(ω 21ω 0 )δ(τ ) 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

/2 -δ = e iω0t +∞ -∞ dω 2π e -iωt) B(ω) i(ω 0 -ω) (A2)
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

iΓ/2 -δ = B(ω 0 ) - m res[B(Ω m )]e i(ω0-Ωm)t ω0 -Ω m (A3)
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 t 0 dτ β(τ )e iω0 τ is valid at every time t ≥ 0 and since t 0 dτ q(τ )e iω0τ vanishes for t = 0 we deduce

B(ω 0 ) = m res[B(Ω m )] ω0 -Ω m + ω 0 -ω 21 (A4) which allows us to rewrite iΓ/2 -δ = lim t→∞ m res[B(Ω m )] 1 -e i(ω0-Ωm)t ω0 -Ω m +ω 0 -ω 21 (A5)
To conclude this evaluation we observe that in the limit where there is a continuum of poles Ω m we can with a good approximation [START_REF] Holstein | Topics in Advanced Quantum Mechanics[END_REF] write

m res[B(Ω m )] e i(ω0-Ωm)t ω0 -Ω m ≃ 2πie -iω0t m res[B(Ω m )]δ(ω 0 -Ω ′ m ) (A6)
Furthermore, in the limit where the poles are near the real axis and where Γ/2 ≫ -Ω ′′ m we have:

B(ω 0 ) ≃ B(ω 0 -i0 + ) ≃ m res[B(Ω m )]P.V.[ 1 ω 0 -Ω ′ m ] +iπ m res[B(Ω m )]δ(ω 0 -Ω ′ m ) + ω 0 -ω 21 (A7)
Inserting Eqs. A6 and A7 in Eq. A5 leads to:

iΓ/2 -δ ≃ m res[B(Ω m )]P.V.[ 1 ω 0 -Ω ′ m ] -iπ m res[B(Ω m )]δ(ω 0 -Ω ′ m ) + ω 0 -ω 21 (A8)
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

B(ω 0 + i0 + ) ≃ m res[B(Ω m )]P.V.[ 1 ω 0 -Ω ′ m ] -iπ m res[B(Ω m )]δ(ω 0 -Ω ′ m ) + ω 0 -ω 21 (A9)
(not the difference of sign with Eq. A7) and we get after comparison wih Eq. A8

iΓ/2 -δ ≃ B(ω 0 + i0 + ) = µ * 1,2 • ω 2 0 c 2 G(x, x ′ , ω 0 + i0 + ) • µ 1,2 (A10)
which is the final result.

Appendix B: The radiated power and the LDOS of a quantum dipole

In order to calculate W Ψ (t) we have to consider the local field in the vicinity of the dipole. Using Eq. 4.19 we get for x → x 0 :

E(x, t) ≃ ω 2 0 c 2 G(x, x 0 , ω 0 + i0 + ) • µ 1,2 e -iω0t (B1)
Moreover, the normed and finite function ∆(x-x 0 ) characterizing the dipole polarization prevents us to obtain divergence in evaluating W Ψ (t). In the case of an homogeneous medium the Green tensor G(x, x 0 , ω 0 + i0 + ) is easily obtained as an integral over polarization states and wavevectors. We get (see [START_REF] Cohen-Tannoudji | Atom-Photon Interactions: Basic Processes and Applications[END_REF]):

G(x 0 , x 0 , ω 0 + i0 + ) = d 3 k (2π) 3 1 k 2 - ω 2 0 ε(ω0) c 2 (I - k ⊗ k ω 2 0 ε(ω0) c 2 ) = +∞ -∞ kdk (2π) 3 2πIf reg (k) k - ω0 √ ε(ω0) c (1 - 2 3 k 2 ω 2 0 ε(ω0) c 2 ) (B2)
In the last line we introduced a regularization function f reg (k) such as f reg (k) ≃ 1 for value near k ≃ 0 but f reg (k) ≃ 0 si |k| → +∞. This trick prevents the divergence and allows us to calculate the integral along a contour in the upper part of the complex plane. Using the residue theorem one get finally

G(x 0 , x 0 , ω 0 + i0 + ) = i ω 0 6πc ε(ω 0 )If reg ( ω 0 ε(ω 0 ) c ). (B3)
At the end we can simplify since f reg ( ω0 √ ε(ω0) c

) ≃ 1. We point out that this result depends on the assumption concerning the convergence of f reg (k) at infinity. It is well known that the Green dyadic propagator is badly defined at the spatial origin and this is clearly another manifestation of this fact. Still the result concerning the imaginary part of G(x 0 , x 0 , ω 0 + i0 + ) is very robust and will keep its absolute meaning since only the real part contains potential divergences. Finally, taking the imaginary part of Eq. B3 leads directly to the result Eq. 5.10. which defines the positive frequency part P (+) diel. (x, t) and the operator M (±) t (for the negative frequency part we simply use P

(-) diel. (x, t) = (P (+) diel. (x, t)) † ).
Finally, if we use the definition P Ψ (x, t) = (µ 1,2 σ(t) + µ * 1,2 σ † (t))∆(xx 0 ) we can obtain a separation between the operators σ(t) and σ † (t) and we thus define P In [60] we showed that we have the dyadic expansion G v (x, x ′ , ω) = G v,⊥ (x, x ′ , ω) + G v,|| (x, x ′ , ω) with for the transverse part

G v,⊥ (x, x ′ , ω) = α,j c 2 Φ α (x)Φ * α (x ′ )ǫ α,j ⊗ ǫα,j ω 2 α -(ω + i0 + ) 2 (D2)
and for the longitudinal part

G v,|| (x, x ′ , ω) = -c 2 α kα ⊗ kα Φ * α (x ′ )Φ α (x) (ω + i0 + ) 2 (D3)
A direct application of L (D6)

Therefore, after integration in the complex plane we get

∆ (+) v (t -t ′ , x, x ′ ) = i α,j E (v) α,j (x) ⊗ E (v) * α,j (x ′ )
×e -iωα(t-t ′ ) Θ(tt ′ ). (D7) and ∆ (-) v (tt ′ , x, x ′ ) = (∆ (+) v (tt ′ , x, x ′ )) * . We emphasize that the longitudinal part G v,|| (x, x ′ , ω) doesn't contribute to ∆ (±) v (tt ′ , 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

1 ω 2 α -(ω + i0 + ) 2 = 1 2ω α [ 1 ω α -ω -i0 + + 1 ω α + ω + i0 + ] = iπ 2ω α (δ(ω -ω α ) -δ(ω + ω α )) + 1 2ω α (P.V.[ 1 ω α -ω ]) -P.V.[ 1 ω α + ω ]) (D8)
we can directly demonstrate the rigorous equivalence

∆ (+) v (τ, x, x ′ ) = i +∞ 0 dω π ω 2 c 2 Imag[G v (x, x ′ , ω)] ×e -iωτ θ(t -t ′ ), (D9)
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 θ(tt ′ ) in Eq. D9). In order to evaluate asymptotically Eq. D9 for long time tt ′ we use the definition We remind that following Feynman [78] we have also the integral δ + (x) := Finally, we deduce the asymptotic result

∆ (+) v (τ, x, x ′ ) ≃ +∞ 0 dω 2π ω 2 c 2 G v (x,
x ′ , ω)e -iωτ , (D13) 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 (±) t [G v (x, x ′ , ω)e -iωτ ] it is thus not difficult to deduce the generality of Eq. 5.19 in the long time limit τ +∞.

2

 2

2

 2

  t 0 dτ β(τ )e iω0 τ = iΓ/2δ = e iω0t γ+i∞ γ-i∞ idp 2π e pt β(p) p + iω 0 (A1) or equivalently using p = γ -iω and defining B(ω+iγ) := β(p): iΓ

Ψ∞ dω 2π ω 2 c 2 L

 2 (x, t) = µ 1,2 σ(t)∆(xx 0 ) † . Now, for the electric field we have by definition E(x, t) = D(x, t) -P diel. (x, t) -P Ψ (x, t). Therefore in the full Hilbert space we can define the operatorE (±) (x, t) = (L (x, t)]. (C6)Finally, from Eq. 3.15 we deduce Eqs. 5.18 and 5.19 which requires only L (±) t . Appendix D: Application of the operator L (±) t in the long time t approximationWe start with the calculation of∆ (±) v (τ, x, x ′ ) = +∞ -(±)τ [e -iωτ G v (x, x ′ , ω)] (D1)

2 α

 2 τ [e -iωτ G v,⊥ (x, x ′ , ω)] = α,j c 2 Φ α (x)Φ * α (x ′ )ǫ α,j ⊗ ǫα,j ω 2 α -(ω + i0 + ) -(ω+i0 + ) 2 = 1 2ωα [ 1 ωα-ω-i0 + + 1 ωα+ω+i0 + ] and 1 x-i0 + = P.V.[ 1 x ] + iπδ(x) for x real thenwe obtainL (±) τ [e -iωτ G v,⊥ (x, x ′ , ω)] = α,j c 2 Φ α (x)Φ * α (x ′ )ǫ α,j ⊗ ǫα,j 2ω α (ωω α ) + P.V.[ 1 ω αω ] -iπδ(ω + ω α ) + P.V.[ 1 ω α + ω ]]e -iωτ(D5)which is actually equivalent toL (±) τ [e -iωτ G v,⊥ (x, x ′ , ω)] = α,j c 2 Φ α (x)Φ * α (x ′ )ǫ α,j ⊗ ǫα,j 2ω α (ω α ∓ ω ∓ i0 + ) e -iωτ .

ω 2 c 2

 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 integrali v (x, x ′ , ω)]e -iωτ = +∞ 0 dω 2π e -iω(τ -|x-x ′ |/c)e -iω(τ +|x-x ′ |/c) 4π|xx ′ | . (D10)

1 x-i0 + = 1 iπ

 11 P.V.[1/x]) + δ(x) for x real. Therefore, we deducei v (x, x ′ , ω)]e -iωτ = δ + (τ -|xx ′ |/c)δ + (τ + |xx ′ |/c) 8π|xx ′ | . (D11)We are interested in the regime τ → +∞ and from its definition δ + (τ + |xx ′ |/c) can be neglected. We thus obtaini +∞ 0 dω π Imag[G v (x, x ′ , ω)]e -iωτ ≃ δ + (τ -|xx ′ |/c) 8π|xx ′ | = +∞ 0 dω 2π G v (x, x ′ , ω)e -iωτ .(D12)

  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

	d dt	S(t) = -iω 21 S(t)

Appendix C: Definition of the positive and negative frequency parts of the electric field operator

As we explained in the Appendix B of [59] (see also [START_REF] Marx | [END_REF]) the positive (respectively negative) frequency part of the electric displacement operator D (±) (x, t) is given by:

where we have defined [START_REF] Marx | [END_REF] the operator L

. With this definition we get Eq. 5.13 which ensures a separation between annihilation and creation operators, i.e., between terms containing c α,j (t) and those containing c † α,j (t). Similarly, for the electric dipole density Pdiel. (t) we introduce a separation between positive and negative frequency part by using the definition

where the linear operator M (±) t 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:


David Schmoranzer
email: david.schmoranzer@neel.cnrs.fr

Rasul Gazizulin

Sébastien Triqueneaux

Eddy Collin

Andrew Fefferman

Development of a sub-mK Continuous Nuclear Demagnetization Refrigerator

Keywords: ultra low temperatures, adiabatic nuclear demagnetization, continuous refrigeration techniques, PrNi 5

come   L'archive ouverte pluridisciplinaire

Introduction

In modern low temperature physics, access to ultra-low temperatures below 1 mK has become not only a hallmark of leading laboratories, but also a practical necessity for the investigation of many physical phenomena including superfluidity of 3 He, quantum states of macroscopic objects or dissipation in amorphous matter [START_REF] Pickett | The European Microkelvin Platform[END_REF]. Adiabatic nuclear demagnetization [START_REF] Lounasmaa | Experimental Principles and Methods below 1 K[END_REF][START_REF] Pobell | Matter and Methods at Low Temperatures[END_REF] of copper or van Vleck paramagnets, combined with 3 He - 4 He dilution refrigeration (DR) has become an indispensable technique for reaching these temperatures, despite the limited amount of experimental time available at ultra-low T on standard demagnetization setups. In practice, the experimental window is often comparable to or less than the spin-lattice relaxation times of the materials of interest, leading to insufficient sample thermalization, e.g., a-Si films show hysteresis between warming and cooling curves in the interval 8-20 mK [START_REF] Schmoranzer | Development of a cryogen-free continuous ADR for the constellation-X mission[END_REF]

SW1 R SW2 R W1A R W1B R W2A R W2B R LINK
Fig. 1 Thermal model of the CNDR setup. MXC -mixing chamber of a dilution regfrigerator, HS -heat switch, NDS -nuclear demagnetization stage. Of the resistive wires shown in the setup, only those connecting the nuclear stages to the heat switches were configured to have a significant thermal resistance. Wire 1a and the sample link, represent copper blocks and have negligible thermal resistance in contrast to the thermal links connecting the nuclear stages to heat switches. Additionally, contact resistances were considered at each connection.

Where available, these were estimated from the literature [START_REF] Mueller | Superconducting aluminum heat switch and plated press-contacts for use at ultralow temperatures[END_REF][START_REF] Blondelle | Electrical Conductance of Bolted Copper Joints for Cryogenic Applications[END_REF], but for contacts of wires to PrNi 5 , crude estimates had to be made.

the lowest temperatures, thermal relaxation time of order 1 hour is observed. At 1 mK, phonon conduction (∝ T 3 ) is expected to drop by a factor of order 1000, leading to much longer relaxation times.

To overcome the limitation given by the time window, a continuous nuclear demagnetization refrigerator (CNDR) has been proposed [5,[START_REF] Toda | Design and expected performance of a compact and continuous nuclear demagnetization refrigerator for sub-mK applications[END_REF], utilizing two demagnetization stages connected in series via superconducting heat switches. In this manuscript, we discuss the simulated performance of a similar refrigerator and comment on some of the details of its design. In the following, we present our thermal model of the CNDR setup, details of the numerical calculation and finally discuss the results of the simulations and their implications for the design of the CNDR.

Thermal Model

We consider a model CNDR system with two identical demagnetization stages, each containing 0.2 mol of PrNi 5 connected to the mixing chamber of a DR and between themselves via linking elements such as Cu or Ag wires and Al heat switches. The schematics of the system are shown in Fig. 1.

In the following, we are assuming that each stage consists of several rods of a fixed diameter, which are all connected to the thermal link leading to the heat switches. For the purposes of the numerical model, this is equivalent to a single rod of said diameter, with its length equal to the sum of the lengths of all constituent rods.

Initially, the two demagnetization stages were thus modeled each as a single entity, characterized only by the respective nuclear and electronic temperatures, T n and T e , with the heat exchange between the nuclear spins and the electrons given by the Korringa law [START_REF] Lounasmaa | Experimental Principles and Methods below 1 K[END_REF]. However, a simple estimate of the thermal relaxation time τ rod of a typical cylindrical rod of PrNi 5 shows that this is not likely to be the case in practice. For a rod of radius R, the thermal relaxation time is τ rod = CρR 2 /κ, where C is the heat capacity in constant magnetic field, ρ the density of PrNi 5 and κ its thermal conductivity [START_REF] Folle | Nuclear Refrigeration Properties of PrNi5[END_REF][START_REF] Meijer | The Thermal Conductivity of PrNi5; Electroplating of PrNi5[END_REF]. For a standard 6 mm rod at ultra-low temperatures, this can easily exceed 1 hour due to the low conductivity of the PrNi 5 alloy. The thermal gradient inside the nuclear refrigerant will limit the available cooling power of the rods and a more detailed model is necessary to take this into account.

We therefore use a discrete layered model of each rod (typically with 5, 10, or 20 co-axial cylindrical shells) where each layer is treated separately and has unique nuclear and electronic temperatures. The layers transfer heat to their neighbours as determined by the electronic heat conductivity of PrNi 5 . Only the outermost layer of the rod is in contact with the thermal link to the heat switches or to the sample volume. Tests showed that 10 layers are sufficient and that the results obtained with 20 do not differ significantly.

The heat flux, Q, between the various parts of the setup is modeled using the thermal conductance G(T ) of the linking element and the temperatures of the end points, T 1 and

T 2 as Q = G(T m )[T 2 -T 1 ], with T m = (T 1 +T 2 )/2
. While this relation is not strictly true for heat switches in superconducting state, it was nonetheless used for its simplicity, as the error is given by a factor of order unity whereas the leakage in the open state of the switch is several orders of magnitude below other considered heat leaks. The thermal conductance, G(T ), of a linking element is defined by the thermal conductivity κ(T ) of the element material, by the cross-section A, length L, and by the contact thermal resistances R 1 (T ), R 2 (T ) at both ends of the element:

1 G(T ) = R 1 (T ) + L κ(T )A + R 2 (T ). (1)
For the Al heat switches, separate relations were employed for their normal and superconducting state, describing heat conduction by electrons (κ ∝ T ) and phonons (κ ∝ T 3 ). The contact thermal resistances R 1 (T ), R 2 (T ) were modeled using a temperature-independent electrical resistance and the Wiedermann-Franz law. The cooling power, Qc , at the mixing chamber plate of the DR is estimated conservatively based on the performance of a Bluefors LD-400 dilution refrigerator:

Qc = αT 2 -q0 , (2)
with the exact values of parameters given by α = 0.05838 W K -2 for the temperature coefficient and q0 = 3.854 µW for the internal heat leak, resulting in a base temperature of 8.125 mK. We note that a base temperature of 6.5 mK has been reached with the LD-400 installed in CNRS Grenoble. The resistances were chosen so that the total resistance between NDS1 and NDS2 is equivalent to 50, 150, or 500 nΩ. Several types of parasitic heat leaks were included in the model. First, a constant heat leak at the sample volume was specified, with values of 5, 10, or 20 nW. Additionally, for each demagnetization stage in changing magnetic field B, three types of heat leaks were considered: (1) eddy current heating, Qeddy = 0.03[WT -2 s 2 ] Ḃ2 , (2) vibrational heating, Qvibr = 10 -8 [WT -1 ] |B|, and (3) constant heat leak, Qconst = 2 nW. For eddy current and vibrational heating, the proportionality constants were estimated from the measurements of Parpia et al. [START_REF] Parpia | Optimization procedure for the cooling of liquid 3He by adiabatic demagnetization of praseodymium nickel[END_REF]. We would like to note that vibrational heating is likely to be higher when the CNDR is implemented on a dry dilution unit. Excess heating such as radiation from the magnet, or conduction along centering rings is included in the 2 nW constant heat leak.

Numerical simulations

The computations were performed by means of integrating the differential equations relating the changes in magnetic fields, changes in temperature, heat fluxes and heat leaks. The code was programmed in Python, with the integration performed using the odeint solver of the SciPy package.

The configuration of the setup was stored in an INI file, and the demagnetization cycles were scripted in a simple custom language providing basic commands for the manipulation of magnetic fields and heat switches. The stability of the computation was governed by the maximum time step permitted for the odeint solver. Introducing the layered model of the demagnetization stages has significantly affected the requirements for the maximum permissible time step to keep the code stable. The initial conditions were chosen so that the entire setup starts at the stable base temperature of the DR unit and both PrNi 5 stages are in zero field. Further evolution of the magnetic fields was fully controlled by the scripts.

Initially, the computation used paramagnetic values of entropy for PrNi 5 , and changes of temperature due to changing magnetic field were calculated using standard expressions [START_REF] Lounasmaa | Experimental Principles and Methods below 1 K[END_REF]. The following parameters were used for PrNi 5 : nuclear spin I = 5/2, nuclear gyromagnetic ratio γ n = 1.71, Knight shift K = 11.2, internal field b = 66 mT [START_REF] Lounasmaa | Experimental Principles and Methods below 1 K[END_REF][START_REF] Pobell | Matter and Methods at Low Temperatures[END_REF][START_REF] Kubota | Nuclear Magnetic Ordering in PrNi5 at 0.4 mK[END_REF], and Korringa constant κ K = 0.001 (an upper estimate from Ref. [START_REF] Pobell | Matter and Methods at Low Temperatures[END_REF]). Even with the internal field b accounted for, this approach led to difficulties at the lowest temperatures, as the spontaneous ordering of PrNi 5 occurring around 0.4 mK is not described within the paramagnetic model.

Therefore, a mathematical model of the entropy of PrNi 5 was constructed based on the measurements of Kubota et al. [START_REF] Kubota | Nuclear Magnetic Ordering in PrNi5 at 0.4 mK[END_REF], with limits corresponding to the paramagnetic behavior at high T and a steeper dependence of entropy on temperature near the ordering temperature, approximating the experimental data on entropy and specific heat. The resulting lines of constant magnetic field B are shown in the temperature-entropy diagram in Fig. 2 and compared with the experimental data of Ref. [START_REF] Kubota | Nuclear Magnetic Ordering in PrNi5 at 0.4 mK[END_REF]. This model was then used to calculate the nuclear demagnetization properties of the two stages: In each time step, first, the adiabatic change in temperature due to the changing magnetic field was evaluated from the model entropy and then, non-adiabatic effects due to heat exchange and heat leaks were added. Fig. 2 Model entropy of PrNi 5 (solid lines) superimposed on the measurements presented in Ref. [START_REF] Kubota | Nuclear Magnetic Ordering in PrNi5 at 0.4 mK[END_REF]. Minor discrepancies exist at low fields (c.f. 37 mT), and might affect the final temperatures below 1 mK, but the principal features related to the ordering are captured qualitatively. The model entropy curves are virtually equivalent to paramagnetic ones down to fields ≈90 mT. Note that the explicit inclusion of the internal field in the paramagnetic model would shift the zero field curve to 66 mT in contrast to experiment and it cannot yield the increased slope (higher heat capacity) at the transition near 0.4 mK.

Results and Discussion

The thermal cycles of the two stages of the refrigerator are illustrated in the entropy-temperature diagram in Fig. 3 for the case of 3 mm rods, 50 nΩ resistance and 5 nW heat leak. The durations of the steps were manually adapted for each value of heat leak to produce near-optimum final temperatures. A sample time trace of the relevant temperatures and magnetic fields is shown in Fig. 4, together with the dependence of the final temperature on the sample heat load. The final temperature is extracted in each case as the maximum of the second stage temperature within a stable cycle.

The main factor limiting the performance of the CNDR is the slow heat exchange between the two stages when thermally coupled at low T. It is therefore necessary to provide a low resistance path for the heat from one PrNi 5 sample to the other. To achieve a value of resistance comparable to 50 nΩ, optimized contacts of high-conductivity wires with the PrNi 5 as well as with the Al heat switch are of particular importance. While suitable solutions for contacting Cu to Al exist [START_REF] Shigematsu | Proceedings of the Sixteenth International Cryogenic Engineering Conference/International Cryogenic Materials Conference[END_REF], our research in this area and practical tests will be reported elsewhere. The more efficient thermalization of the 3 mm rods seems to contribute very little to the performance of the CNDR, as the sample space is only in contact with the outside layer of the rods. Any practical advantage to be gained by replacing the standard 6 mm PrNi 5 rods by 3 mm ones is thus related to the increased surface area of the PrNi 5 stage and to D -h e a t e x c h a n g e w i t h M X C Fig. 3 Thermal cycles of the first (solid, red) and second (dashed, blue) nuclear demagnetization stages (NDS) of the CNDR. NDS 1 cycle: A -demagnetization to 8 mT while thermally decoupled, B -heat exchange with the second stage, C -magnetization to 1.2 T, D -thermal coupling to the mixing chamber. The steps C and D overlap to save time during the magnetization process -the first stage is coupled to the MXC as soon as its temperature exceeds that of the MXC. NDS 2 cycle: 1 -slow demagnetization from 60 to 50 mT while thermally decoupled, 2 -magnetization to 60 mT while connected to the first stage, 3continued cooling by the first stage at constant magnetic field.

the correspondingly reduced thermal boundary resistance, which could help us approach the value of 50 nΩ for the critical thermal link.

Conclusions

We have shown that a continuous nuclear demagnetization refrigerator operating with two PrNi 5 stages connected in series can reach and maintain temperatures below 1 mK at heat loads up to 20 nW. This makes the operation of such a refrigerator feasible on commercial cryogen-free dilution refrigerators. The key element in the setup is the thermal link between the two stages, including a heat switch and it is crucial to optimize thermal resistances of all its components and their contacts. The results of our numerical simulations are in general agreement with those of Ref. [START_REF] Toda | Design and expected performance of a compact and continuous nuclear demagnetization refrigerator for sub-mK applications[END_REF], but we have used a more detailed description of the demagnetization stages, taking into account their geometry and heat conductivity. In addition, we have used a mathematical model of PrNi 5 thermodynamic properties that allow us to describe its refrigeration capabilities in the vicinity of its ferromagnetic ordering occurring near 0.4 mK more precisely than the usual paramagnetic model. 4 Left: Sample time trace of temperatures and magnetic fields of the CNDR setup for 6 mm PrNi 5 rods with a 50 nΩ resistance of the entire link between NDS 1 and NDS 2, with a heat load of 20 nW on the sample. Right: sample heat load vs. the final temperature attainable with the CNDR setup for the indicated rod sizes and resistances of the link between the stages. The final temperature is the maximum temperature encountered during the cycle of the CNDR. As noted above, the exact final temperatures will likely differ from experiment, due to inaccuracies of the model entropy. Considering the maximum error based on Fig. 2, the data obtained for 50 nΩ resistance still remain below 1 mK.

1 m

 1 e a t e x c h a n g e w i t h N D S 2 C -m a g n e t i z a t i o n t o 1 . 2 T M o l a r e n t r o p y ( J K -

  Fig.4Left: Sample time trace of temperatures and magnetic fields of the CNDR setup for 6 mm PrNi 5 rods with a 50 nΩ resistance of the entire link between NDS 1 and NDS 2, with a heat load of 20 nW on the sample. Right: sample heat load vs. the final temperature attainable with the CNDR setup for the indicated rod sizes and resistances of the link between the stages. The final temperature is the maximum temperature encountered during the cycle of the CNDR. As noted above, the exact final temperatures will likely differ from experiment, due to inaccuracies of the model entropy. Considering the maximum error based on Fig.2, the data obtained for 50 nΩ resistance still remain below 1 mK.

  . At

	MXC	wire 1a	HS 1	wire 1b	NDS 1	wire 2a	HS 2	wire 2b	NDS 2	sample link	Sample
			R

Acknowledgements We acknowledge support from the ERC StG grant UNIGLASS No. 714692 and ERC CoG grant ULT-NEMS No. 647917.