Master equation for multilevel interference in a superradiant medium
Aleksei Konovalov Giovanna Morigi
quant-ph

























Theoretische Physik, Saarland University, 66123 Saarbrücken, Germany
















We derive a master equation for a superradiant medium which includes multilevel interference betwen the individual scatterers. The derivation relies on the Born-Markov approximation and implements the coarse graining formalism. The master equation fulfils the Lindblad form and contains terms describing multilevel interference between parallel transitions of a single atom, multi-atom interference between identical transitions, and multi-atom interference between different electronic transitions with parallel dipoles.  This formalism is then applied to determine the excitation spectrum of two emitters using the parameters of the Hydrogen transitions 2S_1/2 →4P_1/2 and 2S_1/2 →4P_3/2, where the gap between the parallel dipoles is of the order of GHz. The distortion of the signal due to the interplay of multilevel and multi-emitter interference is analysed as a function of their distance. These results suggest that interference between parallel dipolar transition can significantly affect the spectroscopic properties of optically dense media.  








Master equation for multilevel interference in a superradiant medium
    Giovanna Morigi
    December 30, 2023
====================================================================







§ INTRODUCTION




Superradiance generally denotes a phenomenon which enhances radiation. In quantum optics, it originates from quantum interference in the light emission by an ensemble of atoms, molecules, or other types of resonant emitters which form an optically-dense medium <cit.>. In free space this requires that the average interparticle distance is smaller than the wavelength of the scattered radiation. Then, the coupling of the individual atomic transitions with the modes of the electromagnetic field can be effectively described in terms of collective dipoles and the radiative properties depend on the collective spin quantum numbers <cit.>. Superradiant (and subradiant) scattering plays a relevant role in the spectroscopy of dense atomic gases <cit.>, it could enhance transport of light in organic semiconductors <cit.>, and it is the key mechanism of recent realizations of ultranarrow lasers <cit.>.

    

Superradiant light scattering is often described by means of a perturbative expansion in the atom-photon interactions and using the Born-Markov approximation <cit.>. Most theoretical treatments focus on two-level dipolar transitions <cit.>, some also including the possible degeneracy of the ground or excited state of the transition <cit.>. These treatments successfully predict experimental measurements at sufficiently low optical densities. Qualitative discrepancies have been found when comparing the predictions of these models with recent experiments with dense atomic media <cit.>. This requires one to assess the effects of terms which are typically discarded or only partially considered.  

In this work we derive a master equation for an optically dense medium and set our focus on vacuum-induced interference <cit.>. Vacuum-induced interference refers to interference phenomena between electronic transitions coupled to common modes of the electromagnetic field. If the transitions are dipolar, they are denoted by parallel dipoles. Interference occurs also when the electromagnetic field modes are in the vacuum, and is qualitatively different from laser-induced interference <cit.>. In closed level structures these effects can be tested by means of quantum beat spectroscopy and are expected to give rise to "steady-state quantum beats" <cit.>. They are also expected to play an important role in high-precision spectroscopy <cit.>. In this work we determine the Born-Markov master equation of multilevel scatterers in an optically dense medium and which includes interference terms between parallel dipoles . For this purpose we derive the master equation by applying the coarse-graining formalism of Ref. <cit.>. The master equation we obtain preserves the Lindblad form and in the limit of one single emitter it reduces to the coarse-grained master equation of Ref. <cit.>. We then apply it to determine the excitation spectrum and the light shift of two identical emitters, each composed by two parallel dipoles sharing the same ground state. 
In this simplified model we show that collective scattering results from the coherent sum of three processes, which we illustrate in Fig. <ref>: (a) the interference between parallel dipoles of the individual atoms, (b) the interference between resonant transitions of different atoms, and (c) the interference between parallel dipoles of different atoms. Here, we argue that in an optically dense medium they can give rise to measurable shifts of the spectroscopic lines.  

This work is organized as follows. In Sec. <ref> we present the derivation of the master equation by eliminating the degrees of freedom of the electromagnetic field within the Born-Markov approximation and by implementing the coarse-graining method developed in Ref. <cit.>. By these means we obtain a superoperator that fulfils the Lindblad form. This superoperator consistently describes interference processes between parallel dipoles of the individual atoms and interference processes of different atoms. In Sec. <ref> we then consider the specific example of two emitters, composed by two parallel dipoles sharing the same ground state, and determine their excitation spectrum using the parameters of the transitions 2S_1/2→4P_1/2 and 2S_1/2→4P_3/2 of Hydrogen atom. By means of a simple fitting function we argue that the interference effects give rise to measurable shifts of the resonance lines. Finally, in Sec. <ref> we draw the conclusions and discuss outlooks of this work. The appendices contain details of the calculations in Sec. <ref> and Sec. <ref>.
   


§ DERIVATION OF THE SUPERRADIANT MASTER EQUATION



In this Section we report the derivation of the Born-Markov master equation for an optically dense atomic or molecular medium. Our derivation follows the lines of textbook derivations <cit.>, and extend it by implementing the coarse-grained method developed in Ref. <cit.>. This allows us to take systematically into account the interference of parallel dipoles and at the same time to preserve the Lindblad form of the master equation. In the single-atom limit our master equation reproduces the one derived in Ref. <cit.>, which includes the interference processes between parallel dipoles in a single atom.

For convenience, in the following we assume an ensemble of emitters with identical electronic transitions. This formalism, nevertheless, can be straightforwardly extended  to ensembles of different particles (which could also be a mixture of atoms and molecules) with quasi-resonant transitions. The relevant assumption is that the emitters are pinned at given positions and are distinguishable particles. Our starting point is the von-Neumann equation governing the coherent dynamics. Below we provide the salient steps leading to the corresponding coarse-grained master equation for the emitters' internal degrees of freedom. 



 §.§ Multilevel emitters interacting with the quantum electromagnetic field


We consider N emitters interacting with the modes of the electromagnetic field (EMF) in the volume V. We assume that the particles are pinned at the positions R⃗_α, with α=1,…,N. We denote by ℋ the Hilbert space of the emitters' internal degrees of freedom and of the EMF's degrees of freedom,  ℋ= ℋ_A⊗ℋ_R. The time evolution of the density matrix χ̂(t), describing the state of photons and emitters, is governed by the von-Neumann equation

    ∂_t χ̂ = [Ĥ,χ̂]/iħ ,
   
where Ĥ is the Hamiltonian determining the dynamics, which we decompose into the sum of the Hamiltonian Ĥ_A for the emitters' (internal) degrees of freedom, the Hamiltonian Ĥ_R for the free EMF, and the emitter-photon interactions V̂:

    Ĥ=Ĥ_A+Ĥ_R+V̂ .
  
We remark here that Ĥ_R≡1̂_A⊗Ĥ_R and  Ĥ_A≡Ĥ_A⊗1̂_R, 
where 1̂_R and by 1̂_A the identity operators in the Hilbert spaces ℋ_R and ℋ_A, respectively. Thus, we use the same notation for the operator Ĥ_j=A,R defined in the extended Hilbert space  ℋ and in the reduced Hilbert space ℋ_j

