Quantum effects in the cooperative scattering of light by atomic clouds

TL;DR

This paper introduces a new simulation method to explore quantum cooperative effects in light scattering by atomic clouds, revealing unique spectral signatures like additional sidebands and asymmetries in the fluorescence spectrum.

Contribution

The study presents a novel simulation technique enabling the investigation of quantum effects in large atomic clouds' light scattering dynamics.

Findings

Observation of additional sidebands at twice the classical Mollow sidebands

Detection of asymmetry in the Mollow triplet

Scaling behavior of quantum signatures with atom number

Abstract

Scattering of classical light by atomic clouds induces photon-mediated effective long-range interactions between the atoms and leads to cooperative effects even at low atomic densities. We introduce a novel simulation technique that allows us to investigate the quantum regime of the dynamics of large clouds of atoms. We show that the fluorescence spectrum of the cloud can be used to probe genuine quantum cooperative effects. Signatures of these effects are the occurrence, and the scaling behavior, of additional sidebands at twice the frequency of the classical Mollow sidebands, as well as an asymmetry of the Mollow triplet.

Figures

Equations

$$\frac{d{\sigma }_j^{-}}{dt}=$$

$$+\frac{\Gamma }{2}\sum _{m\ne j}^N{\sigma }_j^z{\sigma }_m^{-}({\gamma }_{jm}-i{\Delta }_{jm}),$$

$$\frac{d{\sigma }_j^z}{dt}=$$

$$-\Gamma \sum _{m\ne j}^N({\sigma }_m^{+}{\sigma }_j^{-}({\gamma }_{jm}+i{\Delta }_{jm})+\text{h.c.}),$$

$${\gamma }_{jm}=\frac{\sin (k_0|{\mathbf{r}}_j-{\mathbf{r}}_m|)}{k_0|{\mathbf{r}}_j-{\mathbf{r}}_m|},{\Delta }_{jm}=\frac{\cos (k_0|{\mathbf{r}}_j-{\mathbf{r}}_m|)}{k_0|{\mathbf{r}}_j-{\mathbf{r}}_m|}$$

$$S(\omega )=\underset{T\to \infty }{\lim }\underset{t\to \infty }{\lim }{\int }_{-T}^{T}d\tau g^{(1)}(t,\tau )e^{-i\omega \tau },$$

$$g^{(1)}(t,\tau )=\frac{⟨E(\hat{\mathbf{n}},t)E^{†}(\hat{\mathbf{n}},t+\tau )⟩}{⟨E(\hat{\mathbf{n}},t)E^{†}(\hat{\mathbf{n}},t)⟩}$$

$$E^{†}(\hat{\mathbf{n}},t)\propto \sum _{j=1}^N{\sigma }_j^{-}(t)e^{-ik\hat{\mathbf{n}}\cdot {\mathbf{r}}_j}.$$

$$\begin{gathered} S(\omega )=\underset{T\to \infty }{\lim }\underset{t\to +\infty }{\lim }{\int }_{-T}^{T}d\tau e^{-i\omega \tau } \\ \times \left[g^{(1)}(t,\tau )\Theta (\tau )+g^{(1)}(-t,\tau )\Theta (-\tau )\right], \end{gathered}$$

$$⟨{\sigma }_i^a(t){\sigma }_j^b(t+\tau )⟩=\mathrm{Tr}\left\{V(t+\tau ,t)\left[\left(V(t,0)\rho \right){\sigma }_i^a\right]{\sigma }_j^b\right\},$$

$${\rho }_{\text{ss}}=\underset{t\to \infty }{\lim }V(t,0)\rho .$$

$$\underset{t\to \infty }{\lim }⟨{\sigma }_i^a(t){\sigma }_j^b(t+\tau )⟩=\mathrm{Tr}\left[{\sigma }_j^{b}V(\tau ,0)\left({\rho }_{\text{ss}}{\sigma }_i^a\right)\right].$$

$$i{\partial }_t{\mathcal{A}}_{1\dots N}=\mathcal{L}{\mathcal{A}}_{1\dots N},$$

$$\mathcal{L}=\sum _i{\mathcal{L}}_i+\sum _{ij}{\mathcal{L}}_{ij}$$

$$\begin{gathered} {\mathcal{L}}_i{\mathcal{A}}_{1\dots N}=-\frac{{\Delta }_0}{2}\left[{\sigma }_i^z,{\mathcal{A}}_{1\dots N}\right] \\ +\frac{{\Omega }_0}{2}\left[e^{-i{\mathbf{k}}_0\cdot {\mathbf{r}}_i}{\sigma }_i^{-}+e^{i{\mathbf{k}}_0\cdot {\mathbf{r}}_i}{\sigma }_i^{+},{\mathcal{A}}_{1\dots N}\right] \end{gathered}$$

