Theory of multiple quantum coherence signals in dilute thermal gases

Manifestations of dipole-dipole interactions in dilute thermal gases are difficult to sense because of strong inhomogeneous broadening. Recent experiments reported signatures of such interactions in fluorescence detection-based measurements of multiple quantum coherence (MQC) signals, with many characteristic features hitherto unexplained. We develop an original open quantum systems theory of MQC in dilute thermal gases, which allows us to resolve this conundrum. Our theory accounts for the vector character of the atomic dipoles as well as for driving laser pulses of arbitrary strength, includes the far-field coupling between the dipoles, which prevails in dilute ensembles, and effectively incorporates atomic motion via a disorder average. We show that collective decay processes -- which were ignored in previous treatments employing the electrostatic form of dipolar interactions -- play a key role in the emergence of MQC signals.


Introduction
Dipole-dipole interactions enable transport and collective phenomena like dipole blockade in Rydberg gases [1], coherent backscattering (CBS) of light in cold atoms [2], resonant energy transfer in photosynthetic complexes [3,4], interatomic Coulombic decay [5,6], sub- [7] and superradiance [8], and multipartite entanglement [9,10]. A typical manifestation of dipolar interactions are collective displacements of the energy levels of a many-body system [11]. However, in thermal ensembles the observation of these is hindered by motionally induced Doppler shifts. As a result, a mean-field description is often sufficient [12], and the properties of the system can be derived from single-body responses.
Such picture is however not truly complete, since signatures of light-induced dipolar interactions indeed can be observed in so-called multiple quantum coherence (MQC) signals [13] sensing coherent superpositions of the ground and collective excited states of several particles. While originally MQC signals were measured [13] via ultrafast two-dimensional electronic spectroscopy [14], which allows one to probe directly the nonlinear susceptibilities responsible for MQC signals, more recently these signals were extracted from fluorescence emitted by dilute thermal atomic gases [15].
The interpretation of fluorescence-detected MQC signals sparked a debate. On the one hand, it was pointed out [16] that fluorescence as radiated by dilute atomic samples can usually be well understood by considering a system of independent emitters, which corresponds to a mean-field description. On the other hand, the theoretical follow-up [17] to the experiment [15] predicted vanishing double-quantum coherence in the absence of interactions, inferring that the latter are essential to establish collective behaviour. A subsequent experiment [18] further reported significant quantitative deviations from the results based on the model of [17]. Eventually, it was suggested [18] that the double-quantum coherence signals could only be understood by going beyond a two-body model-into the realm of many (N > 2) interacting particles.
The aim of our present article is to settle this debate by elucidating the role of dipole-dipole interactions in MQC signals derived from fluorescence-based measurements. These signals are emitted when atoms are pumped and probed by a series of femtosecond laser pulses, and, until now, their interpretation [13,15,17] relies on two orthodox premises of multi-dimensional nonlinear spectroscopy [19,20]: (a) the laser-atom interactions are sufficiently weak and can be accounted for perturbatively, using double-sided Feynman diagrams, and (b) the dipole-dipole coupling scales as ∼ R −3 with the interatomic distance R, and corresponds to the familiar electrostatic [21], or near-field interaction associated with virtual photons [22].
The theory of MQC signals we present here abandons both these premises. Since fluorescence is an incoherent process, resulting from spontaneous decay of excited emitters (be them independent or interacting), its dynamics can be described with the formalism of open quantum systems [23]. More specifically, we adapt a quantum-optical master equation approach previously used to faithfully model CBS of light by cold atoms [24][25][26][27]. Unlike the diagrammatic approach, our method (i) treats the laser-matter interaction non-perturbatively, and (ii) considers the exact form of dipolar coupling, which, beyond the electrostatic part, includes the exchange of photons via its radiative, or far-field component, scaling as ∼ (k 0 R) −1 , with k 0 = 2π/λ 0 the wave number (λ 0 the wave length) resonant with the individual atomic scatterers' electronic transition frequency ω 0 . This latter extension (ii) in particular ensures the inclusion of irreversible, collective radiative decay processes, which are ignored by the electrostatic modelling [21] of the dipolar coupling.
The radiative part of the dipole-dipole interaction prevails in dilute gases, where the mean interatomic separation is much larger than the resonant optical wave length. Since (k 0 R) −1 1 in this regime, we can treat this interaction perturbatively. Then the average fluorescence intensity can be represented by a series expansion in powers of (k 0 R) −1 , whose subsequent terms correspond (see figure 1) to single, double, etc. scattering contributions. We consequently interpret MQC components distilled from the total fluorescence signal as arising from the multiple scattering of real photons between distant atoms.
Below, we recall the basics of the phase-modulated spectroscopy that was employed to measure MQC signals in dilute thermal gases [15,18,28]. Thereafter, we present our formalism and derive 1QC (single quantum coherence) and 2QC (double quantum coherence) spectra that are in excellent qualitative and reasonable quantitative agreement with experiment [18]. Ultimately, we establish the crucial role of radiative dipole-dipole interactions in the emergence of 2QC signals, which are impossible to observe for independent atoms-in agreement with [17], and we show that a system of two dipole interacting atoms is sufficient for modelling 1QC and 2QC signals-as a consequence of double scattering of real photons-in a dilute thermal gas. We especially highlight collective decay processes as an essential Schematic illustration of the different time scales and of the associated Liouvillians, for a single experimental cycle of duration T cyc . Two ultrashort Gaussian laser pulses with the same width σ arrive at times t 1 and t 2 , separated by a delay τ σ. After the second pulse, fluorescence is monitored during the interval t fl ∼ 10γ −1 τ , until all atoms decay into the ground state. The fluorescence detection time t fl is much longer than the thermal decoherence time τ th > τ, beyond which the atomic motion cannot be ignored. During the ultrashort laser-atom interaction, the dynamics of the atomic operators Q(t) in (12) is governed by the Liouvillian L L . Instead, L int and L γ generate the dynamics during the intervals τ and t fl .
part of the dipolar interactions, without which qualitative agreement with the experiment cannot be achieved.