The emitters' Hamiltonian. The emitters Hamiltonian describes the dynamics of the internal degrees of freedom of N emitters: 
 
    Ĥ_A =∑_α=1^N Ĥ_A_α ,

 where Ĥ_A_α is the Hamiltonian of emitter α=1,…,N at position R⃗_α and we assume that the size of the center-of-mass wavepacket is much smaller than the interparticle distance (in Eq. (<ref>) we omit to explicitly write that Ĥ_A_α is the identity operator in the Hilbert space of the emitters with β≠α).  
 We consider here only the lowest electronic bound states assuming that the system is at room temperature. The spectrum of each emitter is discrete and the Hamiltonian in diagonal form reads

    Ĥ_A_α =∑_n E_n |n ⟩_α⟨ n | ,

with E_n the eigenvalue and |n ⟩_α the corresponding eigenvector for the emitter at the position R⃗_α. In a more general treatment, where the emitters might not be identical and/or in the presence of spatial inhomogeneity, then the energy also depends on the label α. 

The quantum electromagnetic field. We treat the EMF in second quantization and choose the Coulomb gauge. We denote the quantization volume by V and assume periodic boundary conditions. The energy of the field relative to the vacuum energy reads:

    Ĥ_R= ∑_λħω_λâ_λ^†â_λ ,