$$\begin{gathered} {\mathcal{L}}_{ij}{\mathcal{A}}_{1\dots N}={\Delta }_{ij}\left[{\sigma }_i^{+}{\sigma }_j^{-},{\mathcal{A}}_{1\dots N}\right]+ \\ i{\gamma }_{ij}\left({\sigma }_j^{-}{\mathcal{A}}_{1\dots N}{\sigma }_i^{+}-\frac{1}{2}{\sigma }_i^{+}{\sigma }_j^{-}{\mathcal{A}}_{1\dots N}-\frac{1}{2}{\mathcal{A}}_{1\dots N}{\sigma }_i^{+}{\sigma }_j^{-}\right), \end{gathered}$$

$${\Delta }_{ij}$$

$${\gamma }_{ij}$$

$$i{\partial }_t{\mathcal{A}}_i$$

$$i{\partial }_t{\mathcal{A}}_{ij}$$

$$+\sum _{m\ne i,j}{\mathrm{Tr}}_m\left({\mathcal{L}}_{im}+{\mathcal{L}}_{jm}\right){\mathcal{A}}_{ijm},$$

$${\mathcal{A}}_i={\mathrm{Tr}}_{\{k\ne i\}}{\mathcal{A}}_{1\dots N},{\mathcal{A}}_{ij}={\mathrm{Tr}}_{\{k\ne i,j\}}{\mathcal{A}}_{1\dots N},$$

$${\mathcal{A}}_{ij}$$

$${\mathcal{A}}_{ijm}$$

$$i{\partial }_t{\mathcal{A}}_i$$

$$i{\partial }_t{\mathcal{C}}_{ij}$$

$$-{\mathcal{A}}_i{\mathrm{Tr}}_i\left[{\mathcal{L}}_{ij}^S({\mathcal{C}}_{ij}+{\mathcal{A}}_i{\mathcal{A}}_j)\right]$$

$$-{\mathcal{A}}_j{\mathrm{Tr}}_j\left[{\mathcal{L}}_{ij}^S({\mathcal{C}}_{ij}+{\mathcal{A}}_i{\mathcal{A}}_j)\right]$$

$$+\sum _{m\ne i,j}{\mathrm{Tr}}_m\left[{\mathcal{L}}_{im}^S({\mathcal{A}}_i{\mathcal{C}}_{jm}+{\mathcal{A}}_m{\mathcal{C}}_{ij}+{\mathcal{C}}_{ijm})\right]$$

$$+\sum _{m\ne i,j}{\mathrm{Tr}}_m\left[{\mathcal{L}}_{jm}^S({\mathcal{A}}_j{\mathcal{C}}_{im}+{\mathcal{A}}_m{\mathcal{C}}_{ij}+{\mathcal{C}}_{ijm})\right],$$

$$\Gamma =\left\{(0,0),(0,1),(1,0),(1,1)\right\}$$

$$A_{\alpha }=\frac{1}{2}(\text{mathds}1+{\mathbf{r}}_{\alpha }\cdot \mathbf{\sigma }),$$

$${\rho }_{\text{ss}}{\sigma }_i^a=\frac{1}{2}\sum _{{\alpha }_i}{\mathrm{Tr}}_i\left[A_{{\alpha }_i}{\rho }_{\text{ss}}{\sigma }_i^a\right]A_{{\alpha }_i},$$

$$2\operatorname{Tr}_{i}\left[A_{\alpha_{i}}\rho_{\text{ss}}\sigma_{i}^{a}\right]=\operatorname{Tr}_{i}\Bigl{[}\sigma_{i}^{a}({\mathds{1}}_{i}+\sum_{c}r_{\alpha_{i}}^{c}\sigma_{i}^{c})\rho_{\text{ss}}\Bigr{]}\\ =\operatorname{Tr}_{i}\left(\sigma_{i}^{a}\rho_{\text{ss}}\right)+r_{\alpha_{i}}^{a}\operatorname{Tr}_{i}\left(\rho_{\text{ss}}\right)+i\sum_{cd}\varepsilon^{acd}r_{\alpha_{i}}^{c}\operatorname{Tr}_{i}\left(\sigma_{i}^{d}\rho_{\text{ss}}\right).$$

$$s_i^a=\mathrm{Tr}({\sigma }_i^a{\rho }_{\text{ss}}),s_{ij}^{ab}=\mathrm{Tr}({\sigma }_i^a{\sigma }_j^b{\rho }_{\text{ss}}).$$

$$\!\!2\operatorname{Tr}_{i}\left(A_{\alpha_{i}}\rho_{\text{ss}}\sigma_{i}^{a}\right)=(1+s_{i}^{a})\tilde{\rho}^{a}_{\text{ss},\not{\,i}}+i\sum_{cd}\varepsilon^{acd}r_{\alpha_{i}}^{c}(1+s_{i}^{d})\tilde{\rho}^{d}_{\text{ss},\not{\,i}}\\ +\Bigl{[}(r_{\alpha_{i}}^{a}-1)-i\sum_{cd}r_{\alpha_{i}}^{c}\varepsilon^{acd}\Bigr{]}\rho_{\text{ss},\not{\,i}},$$