Phase-modulated spectroscopy
In figure 1 we sketch the spectroscopy experiment used to extract the MQC signals. A thermal ensemble of alkali atoms is excited near-resonantly [15,18] by two subsequent laser pulses injected at times t 1 and t 2 > t 1 , respectively [29]. The fluorescence subsequently emitted by the sample is collected by a photodetector in a direction orthogonal to the exciting field's wave vector k L . This excitation and detection procedure is repeated for M cyc 1 cycles of duration T cyc , with the two pulses garnished with distinct phase tags φ 1 and φ 2 , respectively, which are linearly ramped over subsequent cycles, with distinct rates w 1 = w 2 . These slowly evolving phase tags will turn out to be instrumental to distill the MQC signals from the total fluorescence. T cyc is chosen to significantly exceed the atomic radiative lifetime of a few tens of nanoseconds [30], such that the atoms decay to the ground state [15,18] within a single cycle, and each of the M cyc cycles can be considered in isolation (see figure 2).
To resolve the details of the excitation scheme, we focus on the mth cycle (0 m < M cyc ), which begins at τ m = mT cyc . For positive delays t ∈ [0, T cyc ] with respect to τ m , and in the frame co-rotating with the laser carrier frequency ω L , the injected field has the time-dependent amplitude with the same Gaussian envelope E L (t ) = E 0 exp(−t 2 /2σ 2 ) of both incoming pulses, amplitude E 0 and duration σ. Each phase tag φ j = w j (τ m + t) can be approximated by φ j ≈ w j τ m for modulation frequencies w j ω L [29], 4 and is imprinted on the respective pulse by transmission through an acousto-optical modulator oscillating at frequency w j . Both pulses are collinear, with wave vector k L , but may have different polarizations ε (j) L : in the following we consider the cases ε (1) L = ε (2) L =x, as well as ε (1) L =x and ε (2) L =ŷ. In a typical experiment, σ ∼ 50 fs τ = t 2 − t 1 ∼ 10 ps T cyc ∼ 300 μs. As a consequence of the phase tagging of the driving field, the transient fluorescence intensity Ik(τ , t fl , φ 1 , φ 2 ) in the directionk ⊥ k L exhibits a modulation-in its dependence on τ m through φ 1 and φ 2 -at the frequency difference w 21 = w 2 − w 1 , and at higher harmonics thereof, with the κth harmonic identifying the κth quantum coherence (κQC) [15]. So does the time-integrated (over the mth cycle) fluorescence signal as collected in the far field over times t 2 + t fl , with 0 t fl T cyc − t 2 γ −1 , γ the natural line width of the atomic excited states, faithfully approximated bȳ When incident on a photodetector, the integrated signal generates a photocurrent proportional tō Ik(τ , φ 1 , φ 2 ), which (independently) depends on the pulse delay τ , and-via the phase tags φ 1 , φ 2 (see above)-on the cycling time τ m . Frequency analysis of the photocurrent with respect to τ m , to distill the targeted order κ of MQC, and with respect to τ , to extract the κQC signal strength at the characteristic frequency κω 0 completes the experimental protocol, as elaborated in section 4.4.5 below.