where λ denotes the sum over the EMF modes  and the sum has an upper cutoff given by the energy ħω_cutoff∼ mc^2, with m the electron mass.  The modes are here traveling waves and are fully characterized by the wave vector k⃗_λ and by the transverse polarization e⃗_λ, with the frequency ω_λ=c|k⃗_λ| and c the speed of light in vacuum. Operators â_λ and â_λ^† annihilate and create, respectively, a photon of mode λ, and fulfil the bosonic commutation relations [â_λ,â_λ'^†]=δ_λ,λ' and  [â_λ,â_λ']=0.

The initial state of the EMF field is assumed to be given by the thermal distribution 

    R̂=exp(-Ĥ_R/k_BT)/Z ,

where k_B is Boltzmann's constant, T is the temperature, and Z=Tr{exp(-Ĥ_R/k_BT)} is the partition function. Within the validity of the Born approximation, R̂ gives the state of the EMF at all times. Here we assume room temperatures, T ∼ 300  K.

Emitter-photon interactions. Emitter-photon interactions are here treated in the electric-dipole approximation. Operator V̂ is the sum of the interactions of the fields with each emitter, V̂=∑_α=1^N V̂_α, with

    V̂_α=ħ∑_nΓ̂_n^ασ̂_n^α ,

where the sum is over all pairs of electronic levels n=(n_1,n_2) coupled by an electric dipole transition. 
Here, operator σ̂_n^α describes the transition between |n_1⟩_α and |n_2⟩_α:

    σ̂_n^α≡ |n_1⟩_α⟨ n_2| .

The corresponding coupling strength is determined by the coupling operator Γ̂_n^α, which acts over the degrees of freedom of the electromagnetic field and reads:

    Γ̂_n^α= ∑_λ(g_n^αλâ_λe^ik⃗_λR⃗_α + 
    g̅_n^αλâ^†_λe^-ik⃗_λR⃗_α) .

The coupling strengths g_n^αλ have the dimensions of a frequency and are below given in Gauss units and in the length gauge:

    g_n^αλ=-i√(2πω_λ/ħ V) d⃗_n^α·e⃗_λ ,
          g̅_n^αλ=i√(2πω_λ/ħ V) d⃗_n^α·(e⃗_λ)^* ,

with d⃗_n^α the dipole moment of the transition, which is the matrix element of the dipole operator d̂⃗̂_α and reads d⃗_n^α=_α⟨ n_1|d̂⃗̂_α|n_2⟩_α. We remark that this description applies the long-wave approximation, and thus it is valid when the size of the electronic wave packet is smaller than the optical wavelength. Moreover, in our model we did not include the self-energy which appears in the length gauge (see Refs. <cit.> for an insightful discussion). 

For later convenience we introduce the frequency ω_n:

    ω_n=(E_n_1-E_n_2)/ħ .

By definition it can also take negative values. 

 



 §.§ Master equation for an ensemble of multilevel emitters


We now proceed in deriving the Born-Markov master equation using the coarse-grained formalism. The procedure repeats in the essential steps the one of Ref.<cit.>, with some notable differences due to the many-body nature of the problem.

We first introduce the density matrix ρ̂(t) describing the state of the emitters at time t. Operator ρ̂(t) is defined in the Hilbert space ℋ_ A and is related to the density matrix χ̂(t) by the equation: ρ̂(t)= Tr_R{χ̂(t)}, where Tr_R denotes the partial trace over the degrees of freedom of the EMF. 

We now consider the von-Neumann equation, Eq. (<ref>), and move to the interaction picture with respect to Hamiltonian Ĥ_0=Ĥ_ A+Ĥ_ R. 
We denote the system's density matrix in interaction picture by

    χ̃(t)=Û_0(t)^†χ̂(t)Û_0(t) ,

where we have introduced the unitary operator Û_0(t)=exp(Ĥ_0t/( iħ)). In this representation the reduced density matrix of the system is related to the reduced density matrix in Schrödinger picture by the relation:

    ρ̃(t)= Tr_R{χ̃(t)}= e^-Ĥ_ At/( iħ)ρ̂(t) e^Ĥ_ At/( iħ) .


In interaction picture the unitary operator determining the  time evolution reads:

    Ũ(t,t') = 𝒯exp(- i/ħ∫_t^t' dt_1 Ṽ(t_1)  ) ,
 
where Ṽ(t)=Û_0(t)^†V̂Û_0(t) and 𝒯 denotes the time ordering, such that 

    𝒯Ṽ(t_1) Ṽ(t_2) = Ṽ(t_1) Ṽ(t_2) θ (t_1-t_2) + Ṽ(t_2) Ṽ(t_1) θ (t_2-t_1) ,
 with θ (t) the Heaviside function. Using this formalism, at t'>t the state χ̂̃̂(t) evolves into state 

    χ̃(t')= Ũ(t,t') χ̃(t) Ũ(t,t')^† .




  §.§.§ Dyson equation and Born-Markov approximation
 

Let now Δ t=t'-t>0 denote a finite and sufficiently small time step, which we quantify later. We write the Dyson series of the right-hand side of Eq. (<ref>) till second order in the interaction, but keep the exact form. After tracing out the EMF degrees of freedom we obtain the expression

    ρ̃(t+Δ t) =
      ρ̃(t) +Δ t ∑_αΛ_1^αρ̃(t)
          +Δ t ∑_α,β1/Δ t∫_t-Δ t^t+Δ t dT∫_-Δ t^Δ t dτθ(τ)Λ_2^α,β(T,τ)ρ̃(T-τ) .

The terms Λ_1^α, Λ_2^α on the RHS are linear maps, the subscript indicate the order in the Dyson expansion. In deriving Eq. (<ref>) we have made the Born approximation at the initial time t, namely, we have assumed that there are no quantum correlations at time t between EMF and emitter. This corresponds to writing χ̃(t)=R̂⊗ρ̃(t) where here R̂ is the thermal state of the EMF, Eq. (<ref>). 

 The map Λ_1^α acts over the Hilbert space of the emitter α and is given by:

    Λ_1^αρ̃(t)   =   1/ iħΔ t∫_t^t+Δ t dt_1 Tr_R{[Ṽ_α(t_1) , χ̃(t)] } .
       =   1/ iħ[ ⟨Ṽ_α(t)⟩_R,ρ̃(t)] ,

where between the first and the second line we have applied the Born approximation and introduced the time-averaged operator (here in interaction picture):

    ⟨Ṽ_α(t)⟩_R=1/Δ t∫_t^t+Δ tTr_R{Ṽ_α(t_1) R̂} .

Note that operator V̂_α, Eq. (<ref>), vanishes over the thermal state of the EMF, Eq. (<ref>). 
The second integrand of Eq. (<ref>) contains the Heaviside function θ(τ) and includes also the coupling between different emitters. Its detailed form is reported in Appendix <ref>.

Equation (<ref>) is generally valid for sufficiently short time intervals Δ t, over which one can assume that the Born approximation holds. After some time, in fact, the interactions establish quantum correlations between system and reservoir. These correlations can be neglected when the interactions can be treated perturbatively. 

The master equation becomes local in time when the Markov approximation holds. The Markov approximation consists in approximating ρ̃(T-τ)≈ρ̃(t) in Eq. (<ref>). It is equivalent to the Wigner-Weisskopf approximation for the propagator <cit.>  and is justified when the characteristic time scale τ_R of the correlation function C_αβ(τ), Eq. (<ref>), is orders of magnitude smaller than the system's relaxation time. In a thermal bath the correlation function is composed by a term which decays exponentially with the correlation time τ_R=ħ/k_BT and by power-law tails that can be discarded for typical evolution times <cit.>. At room temperatures, T ∼ 300 K, this time is of the order of τ_R∼ 10^-13 sec. This time shall be compared with the relaxation time of the system. For optical transitions the natural linewidth of a single atom, γ∼ 2π× 10^6-10^8 Hz,  fulfils γτ_R≪ 1. In this limit we can choose the time scale Δ t such that τ_R≪Δ t≪ 1/γ and ignore memory effects in the integral. 

In the presence of dipole-dipole interactions there are some issues to be considered: in first place, superradiance gives rise to an N-fold enhancement of the single atom decay rate, thus when Nγ becomes comparable with 1/τ_R the approximation becomes invalid. This is the regime where one can observe the Dicke phase transition in an ensemble of two-level systems <cit.>, and where the assumptions at the basis of this treatment break down. At the same time, subradiant states can be characterized by extremely small linewidths. Observing their decay requires one to analyse the system's dynamics over long time scales, over which the power-law tails of the correlation function can become important. These considerations suggest that the formalism shall be revisited for media with very high optical dense media. 



  §.§.§ Coarse-grained master equation


In what follows we assume an optically dense medium for which the Born-Markov approximation is valid. Then, from Eq. (<ref>) we derive the Born-Markov master equation (now back in 
Schrödinger picture):

    ∂_tρ̂=1/ iħ[Ĥ_A+Ĥ_S,ρ̂(t)]+ℒ_Dρ̂(t) ,

where Hamiltonian Ĥ_S and superoperator (dissipator) ℒ_D contain both the single-atom as well as the interatomic interference terms between parallel dipoles.  The details of the derivation are  standard and are reported in Appendix <ref>. The master equation is valid for any time t>0 within a grid whose resolution is determined by the coarse-grained time-scale Δ t. As a consequence, the coefficients multiplying the terms of the operator Ĥ_S  and the superoperator ℒ_D are scaled by the function

    Θ_ij^(Δ t)=sin((ω_i+ω_j)Δ t/2)/(ω_i+ω_j)Δ t/2 .

This term selects transitions which are resonant within the resolution set by the coarse-graining time Δ t. For optical transitions, this factor selects a pair of frequencies ω_i and ω_j with opposite signs. Correspondingly, it selects terms in Hamiltonian and dissipator where the pairs of operators σ̂_i^ασ̂_j^β describe an excitation and a de-excitation along two (quasi-)resonant transitions. For convenience, we introduce the operator ζ̂_i^α†≡σ̂_i^α, which describes a transition i_2→ i_1 with ω̅_i=ω_i>0 and dipole moment D⃗_i^α*=d⃗_i^α. Then, the operators appearing in the master equation are of the form ζ̂_i^α†ζ̂_j^β or ζ̂_i^αζ̂_j^β† and the factor (<ref>) now reads

    Θ_ij^(Δ t)=sin((ω̅_i±ω̅_j)Δ t/2)/(ω̅_i±ω̅_j)Δ t/2 .
 
 In what follows we discard the processes where two transitions are simultaneously excited or de-excited, corresponding to the + sign in the argument of Eq. (<ref>).
 
 
Hamilton operator. The Hamiltonian term due to the interaction with the EMF is given by the expression

    Ĥ_S=∑_α⟨V̂_α⟩_R+1/2∑_α,βĤ_αβ^S ,

where ⟨V̂_α⟩_R is given in Eq. (<ref>) and is now reported in Schrödinger picture. This latter term vanishes, since we assume that the EMF is in the thermal state. The Hamilton operator Ĥ_αβ^S contains the frequency shifts and couplings due to the multilevel interference, and is derived from the expression (here given in interaction picture):

    H̃_αβ^S   =   - i/2ħΔ t∫_t^t+Δ t dt_1∫_t^t+Δ t dt_2
          ×θ(t_1-t_2) Tr_R{[Ṽ_α(t_1),Ṽ_β(t_2)]R(t)}+ H.c. .

For α=β it is the Hamilton operator for a single atom and it coincides with the operator derived in Ref.
<cit.>. For  α≠β 
it describes the Hamiltonian terms due to the dipole-dipole interactions, including the interference between all parallel transitions of different atoms. We report it in the form which includes both cases:

    Ĥ_αβ^S = -ħ∑_i,j      [(Δ_ij^αβ-+Δ_ij^αβ(T)) ζ̂_i^α†ζ̂_j^β.
          +.(Δ_ij^αβ+-Δ_ij^αβ(T))^*   ζ̂_i^αζ̂_j^β†]+ H.c. ,

where Δ_ij^αβ(T)=Δ_ij^αβ-(T)-Δ_ij^αβ+(T) and the individual coefficients read (below in Gauss units):

    Δ_ij^αβ±=Θ_ij^(Δ t)D⃗_i^α*·D⃗_j^ β/(2π)^2 ħ c^3𝒫∫_0^ω_ cutdω ω^3/ω±ω_ij F_αβ^ij(R⃗_αβ)
    
          Δ_ij^αβ±(T)=Θ_ij^(Δ t)D⃗_i^α*·D⃗_j^ β/(2π)^2 ħ c^3𝒫∫_0^ω_ cutdω ω^3n(ω,T)/ω±ω_ij F_αβ^ij(R⃗_αβ) .

Here, 𝒫 denotes the Cauchy principal value and ω_ cut is the cutoff frequency. The frequency 

    ω_ij=ω̅_i+ω̅_j/2

is the average between the two transition frequencies, and the coefficient F_αβ^ij(R⃗_αβ) depends also on the distance R⃗_αβ=R⃗_α-R⃗_β between the atoms and on the wave number k=ω/c. It takes the form

    F_αβ^ij(R⃗_αβ) = 
    4π(  j_0(kR_αβ)  [1- (D⃗_i^α·R⃗_αβ)^*(D⃗_j^ β·R⃗_αβ)/D_i^αD_j^βR_αβ^2] .
    
          - . j_1(kR_αβ)/kR_αβ[1- 3(D⃗_i^α·R⃗_αβ)^*(D⃗_j^ β·R⃗_αβ)/D_i^αD_j^βR_αβ^2] ) ,
 
where we used the notation D_i^α=|D⃗_i^α| and R_αβ=|R_αβ|. Here, j_0(x) and j_1(x) are spherical Bessel functions of the first type <cit.>. The dependence on the vector joining the two atoms breaks the spherical symmetry and is at the origin of the anisotropic light emission of superradiance <cit.>. For the case of one atom, N=1, one has F_αα^ij(0)= 8π/3 <cit.>, and Hamiltonian of Eq. (<ref>)
takes the form of the single-atom Hamiltonian of Ref. <cit.>.

Dissipator. The Lindblad term ℒ_D describes the incoherent processes. It can be decomposed into the sum

    ℒ^Dρ̂(t)=∑_α,βℒ_D^αβρ̂(t) ,

where the terms with α=β describe the dissipation of N non-interacting atoms, while the terms with α≠β originate from multiple scattering of resonant photons and vanish when the distance between the atoms exceeds several wavelengths. The individual terms are obtained from the expression in interaction picture

    ℒ̃_D^αβρ̃(t)=1/2ħ^2Δ t∫_t^t+Δ t dt_1∫_t^t+Δ t dt_2 Tr_R{𝒜(t_1,t_2)} ,

where

    𝒜(t_1,t_2)   =   2Ṽ_β(t_1)[ρ̃(t)⊗R̃(t)]Ṽ_α(t_2)
          -[Ṽ_α(t_1)Ṽ_β(t_2),ρ̃(t)⊗R̃(t)]_+

and [,]_+ denotes the anticommutator. After performing the integration and going back to Schrödinger picture, the individual terms take the form

    ℒ_D^αβρ̂(t)   =   ∑_i,j (1+n(ω_ij,T))
          ×(Γ_αβ^ij/2[ζ̂_j^βρ̂(t), ζ̂_i^α†] + Γ_αβ^ij/2[ζ̂_j^β, ρ̂(t) ζ̂_i^α†] ) 
       +   ∑_i,j n(ω_ij,T)
          ×(Γ_αβ^ij*/2[ζ̂_j^β†ρ̂(t), ζ̂_i^α] + Γ_αβ^ij*/2[ζ̂_j^β†, ρ̂(t) ζ̂_i^α] ) ,

with the damping coefficients

    Γ_αβ^ij=Θ_ij^(Δ t)D⃗_i^α*·D⃗_j^ β/2πħ c^3ω_ij^3 F_αβ^ij(k_ij)

and k_ij=ω_ij/c. We note that for i≠ j the damping coefficients are different from zero if the scalar product D⃗_i^α*·D⃗_j^ β≠ 0. 
Master equation (<ref>) fulfils the Lindblad form and take into account the multilevel structure of the quantum emitters. 



  §.§.§ Discussion


We first review the dynamics that the master equation (<ref>) predicts for an very dilute ensemble of emitters (R_αβ→∞), when it is well approximated by N independent experiments with a single atom. In this case the damping coefficients Γ_αα^ii are the Einstein coefficients of spontaneous emission. For i≠ j, instead, the coefficients Γ_αα^ij describe processes where two different transitions with parallel dipoles are simultaneously de-excited. These transitions shall be resonant within the frequency resolution of the coarse graining 1/Δ t.  This process, even though incoherent, is a quantum interference between spectral lines <cit.>. The corresponding terms have been denoted by cross-damping terms in the literature <cit.>. These dynamics have a corresponding hermitian component in the Hamiltonian term Ĥ_S^αα. The coefficients include an energy shift of the electronic states due to the vacuum fluctuations, which for the ground state is the non-relativistic Lamb shift, as well as a shift due to thermal fluctuations of the EMF. Vacuum and thermal fluctuations give also rise to an effective coupling between electronic levels with parallel dipoles and quasi-resonant frequencies, the coupling coefficients are given by Eqs. (<ref>)-(<ref>) after setting α=β.  They can be estimated by using the approximate relation <cit.>

    Δ_ij^αα±≈1/2 (D⃗_i^α*·D⃗_j^ α) Θ_ij^(Δ t)( 1/|D⃗_i^α|^2Δ_ii^αα± + 1/|D⃗_j^α|^2Δ_jj^αα±) .
   

When the interparticle distances are comparable with the wavelength, namely for α≠β, Eq. (<ref>) is the master equation for optically dense media which now includes quantum interference between transitions with parallel dipoles. Keeping only the terms with i=j one obtains the master equation discussed in the literature <cit.>, where the dissipator gives rise to phenomena such as superradiance and subradiance, while the coherent part describes coherent dipole-dipole interaction, including frequency shifts such as the so-called collective Lamb shift <cit.>. Our derivation highlights, in addition, the existence of interference terms between quasi-resonant transitions of different atoms with parallel dipoles both in the incoherent as well as in the coherent part of the master equation. 

We finally remark that, by taking the limit Δ t→ 0, thus for an infinitesimally small coarse-grained time-scale, the function (<ref>) becomes a Dirac-delta function. Then, the coarse-grained master equation reduces to the Born-Markov master equation discussed for instance in Refs. <cit.>. In this limit, however, one discards effects due to the finite time-scale of the reservoir dynamics, thus interference phenomena between parallel transitions which are close in frequency but not exactly resonant. The coarse-graining master equation allows one to include these dynamics in a systematic way. We refer the interested reader to Refs. <cit.> for discussions on the coarse-grained master equation and to the next section for a discussion about the choice of Δ t.
 
 

§ EXCITATION SPECTRUM OF TWO EMITTERS
 
 
  	

We now determine the excitation spectrum of two emitters, which are pinned at the positions R⃗_1=0 and R⃗_2=R⃗ and are uniformly driven by a linearly polarized laser. Their electronic configuration is composed by three electronic levels of Hydrogen, which consists of the ground state |1⟩ and the two excited states |2⟩ and |3⟩. The transitions |1⟩→|2⟩ and |1⟩→|3⟩ are parallel optical dipoles with moments D⃗_12^α and D⃗_13^α, respectively, the transition frequencies are denoted by ω_12 and ω_13 (from now on ω̅_1e=ω_1e>0 with e=2,3). The reduced level structure allows us to highlight the effects of multilevel interference. Despite the fact we consider the parameters of two transitions of the Hydrogen atoms, however, the choice we perform breaks the rotational symmetry of the atoms. This shall be kept in mind when discussing the single-emitter properties.

The dynamics induced by the laser is described by a Hamiltonian term, which is added to the Hamilton operator of Eq. (<ref>). This procedure corresponds to assuming that the laser field is described by a coherent state and to moving to the reference frame where the quantum state of the laser field is in the vacuum <cit.>. We denote by ω_L the laser frequency, and assume that the laser polarization is linear and that the spatial dependence of the laser field wave vector k⃗_L is orthogonal to the vector R⃗ joining the two emitters. The laser-atom Hamiltonian has the form
 
    Ĥ_L=-ħ∑_α=1,2∑_e=2,3 
      	g_1e^α e^- iω_Ltζ̂_1e^α†+ H.c. ,

where we have introduced the Rabi frequency g_1e^α=-d⃗_1e^α·E⃗_L/2ħ, which depends on the electric field amplitude E⃗_L. The master equation takes the form

    ∂_tρ̂=1/ iħ[Ĥ_A+Ĥ_S,ρ̂(t)]+ℒ_Dρ̂(t)+1/ iħ[ Ĥ_L,ρ̂] ,

where now the sums over the atoms run to N=2 and the sums over the internal transitions include just the two transitions with parallel dipolar moments. For simplicity, thus, we can now replace the sum over the transitions i=i_1,i_2 with the sum over the excited state e=2,3. Using the simplified level structure we simplify the Hamiltonian term Ĥ_12^S, Eq. (<ref>), as follows:

    Ĥ_12^S=-  ∑_e,e'=2^3ℱ_c(ω_1e+ω_1e')Ξ_e e'^F(R⃗)ζ̂_1e^1 †ζ̂_1e'^2 + H.c. ,

where ℱ_c(ω_1e+ω_1e') is obtained by means of a smoothening of the fast-oscillating function Θ_ij^(Δ t), see Ref. <cit.> and Sec. <ref>, and

    Ξ_e e'^F(R⃗) =
      D⃗_1e·D⃗_1e'(ω_e e'/c)^3(  y_0(kR) -  y_1(kR)/kR) .

Here we used that the atomic dipole moments are real vectors and introduced the notation ω_e e'=(ω_1e+ω_1e')/2. Moreover, we have used that the dipole moments are orthogonal to the vector connecting the two atoms. When the interference between different transitions is discarded, then ℱ_c(ω_1e+ω_1e')=δ_e,e' and  this term takes  the form of the collective Lamb shift of Ref. <cit.> for the corresponding laser excitation. 

In the dissipator's coefficient we also use the smoothening procedure by replacing Θ_ij^(Δ t) with ℱ_c(ω_i+ω_j). Moreover, we discard the temperature-dependent terms since they give negligibly small contribution at T=300 K and optical frequencies. 



 §.§ Photon-count signal


In order to study the effect of multilevel interference we determine the excitation spectrum S(δ_L) over the whole solid angle and as a function of the laser detuning δ_L=ω_L-ω_12. The excitation spectrum (or photon count signal) is defined as:
 
    S(δ_L) =  ∑_α,β∑_e,e'Γ_αβ^e e' FTr[ζ̂_1e'^βρ̂_ stζ̂_1e^α †] ,
 
and it is calculated for the steady-state density matrix ρ̂_ st, which is the solution of Eq. (<ref>) at eigenvalue zero, ∂_tρ̂_ st=0. In our simulations we take the parameters of the transition 2S→4P of Hydrogen. Specifically the ground state is |1⟩=|2s_1/2, F=0, M_F=0 ⟩, the excited states are
|2⟩= |4p_1/2, F=1, M_F=0 ⟩ and |3⟩= |4p_3/2, F=1, M_F=0 ⟩ and are illustrated in Fig. <ref>. Further details of the parameters are given in Appendix <ref>. The coefficients are calculated taking a coarse grained time scale Δ t=10^-11 sec (see Sec. <ref> for the analysis of the dependence of the results on the choice of the coarse-graining time scale). For further details we refer the reader to the discussion at the end of this section. We note that, for the level scheme which breaks rotational symmetry, the excitation spectrum of a single emitter exhibits non-vanishing shifts even after integration over the whole solid angle <cit.>.

 	

Figure <ref> displays the photon count signal (cyan line) for a given value of the laser intensity and as a function of the laser detuning δ_L for two interatomic distances: (a) R=0.01  μ m  and (b) R=0.1  μ m. These shall be compared with the wavelength λ_12=2π c/ω_12=0.468 μ m such that (a) corresponds to k R≃ 0.13 and (b) to kR≃ 1.3. The orange line gives the signal obtained when one artificially sets the multilevel interference effects to zero (corresponding to setting Θ_ij^(Δ t)→δ(ω̅_i-ω̅_j), namely Δ t→ 0). The mismatch between the cyan and the orange superradiant peaks is caused by the cross-interference terms.

We start with discussing the case R=0.1  μ m, when the interatomic distance is of the order of the wavelength. In this case the photon count signal is dominated by the photon count signal of the individual atoms, the peak maxima are at the frequency of the atomic levels, there are no evident features which could be attributed to superradiance and/or subradiance. Here, the inclusion of cross-interference terms gives rise to a slightly visible discrepancy between the two curves in the frequency interval between the two peaks.  When decreasing the interatomic distance to R=0.01  μ m the spectroscopic lines are splitted into the sub- and superradiant components. The frequency gap between the peaks of the sub- and superradiant components is given by the corresponding diagonal frequency shifts of Eq.  (<ref>). In the next section we determine the line shifts one extracts by analysing these spectra.
 



 §.§ Line shifts due to cross interference


In order to quantify the effect of the cross-interference terms, we determine the line shifts δω_j due to the multilevel interference. We focus on the lines of the superradiant states and extract the shift 
	
    δω_j= 1/2π(x_j^'-x_j) ,

where the quantity x_j^'  (with j=2,3 for |1⟩→ |j⟩) is extracted from the photon-count signal calculated using the master equation (<ref>). The frequency x_j, instead, is obtained by artificially setting all multilevel interference terms to zero, namely, by setting Θ_ij^(Δ t)→δ(ω̅_i-ω̅_j) in the coefficients of Eq. (<ref>). Thus, the frequency x_j includes also the collective Lamb shift.  The line shifts we report are determined  from the photon count signal as a function of the interatomic distance by taking the limit of vanishing Rabi frequencies, and are extracted by fitting the photon count signal using the following function, which is the sum of two Lorentzian curves:
    
    S^LL(x) = a_2/πb_2/2/(x-x_2)^2+(b_2/2)^2 +
        
    
        + a_3/πb_3/2/(x-ω_0-x_3)^2+(b_3/2)^2 ,
 
and ω_0=2πν_0 is the frequency gap between state |2⟩ and |3⟩ and is given in Appendix <ref>. It discards the presence of the subradiant peaks, whose magnitude becomes very small at low Rabi frequencies (for instance, for Rabi frequencies that are 1% the natural linewidth the magnitude is approximately 10^-3,10^-4 smaller than the superradiant ones). Nevertheless, these signals are generally different from zero and give rise to a systematic error in determining the line shift of the superradiant resonance. We remark that the choice of the fitting function is not optimal: In fact, Eq. (<ref>) corresponds to the spectroscopic signal due to the sum of two independent decay processes, and does not properly catch the features due to interference.
Indeed, the data in Fig. <ref> shows that the curves are more similar to Fano-like profiles. Previous studies showed that the excitation spectra of optically dense (homogeneously-broadened) media differ from Lorentz resonances <cit.>. Our choice is thus not going to be a reliable estimate of the shifts induced by multilevel interference. We expect, nevertheless, that it allows us to gain insight into their order of magnitude. 

Figure <ref> shows the line shifts as a function of the interatomic distance: at sufficiently short distances the shifts are significantly larger than the ones predicted for a single emitter and above the systematic error, due to discarding the subradiant peaks and illustrated by the dashed lines. The line shifts tend to increase the frequency gap between the two excited states as R→ 0, while for R→∞ they converge to the values indicated by the dashed lines, that are the shifts we calculate for the case of a single artificial emitter composed by three levels. 
   

We now argue that the observed shifts are due to quantum interference between the processes illustrated in Fig. <ref>. For this purpose we analyse the shifts by considering two artificial cases:  (i) The single-atom cross interference, in which  we only consider the scattering processes displayed in Figures <ref>a) and <ref>b). This corresponds to set Δ_ij^12=0 in  (<ref>) and Γ^ij_12=0 in  (<ref>) for i≠ j. (ii) The interatomic cross interference, in which we discard  scattering processes displayed in Figure <ref>a) and we keep the others. In this case we set Δ_ij^αα=0 and Γ^ij_αα=0 for i≠ j.  We further separately analyse the effect of the cross-damping terms (namely, the terms of the master equations where multilevel interference appears in the dissipator) and of the cross-shift terms (where multilevel interference appears in the Hamiltonian (<ref>)).

We first study the impact of the cross-damping terms versus R and artificially set all terms Δ_ij^αβ=0 with i j in Hamiltonian (<ref>).
Figure <ref>a) represents the results when we include  the cross-damping terms (i) only in the single-atom dissipator (intratomic, α=β), (ii) only in the interatomic dissipator (interatomic, α≠β) and (iii) when we consider both intratomic and interatomic cross-damping terms. In the case (i) the shifts due to the single-atom cross-damping terms at large distance oscillate around a magnitude of ∼ 100 Hz. In the case (ii) the line shifts vanish for R→∞. For vanishing distances the line shifts (i) and (ii) converge to a similar value. The total contribution of the intra- and interatomic cross-damping terms is not additive, as visible when comparing these curves with the ones obtained including both kinds of cross-damping terms.  
Figure <ref>b) displays the impact of the cross-shift terms on the line shifts after artificially setting all terms Γ_ij^12=0 in the dissipator (<ref>). Over  the interval of distances R=[0.1,1] μ m.  The total line shift has some oscillatory behaviour which tends to the single-atom result as R increases. At small R the cross-shift terms become dominant and tend to increase the frequency gap between the spectroscopic lines. 

