Quantum effects in the cooperative scattering of light by atomic clouds
Lorenzo Pucci, Analabha Roy, Tiago Santiago do Espirito Santo, Robin Kaiser, Michael Kastner, Romain Bachelard
QQuantum effects in the cooperative scattering of light by atomic clouds
Lorenzo Pucci, Analabha Roy, Tiago Santiago do Espirito Santo, Robin Kaiser, Michael Kastner,
1, 4 and Romain Bachelard ∗ National Institute for Theoretical Physics (NITheP), Stellenbosch 7600, South Africa Instituto de F´ısica de S˜ao Carlos, Universidade de S˜ao Paulo, 13560-970 S˜ao Carlos, SP, Brazil Universit´e de Nice Sophia Antipolis, CNRS, Institut de Physique de Nice, UMR 7335, F-06560 Valbonne, France Institute of Theoretical Physics, Department of Physics,University of Stellenbosch, Stellenbosch 7600, South Africa (Dated: November 12, 2018)Scattering of classical light by atomic clouds induces photon-mediated effective long-range interac-tions between the atoms and leads to cooperative effects even at low atomic densities. We introducea novel simulation technique that allows us to investigate the quantum regime of the dynamics oflarge clouds of atoms. We show that the fluorescence spectrum of the cloud can be used to probegenuine quantum cooperative effects. Signatures of these effects are the occurrence, and the scalingbehavior, of additional sidebands at twice the frequency of the classical Mollow sidebands, as wellas an asymmetry of the Mollow triplet.
I. INTRODUCTION
Our understanding of the quantum dynamics of many-body systems has benefited from a number of recentachievements. On the experimental side, cold atom sys-tems and ion traps have reached an unprecedented levelof control and allow for the emulation of a large vari-ety of many-body Hamiltonians of interest, including thepossibility of tuning coupling parameters [1–4]. On thenumerical side, progress in the understanding of matrixproduct state has boosted density matrix renormaliza-tion group and other methods, applicable primarily toone-dimensional quantum systems. Reliable simulationsof higher-dimensional systems are in general more diffi-cult, and in many cases impossible.A three-dimensional setting where many exciting dis-coveries have been made, but also many open questionsremain, is light scattering by large clouds of atoms. Theclassical, linear optics regime of such systems has beenextensively studied, and many-body effects, such as su-perradiance [5, 6], modification of the radiation pressureforce [7, 8] or cooperative frequency shifts [9–13], were re-ported. The effective coupling between the atoms in thecloud, mediated by the photon field, turns out to be long-ranged, and as a result cooperative effects occur even indilute clouds, including superradiance [14, 15], Dicke sub-radiance [16], and spectral broadening [17]. The strengthof the cooperative effects depends in these cases on thecloud optical thickness, not the spatial density.Beyond the linear optics regime, saturation effects givethe atoms a nonlinear response to the electric field of thelight. When many atoms are simultaneously in the ex-cited state, the classical coupling of these nonlinear os-cillators leads, for example, to the modification of theline shape of the atomic resonance [18]. Entering thefield of quantum optics, two-time correlations of the ra- ∗ [email protected] diated field provide information on the optical coherenceof the first kind g (1) and on the fluorescence spectrum.For a single atom, saturation leads to the emergenceof the Mollow triplet, a trio of spectral lines of inelas-tically scattered light, one around the laser frequencyand two symmetric ones shifted by the generalized Rabifrequency [19, 20]; see Fig. 1 (left) for an illustration.In a four-wave mixing configuration, these symmetricbands can be amplified as the wave propagates withinthe cloud [21, 22]. Furthermore, this redistribution offrequencies is iterated as two scatterers interact throughtheir radiation [23, 24], which, for dense atomic clouds,results in the presence of additional sidebands at twicethe Rabi frequency for pairs of atoms closer than a wave-length [25]. In the many-body case, quantum correlationsbecome essential when trying to understand the quan-tum features of the collective response of the system.This quantum regime of cooperative light scattering infree space is essentially unexplored. The main obstaclefor a quantum mechanical treatment is, as usual, the ex-ponential growth of the Hilbert space with the numberof atoms. Additionally, and different from, e.g., opti-cal cavities, no suitable collective operators are knownfor describing the dynamics effectively. A recent firststudy of quantum effects in the light scattering by di-lute atomic clouds made use of perturbative techniques,valid in the regime of strong driving, and was able to pre-dict asymmetries in the scattering spectrum, but no ad-ditional spectral lines beyond those of a single atom [26].In this paper we report the discovery of quantum coop-erative effects in light scattering by large dilute clouds ofatoms. We show theoretically that quantum correlationsbuilding up in the cloud induce cooperative sidebands inthe many-body fluorescence spectrum at twice the Rabifrequency, and also lead to a cooperative spectral asym-metry in the Mollow triplet. Both effects scale with theoptical thickness, which is a hallmark of their coopera-tive character. Investigating the angular dependence ofthe scattering spectrum, we find that quantum coopera-tivity is more easily detected at large scattering angles, a r X i v : . [ phy s i c s . op ti c s ] M a r and not in the forward direction. These results provideguidance on where to look for quantum cooperative ef-fects in light scattering experiments on atomic clouds,and are expected to be relevant also for neutral-atom op-tical clocks, Rydberg atoms, and other settings whereeffective long-range interactions play a role.This study of nonperturbative quantum effects in fairlylarge three-dimensional atomic clouds became possibleby a novel simulation technique, combining a discretephase-space representation of spins [27, 28] with higher-order semi-classical evolution equations [29] extendedto driven-dissipative Lindblad dynamics. The methodis highly accurate for higher-dimensional systems withlong-range interactions, and therefore a perfect fit for theproblem at hand. II. MODELING THE ATOMIC CLOUD
For our purposes the cloud of atoms is modelled asan assembly of N two-level systems at fixed, but usuallyrandom, positions r i in three-dimensional space. Using ascalar light model, which is valid for dilute clouds, eachtwo-level system can be described by Pauli spin opera-tors σ ± = ( σ x ± iσ y ) / σ z . This model provides agood description of dilute clouds of atoms cooled belowthe Doppler limit and trapped in a magneto-optical ordipolar trap. Transitions between the two levels of eachatom are driven by a classical planar-wave laser light-field of wave vector k , Rabi frequency Ω , and detuning∆ = ω − ω a from the optical transition. The dynamicsis then described by a set of equations that couple thespin degrees of freedom to the photon field [30]. Perform-ing the rotating-wave, Born and Markov approximationsand eliminating the photon degrees of freedom, equationsof motion can be derived for the atomic internal degreesof freedom in the Heisenberg picture [26, 30], dσ − j dt = (cid:18) i ∆ − Γ2 (cid:19) σ − j + i Ω e i k · r j σ zj + Γ2 N (cid:88) m (cid:54) = j σ zj σ − m ( γ jm − i ∆ jm ) , (1a) dσ zj dt = i Ω (cid:0) e − i k · r j σ − j − h.c. (cid:1) − Γ(1 + σ zj ) − Γ N (cid:88) m (cid:54) = j ( σ + m σ − j ( γ jm + i ∆ jm ) + h.c.) , (1b)where the coefficients γ jm = sin( k | r j − r m | ) k | r j − r m | , ∆ jm = cos( k | r j − r m | ) k | r j − r m | (2)describe the spatial dependence of the light-mediatedlong-range coupling between the atoms, and Γ is the tran-sition linewidth of the two-level atoms. Note that the di-lute regime allows to neglect near-field terms, which scaleas 1 /r , and thus focus on long-range effects. The quantity we want to study, and which is acces-sible in light-scattering experiments, is the fluorescencespectrum S ( ω ) = lim T →∞ lim t →∞ (cid:90) T − T dτ g (1) ( t, τ ) e − iωτ , (3)defined as the Fourier transform of the first-order opticalcoherence g (1) ( t, τ ) = (cid:104) E ( ˆ n , t ) E † ( ˆ n , t + τ ) (cid:105)(cid:104) E ( ˆ n , t ) E † ( ˆ n , t ) (cid:105) (4)in the steady state [realized by the limit t → ∞ in (3)].The optical coherence is defined in terms of the electricfield operator E which, in the far-field limit and emittedin the direction of the normalized vector ˆ n , is given by E † ( ˆ n , t ) ∝ N (cid:88) j =1 σ − j ( t ) e − ik ˆ n · r j . (5)From Eqs. (3)–(5) one can read off that the crucial non-trivial ingredient for calculating the fluorescence spec-trum are the two-time spin–spin correlation functions (cid:104) σ + j ( t ) σ − j ( t + τ ) (cid:105) , time-evolved under the evolution equa-tions (1a) and (1b). III. NUMERICAL METHOD
An analytic calculation of the fluorescence spectrum S was reported in Refs. [23, 24] for the case of two atomspumped at resonance ω a . In that work, sidebands at fre-quencies ω a ±
2Ω were shown to rise in the spectrum,yet the effect is substantial only for atoms closer thana wavelength. A similar effect was predicted for pairsof quantum dots [31]. Treating more than two atomsis much harder. An exact numerical evaluation of theevolution equations (1a) and (1b) is possible for about adozen atoms or so, but these numbers are too small toinvestigate scaling behavior or extrapolate to experimen-tally relevant system sizes.Here we deal with these difficulties by developing asimulation method for driven-dissipative quantum me-chanical evolution equations. Our method builds upon asimulation technique that makes use of a discrete phase-space representation of the Pauli operators [27, 28], andtime-evolves phase space points as well as correlation co-efficients according to semiclassical evolution equations[29]. We developed these techniques further by extend-ing them to nonunitarily-evolving driven-dissipative sys-tems, and made the method applicable to the compu-tation of two-time correlation functions by making useof the quantum regression theorem. The main featuresof the method are: (i) The use of semiclassical time-evolution equations makes this method particularly suit-able for systems with long-range interactions and systemsin higher spatial dimension. (ii) The explicit incorpora-tion of correlation coefficients in the quasiclassical equa-tions of motion allows for the accurate computation ofcorrelation functions; for N = 2 the exact results are re-covered. (iii) The numerical cost scales polynomially, notexponentially, with the system size N ; hundred and moreatoms can be treated.The novel numerical method, described in the follow-ing, is developed for calculating the fluorescence spec-trum (3) of a cloud of two-level atoms as described by thetime-evolution equations (1a) and (1b), but is also appli-cable, more generally, for the calculation of two-time cor-relation functions in long-range spin models. The atomiccloud we consider is three-dimensional, and the atom-atom interactions are long-ranged, which makes the sys-tem particularly suitable for simulation methods basedon the quasiclassical time-evolution of discrete phase-space points [28, 29]. However, these methods are for-mulated for unitarily evolving quantum spin systems andgive access to equal-time correlation functions, but not totwo-time correlations. In Sec. III A we rewrite the quan-tities needed for calculating the fluorescence spectrumin a form that is more amenable to a semi-classical time-evolution. In Sec. III B we derive semi-classical equationsof motion for driven-dissipative, non-unitarily evolvingquantum spin systems. In Sec. III C we introduce a dis-crete phase space representation and apply it to the cal-culation of the two-time correlation functions that arerequired for obtaining the fluorescence spectrum. The re-sulting novel simulation scheme is benchmarked againstexact results in Appendix B. A. Quantum regression
We extract the spectral properties of the radiated lightin direction ˆ n from the first order optical coherence g (1) ( t, τ ) (4) by evaluating the expression for the spec-trum (3), which can be rewritten as S ( ω ) = lim T →∞ lim t → + ∞ (cid:90) T − T dτ e − iωτ × (cid:104) g (1) ( t, τ )Θ( τ ) + g (1) ( − t, τ )Θ( − τ ) (cid:105) , (6)where Θ denotes the Heaviside step function. g ( − t, τ )is a short-hand notation for flipping the sign of t of theunitary terms in (1b) while keeping the dissipative onesunchanged, which corresponds to time-reversing the dy-namics of the total system (two-level atoms plus photons)before eliminating the photonic field. Taking the limit t → + ∞ inside the integral, we see that lim t →∞ g (1) ( t, τ )and lim t →∞ g (1) ( − t, τ ) are the two limits required forcalculating S . Here, for being definite, we discuss thefirst case only, as the second can be treated analogously.According to Eqs. (4) and (5), the optical coherence g (1) ( t, τ ) can be expressed in terms of two-time correla-tions of Pauli spin operators, and hence our main ob- ject of interest will be lim t →∞ (cid:10) σ ai ( t ) σ bj ( t + τ ) (cid:11) for all a, b ∈ { x, y, z } .Making use of the quantum regression theorem [32, 33],we can write the two-time spin–spin correlation as (cid:10) σ ai ( t ) σ bj ( t + τ ) (cid:11) = Tr (cid:8) V ( t + τ, t ) [( V ( t, ρ ) σ ai ] σ bj (cid:9) , (7)where ρ is the initial state and V ( t, t ) denotes the time-evolution operator corresponding to the equations of mo-tion (1a) and (1b). The coefficients of these differentialequations do not explicitly depend on time, which im-plies V ( t + t , t ) = V ( t,
0) for all t and t . Under thisevolution and in the limit t → ∞ , we expect ρ to evolvetowards a steady state, which we denote by ρ ss = lim t →∞ V ( t, ρ. (8)This allows us to write the long-time limit of the two-timecorrelation function aslim t →∞ (cid:10) σ ai ( t ) σ bj ( t + τ ) (cid:11) = Tr (cid:2) σ bj V ( τ,
0) ( ρ ss σ ai ) (cid:3) . (9)The nontrivial object in this expression is the time-evolved operator V ( τ,
0) ( ρ ss σ ai ), which we are going tocalculate in Sec. III C. B. Time-evolution equations
As a starting point for calculating V ( τ,
0) ( ρ ss σ ai ) weneed the steady state density operator ρ ss , which we cal-culate approximately by a method described in the fol-lowing. This method is applicable to arbitrary N -spintrace-1 operators A ...N under the dynamics generatedby a Lindblad operator. Starting from the initial densityoperator, A ...N = ρ , will allow us to obtain ρ ss after suf-ficiently long evolution times. Other choices of A ...N willbe used to calculate the time-evolution of other trace-1operators that appear when calculating V ( τ,
0) ( ρ ss σ ai ) inSec. III C.The propagator V that induces the evolution equations(1a) and (1b) can be written as a Lindblad differentialequation i∂ t A ...N = L A ...N , (10)where the Lindblad operator L = (cid:88) i L i + (cid:88) ij L ij (11)consists of on-site terms L i A ...N = − ∆ σ zi , A ...N ]+ Ω (cid:2) e − i k · r i σ − i + e i k · r i σ + i , A ...N (cid:3) (12)and pair interactions L ij A ...N = ∆ ij (cid:2) σ + i σ − j , A ...N (cid:3) + iγ ij (cid:0) σ − j A ...N σ + i − σ + i σ − j A ...N − A ...N σ + i σ − j (cid:1) , (13)with ∆ ij = − Γ2 cos ( k | r i − r j | ) k | r i − r j | for i (cid:54) = j, i = j, (14) γ ij = Γ sin ( k | r i − r j | ) k | r i − r j | for i (cid:54) = j, i = j, (15)and Γ = d k / (2 π (cid:126) (cid:15) ).Taking partial traces on both sides of (11) one obtains,in the spirit of the Bogoliubov-Born-Green-Kirkwood-Yvon (BBGKY) hierarchy [34], a set of coupled evolutionequations, i∂ t A i = L i A i + (cid:88) m (cid:54) = i Tr m L im A im , (16a) i∂ t A ij = ( L i + L j + L ij ) A ij + (cid:88) m (cid:54) = i,j Tr m ( L im + L jm ) A ijm , (16b)where the reduced operators A i = Tr { k (cid:54) = i } A ...N , A ij = Tr { k (cid:54) = i,j } A ...N , (17)are defined as partial traces over all sites except the in-dexed ones. It is expected for the dilute regime that thetwo-atom coupling dominates over three- or more-atomcoupling, which makes a cluster expansion suitable forthe truncation of the hierarchy of equations. By meansof a cluster expansion we separate the reduced A opera-tors into product and connected parts, A ij = A i A j + C ij , (18a) A ijm = A i A j A m + A i C jm + A j C im + A m C ij + C ijm , (18b)which implicitly defines the connected operators C . Sub-stituting these definitions into (16a) and (16b), werewrite the first two equations of the BBGKY hierarchyas i∂ t A i = L i A i + (cid:88) m (cid:54) = i Tr (cid:2) L Sim ( C im + A i A m ) (cid:3) , (19a) i∂ t C ij = ( L i + L j ) C ij + L Sij ( C ij + A i A j ) (19b) − A i Tr i (cid:2) L Sij ( C ij + A i A j ) (cid:3) − A j Tr j (cid:2) L Sij ( C ij + A i A j ) (cid:3) + (cid:88) m (cid:54) = i,j Tr m (cid:2) L Sim ( A i C jm + A m C ij + C ijm ) (cid:3) + (cid:88) m (cid:54) = i,j Tr m (cid:2) L Sjm ( A j C im + A m C ij + C ijm ) (cid:3) , where L Sij = L ij + L ji . Equation (19b) contains three-spin connected contributions C ijm , the time-evolution ofwhich depends on four-spin terms, and so on. To turnthis into a numerically tractable problem, we truncatethe BBGKY hierarchy at second order by neglecting theterms C ijm in (19b). As stated in [29], if we also ne-glected the C ij terms, we would recover the classical time-evolution equation presented in [28], in that sense, moreterms of the truncated hierarchy means more quantumcorrections for the spins dynamics. The two-spin con-nected contribution C ij are related to the two-spin quan-tum correlations, which implies that, unlike in the classi-cal case, the two-spin connected correlations do not van-ish, (cid:104) σ ± ,zj σ ± ,zm (cid:105) − (cid:104) σ ± ,zj (cid:105)(cid:104) σ ± ,zm (cid:105) (cid:54) = 0. This is the main ap-proximation made in the numerical scheme, and it givesgood results whenever genuine three- and more-spin con-nected contributions are negligible.To bring the resulting truncated operator equationsinto a numerically tractable form, we expand all oper-ators in terms of Pauli spin operators and we obtain aset of coupled ordinary differential equations (see Ap-pendix A), which can be integrated by standard numer-ical methods. C. Calculation of V ( τ,
0) ( ρ ss σ ai ) We calculate the steady-state density operator ρ ss bysetting to zero the left-hand sides of (A2) and (A3), andnumerically solving the resulting algebraic equations bya standard Newton-Krylov solver. The stationary valuesof the a- and c-coefficients in those equations encode therequired information on ρ ss .Starting from the thus obtained steady-state densityoperator ρ ss , we expand ρ ss σ ai in terms of so-called phasepoint operators, rewrite the result in terms of trace-1 op-erators, and then use again the time-evolution equationsof Secs. B and C in order to obtain V ( τ,
0) ( ρ ss σ ai ).The discrete phase-space representation of a singlespin-1 / { (0 , , (0 , , (1 , , (1 , } (20)consisting of four points, each of which has an associ-ated three-vector, r (0 , = (1 , , r (0 , = ( − , − , r (1 , = (1 , − , − r (1 , = ( − , , − α ∈ Γ one assigns a so-called phasepoint operator A α = ( + r α · σ ) , (21)where σ = ( σ x , σ y , σ z ) is the vector of Pauli operators.The phase point operators form a basis, and any operatoron C can be expressed as a linear combination of the fouroperators A α . Similarly one could expand an operator onthe tensor product Hilbert space ( C ) ⊗ N of N spin-1 / A α . Here we follow a different approach anddo this expansion only for the i th factor of the productspace, which yields ρ ss σ ai = 12 (cid:88) α i Tr i [ A α i ρ ss σ ai ] A α i , (22)where Tr i denotes a partial trace over the i th factor of thetensor product Hilbert space. By elementary spin alge-bra, the coefficients of this expansion, which are operator-valued and act on ( C ) ⊗ ( N − , can be written as2 Tr i [ A α i ρ ss σ ai ] = Tr i (cid:104) σ ai ( i + (cid:88) c r cα i σ ci ) ρ ss (cid:105) = Tr i ( σ ai ρ ss ) + r aα i Tr i ( ρ ss ) + i (cid:88) cd ε acd r cα i Tr i (cid:0) σ di ρ ss (cid:1) . (23)Next we rewrite (23) as a linear combination of trace-1operators, such that the time-evolution scheme of Sec.III B can be applied to each of those operators. To thispurpose it is convenient to (partially) expand operatorsin the tensor product basis of Pauli spin operators, wherewe denote the expansion coefficients by s ai = Tr( σ ai ρ ss ) , s abij = Tr( σ ai σ bj ρ ss ) . (24) Starting from (23) we can write i ( A α i ρ ss σ ai ) = (1+ s ai )˜ ρ a ss , (cid:54) i + i (cid:88) cd ε acd r cα i (1+ s di )˜ ρ d ss , (cid:54) i + (cid:104) ( r aα i − − i (cid:88) cd r cα i ε acd (cid:105) ρ ss , (cid:54) i , (25)where we have defined ρ ss , (cid:54) i = Tr i ρ ss , ˜ ρ a ss , (cid:54) i = Tr i [( + σ ai ) ρ ss ]1 + s ai , (26)both of which are trace-1 operators on ( C ) ⊗ ( N − .Those operators can also be expanded in terms of Paulispin operators, and the corresponding expansion coeffi-cients can be expressed in terms of the coefficients (24)of ρ ss ,˜ s a k ,ak = Tr (cid:0) σ a k k ˜ ρ ass, (cid:54) i (cid:1) = s a k k + s aa k ik s ai , (27a)˜ s a j a k ,ajk = Tr (cid:0) σ a j j σ a k k ˜ ρ ass, (cid:54) i (cid:1) = s a j a k jk + s aa j a k ijk s ai , (27b)and so on. Inserting (22) and (25) into (9) we obtainlim t →∞ (cid:10) σ + i ( t ) σ − j ( t + τ ) (cid:11) = 14 (cid:88) α i (cid:8) (1 + s xi )(1 − r zα i )( a x ; xj ; α i ( τ ) − ia y ; xj ; α i ( τ )) + (1 + s yi )(1 − r zα i )( a y ; yj ; α i ( τ ) + ia x ; yj ; α i ( τ ))+ (1 + s zi ) (cid:2) r xα i a x ; zj ; α i ( τ ) + r yα i a y ; zj ; α i ( τ ) + i ( r yα i a x ; zj ; α i ( τ ) − r xα i a y ; zj ; α i ( τ )) (cid:3) + ( r zα i − (cid:2) a xj ; α i ( τ ) + a yj ; α i ( τ ) + i ( a xj ; α i ( τ ) − a yj ; α i ( τ )) (cid:3)(cid:9) (28)with a b ; aj ; α i ( τ ) = Tr (cid:2) σ bj V ( τ, ρ ass, (cid:54) i A α i ) (cid:3) , (29a) a bj ; α i ( τ ) = Tr (cid:2) σ bj V ( τ, ρ ss, (cid:54) i A α i ) (cid:3) . (29b)According to (29a) and (29b), in order to obtain thedesired two-time correlation function (28), we need to It turns out to be advantageous to expand in an overcom-plete basis, using the phase point operators corresponding to r (cid:48) (0 , = (1 , − , r (cid:48) (0 , = ( − , , r (cid:48) (1 , = (1 , , − r (cid:48) (1 , = ( − , − , − ρ ss σ ai directly in terms oftrace-1 operators, is also feasible, but again turned out to per-form worse than the scheme described here. We tested other ways of expressing (23) in terms of trace-1 op-erators, but for our purposes none of them turned out to beadvantageous in terms of accuracy or computational cost. calculate for each i and each phase space operator A α i the time-evolution of four operators, namely ˜ ρ ass, (cid:54) i A α i for a ∈ { x, y, z } , and ρ ss, (cid:54) i A α i . We do this by making useof the method developed in Sec. III B, letting A ...N in(10) take the role of each of the four mentioned opera-tors. The computational cost of the method scales like3 N (3 N − / N . Applying it tothe four trace-1 operators, the eight phase point opera-tors, and the N lattice sites required in (28), (29a), and(29b), results in an overall computational cost that scalesasymptotically like N . IV. RESULTS
We performed simulations for system sizes N = 14,28, 48, 72 and 96 atoms at a fixed low density ρ . Theatoms are placed at random in a spherical volume, butwith the constraint of a minimal distance of one fourth FIG. 1. (a) Fluorescence spectrum for a cloud of density ρ/k = 0 .
1, driven at Ω = 20Γ at resonance (∆ = 0), for N = 14 , , ,
72 and 96 atoms, with the inset showingthe behaviour of the peaks at ω ≈ . Amplitude of theadditional Mollow sidebands as a function of (b) the opticalthickness b for different densities (Ω = 20Γ, N = 72 and∆ = 0) and (c) of the Rabi frequency ( ρ/k = 0 . N = 72and ∆ = 0). The amplitude is defined as (cid:82) + δω − δω | S − S | dω ,where S is the single-atom spectrum, and δω a suitably cho-sen integration range. The dash-dotted line in (c) refers to apower-law fit ( A ≈ . / Γ) − . ). of the mean distance between neighbours, such as torule out unwanted noncooperative effects due to acciden-tally close pairs of atoms. Fig. 1(a) shows the numer-ically computed fluorescence spectra for various N , atfixed density and laser parameters. The three prominentpeaks in the plot at ω = 0 and ± Ω form the Mollowtriplet [19, 20]. The first main result is the observationof additional sidebands in the fluorescence spectrum at ω = ± . These sidebands are genuine quantum effects,as they require the presence of quantum pair correlations.Indeed if connected correlations between different siteswere absent and two-time correlations would factorize, (cid:104) σ ± ,zj σ ± ,zm (cid:105) = (cid:104) σ ± ,zj (cid:105)(cid:104) σ ± ,zm (cid:105) for j (cid:54) = m , one would have (cid:104) σ + j ( t ) σ − m ( t + τ ) (cid:105) = (cid:104) σ + j ( t ) (cid:105)(cid:104) σ − m ( t ) (cid:105) in the steady-stateregime, and in that case the inelastic ( ω (cid:54) = 0) spectrumwould only depend on single-site two-time correlations (cid:104) σ + j ( t ) σ − j ( t + τ ) (cid:105) . The factorizing terms (cid:104) σ ± ,zj (cid:105)(cid:104) σ ± ,zm (cid:105) may modify the local Rabi frequency experienced byeach atom and inhomogeneously broaden the single-atomMollow triplet, but cannot give rise to higher-order side-bands. In other words, a classical treatment of the Hamil-tonian (or, as in this case, Lindbladian) results in the ab-sence of connected correlations between different spins,and no additional sidebands are observed.The novel sidebands are true cooperative effects. Ifthe sidebands were two-atom or few-atom effects, theirpeak height would depend only on the spatial density.In Figs. 1(a) and (b), however, we observe that thesidebands grow with the number of atoms N even atfixed density, and scale linearly with the optical thick-ness b = 2 N/ ( kR ) rather than with the spatial density.This effect can be attributed to the long-range nature ofthe effective interactions (2) between the atomic internal FIG. 2. Left: Fluorescence spectrum for an atomic cloudof density ρ/k = 0 .
1, driven at Ω = 20Γ out of resonance(∆ = − Ω / θ = 0, π/
10, 3 π/ π/
5, 4 π/
5. The asymmetry of the sidebands is clearly visible.Right: The asymmetry of the spectrum in the forward direc-tion ( θ = 0), quantified by ( A − − A + ) / ( A + A + ) where A ± isthe amplitude of the sideband at ± Ω , plotted as a function ofthe optical thickness for different densities and system sizes. degrees of freedom. The scaling with b is reminiscent ofcooperative phenomena in the linear optics regime [35],but is here observed, for the first time to the best ofour knowledge, for a quantum cooperative phenomenonin free space. Furthermore, although single scatteringprocesses may exhibit quantum optics interferences phe-nomena [36], they cannot capture the additional side-bands. These sidebands at ± can be understood asthe first step of a higher order harmonic generation pro-cess, where the next orders could be studied by includ-ing higher-order quantum correlations. However, the lowrelative intensity of approximately 10 makes it hard todetect the peaks at ± experimentally. With that inmind, we searched for traces of quantum cooperativity inthe Mollow triplet bands ( ± Ω ).The second result is the observation that quantum co-operativity breaks the symmetry of the spectrum. Thesingle-atom fluorescence spectrum is always symmetricwith respect to the frequency of the driving light, inde-pendently of the detuning of the driving from the atomicresonance [19]. For large atomic clouds and in the pres-ence of detuning (∆ (cid:54) = 0), it was predicted that coherenceeffects may induce an asymmetry of the Mollow side-bands in the forward scattering direction [26]. Our simu-lations show a similar effect for the scattering of detunedlight, where the Mollow sidebands at ω = ± Ω exhibita significant asymmetry (Fig. 2 left). This asymmetryscales with the optical thickness b (Fig. 2 right), whichconfirms the cooperative nature of this effect. In the ab-sence of quantum correlations the spectrum, being com-posed of the sum of N symmetric spectra, is necessar-ily symmetric, which confirms the genuine quantumnessof the observed asymmetry. However, going beyond theprediction of Ref. [26], we here observe that the asym-metry is also present outside the forward lobe, i.e., forscattering angles θ ≥ /kR where diffuse light dominates(Fig. 2 left). Surprisingly, the asymmetry is inverted inthe forward direction ( θ < /kR ) in comparison with θ > /kR . Experimentally the asymmetry of the stan-dard Mollow sidebands, which reaches ∼
30% in our sim-
FIG. 3. Angular dependence of the fluorescence spectrum ofa cloud of N = 72 atoms with density ρ = 0 . k , driven bya field with Ω = 20Γ at resonance, which corresponds to asaturation parameter s = 3200. The inset shows the elasticand integrated inelastic spectra ( S el ( θ ) = S ( ω = 0 , θ ) and S inel ( θ ) = (cid:82) ω (cid:54) =0 S ( ω, θ ) dω ), respectively, as discussed in thetext. ulations, should be relatively easy to detect.In Fig. 3 we show the fluoresence spectrum as a func-tion of the scattering angle θ and the frequency ω inthe regime of deep saturation, where most of the light isexpected to be scattered inelastically ( ω (cid:54) = 0). In thisregime the portion of elastically scattered light for a sin-gle atom goes as 1 /s = (∆ + Γ / / at large Ω , s being referred as the saturation parameter, so mostof the light is scattered inelastically. We clearly observethe quasi-isotropic inelastic Mollow triplet and higher or-der sidebands. For the parameters considered, a strongelastic component is particularly visible in the forwarddirection (see inset), which we attribute to the construc-tive interference of the elastic component of the electricfield in the forward direction (which, according to linearoptics, is expected to scale like N with the system size).This indicates that signatures of quantum cooperativity,which are intimately connected to inelastic scattering,may be more easily detected at larger scattering angles,and not in the forward direction. We note however thatin the forward direction the inelastic component exhibitsa small dip. The physical origin of this feature remainsto be understood. V. CONCLUSIONS
We have reported the discovery of signatures of quan-tum cooperativity in the fluorescence spectrum of largedilute atomic clouds. The rise of additional sidebandsat frequencies ± from the central line, as well as the asymmetry in the spectrum of the cloud driven outof resonance, are identified as proper quantum effectsthat cannot occur in the absence of genuine quantumcorrelations. Moreover, by analyzing parameter depen-dences and scaling properties, the cooperative nature ofthe observed phenomena is revealed. Cooperativity is ul-timately related to the long-range character of the effec-tive atom–atom interactions induced by the photon field.The deficiency in inelastically scattered photons in theforward cone ( | θ | ≤ /kR ) is particularly interesting, be-cause it implies that the forward direction, which has longbeen considered a natural candidate for probing cooper-ative phenomena in the linear-optics regime [37], maybe less suitable for probing quantum cooperative effects.Furthermore, while the second-order optical coherence isusually considered the ideal candidate for revealing thequantum nature of the light scattering by atoms, withphoton bunching [38] and anti-bunching [39] as paradig-matic signatures, we show in this paper that the first-order optical coherence g (1) , which witnesses the quan-tum nature of the atom-light coupling, already showsclear signatures of quantum cooperativity. More gener-ally, our results suggest that the quantum optics regimeof an optically deep system is substantially richer than itssingle-atom physics, and holds much promise for furtherstudies of cooperative effects. This may become relevantfor neutral atom optical clocks or other long-range quan-tum systems such as Rydberg atoms. To gain access tothis regime on the computational side, the simulationtechnique developed in this paper, based on a truncationof the hierarchy of correlations, proves to be a powerfultool.The cooperative nature of the observed effects suggeststhat dilute atomic clouds might be used as experimentalplatforms for quantum-simulating plasmas, free electronlasers, and other quantum long-range interacting systemsin which cooperativity plays an essential role. ACKNOWLEDGMENTS
We thank N. Piovella for stimulating discussions. A.R.acknowledges financial support from the National Re-search Foundation of South Africa. This work has re-ceived financial support from FAPESP and CNPq. M.K.acknowledges financial support from the National Re-search Foundation of South Africa via the IncentiveFunding and the Competitive Programme for Rated Re-searchers. M.K. and R.B. acknowledge financial supportby a FAPESP/SU collaboration grant.
Appendix A: Time-evolution of the Pauli expansioncoefficients
We expand all operators in (18a) and (18b), except C ijm that is neglected, in terms of Pauli spin operators, A i = 12 ( + a i · σ i ) , C ij = 14 (cid:88) a,b ∈{ x,y,z } c abij σ ai σ bj . (A1)Inserting these expansions into (19a) and the truncatedversion of (19b), and making use of Lindblad equations(10)–(13), we obtain time-evolution equations for thePauli expansion coefficients, ∂ t a ai = (cid:88) b a bi (cid:8) − ∆ ε zba + Ω (cid:2) ε xba cos( − k · r i ) + ε yba sin( − k · r i ) (cid:3)(cid:9) − Γ (cid:2) a ai (1 + δ za ) + δ za (cid:3) + (cid:88) b (cid:88) m (cid:54) = i (cid:110) ε xba (cid:104) ∆ im (cid:0) a bi a xm + c bxim (cid:1) − γ im (cid:16) a bi a ym + c byim (cid:17)(cid:105) + ε yba (cid:104) ∆ im (cid:16) a bi a ym + c byim (cid:17) + γ im (cid:0) a bi a xm + c bxim (cid:1)(cid:105)(cid:111) (A2)and ∂ t c abij = (cid:88) c c cbij {− ∆ ε zca + Ω [cos( − k · r i ) ε xca + sin( − k · r i ) ε yca ] } (A3)+ (cid:88) c c acij (cid:8) − ∆ ε zcb + Ω (cid:2) cos( − k · r j ) ε xcb + sin( − k · r j ) ε ycb (cid:3)(cid:9) − Γ c abij (cid:18) δ az + δ bz (cid:19) − γ ij (cid:88) c,d (cid:0) c cdij + a ci a dj (cid:1) (cid:0) ε xca ε dxb + ε yca ε dyb (cid:1) + (cid:88) c c cbij (cid:88) m (cid:54) = ij (cid:104) a xm (cid:16) ∆ im ε xca + γ im ε yca (cid:17) + a ym (cid:16) ∆ im ε yca − γ im ε xca (cid:17)(cid:105) + (cid:88) c c acij (cid:88) m (cid:54) = ij (cid:104) a xm (cid:16) ∆ jm ε xcb + γ jm ε ycb (cid:17) + a ym (cid:16) ∆ jm ε ycb − γ jm ε xcb (cid:17)(cid:105) + (cid:88) c a cj (cid:104) δ ax (cid:16) ∆ ij ε xcb + γ ij ε ycb (cid:17) + δ ay (cid:16) ∆ ij ε ycb − γ ij ε xcb (cid:17)(cid:105) + (cid:88) c a ci (cid:104) δ bx (cid:16) ∆ ij ε xca + γ ij ε yca (cid:17) + δ by (cid:16) ∆ ij ε yca − γ ij ε xca (cid:17)(cid:105) − (cid:88) c a bj (cid:104)(cid:0) c cxij + a ci a xj (cid:1) (cid:16) ∆ ij ε xca + γ ij ε yca (cid:17) + (cid:0) c cyij + a ci a yj (cid:1) (cid:16) ∆ ij ε yca − γ ij ε xca (cid:17)(cid:105) − (cid:88) c a ai (cid:104)(cid:0) c xcij + a xi a cj (cid:1) (cid:16) ∆ ij ε xcb + γ ij ε ycb (cid:17) + (cid:0) c ycij + a yi a cj (cid:1) (cid:16) ∆ ij ε ycb − γ ij ε xcb (cid:17)(cid:105) + (cid:88) c a ci (cid:88) m (cid:54) = ij (cid:104)(cid:16) c xbmj ∆ im − c ybmj γ im (cid:17) ε xca + (cid:16) c ybmj ∆ im + c xbmj γ im (cid:17) ε yca (cid:105) + (cid:88) c a cj (cid:88) m (cid:54) = ij (cid:104)(cid:16) c axim ∆ jm − c ayim γ jm (cid:17) ε xcb + (cid:16) c ayim ∆ jm + c axim γ jm (cid:17) ε ycb (cid:105) + (cid:88) c (cid:88) m (cid:54) = ij (cid:104)(cid:16) c cbxijm ∆ im − c cbyijm γ im (cid:17) ε xca + (cid:16) c cbyijm ∆ im + c cbxijm γ im (cid:17) ε yca (cid:105) + (cid:88) c (cid:88) m (cid:54) = ij (cid:104)(cid:16) c acxijm ∆ jm − c acyijm γ jm (cid:17) ε xcb + (cid:16) c acyijm ∆ jm + c acyijm γ jm (cid:17) ε ycb (cid:105) . These equations form a set of coupled ordinary differ-ential equations, which can be integrated by standard numerical methods.We calculate the steady-state density operator ρ ss by FIG. 4. Fluorescence spectrum for a cloud of N = 7 atoms,driven at Ω = 20Γ at resonance (∆ = 0), and for densities ρ/k = 0 .
03 (black), 0 . . setting to zero the left-hand sides of (A2) and (A3), andnumerically solving the resulting algebraic equations by astandard Newton-Krylov solver. The stationary values ofthe a - and c -coefficients encode the required information on ρ ss . Appendix B: Benchmarking against exact results
The accuracy of the proposed simulation method istested by benchmarking the fluorescence spectrum S against exact results. Up to N = 7 spins (atoms) couldbe dealt with exactly by using the “Quantum Toolboxin Python” [40, 41], a module tailored to simulate thedynamics of open quantum systems and especially thoseof quantum optics.As shown in Fig. 4 for Rabi frequency Ω = 20Γ, fordensities up to ρ/k = 0 . ρ/k = 0 . leadto a similar degree of agreement (not shown).Besides the spectra, we also benchmarked other rele-vant quantities, including the steady state ρ ss calculatedaccording to Sec. C, as well as the the two-time correla-tions evolved from the latter. All show very good agree-ment for densities up to ρ/k ∼ . [1] J. Simon, W. S. Bakr, R. Ma, M. E. Tai, P. M. Preiss,and M. Greiner, “Quantum simulation of antiferromag-netic spin chains in an optical lattice,” Nature , 307–312 (2011).[2] J. W. Britton, B. C. Sawyer, A. C. Keith, C.-C. J.Wang, J. K. Freericks, H. Uys, M. J. Biercuk, andJ. J. Bollinger, “Engineered two-dimensional Ising inter-actions in a trapped-ion quantum simulator with hun-dreds of spins,” Nature , 489–492 (2012).[3] R. Islam, C. Senko, W. C. Campbell, S. Korenblit,J. Smith, A. Lee, E. E. Edwards, C.-C. J. Wang, J. K.Freericks, and C. Monroe, “Emergence and frustrationof magnetism with variable-range interactions in a quan-tum simulator,” Science , 583–587 (2013).[4] I. M. Georgescu, S. Ashhab, and F. Nori, “Quantumsimulation,” Rev. Mod. Phys. , 153–185 (2014).[5] R. H. Dicke, “Coherence in spontaneous radiation pro-cesses,” Phys. Rev. , 99–110 (1954).[6] R. A. de Oliveira, M. S. Mendes, W. S. Martins, P. L.Saldanha, J. W. R. Tabosa, and D. Felinto, “Single-photon superradiance in cold atoms,” Phys. Rev. A ,023848 (2014).[7] Ph. W. Courteille, S. Bux, E. Lucioni, K. Lauber, T. Bi-enaim´e, R. Kaiser, and N. Piovella, “Modification of ra-diation pressure due to cooperative scattering of light,”Eur. Phys. J. D , 69–73 (2010).[8] H. Bender, C. Stehle, S. Slama, R. Kaiser, N. Piovella,C. Zimmermann, and Ph. W. Courteille, “Observationof cooperative Mie scattering from an ultracold atomiccloud,” Phys. Rev. A , 011404 (2010). [9] R. Friedberg, S. R. Hartmann, and J. T. Manassah,“Frequency shifts in emission and absorption by reso-nant systems ot two-level atoms,” Phys. Rep. , 101–179(1973).[10] R. R¨ohlsberger, K. Schlage, B. Sahoo, S. Couet, andR. R¨uffer, “Collective Lamb shift in single-photon super-radiance,” Science , 1248–1251 (2010).[11] J. Keaveney, A. Sargsyan, U. Krohn, I. G. Hughes,D. Sarkisyan, and C. S. Adams, “Cooperative Lamb shiftin an atomic vapor layer of nanometer thickness,” Phys.Rev. Lett. , 173601 (2012).[12] S. D. Jenkins, J. Ruostekoski, J. Javanainen, R. Bour-gain, S. Jennewein, Y. R. P. Sortais, and A. Browaeys,“Optical resonance shifts in the fluorescence of thermaland cold atomic gases,” Phys. Rev. Lett. , 183601(2016).[13] S. L. Bromley, B. Zhu, M. Bishof, X. Zhang, T. Bothwell,J. Schachenmayer, J. T. L. Nicholson, R. Kaiser, S. F.Yelin, M. D. Lukin, A.M. Rey, and J. Ye, “Collectiveatomic scattering and motional effects in a dense coherentmedium,” Nat. Comm. , 11039 (2016).[14] S. J. Roof, K. J. Kemp, M. D. Havey, and I. M. Sokolov,“Observation of single-photon superradiance and the co-operative Lamb shift in an extended sample of coldatoms,” Phys. Rev. Lett. , 073003 (2016).[15] M. O. Ara´ujo, I. Kreˇsi´c, R. Kaiser, and W. Guerin, “Su-perradiance in a large and dilute cloud of cold atoms inthe linear-optics regime,” Phys. Rev. Lett. , 073002(2016).[16] W. Guerin, M. O. Ara´ujo, and R. Kaiser, “Subradiance in a large cloud of cold atoms,” Phys. Rev. Lett. ,083601 (2016).[17] R. T. Sutherland and F. Robicheaux, “Coherent for-ward broadening in cold atom clouds,” Phys. Rev. A ,023407 (2016).[18] B. Zhu, J. Cooper, J. Ye, and A. M. Rey, “Light scat-tering from dense cold atomic media,” Phys. Rev. A ,023612 (2016).[19] B. R. Mollow, “Power spectrum of light scattered by two-level systems,” Phys. Rev. , 1969–1975 (1969).[20] F Schuda, C R Stroud Jr., and M Hercher, “Observationof the resonant Stark effect at optical frequencies,” J.Phys. B , L198–L202 (1974).[21] G. S. Agarwal and Robert W. Boyd, “Quantum theory ofrabi sideband generation by forward four-wave mixing,”Physical Review A , 4019–4027 (1988).[22] A. J. Traverso, C. O’ Brien, B. H. Hokr, J. V. Thomp-son, L. Yuan, C. W. Ballmann, A. A. Svidzinsky, G. I.Petrov, M. O. Scully, and V. V. Yakovlev, “Directionalcoherent light via intensity-induced sideband emission,”Light: Science & Applications , e16262 (2017).[23] I. R. Senitzky, “Sidebands in strong-field resonance fluo-rescence,” Phys. Rev. Lett. , 1334–1337 (1978).[24] G. S. Agarwal, R. Saxena, L. M. Narducci, D. H. Feng,and R. Gilmore, “Analytical solution for the spectrumof resonance fluorescence of a cooperative system of twoatoms and the existence of additional sidebands,” Phys.Rev. A , 257–259 (1980).[25] Y. Ben-Aryeh and C. M. Bowden, “Resonance fluores-cence spectra of two driven two-level atoms,” IEEE J.Quant. Electron. , 1376–1382 (1988).[26] J. R. Ott, M. Wubs, P. Lodahl, N. A. Mortensen, andR. Kaiser, “Cooperative fluorescence from a stronglydriven dilute cloud of atoms,” Phys. Rev. A , 061801(2013).[27] W. K. Wootters, “A Wigner-function formulation offinite-state quantum mechanics,” Ann. Phys. , 1–21(1987).[28] J. Schachenmayer, A. Pikovski, and A. M. Rey, “Many-body quantum spin dynamics with Monte Carlo trajec-tories on a discrete phase space,” Phys. Rev. X , 011022(2015).[29] L. Pucci, A. Roy, and M. Kastner, “Simulation of quan- tum spin dynamics by phase space sampling of Bogo-liubov-Born-Green-Kirkwood-Yvon trajectories,” Phys.Rev. B , 174302 (2016).[30] R. H. Lehmberg, “Radiation from an N -atom system. I.General formalism,” Phys. Rev. A , 883–888 (1970).[31] G. Angelatos and S. Hughes, “Entanglement dynamicsand Mollow nonuplets between two coupled quantumdots in a nanowire photonic-crystal system,” Phys. Rev.A , 051803 (2015).[32] C. W. Gardiner and P. Zoller, The Quantum World ofUltra-Cold Atoms and Light: Foundations of QuantumOptics (Imperial College Press, London, 2015).[33] D. A. Steck, “Quantum and atom optics,” Sec. 5.7.3,available online at http://steck.us/teaching , (revision0.11.5, 27 November 2016).[34] M. Bonitz, Quantum Kinetic Theory (Springer, Cham,2015).[35] W. Guerin, M. T. Rouabah, and R. Kaiser, “Light in-teracting with atomic ensembles: collective, cooperativeand mesoscopic effects,” J. Mod. Opt. , 1–13 (2016).[36] Philippe Grangier, G´erard Roger, Alain Aspect, AntoineHeidmann, and Serge Reynaud, “Observation of pho-ton antibunching in phase-matched multiatom resonancefluorescence,” Phys. Rev. Lett. , 687–690 (1986).[37] M. O. Scully, E. S. Fry, C. H. R. Ooi, and K. W´odkiewicz,“Directed spontaneous emission from an extended ensem-ble of N atoms: Timing is everything,” Phys. Rev. Lett. , 010501 (2006).[38] B. L. Morgan and L. Mandel, “Measurement of photonbunching in a thermal light beam,” Phys. Rev. Lett. ,1012–1015 (1966).[39] H. J. Kimble, M. Dagenais, and L. Mandel, “Photonantibunching in resonance fluorescence,” Phys. Rev. Lett. , 691–695 (1977).[40] J. R. Johansson, P. D. Nation, and F. Nori, “QuTiP:An open-source Python framework for the dynamics ofopen quantum systems,” Comput. Phys. Commun. ,1760–1772 (2012).[41] J. R. Johansson, P. D. Nation, and F. Nori, “QuTiP 2:A Python framework for the dynamics of open quan-tum systems,” Comput. Phys. Commun.184