Hamiltonian
Typical experiments to sense MQC signals are carried out at room temperatures with atomic clouds of low densities ρ ∼ 10 7 -10 11 cm −3 [15,18]. These are very dilute gases for which ρk −3 0 1, with the atoms typically in the far-field of each other. Given that collisional broadening (see, e.g., [31]) in such ensembles is at least one order of magnitude narrower than the natural line width [28], we can neglect collisions between the atoms to a good approximation.
Our model considers N identical, initially unexcited atoms at random positions r α (t) = r α0 + v α t, where r α0 and v α are, respectively, the initial coordinate and velocity of atom α = 1, . . . , N. We assume a homogeneous and isotropic distribution of the initial coordinates r α0 (which implements the above far-field configuration), while the Cartesian components of the velocities v α can be drawn from the Maxwell-Boltzmann distribution corresponding to the temperature of the cloud. At typical room temperatures the atomic velocity is ∼ 10 2 m/s; hence, the atomic momentum is much larger (for 87 Rb atoms probed in [18], by about three orders of magnitude) than the photon momentum. Under this condition we can neglect the recoil effect and treat the atomic motion classically.
Then, the total Hamiltonian of the atomic ensemble interacting with an electromagnetic bath represented by an infinite set of quantized harmonic oscillators in their vacuum state, and with the classical, near-resonant pump and probe pulses (1) reads where defines the free atoms' Hamiltonian, with ω 0 ≈ ω L , and D † α (D α ) the atomic raising (lowering) operators incorporating arbitrary degeneracies of the ground and excited state sublevels.
• H F = ks ω k a † ks a ks describes the contribution from the quantized, free electromagnetic environment, with a † ks (a ks ) the creation (annihilation) operator of a photon with wave vector k and polarization s. • H AL and H AF mediate the coupling between the atomic degrees of freedom and, respectively, the injected laser field and the quantized field. We employ the electric dipole approximation to write the Hamiltonians H AL and H AF as where d is the atomic dipole matrix element [31], which is assumed to be real, without loss of generality, E L (r α (t), t) is given by (1) and E (+/−) (r α (t)) is the positive/negative-frequency part of the quantized field given by with 0 the free space permittivity, V the quantization volume, and ε ks the polarization vector. In equation (4) for H AL , we have additionally applied the rotating-wave (resonant) approximation, whereas counter-rotating terms must be kept in equation (5) for H AF to obtain correct collective level shifts [32][33][34]. MQC signals were observed [18] in 87 Rb vapours, with the probed transitions endowed with fine and hyperfine structures, leading to rich multi-component features of the signals. We here ignore these features and specialize to atoms equipped with a J g = 0 → J e = 1 transition, which involves a non-degenerate ground state and three sublevels in the excited state (see figure 1). While this simplifies our model with respect to the 87 Rb experiment [18], it still incorporates the experimentally important excited state degeneracy (a feature already employed earlier [24][25][26] to account for the relevant optical transition 1 S 0 → 1 P 1 in experiments on the CBS of light by cold strontium atoms [35]). A crucial consequence thereof is the presence of levels that are not driven by the polarized laser fields, but still can be populated by the dipole-dipole interaction. Such processes can be identified using different pump-probe polarizations and observation directions [18]. Unlike previous treatments based on two-level atoms without any internal degeneracy [13,16,17,[36][37][38], our model is adequate to capture the consequences of this redistribution of energy between angular momentum states.
The atomic lowering operator associated with a J g = 0 → J e = 1 transition is given by where σ α lk = |l α k| α , the excited state sublevels |2 , |3 , |4 (see figure 1) have magnetic quantum numbers −1, 0, +1, and couple to field Note that, in equation (3), we omit the kinetic energy of the atomic centre-of-mass motion, as well as atomic contact interaction terms [39], which must be taken into account in dense atomic systems, where ρk −3 0 1, but can be safely neglected in dilute ensembles [40].

Fluorescence intensity
The experimentally detected total signal is expressed through the quantum-mechanical expectation value of the intensity operator, averaged over atomic configurations, where · · · = Tr{. . . (0)}. The initial density operator (0) of the atom-field system factorises into the ground state of the atoms and the vacuum state of the field. The positive frequency part E (+) k (t) of the operator for the electric far-field can be expressed through the atomic lowering operator of the source atoms [41,42], 6 where f = ω 2 0 d/(4π 0 c 2 R d ), c is the speed of light, R d is the distance from the center of mass of the cloud to the detector, and k is the wave vector of the scattered light. In writing (9), we assumed that R d r αβ (t) = |r α (t) − r β (t)| for any α, β. Finally, · · · conf stands for the configuration (or disorder) average. It results from the random initial positions of the atoms within the cloud, and from their thermal motion during the fluorescence detection time t fl . We will treat the thermal motion effectively, through configuration averaging which involves an integration over the length and direction of the vectors r αβ connecting pairs of atoms (see section 4.4.4 for more details), with mean interatomic distancer. 7 Plugging (9) into (8), with atomic levels from (7), and k = kx, we obtain e ikx·r αβ σ α z1 σ β 1z + σ α y1 σ β 1y conf (10) and Iŷ follows upon replacingx →ŷ, y → x in (10). Note that, in the above equation, we dropped the time argument and introduced a new (degenerate) atomic excited manifold basis,

Master equation
Evaluation of the fluorescence intensity (10) includes certain time-dependent atomic dipole correlators like σ α z1 σ β 1z . For immobile atoms it is technically involved, though fundamentally straightforward, to derive a Lehmberg-type master equation for an arbitrary atomic operator [43] (or, equivalently, for the atomic density operator [41]) under the standard Born and Markov approximations [23]. Such a master equation contains a Liouvillian describing the dipole-dipole coupling, which is a function of the distance r αβ between the interacting particles. In a thermal ensemble, the interatomic distance is time-dependent, such 5 Equivalently, we denote the unit polarization vectors in the Cartesian basis byê x ,ê y , andê z =ê 0 =ẑ. 6 Strictly speaking, due to the coupling described by the Hamiltonian H AF (see (5)), any atomic and field operators are atom-field operators at t > 0. 7r 0.554ρ −1/3 [20], which for a particle density ρ 10 7 -10 11 cm −3 [15,18] is about 100-1 μm. that employing the standard approximations in this regime requires justification. As it turns out, the Born and Markov approximations may also prevail for moving atoms-provided that the Doppler shift k 0 v α (v α = |v α | and k 0 = ω 0 /c) amounts to a negligibly small phase shift k 0 v α r αβ (t)/c 1 during the propagation time r αβ (t)/c of the interaction between the atoms [34]. In other words, the atoms can be considered 'frozen' on the time scale r αβ (t)/c. This condition is well satisfied, for example, for r αβ (t) =r ∼ 10 μm and a Doppler shift k 0 v α ∼ 600 MHz, since k 0 v α r αβ (t)/c ∼ 2 × 10 −5 . We therefore obtain the same form of the dipole-dipole interaction as for fixed atoms, where r αβ is substituted by r αβ (t) [34].
Thus, an arbitrary time-dependent atomic dipole correlator can be obtained from the quantum master equation [42,43] (written in the interaction picture with respect to the free-atom Hamiltonian H 0 ) 8 : where Q is an arbitrary atomic operator (on the space of N atoms) averaged over the free field subsystem, with L α γ the Lindblad superoperator describing the atomic relaxation due to the electromagnetic bath; and L int = N α =β=1 L αβ , with L αβ the Liouvillian describing the dipole-dipole interactions mediated by the (common) electromagnetic bath. The Liouvillians read [24,25,33,42] with the laser-atom detuning incorporating the Doppler shift k L · v α of the atomic transition frequency factored out from the laser field (1), the spontaneous decay rate γ (identical for all excited state levels), and with the dipole-dipole interaction tensor ← → T (k 0 r αβ (t),n) which, withn = r αβ /r αβ = (sin θ cos φ, sin θ sin φ, cos θ) and ξ(t) = k 0 r αβ (t), reads In this equation, 1 is a 3 × 3 unit matrix andnn is a 3 × 3 dyadic. Since we assume linear polarization of the incoming laser pulses, it is most convenient to express the matrix elements of ← → T (ξ(t),n) in the Cartesian basis,

Near-versus far-field form of the dipole-dipole interaction
Note that the near-field term in equation (17), scaling as ξ −3 (we drop for brevity the time argument), is the ξ 1 asymptote of the dipole-dipole coupling tensor (17) [43], and is associated with the electrostatic interaction. Because of its (slow) power-law decay with the distance, this interaction is considered to be pertinent at any interatomic separations [20], while the far-field contribution in equation (17) (see below) complicates the interpretation of nonlinear signals [44]. It is therefore the electrostatic interaction that was employed in all previous theoretical treatments dealing with MQC signals in dilute atomic gases [13,17,[36][37][38], where atomic angular momentum was also ignored.
To identify the associated form of dipolar interactions, we consider a simple two-level atomic structure rather than our J g = 0 → J e = 1 transition, i.e. we keep only level |3 in the excited state manifold. The tensor (17) then reduces to its matrix element ← → T zz corresponding to the magnetic number zero, which, in the near-field limit ∝ ξ −3 , reads Recalling the definition γ = 4ω 3 0 d 2 /(4π 0 3 c 2 ) [31], and introducing the dipole moment operators d(σ 13 + σ 31 ) of the two-level systems, with d = dê z , equations (12), (16) and (18) lead to the following Heisenberg equation of motion for two atoms α and β: where we have identified the electrostatic dipole-dipole interaction V d-d . We can see from equation (19) that the interatomic potential V d-d generates purely Hamiltonian dynamics, which does not affect the ground and doubly excited states |1 α |1 β and |3 α |3 β , respectively, of the uncoupled system. It does, however, shift both uncoupled single-excitation eigenstates |3 α |1 β and |1 α |3 β . In quantum electrodynamics, these shifts are interpreted as originating from the exchange of virtual photons [39].
In contrast to the electrostatic part (19) of the interaction, its radiative part accounts for the exchange of real, or transverse, photons [39], and generates both, collective shifts and collective decay of the atomic excited states, associated, respectively, with the tensors (19), we here are dealing with dilute gases characterized by interatomic distances that are much larger than the optical wavelength. Consequently, in the far-field limit ξ 1, the tensor ← → T retains its real and imaginary parts, scaling as ξ −1 , by equation (17). In that limit, 1, and decay rates and level shifts can be deduced from which are both oscillating functions of ξ, with identical amplitude. This means that both terms contribute equally to the dynamics and, hence, the non-Hermitian dissipative terms generated by ← → Γ (ξ,n) cannot be neglected.
With the definition (20) of ← → Γ and ← → Ω , the interaction Liouvillian from (16) can be decomposed as This distinction between contributions to L αβ will be essential for our interpretation of MQC spectra in section 5 below. Let us therefore describe their action on specific operators Q in detail: (21)) transfers one excitation from atom β to atom α. When (D † α ) r (D β ) j acts on the projector Q = σ α 11 σ β ff from the left (where (f ∈ {x, y, z}), it transforms the latter into the operator Q = δ jf σ α r1 σ β 1f of electronic coherences of both atoms. This transformation amounts to the dynamical equatioṅ which couples the exited state populations to the electronic coherences.
The action of σ α r1 σ β 1j from the right can be described analogously. • Since the raising and lowering operators of different atoms in L (2) αβ (see (22)) act from the left and right side on an atomic operator, this Liouvillian performs the transformation to a two-atom subspace with an extra excitation. As a result, an operator Q = σ α 1j σ β r1 describing electronic coherences of both atoms is transformed into that of Zeeman coherences and excited state populations, . We can write this as a dynamical equation, Q = 2 fl ← → Γ fl Q , stating that Zeeman coherences and excited state populations of two atoms (embodied in Q ) are transformed to their electronic coherences (embodied in Q ) via a collective decay process.

Separation of time scales
Our approach to solve (12), with ingredients (13)-(17), analytically is based on the observation that the atomic dynamics during each experimental cycle (see figure 2) is characterized by four different time scales (recall section 2). The first, ultrashort time scale is set by the laser pulse length σ ∼ 50 fs (which enters through (1)). The interpulse delay τ ∼ 10 ps, which parametrizes (2), determines the second, short time scale. After the second pulse, atomic fluorescence is collected until the end of the cycle, at time T cyc (entering (2), via φ 1 and φ 2 , through τ m ), until all atoms relax into their ground states by emission of a photon. This stage of the fluorescence detection establishes the third time scale, which typically lasts for about four orders of magnitude longer than the short time scale-for several microseconds. Finally, within the long time scale, the typical thermal coherence time τ th [45] of the order of 10 ns defines the fourth time scale, 9 which demarcates the regime of the atomic dynamics beyond which the Doppler effect is not negligible.
Such separation of time scales motivates a simplified version of (12), which decomposes into terms relevant at distinct time scales. We thus completely ignore the incoherent dynamics generated by L γ and L int during that ultrashort time scale, while, on longer time scales, we solve (12) including all terms from (13) and (16) except for the Liouvillians L L which generate the dynamics induced by the (on these long time scales) 'kick'-like laser-atom interaction. The thus defined propagation over ultrashort and longer time scales is then interchanged twice, to cover a single experimental cycle of duration T cyc . This approach resembles a stroboscopic description of the dynamics familiar, e.g., from stroboscopic maps [46] or also from micromaser theory [47].
Furthermore, we need to blend in the atoms' thermal motion, which can be safely neglected on the ultrashort and short time scales. The atomic configuration set by r α0 is therefore considered frozen on these time scales. On the long time scale, however, atomic motion does affect the dipole-dipole coupling L int via the then time-dependent interatomic vectors r αβ (t) = r α0 − r β0 + (v α − v β )t. This leads to rapidly oscillating, r αβ (t)-dependent phase factors in certain terms of our solutions of (12), which average out for many repetitions of experimental cycles of duration T cyc , each of them being much longer than the thermal coherence time τ th . Therefore, we effectively account for the thermal atomic motion by averaging over the atomic configurations. As a result, atomic dipole correlators in (10) with position-dependent phases vanish, regardless of the observation direction. In contrast, so-called incoherent terms [25], which do not exhibit such phases and correspond to α = β in equation (10), prevail. We note that, although this reduces the observable to a sum of single-atom correlators, these are obtained from the solution of the master equation (12) for N dipole-dipole interacting atoms. Therefore, through these interactions, the fluorescence intensity (10) still depends on the surrounding atoms. An analogous reasoning applies for the position-dependent phases in the coupling (17) and in the pulses' envelope (see equation (1)), as will be detailed in section 4.4.3.

Dynamics of independent atoms: single-atom case
In the following, we present our derivation of Ik(t) from (8). For noninteracting atoms (that is, L int = 0), it suffices to find a solution of a single-atom equation, as we show in this subsection. The dynamical behaviour of N atoms is then given by the tensor product of evolution laws for individual atoms, as we discuss in section 4.3.
By equation (12), an arbitrary single-atom operator Q α obeys the equation of motioṅ where the Liouvillians on the right-hand side of (23) are given by (13) and (14).

. Ultrashort time scale-atom-laser interaction
As argued in section 4.1, during the ultrashort intervals of non-vanishing laser intensity around the pulses' arrival times t j , dissipation processes generated by L α γ are insignificant. Hence, we need to solve the following equation of motion generated by the Liouvillian L α L : where Ω L (t) = 2dE L (t − t j )/ is the time-dependent (real-valued) Rabi frequency, and This is a generalized atomic dipole operator (whose index α we drop for simplicity), with S † = D † α · ε L , S = D α · ε * L atomic raising and lowering operators, according to (7), and the time-independent phase which includes the phase tag w j τ m proper, as well as the phases imprinted by the carrier frequency according to (1). The laser-atom detuning δ α is given by equation (15). Since the interaction with the Gaussian laser pulses occurs instantaneously in comparison to all other time scales (see section 4.1 and figure 2), these can be modelled by delta pulses carrying the same energy E pulse . Such pulses are characterised solely by their area which can be found using the relation (see, for instance, [48]) between the energy and the intensity of a pulse, with w 0 the laser beam's cross-section. Solving (28) for the field amplitude E 0 , we obtain ϑ = 2d(σcμ 0 E pulse √ π/w 0 ) 1/2 / , with μ 0 the free space permeability. Furthermore, under the delta pulse approximation, equation (24) is simplified by replacing Ω L (t) by a properly weighted Dirac delta-function, and using the well-known formula (see, e.g. [49]) where and Δ α = k L · v α is the Doppler shift of atom α. Now, the action of a delta kick at time t j on an atom can be expressed through a unitary superoperator R α j as where Q α (t − j ) (Q α (t + j )) is the single-atom operator before (after) the kick. We can deduce the transformation performed by R α j in an algebraic form. To that end, we notice that the operators S † and S on the right-hand side of equation (25) are generalizations of the Pauli pseudo-spin operators σ + and σ − to the case of degenerate excited state sublevels. Therefore, the algebraic properties of S † , S are similar to those of σ + , σ − . To illustrate this by an example, let us consider laser polarizations ε L =x and ε L =ŷ, for which the operator S reads, respectively (see also (11)), and the corresponding raising operators follow by Hermitian conjugation of (32). Using equations (25) and (32), we obtain which is a projector onto the subspace of the laser-driven levels |1 , |2 , and |4 . This implies and, hence, The last two terms of the obtained result are formally equivalent to the known operator identity for two-level systems (see, e.g. [50]), and describe an azimuthal rotation on the Bloch sphere by an angle ϑ. In addition, the (trivial) evolution of laser-undriven levels (in the basis (11), and for ε L =x (ε L =ŷ), these are the levels |y and |z (|x and |z )) is included in (35) via the projector 1 − M 2 j . With the aid of equations (35) and (25), we obtain [51] the following expansion for the superoperator R α j from (31): with superoperators R [p] j independent of ϕ j , and defined by their action on a single-atom operator Q α : where Equations (36) and (37) are crucial for the subsequent analysis of MQC signals in section 4.4.5, since they depend on the modulated phase ϕ j given by (30). We note that the simple structure of the expansion (36) is a direct consequence of the RWA [52] applied in the atom-laser Hamiltonian H AL as given in (4). While the decomposition (36) is general, in our scenario the atoms are initially unexcited. In this case, the atomic density operator after the first laser pulse does not feature factors e ±2iϕ 1 , which is important for the disorder average and signal demodulation (see sections 4.4.4 and 4.4.5). For thex-polarized pump pulse, using (11), (32), (36) and (37), we deduce where we additionally introduced the indexx, to make the polarisation explicit, and |x = S † |1 . The single-atom density operator after the pump pulse thus exhibits electronic coherences, modulated as e ±iϕ 1 , and populations of the ground and excited state sublevels dependent on the pulse area. In section 4.4.3 below, by combining delta-pulse excitations with the subsequent time evolution due to L γ + L int , we obtain an analytical formula for the two-atom state after the second pulse. In fact, for typical interpulse delays τ ∼ 10 ps, the dynamics generated by the Liouvillians L γ and L int can be neglected to a good approximation. With this in mind, we can derive explicit expressions for the states |ψ 1 ψ 1 | := R α 2,x R α 1,x |1 1| and |ψ 2 ψ 2 | := R α 1,x R α 2,ŷ |1 1| of independent (since we ignore L int ) atoms after a sequence of twox-x-andx-ŷ-polarized pulses, respectively. From equations (31), (35) and (39) we obtain = cos 2 ϑ 2 + e i(ϕ 1 −ϕ 2 ) sin 2 ϑ 2 |1 − i(e iϕ 1 + e iϕ 2 ) sin ϑ 2 cos Let us now highlight those matrix elements of |ψ 1 ψ 1 | and |ψ 2 ψ 2 |, respectively, which can contribute to MQC signals upon demodulation. Specifically, these are affected by the phase modulation only through the difference ϕ 1 − ϕ 2 . It follows from (40) that if both pulses have identical (x-x) linear polarization, the second pulse transforms electronic coherences exhibited by (39) into ground and excited state populations, σ 11 and σ xx , respectively, (see figure 3(a)), both modulated as cos(ϕ 1 − ϕ 2 ). In contrast, equation (41) shows that, if the pulses are perpendicularly polarized (x-ŷ), the Zeeman coherences (off-diagonal density matrix elements in the excited state manifold) σ xy and σ yx are created (see figure 3(b)), and modulated as e ±i(ϕ 1 −ϕ 2 ) . Note that this phase modulation at short delays τ is independent of the fact that L γ and L int are dropped. Before we proceed, let us notice that, so far, we described the action of the laser pulses on the atoms non-perturbatively in the pulse area ϑ. Strong pulses (e.g., ϑ ≈ π) induce a nonlinear response of the atoms [19,53] and can be used as a control and diagnostic tool in electronic spectroscopy [54,55]. Furthermore, we have shown that strong pulses substantially modify the MQC signals [42]. However, in our present contribution we seek to elucidate the impact of the dipolar interactions on these signals. Thus, we adhere to the experimental conditions of [18] and evaluate (37) for small ϑ, in the perturbative regime.

Beyond the ultrashort time scale
We now derive from (23) the time evolution generated by the Liouvillian L α γ as given by (14): with the formal solution (for simplicity, we take t 0 = 0) The explicit solution is derived in appendix A, with the result Q α (t) = P e Q α (0)P g e −γt/2 + P g Q α (0)P e e −γt/2 where P e = D † α · D α is the projector on the subspace of excited states, and P g = 1 − P e . Equation (44) has a transparent interpretation: the first line governs the exponential decay of electronic coherences with rate γ/2, while the remaining terms ensure trace conservation and the decay of excited state populations and Zeeman coherences with rate γ.

Dynamics of independent atoms: multi-atom case
An arbitrary operator Q(t) of N atoms can now be obtained by combining single-atom operators Q α in a tensor product, Q(t) = ⊗ N α=1 Q α (t). Consistently, the evolution of a composite operator Q(t) in the case of non-interacting atoms is governed by superoperators where L γ = N α=1 L γ α , the action of R α j is defined by equation (37), and that of L γ α by (44). The superoperator R L j describes the coherent interaction of the individual atoms with the ultrashort incoming pulses, whereas exp(L γ t) governs their incoherent dynamics when the laser field is switched off.
We stress that, by (44), this dynamics amounts to an exponential decay of all atomic variables in the multi-atom excited state manifold (that is, all excited state populations, as well as Zeeman and electronic coherences), while the collective ground state of all atoms is unaffected by exp(L γ t).
Finally, we introduce the resolvent superoperator G γ (z), z ∈ C, which is the Laplace image of the evolution superoperator exp(L γ t), and will be used in the next subsection, where we present the case of interacting atoms.

Dynamics of the dipole-dipole interacting atoms 4.4.1. Incorporation of L int
As discussed in section 4.1, the dipole-dipole interaction is relevant on short as well as on long time scales. Due to the weakness of the coupling parameter |g(ξ)| in the far field (see section 3.4), the coupling L int can be readily included perturbatively into our solutions (45) for independent atoms. Note that this treatment is distinct from the traditional one [13,17,[36][37][38], where the dipole-dipole interactions are accounted for non-perturbatively, via transformation to the exciton basis [19].
The starting point for the perturbative expansion in L int is the equation of motion for an N-atom operator Q, recall (12) and figure 2, valid on short and long time scales. Solving (47) perturbatively in L int , by the Laplace transform method [56], we obtain a series expansionQ whereQ(z) = ∞ t 0 dt exp(−zt )Q(t ), and Q(t 0 ) is an operator which encodes the initial condition, to be specified below. The time evolution of Q(t) then follows from (48) via the convolution theorem for Laplace transforms [56].

Combining piecewise solutions
Equation (47) holds during the (short) interpulse delay τ , and during the (long) fluorescence collection time t fl (see figure 2). Therefore, when solving (47) with the aid of the Laplace transform, we introduce distinct variables for each interval, z 1 and z 2 , respectively. During the short interpulse delay, the time evolution is given by the inverse Laplace transform with z 0 ∈ R chosen such that z 1 lies in the half-plane of convergence [56], andQ(z 1 ) is given by (48). The evolution during the long time interval is once again defined by (49), with the same z 0 , and upon the replacements The dynamics of the atomic operators during the ultrashort intervals can be absorbed in the Laplace transform solutionsQ(z 1 ) andQ(z 2 ), through the initial conditions specified by Q(t 0 ) in (48). The short interval τ starts after the pump pulse at time t 1 . Hence, by (31), we must set , it describes atomic operators before any interaction took place. Likewise, for the long time interval starting after the probe pulse, we set Q(t 0 ) ≡ Q(t + 2 ) = R L 2 Q(t − 2 ) in (48), with Q(t − 2 ) given by the inverse Laplace transform (49) with t − 2 = τ . Thus, through the initial conditions Q(t + 2 ), a solution during the long time interval includes the action of two laser pulses and the evolution during the short time interval. Combination of piecewise solutions, each of which is valid on distinct time intervals, finally provides the evolution of an arbitrary atomic operator throughout the entire excitation-emission cycle.

Single and double scattering contributions
The formal solution (48) applies for an arbitrary number N of atoms, though we focus on the case N = 2 in the following. At a first glance, this may appear a gross oversimplification of the experimentally studied physical system, including many atomic scatterers. Note, however, that we are dealing with a dilute thermal gas, whose fluorescence signal can be oftentimes inferred from the results for N = 1 [16]. Yet, as we show below, one cannot measure a double quantum coherence signal by monitoring fluorescence of one, nor of two non-interacting atoms. Nor need we consider more than two atoms exchanging photons. Indeed, the atoms are in the far field of, and detuned from each other. The latter condition stems from the motion-induced Doppler shifts being much larger than the natural linewidth. The photon scattering mean free path [57] in such medium can be estimated as 10 withσ sc the mean scattering cross-section (see appendix B). Therefore, a thermal atomic vapour of a reasonable size (say, a spherical cloud with a diameter D of a few centimetres) is characterized by a very low optical thickness [40] b = D/ sc 1, such that long scattering sequences of photons, mediating interactions between many particles, are highly improbable. Consistently, single and double scattering contributions do capture the essential physics in optically thin ensembles of scatterers [58], and we restrict our perturbative expansion (48) to second order in L int . This is the lowest-order contribution which survives the disorder average, and it accounts for the double scattering 11 of a single photon [26]. Thereby, we ignore recurrent scattering [59], where a photon bounces back and forth between the atoms, and which corresponds to a higher-order expansion in L int . Neglecting such processes is fully justified in optically thin media [40].
Indeed, as we will see below, the above framework enables deducing 1QC and 2QC signals in good qualitative agreement with the experiment [18], and establishes the relevance of few-body systems in the characterization of MQC signals emitted by dilute thermal atomic ensembles.
The single and double scattering contributions can be derived from the two-dimensional Laplace transformsQ [n] (z 1 , z 2 , ϕ 1 , ϕ 2 with n = 0, 2, which are obtained from equation (48) and the contribution of piecewise solutions as described in section 4.4.2. Note that the phase tags (from equation (30)) in the argument of (51) are absorbed in R L 1 and R L 2 , via (36) and (45). To infer the detected fluorescence intensity from the operator equation (51), we need specific atomic correlation functions which correspond to α = β in equation (10) (see section 4.1). These are obtained by representing the operator equations in matrix form, via the complete basis of the operator space for two atoms (see appendix C), and subsequently evaluating the quantum mechanical expectation value, followed by the configuration average: where G γ (z), R L j , and V are 256 × 256 matrices representing operators G γ (z), and R L j , L int , respectively, with elements given by while [Q(0)] k = Tr Q † k |1 1 1| 1 ⊗ |1 2 1| 2 is a vector with 256 elements. Let us finally represent the Liouvillians L (1) int and L (2) int , see (21) and (22), with the matrices

Configurational average
The configuration or disorder average in (52) is necessary to capture the impact of changes of the random atomic configuration during the fluorescence detection (section 4.1). As a result of averaging, all terms carrying position-dependent phases vanish. In addition to interference contributions with α = β in equation (10) which are already excluded in (52), the disorder average obliterates all terms with r α0 -dependent phases-which are due to the action of the two incoming laser pulses (see equations (30) and (36)). What remains are terms carrying position-independent phases that originate in single scattering processes, independent of V, and double scattering processes, quadratic in V. These single and double scattering contributions have weights |g(ξ)| 0 = 1 and |g(ξ)| 2 = 1/ξ 2 , respectively. Furthermore, only particular products of the matrix elements of ← → T (ξ,n) (which enter (52) at second order in L int , through (16)) contribute to the average signal upon integration over the isotropic distribution of the vectorn [25,26]: ← → T kl ← → T * mn conf = 0 if and only if the indices k, l, m, n are pairwise identical, cf appendix D. This amounts to the suppression, in the average fluorescence intensity, of all contributions stemming from the interference of atomic excited state sublevels with different magnetic quantum numbers.
Finally, note that through the matrix R L j the atomic correlation functions carry phases Δ α t j proportional to the Doppler shifts, by virtue of (30), (36). Unlike the uniform spatial distribution of the atoms, the Doppler shifts Δ α are non-uniformly Maxwell-Boltzmann distributed. Therefore, the result of averaging over Δ α is different from that of averaging over random atomic positions; it would not eliminate the terms including velocity-dependent phases, but instead result in inhomogeneous line broadening [31] of MQC spectra [42]. In our present contribution, we ignore this effect and set the Doppler shifts in equation (30) to zero, while focussing on the general structure and physical origin of the resonances of 1QC and 2QC spectra.
The summation range in (55) is a consequence of equation (39): since, for n = 0, the time evolution (52) is local in the Hilbert space of each atom, individual atomic operators can depend only on the phase tags e ±iϕ 21 , and not on e ±2iϕ 21 . In contrast, fluorescence intensity emerging due to double scattering can carry phases ±2ϕ 21 inĨ [2] k (z 1 , z 2 , ϕ 21 ), giving rise to 2QC signals. Taking the inverse Laplace transforms (see (49)) of the functions defined by (55) and (56) with respect to z 1 and z 2 , we obtain the associated time-domain expressions I [0] l,k (τ , t fl , ϕ 21 ) and I [2] l,k (τ , t fl , ϕ 21 ), I [2] k (τ , t fl , ϕ 21 ) = 2 l=−2 I [2] l,k (τ , t fl )e iϕ 21 l , where I [0] l,k (τ , t fl ) and I [2] l,k (τ , t fl ) are exponentially decaying functions of τ and t fl (recall our discussion in section 4.3), related toĨ [n] l,k (z 1 , z 2 ) by the two-dimensional Laplace transform: The total transient intensity is the sum of (57) and (58), which can be used to deduce 1QC and 2QC spectra, as outlined in section 2, and discussed below.

Signal demodulation
The goal of the demodulation procedure is to extract only those specific components from the fluorescence signal, which oscillate at defined integer multiples κ of the modulation frequency w 21 . The spectra of these signals, generally referred to as MQCs, are extracted from (2), after replacement of φ 1 and φ 2 by the phase difference ϕ 21 , via the following two-fold, half-sided Fourier transforms, with κ = 1, 2 the MQC orders here considered. (60), with I [0] k (τ , t fl , ϕ 21 ) and I [2] k (τ , t fl , ϕ 21 ) from (57) and (58), into (61), and integration over τ m yields the result: 1,k ,Ĩ [2] 1,k , andĨ [2] 2,k , respectively (see (63) and (64)). Horizontal red lines indicate the incoming laser pulses, vertical yellow lines the emitted fluorescence, oblique blue lines the scattered fields via dipole-dipole interactions. Solid and dashed lines indicate single field amplitudes and their complex conjugates (see also figure 1). The difference between the double scattering processes (b) and (c) is that, in the former case, fluorescence stems from the excited state sublevels of the second atom that are driven by the scattered, but not by the incident field, while, in the latter case, by both the scattered and the incident field. As discussed in section 5, processes (a), (b) and (c) are associated with 1QC and 2QC signals, respectively.
To summarise, through (63) and (64) we connect 1QC and 2QC spectra as measured by methods of phase-modulated spectroscopy to specific components of the fluorescence intensity emitted by two atoms. The intensity is expressed via (10) (with N = 2 and, as a consequence of the configuration average over the atomic random positions, α = β) through the excited state populations of single atoms that may be independent or interacting with each other by exchanging a single photon. Furthermore, while equation (63) indicates that the 1QC signal is fed by both, single and double scattering processes, equation (64) shows that 2QC signal can stem solely from double scattering, see figure 4. In other words, whereas 1QC signals may be observed on independent particles, 2QC signals are always an unambiguous proof of interatomic interactions. Thus the origin of 2QC signals in phase-modulated spectroscopy on disordered samples [15,18] is substantially different from the one outlined in [16], where 2QC signals stem from two-atom observables of non-interacting atoms.
The dipolar interactions are represented by excitation exchange and collective decay processes encapsulated in the Liouvillians L (1) αβ and L (2) αβ defined by (21) and (22). In section 5, we present the spectra (63) and (64) and elucidate the significance of these Liouvillians for the emergence of 1QC and 2QC signals, which explains a number of experimentally observed features of MQC spectra [18].

Spectra of single and double quantum coherence signals
In figure 5 we plot the single and double coherence spectra Sk(ω; 1) and Sk(ω; 2), given by (63) and (64), for k =x,ŷ, and for parallel (x-x) and perpendicular (x-ŷ) pump-probe polarization channels (in brief, and ⊥ channels). As discussed in section 4.4.4, these signals emerge from the configuration average. For illustration, the background of each panel in figure 5 shows MQC spectra stemming from pairs of atoms at fixed locations. These were calculated by omitting the average · · · conf in (52) and, instead, by considering randomly drawn configurations (ξ,n), which determine the tensor ← → T and the phase factor exp(ik · r 12 )-as detailed in the caption of figure 5. The deviations between background and averaged  (10) and (61)], for detection directionsk =x (left) andk =ŷ (right), and for quantum coherences of order κ = 1 (top) and κ = 2 (bottom). Symbols and ⊥ indicate parallel (x-x) and perpendicular (x-ŷ) pump-probe polarizations, respectively. Parameter values match the experimental ones [18]: average interatomic distancer ≈ 10 μm (which corresponds to a particle density ρ = 10 8 cm −3 ), laser pulse areas ϑ = 0.14π and durations σ = 21 fs (tuned to resonance with the D2-line of 87 Rb atoms [30]), transition wavelength λ 0 = 790 nm, and spontaneous decay rate γ ≈ 2π× 6.067 MHz. This choice of parameters results in a mean scaled interatomic distancē ξ = k 0r ≈ 80, with k 0 the atomic transition's wave number. Thin lines in the background show MQC spectra (with the same colour coding for the real and imaginary parts as the average spectra) for 100 random configurations (ξ,n), with ξ distributed uniformly in the interval 67.2 < ξ < 92.8, andn distributed isotropically. The insets show the signals obtained by substituting ← → T → i ← → Ω , which amounts to ignoring the collective decay processes. spectra are striking for both parallel and perpendicular pulses: the single realizations have little resemblance with the true (configuration averaged) line shapes, from which they oftentimes differ by orders of magnitude. The absence of the configuration average constitutes one of several reasons for a substantial qualitative and quantitative mismatch between the model including two fixed molecules [17] and the observations reported in [18]. Further reasons will be identified below. Henceforth, we will discuss our configuration-averaged spectra. We stress that, while our general analytical solutions, derived from (51) and (52), assume that dipolar interactions act during both, the interpulse delay and fluorescence harvesting, we have checked numerically that the 1QC and 2QC spectra can be fully understood if we include the interatomic interactions only during the fluorescence detection stage. This is reasonable, given the typically short duration of ∼ 10 ps of the interpulse delay and the weakness of the dipolar coupling. Therefore, throughout this section we interpret the emergence of 2QC spectra (as much as contributions to 1QC spectra) as the result of dipole-dipole interactions during the long time scale, i.e. after the atoms have interacted with the two incoming laser pulses.
In the channel, the real and imaginary parts of the 1QC spectra exhibit, respectively, absorptive and dispersive resonances at ω = ω 0 . In the here realized regime of relatively weak (perturbative) driving fields, these line shapes are consistently proportional to the real and imaginary parts of the linear susceptibility [29]. For both choices of pump-probe polarizations and observation directions, the 2QC spectra feature qualitatively identical line shapes centred at ω = 2ω 0 , albeit with opposite sign to the 1QC spectra. The sign flip of the 2QC spectra can be understood by recalling that they originate from double scattering. Thereby, each atom effectively interacts with four fields: two laser pulses and two far-fields scattered by the other atom, see figure 4 (c). This nonlinear four-wave mixing process is described by the third-order susceptibility [19,60], whose sign is opposite to that of the linear susceptibility [60]. We would like to stress that the sign flip of the 2QC spectra was experimentally observed [28], but its interpretation has been lacking.
In the channel, the peak magnitude of the 1QC signal fork =x is four orders of magnitude smaller than fork =ŷ, see figure 5. This significant difference arises because the signal fork =ŷ is dominated by single scattering from individual atoms, as illustrated by figure 4(a) (i.e., by theĨ [0] 1,k component in (63)). This contribution to the single coherence signal is independent of the interatomic distance. Fork =x, on the contrary, the atoms in their |x state cannot directly emit in thex-direction. Hence, the signal is due to double scattering contributionĨ [2] 1,k (see figure 4(b)), wherein the excitation is transferred from the state |x of one atom to the state |y or |z of the other-which was in its ground state prior to the interaction. Such an interaction process, involving the exchange of an excitation between the atoms, is mediated solely by the  1 and figure 3) to the excited state populations of the atoms which give rise to fluorescence signals. Only those variables that contribute to 2QC signals are indicated. (a) In the ⊥ channel and for the observation directionk =x, Zeeman coherences σ α xy σ β xy (red arcs) of independent atoms are transformed, by the collective decay process mediated by the matrix V (2) , and accompanied by the emission of a single photon, into electronic coherences (orange arcs), e.g., σ α 1y σ β x1 . By the interaction represented by the matrix V (1) , the atomic state |y α is excited, while atom β is de-excited, which is described by the correlator σ α yy σ β 11 (black dot). The fluorescence emitted by atom α contributes to the 2QC signal. (b) In the channel and for the observation directionk =ŷ, the pathway indicates the following sequence of processes: In the channel and for the observation directionsk =ŷ orx, additional pathways are mediated solely by the matrices V (1) Furthermore, the magnitudes of the spectra Sx(ω; 1) and Sx(ω; 2) are reduced by a factor of two, in accordance with the facts that these signals are mediated solely by V (1) , hence, ∝ ← → Γ xr ← → Γ xr conf + ← → Ω xr ← → Ω xr conf (see the above paragraph), and that ← → Γ xr ← → Γ xr conf = ← → Ω xr ← → Ω xr conf (r = y, z). As for the 2QC spectrum Sŷ(ω; 2), its magnitude is diminished and its sign is flipped when the collective decay and excitation exchange processes proportional to the tensor ← → Γ are switched off. Thus, using the static form of the dipole-dipole interaction [13,16,17,36,37] (that is, completely ignoring ← → Γ ) results in overlooking quantitative and/or qualitative features of 1QC and 2QC signals which were measured experimentally [18,28], in both ⊥ and channels, and forx andŷ observation directions.

Conclusions
We developed an original, microscopic, open quantum systems theory of MQC signals in fluorescence-detection based measurements on dilute thermal gases. In general, such signals are regarded as sensitive probes of dipole-dipole interactions in many-particle systems [13], but recent experiments [15] initiated some debate on their true nature [16][17][18]62]. One of the purposes of this work was to resolve the persisting controversy. Our theory is characterised by several distinctive features. First, we incorporated the radiative part of the dipole-dipole interactions, including the concomitant collective decay processes which are ignored in other treatments [13,15,17,[36][37][38]. Due to the large interatomic distance in dilute atomic vapours, this interaction is here treated perturbatively. In this perspective, 2QC signals are sensing the exchange of a single photon within pairs of atoms, and the two-atom correlations resulting from this fundamental interaction process. Furthermore, we considered a realistic atomic level structure, equipping the atoms with angular momentum, and treated the coupling of such atoms to polarized laser pulses of arbitrary strength non-perturbatively. In addition, we effectively accounted for the thermal atomic motion with a configuration average. While it is the combination of all these features that allowed us to obtain excellent qualitative, and reasonable quantitative agreement with hitherto partially unexplained experimental observations [18,28], here we would like to highlight the role of collective decay processes in the emergence of MQC signals. Such processes are completely ignored in treatments which employ the electrostatic form of the dipole-dipole interactions [13,15,17,[36][37][38].
In the present contribution, we focussed on 1QC and 2QC signals, which we obtained by considering a two-atom model. On the one hand, in contrast to previous claims [16], we have shown that the detection of 2QC signals unambiguously indicates the presence of dipolar interactions-mediated by the exchange of real photons. On the other hand, our results show that substantial deviations between theory [17] and experiment [15,18] do not necessitate the consideration of many (N > 2) interacting particles to predict the behaviour of 2QC signals; two randomly arranged interacting particles encapsulate the relevant physics.
Our method defines a general and versatile framework for the systematic assessment of higher-order contributions-which were observed experimentally [15,18,28]-by combining diagrammatic expansions in the far-field dipole-dipole coupling with single-atom responses to laser pulses of arbitrary strength [26,27]. In the future, we will explore the nature of the correlations-quantum or classical-arising from the exchange of a photon between the atoms. Finally, we will generalize our method to resolve more subtle features of MQC spectra [28] due to fine and hyperfine structures of the involved dipole transitions.
where A is a generalized superoperator. Master equation (12) is an example of such equations. Equation (C.1) has the formal solution

Q(t) = exp(A t)Q(0). (C.2)
For individual atoms, such solutions can be obtained in a relatively compact analytical form (see section 4.2), which simplifies their physical interpretation. The superoperator exp(A t) in equation (C.2) generates coupled dynamics of arbitrary atomic operators lumped together in a vector Q. We choose the latter's elements Q n to form a complete orthonormal basis, which for a single atom reads: Q = (1/2, μ 1 /2, μ 2 /2, μ 3 /2, σ 14 , σ 41 , σ 13 , σ 31 , σ 12 , σ 21 , σ 34 , σ 43  The elements of Q satisfy the properties of orthonormality, Tr{Q † m Q n } = δ nm and completeness, O = 16 n=1 c n Q n , with O an arbitrary atomic operator and expansion coefficients c n . For N atoms, the complete basis is given by the tensor product Q 1 ⊗ Q 2 ⊗ · · · ⊗ Q N . In this case, vector Q has dimension 16 N , and it evolves according to where A(t) = [a nm (t)] is a (16 × 16) N matrix, whose elements a nm (t) = Tr{Q † m exp(A t)Q n }, with Q n the elements of the N-atom basis set. Performing a quantum mechanical average on both sides of equation (C.8) with respect to the initial atomic density operator, we obtain Q(t) = A(t) Q(0) , ( C . 9 ) wherefrom we can express the expectation value of any N-atom observable.