The  behaviour at short distances is diplayed in Fig. <ref>. Here it is evident that the cross-shift terms are responsible for large shifts of the lines. Below R=48 nm (which corresponds to R∼λ/10) the shift of the line |1⟩→|3⟩ increase rapidly to the magnitude of 0.6 MHz, which starts to be comparable with the natural linewidth for optical transitions. 
 
 

       




 §.§ About the coarse graining time scale



The use of the coarse graining master equation allows one to derive ab initio a master equation fulfilling the Lindlad form and yet systematically including the cross-interference terms. The drawback is the explicit dependence on the coarse graning time, which becomes visible in the  functional form of Θ_ij^(Δ t), Eq. (<ref>), and which multiplies all coefficients for i≠ j. 
We note that this function determines the frequency window, for which the interference of two parallel dipolar transitions give rise to relevant contributions to the dynamics. 

One striking property is that Θ_ij^(Δ t) gives rise to strong oscillations of the coefficients with Δ t. The oscillations are majorly due to the sharp time intervals over which the dynamics has been divided and could be eliminated by introducing a smoothening, for instance by taking a Gaussian function of width Δ t and calculating the convolution <cit.>

    Θ_ij^(Δ t)→ℱ_c  (ω_i+ω_j)= ∫_0^∞ dx Θ_ij(x) e^-x^2/Δ t^2/√(πΔ t/2) .
 