Full text

I Introduction

Our understanding of the quantum dynamics of many-body systems has benefited from a number of recent achievements. On the experimental side, cold atom systems and ion traps have reached an unprecedented level of control and allow for the emulation of a large variety of many-body Hamiltonians of interest, including the possibility of tuning coupling parameters Simon et al. (2011); Britton et al. (2012); Islam et al. (2013); Georgescu et al. (2014). On the numerical side, progress in the understanding of matrix product state has boosted density matrix renormalization group and other methods, applicable primarily to one-dimensional quantum systems. Reliable simulations of higher-dimensional systems are in general more difficult, and in many cases impossible.

A three-dimensional setting where many exciting discoveries have been made, but also many open questions remain, is light scattering by large clouds of atoms. The classical, linear optics regime of such systems has been extensively studied, and many-body effects, such as superradiance Dicke (1954); de Oliveira et al. (2014), modification of the radiation pressure force Courteille et al. (2010); Bender et al. (2010) or cooperative frequency shifts Friedberg et al. (1973); Röhlsberger et al. (2010); Keaveney et al. (2012); Jenkins et al. (2016); Bromley et al. (2016), were reported. The effective coupling between the atoms in the cloud, mediated by the photon field, turns out to be long-ranged, and as a result cooperative effects occur even in dilute clouds, including superradiance Roof et al. (2016); Araújo et al. (2016), Dicke subradiance Guerin et al. (2016a), and spectral broadening Sutherland and Robicheaux (2016). The strength of the cooperative effects depends in these cases on the cloud optical thickness, not the spatial density.

Beyond the linear optics regime, saturation effects give the atoms a nonlinear response to the electric field of the light. When many atoms are simultaneously in the excited state, the classical coupling of these nonlinear oscillators leads, for example, to the modification of the line shape of the atomic resonance Zhu et al. (2016). Entering the field of quantum optics, two-time correlations of the radiated field provide information on the optical coherence of the first kind $g^{(1)}$ and on the fluorescence spectrum. For a single atom, saturation leads to the emergence of the Mollow triplet, a trio of spectral lines of inelastically scattered light, one around the laser frequency and two symmetric ones shifted by the generalized Rabi frequency Mollow (1969); Schuda et al. (1974); see Fig. 1 (left) for an illustration. In a four-wave mixing configuration, these symmetric bands can be amplified as the wave propagates within the cloud Agarwal and Boyd (1988); Traverso et al. (2017). Furthermore, this redistribution of frequencies is iterated as two scatterers interact through their radiation Senitzky (1978); Agarwal et al. (1980), which, for dense atomic clouds, results in the presence of additional sidebands at twice the Rabi frequency for pairs of atoms closer than a wavelength Ben-Aryeh and Bowden (1988). In the many-body case, quantum correlations become essential when trying to understand the quantum features of the collective response of the system. This quantum regime of cooperative light scattering in free space is essentially unexplored. The main obstacle for a quantum mechanical treatment is, as usual, the exponential growth of the Hilbert space with the number of atoms. Additionally, and different from, e.g., optical cavities, no suitable collective operators are known for describing the dynamics effectively. A recent first study of quantum effects in the light scattering by dilute atomic clouds made use of perturbative techniques, valid in the regime of strong driving, and was able to predict asymmetries in the scattering spectrum, but no additional spectral lines beyond those of a single atom Ott et al. (2013).

In this paper we report the discovery of quantum cooperative effects in light scattering by large dilute clouds of atoms. We show theoretically that quantum correlations building up in the cloud induce cooperative sidebands in the many-body fluorescence spectrum at twice the Rabi frequency, and also lead to a cooperative spectral asymmetry in the Mollow triplet. Both effects scale with the optical thickness, which is a hallmark of their cooperative character. Investigating the angular dependence of the scattering spectrum, we find that quantum cooperativity is more easily detected at large scattering angles, and not in the forward direction. These results provide guidance on where to look for quantum cooperative effects in light scattering experiments on atomic clouds, and are expected to be relevant also for neutral-atom optical clocks, Rydberg atoms, and other settings where effective long-range interactions play a role.

This study of nonperturbative quantum effects in fairly large three-dimensional atomic clouds became possible by a novel simulation technique, combining a discrete phase-space representation of spins Wootters (1987); Schachenmayer et al. (2015) with higher-order semi-classical evolution equations Pucci et al. (2016) extended to driven-dissipative Lindblad dynamics. The method is highly accurate for higher-dimensional systems with long-range interactions, and therefore a perfect fit for the problem at hand.

II Modeling the atomic cloud