This smoothening procedure delivers the new damping coefficients

    Γ_αβ^ij  (F)=ℱ_c   (ω_i+ω_j)D⃗_i^α*·D⃗_j^ β/2πħ c^3ω_ij^3  F_αβ^ij(k_ij) ,

which preserves the Lindblad form of the density matrix. Similarly we obtain the cross-coupling Hermitian terms after the smoothening.

Even after this smoothening, the coefficients of the master equation still depend on the choice of Δ t. For the master equation to be valid, their value shall be independent on the specific choice of Δ t over an interval of value. A rigorous lower bound for Δ t can be found by imposing the positivity of the Lindblad equation, as discussed in Ref. <cit.>. An heuristic approach is based on identifying the coarse-grained time for which the scattering properties are stable over several orders of magnitude, such that τ_R≪Δ t and Δ t is smaller than the smallest rate of the system dynamics. Figure <ref> shows the line shifts for different values of the coarse-graining time. 
The results do not vary over the interval of values of Δ t, over which we expect that the time-scale separation ansatz holds. They start to appreciably vary for Δ t > 10^-10 sec, and thus when Δ t becomes comparable with the natural lifetime of the excited states, which is here of the order of 10^-8 sec.








   	



§ CONCLUSIONS
 

In this work we have presented the systematic derivation of a master equation for an optically dense medium, which is composed by multilevel emitters. The master equation fulfils the Lindblad theorem <cit.> and includes the effect of interference between transitions which have parallel dipoles. This interference is induced by vacuum effects and gives rise to additional terms in the dissipator and Hamiltonian which can mutually interfere and whose strengths depend on the mean interparticle distance. 

We have provided a numerical example where we have applied our master equation to two identical emitters each consisting of  two parallel dipoles with a common ground state. We have shown that, even if the dipoles are not resonant, vacuum induced interference gives rise to measurable effects in the excitation spectrum. We have verified that the magnitude of the shifts depend on the ratio between the frequency gap between the interfering dipoles and their average linewidth and increases as this ratio decreases <cit.>, they become more evident when the interparticle distance decreases and emerge from the interplay of the interference between parallel dipoles of a single emitter and of the two emitters. Moreover, for realistic configurations the photodetection signal depends on the angle of emission and can be larger for certain directions <cit.>. 

Future work shall focus on alkali or alkali-earth metal atoms, consider the full sublevel structure and analyse the spectrum at different detection angles. A more accurate choice of the fitting functions shall provide a better estimate of the line shift due to multi-level interference <cit.>. This model, moreover, can be extended to Rydberg transitions <cit.>, where the multilevel interference is expected to be more prominent <cit.>, and to molecules <cit.>. 

The master equation here derived  can be extended and applied to studying propagation of quantum light in superradiant media and confined geometries <cit.>. By means of the input-output formalism <cit.> one can extract from our model the coherence properties of the scattered light and analyse the effect of vacuum induced interference on field- and intensity-intensity correlation functions.  Future studies will analyse its prediction  on light transport in a disordered medium <cit.> and in an ordered array of emitters <cit.> for level configurations where vacuum-induced interference is expected to be relevant.

The authos are grateful to Guido Pupillo, Johannes Schachenmeyer, and to the ITN Network ColOpt members for scientific discussions. We thank Andreas Buchheit for help in the first stages of this project and Anette Messinger for discussions and for careful reading of this manuscript. Funding by the EU ITN Network ColOpt and by the German Research Foundation (DFG, Priority Programme No. 1929, GiRyd) is gratefully acknowledged. 
    
   



§ DERIVATION OF THE BORN-MARKOV MASTER EQUATION IN THE COARSE GRAINING FORMALISM



The second integrand on the right-hand side of Eq. (<ref>) is here reported after applying the Born approximation:

    Λ_2^α,β(T,τ)ρ̃(τ_-)   =   ∑_i,jC̅_ij^αβ(τ)([σ̃_j^β(τ_-) ρ̃(τ_-), σ̃_i^α(τ_+)].
    
     
          .+[σ̃_i^α(τ_+),ρ̃(τ_-)σ̃_j^β(τ_-)]) +H.c. ,


where τ_±=T±τ. Subscript i labels a pair of level coupled by a non-vanishing dipole moment: i≡ i_1,i_2 with dipole moment d⃗_i^α=  _α⟨ i_1|d⃗|i_2⟩_α. The function C̅_ij^αβ(τ) specifically reads

    C̅_ij^αβ(τ)   =   ∑_λ(g_i^λg̅_j^λ (n(ω_λ,T)+1) e^- iω_λτ e^ ik⃗_λ·(R⃗_α-R⃗_β).
          .+g̅_i^λ  g_j^λ n(ω_λ,T) e^ iω_λτ e^- ik⃗_λ·(R⃗_α-R⃗_β)) ,