For our purposes the cloud of atoms is modelled as an assembly of $N$ two-level systems at fixed, but usually random, positions $\mathbf{r}_{i}$ in three-dimensional space. Using a scalar light model, which is valid for dilute clouds, each two-level system can be described by Pauli spin operators $\sigma^{\pm}=(\sigma^{x}\pm i\sigma^{y})/2$ and $\sigma^{z}$. This model provides a good description of dilute clouds of atoms cooled below the Doppler limit and trapped in a magneto-optical or dipolar trap. Transitions between the two levels of each atom are driven by a classical planar-wave laser light-field of wave vector $\mathbf{k}_{0}$, Rabi frequency $\Omega_{0}$, and detuning $\Delta_{0}=\omega_{0}-\omega_{a}$ from the optical transition. The dynamics is then described by a set of equations that couple the spin degrees of freedom to the photon field Lehmberg (1970). Performing the rotating-wave, Born and Markov approximations and eliminating the photon degrees of freedom, equations of motion can be derived for the atomic internal degrees of freedom in the Heisenberg picture Lehmberg (1970); Ott et al. (2013),

[TABLE]

where the coefficients

[TABLE]

describe the spatial dependence of the light-mediated long-range coupling between the atoms, and $\Gamma$ is the transition linewidth of the two-level atoms. Note that the dilute regime allows to neglect near-field terms, which scale as $1/r^{3}$, and thus focus on long-range effects.

The quantity we want to study, and which is accessible in light-scattering experiments, is the fluorescence spectrum

[TABLE]

defined as the Fourier transform of the first-order optical coherence

[TABLE]

in the steady state [realized by the limit $t\to\infty$ in (3)]. The optical coherence is defined in terms of the electric field operator $E$ which, in the far-field limit and emitted in the direction of the normalized vector $\hat{\bm{n}}$, is given by

[TABLE]

From Eqs. (3)–(5) one can read off that the crucial nontrivial ingredient for calculating the fluorescence spectrum are the two-time spin–spin correlation functions $\langle\sigma_{j}^{+}(t)\sigma_{j}^{-}(t+\tau)\rangle$, time-evolved under the evolution equations (1a) and (1b).

III Numerical method

An analytic calculation of the fluorescence spectrum $S$ was reported in Refs. Senitzky (1978); Agarwal et al. (1980) for the case of two atoms pumped at resonance $\omega_{a}$. In that work, sidebands at frequencies $\omega_{a}\pm 2\Omega$ were shown to rise in the spectrum, yet the effect is substantial only for atoms closer than a wavelength. A similar effect was predicted for pairs of quantum dots Angelatos and Hughes (2015). Treating more than two atoms is much harder. An exact numerical evaluation of the evolution equations (1a) and (1b) is possible for about a dozen atoms or so, but these numbers are too small to investigate scaling behavior or extrapolate to experimentally relevant system sizes.

Here we deal with these difficulties by developing a simulation method for driven-dissipative quantum mechanical evolution equations. Our method builds upon a simulation technique that makes use of a discrete phase-space representation of the Pauli operators Wootters (1987); Schachenmayer et al. (2015), and time-evolves phase space points as well as correlation coefficients according to semiclassical evolution equations Pucci et al. (2016). We developed these techniques further by extending them to nonunitarily-evolving driven-dissipative systems, and made the method applicable to the computation of two-time correlation functions by making use of the quantum regression theorem. The main features of the method are: (i) The use of semiclassical time-evolution equations makes this method particularly suitable for systems with long-range interactions and systems in higher spatial dimension. (ii) The explicit incorporation of correlation coefficients in the quasiclassical equations of motion allows for the accurate computation of correlation functions; for $N=2$ the exact results are recovered. (iii) The numerical cost scales polynomially, not exponentially, with the system size $N$; hundred and more atoms can be treated.

The novel numerical method, described in the following, is developed for calculating the fluorescence spectrum (3) of a cloud of two-level atoms as described by the time-evolution equations (1a) and (1b), but is also applicable, more generally, for the calculation of two-time correlation functions in long-range spin models. The atomic cloud we consider is three-dimensional, and the atom-–atom interactions are long-ranged, which makes the system particularly suitable for simulation methods based on the quasiclassical time-evolution of discrete phase-space points Schachenmayer et al. (2015); Pucci et al. (2016). However, these methods are formulated for unitarily evolving quantum spin systems and give access to equal-time correlation functions, but not to two-time correlations. In Sec. III.1 we rewrite the quantities needed for calculating the fluorescence spectrum in a form that is more amenable to a semi-classical time-evolution. In Sec. III.2 we derive semi-classical equations of motion for driven-dissipative, non-unitarily evolving quantum spin systems. In Sec. III.3 we introduce a discrete phase space representation and apply it to the calculation of the two-time correlation functions that are required for obtaining the fluorescence spectrum. The resulting novel simulation scheme is benchmarked against exact results in Appendix B.

III.1 Quantum regression

We extract the spectral properties of the radiated light in direction $\hat{\bm{n}}$ from the first order optical coherence $g^{(1)}(t,\tau)$ (4) by evaluating the expression for the spectrum (3), which can be rewritten as

[TABLE]