where n(ω,T)=1/[exp(ħω/k_BT)-1] is the mean photon number at frequency ω and temperature T and the sum over the modes is bounded by the cutoff frequency ω_cut. In the continuum limit it is given by the expression

    C̅_ij^αβ(τ)   →   ∫_0^ω_cutd ω/(2π)^2ħ c^3 ω^3([1+ n(ω,T) ]e^-iωτ. 
          .+ n(ω,T) e^iωτ)  F̅^ij(k,R⃗_αβ) .



Assuming the Born-Markov approximation, we can write ρ̃(T-τ)≈ρ̃(T) in Eq. (<ref>) <cit.>. We also note that, consistently with the Markov approximation, ρ̃(T) is essentially constant over the interval of integration [t,t+Δ t] of the variable T. We then set ρ̃(T)= ρ̃(t̅) with t̅=t+Δ t/2. Using that σ̃_j^β(t_1)= e^ iω_j(t_1-t̅)σ̃_j^β(t̅), we first rewrite Eq. (<ref>) as

    Λ_2^α,β(T,τ)ρ̃(τ_-)   ≈   ∑_i,j𝒞_ij^αβ(T,τ)
    ([σ̃_j^β(t̅) ρ̃(t̅), σ̃_i^α(t̅)].
    
     
          .+[σ̃_i^α(t̅),ρ̃(t̅)σ̃_j^β(t̅)]) +H.c. ,

where 

    𝒞_ij^αβ(T,τ)=C̅_ij^αβ(τ) e^ i(ω_i+ω_j)T e^ i(ω_i-ω_j)τ.

Using now that σ̃^α_j(t)=exp( iω_jt)σ_j^α in Eq. (<ref>), the time integrals take the form:

    1/2Δ t∫_-Δ t^Δ t dT e^± i(ω_i-ω_j) T/2∫_-Δ t^Δ t  dτ θ(τ)C_αβ(τ)  e^± i(ω_i+ω_j) τ/2
          =Θ_ij^(Δ t)∫_-Δ t^Δ t  dτ θ(τ)C_αβ(τ)  e^± i(ω_i-ω_j) τ/2 ,

where 

    Θ_ij^(Δ t)=sin((ω_i+ω_j)Δ t/2)/(ω_i+ω_j)Δ t/2 .

When the transition are in the optical range, this function selects secular terms. For this reason, in the following we restrict the sum to all pairs such that ω_i>0. 
The second integral is evaluated after approximating the extrema of integration by [-Δ t,Δ t]→[-∞,∞], which is consistent with the assumption that C(τ) decays to zero over time scales much shorter than Δ t.



§ PARAMETERS OF THE SIMULATION



The magnitude of the fine structure splitting for 4p state is taken to be ν_0≈1.367 GHz <cit.> and includes also the Hyperfine structure splitting and QED corrections. We neglect thermal effects: we set n(ω_1e)=0, which is a good approximation at room temperature T=300 K. Moreover, we take the following values for the radiative shifts: Δ_22^S=-2π× 1401.52 kHz for the state 4p_1/2 and Δ_33^S=2π× 1767.30 kHz for the state 4p_3/2 <cit.>. We then construct the atomic cross-shift term between the excited states using relation (<ref>): Δ_23^S= Δ_32^S=2π× 366.2 kHz using the relation between dipole moments of the corresponding transitions: d_12=(1/3)d_R, d_13=-(√(2)/3)d_R, where d_R is the radial integral d_R=⟨ 2s|r|4p ⟩ = 1.28 [a.u.]. The values for the natural line width are  γ_2=Γ_22 = 2π× 511 kHz for the state 4p_1/2 and γ_3≡Γ_33 =2π× 1022 kHz for the state 4p_3/2. All the cross-interference terms both for the dissipator and the Lamb shift were computed with the coarse-graining time Δ t = 10^-11 sec. The computational checks showed that the solutions of the master equation for the chosen system remains stable in this coarse-graining time region, see the Figure <ref> in Sec. <ref>.    



§ DETERMINATION OF THE LINE SHIFTS



Line shifts due to cross-interference for a single emitter
In this appendix we illustrate the procedure we apply in order to determine the line shifts due to the multilevel quantum interference terms. We provide the example of a single emitter, and we refer to it using the wording "single atom". However, due to the special structure we assume the emitter is not rotationally invariant, which changes the spectroscopic properties and gives rise to a global line shift when integrating the photon count signal over the whole solid angle. 

The line shifts for a single atom are defined as  
 
    Δ_j(g) = 1/2π(x_j^g - x_j^eigen),
    
where x_j^eigen is the eigenfrequency of the j^th transition  (j=1 for |1⟩→ |2⟩ and j=2 for |1⟩→ |3⟩) which we extract from the master equation when we set all multilevel interference terms to zero, and x_j^g is the line position obtained by fitting the photon count signal using the fitting function  (<ref>). The line shifts Δ_j(g) depends on the laser intensity and thus on the Rabi frequency g.  The line shift we identify corresponds to the limit:
 
    Δ_j =lim_g→ 0Δ_j(g).
    
  
 The limit is extracted from our numerical analysis: We evaluate it for decreasing values of g. We report the behaviour in Fig. <ref> a) when we set to zero the multilevel interference terms and b) for the full master equation.  The presence of the interference terms shifts both peaks to the magnitudes ± 195 Hz for vanishing laser intensity. 
 






Dicke:1954

R. H. Dicke, Phys. Rev. 93, 99 (1954).

Friedberg:1972
R. Friedberg, S: Hartmann, and J. Massanah, Phys. Rep. 7, 101 (1973).

GrossHaroche:1982  
	
	M. Gross and S. Haroche, Phys. Rep. 93, 301-396(1982).  

Bromley:2016
	 
Sarah L. Bromley, Bihui Zhu, Michael Bishof, Xibo Zhang, Tobias Bothwell, Johannes Schachenmayer, Travis L. Nicholson, Robin Kaiser, Susanne F. Yelin, Mikhail D. Lukin, Ana Maria Rey, and Jun Ye,
Nature Communications 7, 11039 (2016).

Grimes:2016  

David D.Grimes, Stephen L. Coy, Timothy J. Barnum, Yan Zhou, Susanne F. Yelin, and Robert W. Field,

Phys. Rev. A 95, 043818 (2017). 

Araujo:2016

M. O. Araujo, I. Kresic, R. Kaiser, and W. Guerin, Phys.Rev.Lett. 117, 073002 (2016).

Corman:2017

L. Corman, J.-L. Ville, R. Saint-Jalm, M. Aidelsburger, T. Bienaimé, S. Nascimbène, J. Dalibard, and J. Beugnon,
Phys. Rev. A 96, 053629 (2017).

Peyrot:2018

T. Peyrot, Y. R. P. Sortais, A. Browaeys, A. Sargsyan, D. Sarkisyan, J. Keaveney, I. G. Hughes, and C. S. Adams,
Phys. Rev. Lett. 120, 243401 (2018).

Jennewein:2018

Stephan Jennewein, Ludovic Brossard, Yvan R. P. Sortais, Antoine Browaeys, Patrick Cheinet, Jacques Robert, and Pierre Pillet, Phys. Rev. A 97, 053816 (2018).

Pellegrino:2014
J. Pellegrino, R. Bourgain, S. Jennewein, Y.R.P. Sortais, A. Browaeys, S.D. Jenkins and J. Ruostekoski, Phys. Rev. Lett. 113, 133602 (2014).

Orgiu:2015

E. Orgiu, J. George, J. A. Hutchison, E. Devaux, J. F. Dayen, B. Doudin, F. Stellacci, C. Genet, J. Schachenmayer, C. Genes, G. Pupillo, P. Samori, and T. W. Ebbesen,

Nature Materials 14, 1123-1129 (2015).	 

 Meiser:2009
D. Meiser, Jun Ye, D. R. Carlson, and M. J. Holland,

Phys. Rev. Lett. 102, 163601 (2009).

Norcia:2016
M. A. Norcia, M. N. Winchester, J. R. K. Cline, and J. K. Thompson, 
Science advances 2, e1601231 (2016);

M. A. Norcia, J. R. K. Cline, J. A. Muniz, J. M. Robinson, R. B. Hutson, A. Goban, G. E. Marti, J. Ye, and J. K. Thompson,

Phys. Rev. X 8, 021036 (2018).

Agarwal

G. S. Agarwal, Phys. Rev. A 4, 1791 (1971).

Lehmberg:1970 
	
R. H. Lehmberg, Phys. Rev. A 2, 883 (1970); 


Phys. Rev. A 2, 889 (1970).

Milonni:1974

P. W. Milonni and P. L. Knight,  Phys. Rev. A 10, 1096 (1974).

James:1993 Daniel F. V. James, Phys. Rev. A 47, 1336 (1993).   

Carmichael H. J. Carmichael,  An Open Systems Approach to Quantum Optics
(Springer-Verlag, Berlin, 1993).

Fleischhauer:1999

M.Fleischhauer and S. F. Yelin, Phys. Rev. A 59, 2427 (1999).

Yelin:Review

G. D. Lin and S. F. Yelin, Advances in Atomic, Molecular, and Optical Physics 61, 295-329 (2012).

Breuer  H.-P. Breuer and F. Petruccione, The theory of open quantum systems (Oxford University Press, Oxford, 2002).

Zhu:2016 
B. Zhu, J. Cooper, J. Ye, and A. M. Rey, Phys. Rev. A 94, 023612 (2016).    

Lin:2012 G.-D. Lin and S. F. Yelin, Phys. Rev. A 85, 033831 (2012).


Milonni:1976    P. Milonni, Phys. Rep. 25, 1 (1976).

Cardimona:1983
D. A. Cardimona, M. G. Raymer, and C. R. Stroud, Jr., J. Phys. B 15, 55 (1982); D. A. Cardimona and C. R. Stroud, Phys. Rev. A 27, 2456 (1983).

Ficek Z. Ficek and S. Swain, Quantum Interference and Coherence: Theory and Experiments (Springer, New York, 2005).

Kiffner:2010
M. Kiffner, M. Macovei, J. Evers, and C. H. Keitel, Progress in Optics 55, 85 (2010). 

Berman:1998

P. R. Berman, Phys. Rev. A 58, 4886 (1998).

Buchheit:2016  A. A. Buchheit and G. Morigi, Phys. Rev. A 94, 042111 (2016).

Horbatsch

A. Marsman, M. Horbatsch, and E. A. Hessels,
Phys. Rev. A 96, 062111 (2017).   

Udem:2019

T. Udem, L. Maisenbacher, A. Matveev,  V. Andreev,  A. Grinin,  A. Beyer, N. Kolachevsky, R. Pohl,  D. C. Yost, and  T. W. Hansch, Ann. Phys. (Berlin) 531,  1900044 (2019). 
 
Yost:2014

D. C. Yost, A. Matveev, E. Peters, A. Beyer, T. W. Hansch, and Th. Udem,
Phys. Rev. A 90, 012512 (2014).

Lidar:2001  D. A. Lidar, Z. Bihary, and K. B. Whaley, Chem. Phys. 268, 35 (2001).

Majenz:2013
C. Majenz, T. Albash, H.-P. Breuer, and D. A. Lidar, Phys. Rev. A 88, 012103 (2013).

Cohen-Tannoudij
C. Cohen-Tannoudji, J. Dupont-Roc, and G. Grynberg, Atom-Photon Interactions: Basic Processes and Applications (Wiley-VCH, 2004).

Rubio:2018
V. Rokaj, D. M. Welakuh, M. Ruggenthaler,  and A. Rubio, J. Phys. B: At. Mol. Opt. Phys. 51, 034005 (2018).

Ingold
G. L. Ingold, in Coherent Evolution in Noisy Environments, edited by A. Buchleitner and K. Hornberger (Springer, 2002), pp. 1-53.

Hepp
K. Hepp and E. H. Lieb, Ann. Phys. 76, 360 (1973). 

Abramowitz-Stegun
M. Abramowitz and I. A. Stegun: Handbook of Mathematical Functions (National bureau of standards and technologies, June 1964).

Scully:2009
M. O. Scully, Phys. Rev. Lett. 102, 143601 (2009).

Roehlsberger:2011
R. Rölhsberger, K. Schlage, B. Sahoo, S. Couet, and R. Roeffer, Science 328, 1248 (2010).

Scully:book
M. O. Scully and M. S. Zubairy, Quantum Optics (Cambridge University Press, 1997).

Farina:2019

D. Farina and V. Giovannetti,
Phys. Rev. A 100, 012107 (2019).

Giannelli:2019 L.  Giannelli,  T.  Schmit,  and  G.  Morigi  Physica Scripta 94, 014012  (2018).

Putnam:2016
G. Putnam, G. D. Lin, and S. F. Yelin, preprint arXiv:1612.04477 (2016).


Jentschura:1997

U. D. Jentschura, G. Soff, and P. J. Mohr, 
Phys. Rev. A 56, 1739 (1997).


Lahaye:2017

S. de Léséleuc, D. Barredo, V. Lienhard, A. Browaeys, and T. Lahaye,

Phys. Rev. Lett. 119, 053202 (2017).

Zelevinsky:2015 

B. H. McGuyer, M. McDonald, G. Z. Iwata, M. G. Tarallo, W. Skomorowski, R. Moszynski, and T. Zelevinsky,
Nat. Phys. 11, 32 (2015).

Zhou:2017
Yao Zhou, Zihao Chen, and Jung-Tsung Shen, Phys. Rev. A 95, 043832 (2017).

Chen:2018

Zihao Chen, Yao Zhou, and Jung-Tsung Shen
Phys. Rev. A 98, 053830 (2018).

Zhou:2020

Yao Zhou, Zihao Chen, and Jung-Tsung Shen
Phys. Rev. A 101, 043831 (2020).

Collett:1985
C.W. Gardiner and M. J. Collett, Phys. Rev. A 31, 3761 (1985).

Javanainen:2014
J. Javanainen, J. Ruostekoski, Y. Li, and S.-M. Yoo,
Phys. Rev. Lett. 112, 113603 (2014).

Cottier:2018

F. Cottier, R. Kaiser, and R. Bachelard,
Phys. Rev. A 98, 013622 (2018).

Longo:2015

P. Longo, C. H. Keitel, and J. Evers,
Scientific Reports 6, 23628 (2016).

Facchinetti:2016

G. Facchinetti, S. D. Jenkins, and J. Ruostekoski,
Phys. Rev. Lett. 117, 243601 (2016).

Kolachevsky  N. Kolachevsky, A. Beyer, L. Maisenbacher, A. Matveev, R. Pohl, K. Khabarova, A. Grinin, T. Lamour, D. C.
	Yost, T. W. Haensch, and Th. Udem,  AIP Conference Proceedings 1936, 020015 (2018). 