where $\Theta$ denotes the Heaviside step function. $g(-t,\tau)$ is a short-hand notation for flipping the sign of $t$ of the unitary terms in (1b) while keeping the dissipative ones unchanged, which corresponds to time-reversing the dynamics of the total system (two-level atoms plus photons) before eliminating the photonic field. Taking the limit $t\rightarrow+\infty$ inside the integral, we see that $\lim_{t\to\infty}g^{(1)}(t,\tau)$ and $\lim_{t\to\infty}g^{(1)}(-t,\tau)$ are the two limits required for calculating $S$. Here, for being definite, we discuss the first case only, as the second can be treated analogously. According to Eqs. (4) and (5), the optical coherence $g^{(1)}(t,\tau)$ can be expressed in terms of two-time correlations of Pauli spin operators, and hence our main object of interest will be $\lim_{t\to\infty}\left<{\sigma_{i}^{a}(t)\sigma_{j}^{b}(t+\tau)}\right>$ for all $a,b\in\{x,y,z\}$.

Making use of the quantum regression theorem Gardiner and Zoller (2015); Steck , we can write the two-time spin–spin correlation as

[TABLE]

where $\rho$ is the initial state and $V(t,t_{0})$ denotes the time-evolution operator corresponding to the equations of motion (1a) and (1b). The coefficients of these differential equations do not explicitly depend on time, which implies $V(t+t_{0},t_{0})=V(t,0)$ for all $t$ and $t_{0}$. Under this evolution and in the limit $t\to\infty$, we expect $\rho$ to evolve towards a steady state, which we denote by

[TABLE]

This allows us to write the long-time limit of the two-time correlation function as

[TABLE]

The nontrivial object in this expression is the time-evolved operator $V(\tau,0)\left(\rho_{\text{ss}}\sigma_{i}^{a}\right)$, which we are going to calculate in Sec. III.3.

III.2 Time-evolution equations

As a starting point for calculating $V(\tau,0)\left(\rho_{\text{ss}}\sigma_{i}^{a}\right)$ we need the steady state density operator $\rho_{\text{ss}}$, which we calculate approximately by a method described in the following. This method is applicable to arbitrary $N$-spin trace-1 operators $\mathscr{A}_{1\ldots N}$ under the dynamics generated by a Lindblad operator. Starting from the initial density operator, $\mathscr{A}_{1\ldots N}=\rho$, will allow us to obtain $\rho_{\text{ss}}$ after sufficiently long evolution times. Other choices of $\mathscr{A}_{1\ldots N}$ will be used to calculate the time-evolution of other trace-1 operators that appear when calculating $V(\tau,0)\left(\rho_{\text{ss}}\sigma_{i}^{a}\right)$ in Sec. III.3.

The propagator $V$ that induces the evolution equations (1a) and (1b) can be written as a Lindblad differential equation

[TABLE]

where the Lindblad operator

[TABLE]

consists of on-site terms

[TABLE]

and pair interactions

[TABLE]

with

[TABLE]

and $\Gamma=d^{2}k_{0}^{3}/(2\pi\hbar\epsilon_{0})$.

Taking partial traces on both sides of (11) one obtains, in the spirit of the Bogoliubov-Born-Green-Kirkwood-Yvon (BBGKY) hierarchy Bonitz (2015), a set of coupled evolution equations,

[TABLE]

where the reduced operators

[TABLE]

are defined as partial traces over all sites except the indexed ones. It is expected for the dilute regime that the two-atom coupling dominates over three- or more-atom coupling, which makes a cluster expansion suitable for the truncation of the hierarchy of equations. By means of a cluster expansion we separate the reduced $\mathscr{A}$ operators into product and connected parts,

[TABLE]

which implicitly defines the connected operators $\mathscr{C}$. Substituting these definitions into (16a) and (16b), we rewrite the first two equations of the BBGKY hierarchy as

[TABLE]

where $\mathcal{L}^{S}_{ij}=\mathcal{L}_{ij}+\mathcal{L}_{ji}$. Equation (19b) contains three-spin connected contributions $\mathscr{C}_{ijm}$, the time-evolution of which depends on four-spin terms, and so on. To turn this into a numerically tractable problem, we truncate the BBGKY hierarchy at second order by neglecting the terms $\mathscr{C}_{ijm}$ in (19b). As stated in Pucci et al. (2016), if we also neglected the $\mathscr{C}_{ij}$ terms, we would recover the classical time-evolution equation presented in Schachenmayer et al. (2015), in that sense, more terms of the truncated hierarchy means more quantum corrections for the spins dynamics. The two-spin connected contribution $\mathscr{C}_{ij}$ are related to the two-spin quantum correlations, which implies that, unlike in the classical case, the two-spin connected correlations do not vanish, $\langle\sigma_{j}^{\pm,z}\sigma_{m}^{\pm,z}\rangle-\langle\sigma_{j}^{\pm,z}\rangle\langle\sigma_{m}^{\pm,z}\rangle\neq 0$. This is the main approximation made in the numerical scheme, and it gives good results whenever genuine three- and more-spin connected contributions are negligible.

To bring the resulting truncated operator equations into a numerically tractable form, we expand all operators in terms of Pauli spin operators and we obtain a set of coupled ordinary differential equations (see Appendix A), which can be integrated by standard numerical methods.

III.3 Calculation of $V(\tau,0)\left(\rho_{\text{ss}}\sigma_{i}^{a}\right)$

We calculate the steady-state density operator $\rho_{\text{ss}}$ by setting to zero the left-hand sides of (31) and (32), and numerically solving the resulting algebraic equations by a standard Newton-Krylov solver. The stationary values of the a- and c-coefficients in those equations encode the required information on $\rho_{\text{ss}}$.

Starting from the thus obtained steady-state density operator $\rho_{\text{ss}}$, we expand $\rho_{\text{ss}}\sigma_{i}^{a}$ in terms of so-called phase point operators, rewrite the result in terms of trace-1 operators, and then use again the time-evolution equations of Secs. B and C in order to obtain $V(\tau,0)\left(\rho_{\text{ss}}\sigma_{i}^{a}\right)$.

The discrete phase-space representation of a single spin-$1/2$ degree of freedom as introduced by Wootters Wootters (1987) is based on a discrete phase space

[TABLE]

consisting of four points, each of which has an associated three-vector, $\bm{r}_{(0,0)}=(1,1,1)$, $\bm{r}_{(0,1)}=(-1,-1,1)$, $\bm{r}_{(1,0)}=(1,-1,-1)$, and $\bm{r}_{(1,1)}=(-1,1,-1)$. To each phase space point $\alpha\in\Gamma$ one assigns a so-called phase point operator

[TABLE]

where $\bm{\sigma}=\left(\sigma^{x},\sigma^{y},\sigma^{z}\right)$ is the vector of Pauli operators. The phase point operators form a basis, and any operator on ${\mathds{C}}^{2}$ can be expressed as a linear combination of the four operators $A_{\alpha}$. Similarly one could expand an operator on the tensor product Hilbert space $({\mathds{C}}^{2})^{\otimes N}$ of $N$ spin-$1/2$ degrees of freedom in the corresponding tensor product basis of $A_{\alpha}$. Here we follow a different approach and do this expansion only for the $i$th factor of the product space, which yields 111It turns out to be advantageous to expand in an overcomplete basis, using the phase point operators corresponding to $\bm{r}_{(0,0)}^{\prime}=(1,-1,1)$, $\bm{r}_{(0,1)}^{\prime}=(-1,1,1)$, $\bm{r}_{(1,0)}^{\prime}=(1,1,-1)$, $\bm{r}_{(1,1)}^{\prime}=(-1,-1,-1)$, in addition to those defined above. While such an expansion gives identical results on an exact level, differences arise when approximating the time-evolution. Avoiding the expansion altogether, and writing $\rho_{\text{ss}}\sigma_{i}^{a}$ directly in terms of trace-$1$ operators, is also feasible, but again turned out to perform worse than the scheme described here.

[TABLE]

where $\operatorname{Tr}_{i}$ denotes a partial trace over the $i$th factor of the tensor product Hilbert space. By elementary spin algebra, the coefficients of this expansion, which are operator-valued and act on $({\mathds{C}}^{2})^{\otimes(N-1)}$, can be written as

[TABLE]

Next we rewrite (23) as a linear combination of trace-$1$ operators, such that the time-evolution scheme of Sec. III.2 can be applied to each of those operators. To this purpose it is convenient to (partially) expand operators in the tensor product basis of Pauli spin operators, where we denote the expansion coefficients by

[TABLE]

Starting from (23) we can write 222We tested other ways of expressing (23) in terms of trace-1 operators, but for our purposes none of them turned out to be advantageous in terms of accuracy or computational cost.

[TABLE]

where we have defined

[TABLE]

both of which are trace-$1$ operators on $({\mathds{C}}^{2})^{\otimes(N-1)}$. Those operators can also be expanded in terms of Pauli spin operators, and the corresponding expansion coefficients can be expressed in terms of the coefficients (24) of $\rho_{\text{ss}}$,

[TABLE]

and so on. Inserting (22) and (25) into (9) we obtain

[TABLE]

with

[TABLE]

According to (29a) and (29b), in order to obtain the desired two-time correlation function (28), we need to calculate for each $i$ and each phase space operator $A_{\alpha_{i}}$ the time-evolution of four operators, namely $\tilde{\rho}^{a}_{ss,\not{\,i}}A_{\alpha_{i}}$ for $a\in\{x,y,z\}$, and $\rho_{ss,\not{\,i}}A_{\alpha_{i}}$. We do this by making use of the method developed in Sec. III.2, letting $\mathscr{A}_{1\ldots N}$ in (10) take the role of each of the four mentioned operators. The computational cost of the method scales like $3N(3N-1)/2$ with the system size $N$. Applying it to the four trace-$1$ operators, the eight phase point operators, and the $N$ lattice sites required in (28), (29a), and (29b), results in an overall computational cost that scales asymptotically like $N^{3}$.

IV Results

We performed simulations for system sizes $N=14$, $28$, $48$, $72$ and $96$ atoms at a fixed low density $\rho$. The atoms are placed at random in a spherical volume, but with the constraint of a minimal distance of one fourth of the mean distance between neighbours, such as to rule out unwanted noncooperative effects due to accidentally close pairs of atoms. Fig. 1(a) shows the numerically computed fluorescence spectra for various $N$, at fixed density and laser parameters. The three prominent peaks in the plot at $\omega=0$ and $\pm\Omega_{0}$ form the Mollow triplet Mollow (1969); Schuda et al. (1974). The first main result is the observation of additional sidebands in the fluorescence spectrum at $\omega=\pm 2\Omega_{0}$. These sidebands are genuine quantum effects, as they require the presence of quantum pair correlations. Indeed if connected correlations between different sites were absent and two-time correlations would factorize, $\langle\sigma_{j}^{\pm,z}\sigma_{m}^{\pm,z}\rangle=\langle\sigma_{j}^{\pm,z}\rangle\langle\sigma_{m}^{\pm,z}\rangle$ for $j\neq m$, one would have $\langle\sigma_{j}^{+}(t)\sigma_{m}^{-}(t+\tau)\rangle=\langle\sigma_{j}^{+}(t)\rangle\langle\sigma_{m}^{-}(t)\rangle$ in the steady-state regime, and in that case the inelastic ($\omega\neq 0$) spectrum would only depend on single-site two-time correlations $\langle\sigma_{j}^{+}(t)\sigma_{j}^{-}(t+\tau)\rangle$. The factorizing terms $\langle\sigma_{j}^{\pm,z}\rangle\langle\sigma_{m}^{\pm,z}\rangle$ may modify the local Rabi frequency experienced by each atom and inhomogeneously broaden the single-atom Mollow triplet, but cannot give rise to higher-order sidebands. In other words, a classical treatment of the Hamiltonian (or, as in this case, Lindbladian) results in the absence of connected correlations between different spins, and no additional sidebands are observed.

The novel sidebands are true cooperative effects. If the sidebands were two-atom or few-atom effects, their peak height would depend only on the spatial density. In Figs. 1(a) and (b), however, we observe that the sidebands grow with the number of atoms $N$ even at fixed density, and scale linearly with the optical thickness $b_{0}=2N/(kR)^{2}$ rather than with the spatial density. This effect can be attributed to the long-range nature of the effective interactions (2) between the atomic internal degrees of freedom. The scaling with $b_{0}$ is reminiscent of cooperative phenomena in the linear optics regime Guerin et al. (2016b), but is here observed, for the first time to the best of our knowledge, for a quantum cooperative phenomenon in free space. Furthermore, although single scattering processes may exhibit quantum optics interferences phenomena Grangier et al. (1986), they cannot capture the additional sidebands. These sidebands at $\pm 2\Omega_{0}$ can be understood as the first step of a higher order harmonic generation process, where the next orders could be studied by including higher-order quantum correlations. However, the low relative intensity of approximately $10^{−3}$ makes it hard to detect the peaks at $\pm 2\Omega_{0}$ experimentally. With that in mind, we searched for traces of quantum cooperativity in the Mollow triplet bands ($\pm\Omega_{0}$).

The second result is the observation that quantum cooperativity breaks the symmetry of the spectrum. The single-atom fluorescence spectrum is always symmetric with respect to the frequency of the driving light, independently of the detuning of the driving from the atomic resonance Mollow (1969). For large atomic clouds and in the presence of detuning ($\Delta\neq 0$), it was predicted that coherence effects may induce an asymmetry of the Mollow sidebands in the forward scattering direction Ott et al. (2013). Our simulations show a similar effect for the scattering of detuned light, where the Mollow sidebands at $\omega=\pm\Omega$ exhibit a significant asymmetry (Fig. 2 left). This asymmetry scales with the optical thickness $b_{0}$ (Fig. 2 right), which confirms the cooperative nature of this effect. In the absence of quantum correlations the spectrum, being composed of the sum of $N$ symmetric spectra, is necessarily symmetric, which confirms the genuine quantumness of the observed asymmetry. However, going beyond the prediction of Ref. Ott et al. (2013), we here observe that the asymmetry is also present outside the forward lobe, i.e., for scattering angles $\theta\geq 1/kR$ where diffuse light dominates (Fig. 2 left). Surprisingly, the asymmetry is inverted in the forward direction ($\theta<1/kR$) in comparison with $\theta>1/kR$. Experimentally the asymmetry of the standard Mollow sidebands, which reaches $\sim 30\%$ in our simulations, should be relatively easy to detect.

In Fig. 3 we show the fluoresence spectrum as a function of the scattering angle $\theta$ and the frequency $\omega$ in the regime of deep saturation, where most of the light is expected to be scattered inelastically ($\omega\neq 0$). In this regime the portion of elastically scattered light for a single atom goes as $1/s=(\Delta_{0}^{2}+\Gamma^{2}/4)/2\Omega_{0}^{2}$ at large $\Omega_{0}$, $s$ being referred as the saturation parameter, so most of the light is scattered inelastically. We clearly observe the quasi-isotropic inelastic Mollow triplet and higher order sidebands. For the parameters considered, a strong elastic component is particularly visible in the forward direction (see inset), which we attribute to the constructive interference of the elastic component of the electric field in the forward direction (which, according to linear optics, is expected to scale like $N^{2}$ 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 that in the forward direction the inelastic component exhibits a small dip. The physical origin of this feature remains to be understood.

V Conclusions

We have reported the discovery of signatures of quantum cooperativity in the fluorescence spectrum of large dilute atomic clouds. The rise of additional sidebands at frequencies $\pm 2\Omega_{0}$ from the central line, as well as the asymmetry in the spectrum of the cloud driven out of resonance, are identified as proper quantum effects that cannot occur in the absence of genuine quantum correlations. Moreover, by analyzing parameter dependences and scaling properties, the cooperative nature of the observed phenomena is revealed. Cooperativity is ultimately related to the long-range character of the effective atom–atom interactions induced by the photon field. The deficiency in inelastically scattered photons in the forward cone ($\lvert\theta\rvert\leq 1/kR$) is particularly interesting, because it implies that the forward direction, which has long been considered a natural candidate for probing cooperative phenomena in the linear-optics regime Scully et al. (2006), may be less suitable for probing quantum cooperative effects. Furthermore, while the second-order optical coherence is usually considered the ideal candidate for revealing the quantum nature of the light scattering by atoms, with photon bunching Morgan and Mandel (1966) and anti-bunching Kimble et al. (1977) as paradigmatic signatures, we show in this paper that the first-order optical coherence $g^{(1)}$, which witnesses the quantum nature of the atom-light coupling, already shows clear signatures of quantum cooperativity. More generally, our results suggest that the quantum optics regime of an optically deep system is substantially richer than its single-atom physics, and holds much promise for further studies of cooperative effects. This may become relevant for neutral atom optical clocks or other long-range quantum systems such as Rydberg atoms. To gain access to this regime on the computational side, the simulation technique developed in this paper, based on a truncation of the hierarchy of correlations, proves to be a powerful tool.

The cooperative nature of the observed effects suggests that dilute atomic clouds might be used as experimental platforms for quantum-simulating plasmas, free electron lasers, and other quantum long-range interacting systems in which cooperativity plays an essential role.

Acknowledgements.

We thank N. Piovella for stimulating discussions. A.R. acknowledges financial support from the National Research Foundation of South Africa. This work has received financial support from FAPESP and CNPq. M.K. acknowledges financial support from the National Research Foundation of South Africa via the Incentive Funding and the Competitive Programme for Rated Researchers. M.K. and R.B. acknowledge financial support by a FAPESP/SU collaboration grant.

Appendix A Time-evolution of the Pauli expansion coefficients

We expand all operators in (18a) and (18b), except $\mathscr{C}_{ijm}$ that is neglected, in terms of Pauli spin operators,

[TABLE]

Inserting these expansions into (19a) and the truncated version of (19b), and making use of Lindblad equations (10)–(13), we obtain time-evolution equations for the Pauli expansion coefficients,

[TABLE]

and

[TABLE]

These equations form a set of coupled ordinary differential equations, which can be integrated by standard numerical methods.

We calculate the steady-state density operator $\rho_{\text{ss}}$ by setting to zero the left-hand sides of (31) and (32), and numerically solving the resulting algebraic equations by a standard Newton-Krylov solver. The stationary values of the $a$- and $c$-coefficients encode the required information on $\rho_{\text{ss}}$.

Appendix B Benchmarking against exact results

The accuracy of the proposed simulation method is tested by benchmarking the fluorescence spectrum $S$ against exact results. Up to $N=7$ spins (atoms) could be dealt with exactly by using the “Quantum Toolbox in Python” Johansson et al. (2012, 2013), a module tailored to simulate the dynamics of open quantum systems and especially those of quantum optics.

As shown in Fig. 4 for Rabi frequency $\Omega_{0}=20\Gamma$, for densities up to $\rho/k^{3}=0.1$ our simulation results are in very good agreement with exact results, for the main as well as the secondary Mollow sidebands. For larger densities ($\rho/k^{3}=0.3$ in that same figure), when the coupling gets stronger and higher order correlations are expected to become more relevant, the two spectra exhibit more substantial deviations, although the agreement is still acceptable. Other values of the driving frequency $\Omega_{0}$ lead to a similar degree of agreement (not shown).

Besides the spectra, we also benchmarked other relevant quantities, including the steady state $\rho_{\text{ss}}$ calculated according to Sec. C, as well as the the two-time correlations evolved from the latter. All show very good agreement for densities up to $\rho/k^{3}\sim 0.3$.