Stochastic electrodynamics simulations for collective atom response in optical cavities

TL;DR

This paper uses stochastic electrodynamics simulations to explore the collective optical behavior of atoms in a cavity, revealing how subradiant modes and atomic correlations influence the cavity response and quantum effects.

Contribution

It introduces a semiclassical approximation for atomic saturation in stochastic electrodynamics simulations, accurately capturing quantum effects near subradiant resonances.

Findings

Subradiant modes can be coupled and resolved with spatially varying atomic transition frequencies.

Cavity response reveals atomic position correlations and quantum many-body phases.

Semiclassical models qualitatively match quantum responses outside subradiant resonances.

Abstract

We study the collective optical response of an atomic ensemble confined within a single-mode optical cavity by stochastic electrodynamics simulations that include the effects of atomic position correlations, internal level structure, and spatial variations in cavity coupling strength and atom density. In the limit of low light intensity the simulations exactly reproduce the full quantum field-theoretical description for cold stationary atoms and at higher light intensities we introduce semiclassical approximations to atomic saturation that we compare with the exact solution in the case of two atoms. We find that collective subradiant modes of the atoms, with very narrow linewidths, can be coupled to the cavity field by spatial variation of the atomic transition frequency and resolved at low intensities, and show that they can be specifically driven by tailored transverse pumping beams. We show that the cavity optical response, in particular both the subradiant mode profile and the resonance shift of the cavity mode, can be used as a diagnostic tool for the position correlations of the atoms and hence the atomic quantum many-body phase. The quantum effects are found to be most prominent close to the narrow subradiant mode resonances at high light intensities. Although an optical cavity can generally strongly enhance quantum fluctuations via light confinement, we show that the semiclassical approximation to the stochastic electrodynamics model provides at least a qualitative agreement with the exact optical response outside the subradiant mode resonances even in the presence of significant saturation of the atoms.

Figures

Equations

$${\mathcal{H}}_{int}=-\frac{1}{{ϵ}_0}\hat{\mathbf{P}}(x)\cdot \hat{\mathbf{D}}(x).$$

$${\hat{\mathbf{D}}}^{+}(x)=\sum _q{\zeta }_q{\hat{\mathbf{e}}}_q{\hat{a}}_q(t)f_q(x),{\zeta }_q=\sqrt{\frac{\hslash {ϵ}_0{\omega }_q}{2L}},$$

$${\mathcal{H}}_{int}\approx -\hslash \left[\mathbf{g}(x)\cdot {\hat{\mathbf{P}}}^{-}(x)\hat{a}(t)+{\mathbf{g}}^{\ast }(x)\cdot {\hat{\mathbf{P}}}^{+}(x){\hat{a}}^{†}(t)\right].$$

$$H$$

$$-\hslash \int \mathrm{d}x\mathbf{h}(x)\cdot {\hat{\mathbf{P}}}^{-}(x)+{\mathbf{h}}^{\ast }(x)\cdot {\hat{\mathbf{P}}}^{+}(x)$$

$$-\hslash \int \mathrm{d}x\left[\mathbf{g}(x)\cdot {\hat{\mathbf{P}}}^{-}(x)\hat{a}+{\mathbf{g}}^{\ast }(x)\cdot {\hat{\mathbf{P}}}^{+}(x){\hat{a}}^{†}\right]$$

$$-\hslash {\Delta }_c{\hat{a}}^{†}\hat{a}-\hslash \eta \left(\hat{a}+{\hat{a}}^{†}\right),$$

$${\mathcal{H}}_g$$

$${\mathcal{H}}_e$$

$${\hat{\mathbf{P}}}^{+}(x)$$

$${\hat{\mathbf{P}}}_{\nu \eta }^{+}(x)$$

$${\mathbf{d}}_{g\nu e\eta }\equiv \mathcal{D}\sum _{\sigma }{\hat{\mathbf{e}}}_{\sigma }⟨e\eta ;1g|1\sigma ;g\nu ⟩\equiv \mathcal{D}\sum _{\sigma }{\hat{\mathbf{e}}}_{\sigma }{\mathcal{C}}_{\nu ,\eta }^{(\sigma )}.$$

$$\frac{d\hat{\varrho }}{dt}=\frac{1}{i\hslash }\left[H,\hat{\varrho }\right]+\kappa \left(2\hat{a}\hat{\varrho }{\hat{a}}^{†}-{\hat{a}}^{†}\hat{a}\hat{\varrho }-\hat{\varrho }{\hat{a}}^{†}\hat{a}\right).$$

$$\hat{a}={\hat{a}}_F+\frac{i}{\bar{\kappa }}\int \mathrm{d}x{\hat{\mathbf{P}}}^{+}(x)\cdot {\mathbf{g}}^{\ast }(x),$$

$$\frac{1}{\bar{\kappa }}=\frac{\kappa +i{\Delta }_c}{{\kappa }^2+{\Delta }_c^2},$$

$${\dot{\hat{\psi }}}_{g\nu }(x,t)$$

$$+\frac{1}{{\bar{\kappa }}^{\ast }}\int \mathrm{d}{x}^{\prime }{\mathbf{d}}_{g\nu e\eta }\cdot {\mathsf{G}}_c^{\ast }(x,x^{\prime }){\hat{\mathbf{P}}}^{-}(x^{\prime }){\hat{\psi }}_{e\eta }(x)$$

$$+i{\mathbf{g}}^{\ast }(x)\cdot {\mathbf{d}}_{g\nu e\eta }{\hat{a}}_F^{†}{\hat{\psi }}_{e\eta }(x),$$

$${\dot{\hat{\psi }}}_{e\eta }(x,t)$$

$$-\frac{1}{\bar{\kappa }}\int \mathrm{d}{x}^{\prime }{\mathbf{d}}_{e\eta g\nu }\cdot {\mathsf{G}}_c(x,x^{\prime }){\hat{\mathbf{P}}}^{+}(x^{\prime }){\hat{\psi }}_{g\nu }(x)$$

$$+i\mathbf{g}(x)\cdot {\mathbf{d}}_{e\eta g\nu }{\hat{\psi }}_{g\nu }(x){\hat{a}}_F$$

$$-\frac{1}{\bar{\kappa }}{\mathbf{d}}_{e\eta g\tau }\cdot {\mathsf{G}}_c(x,x){\mathbf{d}}_{g\tau e\zeta }{\hat{\psi }}_{e\zeta },$$

$${\Delta }_{e\eta }=$$

$${\Delta }_{g\nu }=$$

$${\mathsf{G}}_c(x,x^{\prime })=\mathbf{g}(x){\mathbf{g}}^{\ast }(x^{\prime }).$$

$$\left[{\hat{\psi }}_{g\nu }(x),{\hat{a}}_F\right]=-\frac{i}{\bar{\kappa }}{\mathbf{g}}^{\ast }(x)\cdot {\mathbf{d}}_{g\nu e\eta }{\hat{\psi }}_{e\eta }(x),$$

$$\frac{d}{dt}{\hat{P}}^{+}(x)$$

$$+\frac{{\mathcal{D}}^2}{\bar{\kappa }}\int \mathrm{d}{x}^{\prime }{G}_c(x,x^{\prime }){\hat{\psi }}_e^{†}(x){\hat{P}}^{+}(x^{\prime }){\hat{\psi }}_e(x)$$

$$-\frac{{\mathcal{D}}^2}{\bar{\kappa }}\int \mathrm{d}{x}^{\prime }{G}_c(x,x^{\prime }){\hat{\psi }}_g^{†}(x){\hat{P}}^{+}(x^{\prime }){\hat{\psi }}_g(x)$$

$$-i{\mathcal{D}}^2\left[{\hat{\psi }}_e^{†}{\hat{\psi }}_e-{\hat{\psi }}_g^{†}{\hat{\psi }}_g\right]\left[g(x){\hat{a}}_F+h(x)\right],$$

$$\frac{d}{dt}{\hat{\psi }}_e^{†}{\hat{\psi }}_e$$

$$=-2{\mathcal{D}}^2\text{mbox}Re\left[\frac{1}{\bar{\kappa }}\right]G_c(x,x){\hat{\psi }}_e^{†}{\hat{\psi }}_e$$

$$-\frac{\mathcal{D}}{\bar{\kappa }}\int \mathrm{d}{x}^{\prime }{G}_c(x,x^{\prime }){\hat{\psi }}_e^{†}(x){\hat{P}}^{+}(x^{\prime }){\hat{\psi }}_g(x)$$

$$-\frac{\mathcal{D}}{{\bar{\kappa }}^{\ast }}\int \mathrm{d}{x}^{\prime }{G}_c^{\ast }(x,x^{\prime }){\hat{\psi }}_g^{†}(x){\hat{P}}^{-}(x^{\prime }){\hat{\psi }}_e(x)$$

$$+ig(x){\hat{P}}^{-}{\hat{a}}_F-i{g}^{\ast }(x){\hat{a}}_F^{†}{\hat{P}}^{+}$$

$$+i\left[h(x){\hat{P}}^{-}-h^{\ast }(x){\hat{P}}^{+}\right],$$

$$\frac{d}{dt}{\hat{P}}^{+}(x)$$

$$-\frac{{\mathcal{D}}^2}{\bar{\kappa }}\int \mathrm{d}{x}^{\prime }{G}_c(x,x^{\prime }){\hat{\psi }}_g^{†}(x){\hat{P}}^{+}(x^{\prime }){\hat{\psi }}_g(x)$$

$$+i{\mathcal{D}}^2{\hat{\psi }}_g^{†}{\hat{\psi }}_g\left[g(x){\hat{a}}_F+h(x)\right].$$

$$P_2(x;x^{\prime })\equiv ⟨{\hat{\psi }}_g^{†}(x){\hat{P}}^{+}(x^{\prime }){\hat{\psi }}_g(x)⟩.$$

$$\frac{d}{dt}$$

$$-\frac{{\mathcal{D}}^2}{\bar{\kappa }}{G}_c(x^{\prime },x)P_2(x^{\prime };x)$$

Full text

I Introduction

The advent of optical cavities revolutionized the study of quantum interactions of light and atoms by enabling the study of a single quantized light field with matter Carmichael (2007); Ritsch et al. (2013). More recently, advances in trapping cold atoms Kruse et al. (2003); Gupta et al. (2007) in single-mode high-finesse optical cavities have greatly enhanced the control and tunability of such systems. Current experiments are now able to work with quantum degenerate atoms Brennecke et al. (2007); Slama et al. (2007); Colombe et al. (2007); Murch et al. (2008); Brahms et al. (2012); Botter et al. (2013); Schmidt et al. (2014); Kessler et al. (2014) coupling the study of quantum optics with quantum many-body physics.

Confining many atoms and light inside a cavity can lead to a response that is very different from that of a single atom. The typical signature of such collective behavior occurs when the emission rates of collective modes are enhanced (superradiant) or suppressed (subradiant) compared to the emission from a single atom, and cannot be accounted for by a model of atoms which scatter light independently. Within an optical cavity these collective effects can lead to phenomena such as superradiant lasers Bohnet et al. (2012) or realizations of the Dicke model Baumann et al. (2010). In multimode cavities, approximate descriptions of many-atom systems have been shown to result, e.g., in a rich phenomenology of different quantum phases Gopalakrishnan et al. (2009); Strack and Sachdev (2011). The wide variety of physics accessible by cavity-atom systems has necessitated a number of different treatments for the studies of many-particle interactions with the cavity field Moore et al. (1999); Horak and Ritsch (2001); Gardiner et al. (2001); Nagy et al. (2008); Szirmai et al. (2009); Venkatesh and O’Dell (2013); Mekhov and Ritsch (2009); Niedenzu et al. (2013); Lee and Ruostekoski (2014, 2015); Kuraptsev and Sokolov (2016); Goldwin et al. (2014); Samoylova et al. (2015); Elliott et al. (2015a, b); Bhattacherjee (2015); Joshi and Larson (2015); Mazzucchi et al. (2016); Elliott and Mekhov (2016), including studies of cavity optomechanics Kippenberg and Vahala (2007, 2008); Meystre (2013); Aspelmeyer et al. (2013). Cold atomic ensembles trapped inside or close to the waveguides and fibers form a closely related system where the light field mediating interactions between atoms is confined in one dimension Goban et al. (2015); Vetsch et al. (2010); Goban et al. (2012); Yalla et al. (2014); Chang et al. (2013); Shen and Fan (2005); Ruostekoski and Javanainen (2016).

Strong collective optical response of cold and dense atomic gases in free space has attracted considerable experimental interest Balik et al. (2013); Bienaimé et al. (2010); Pellegrino et al. (2014); Kwong et al. (2014); Jenkins et al. (2016a); Jennewein et al. (2016); Bromley et al. (2016); Guerin et al. (2016); Bons et al. (2016). The atoms can exhibit collective radiative resonance linewidths and line shifts, and experience recurrent scattering where light is scattered more than once by the same atom – related effects of which have been actively investigated theoretically in cold atom vapors Morice et al. (1995); Ruostekoski and Javanainen (1997a, b); Javanainen et al. (1999); Ruostekoski and Javanainen (1999); Berhane and Kennedy (2000); Clemens et al. (2003); Chomaz et al. (2012); Jenkins and Ruostekoski (2012); Antezza and Castin (2013); Javanainen et al. (2014); Javanainen and Ruostekoski (2016); Kuraptsev and Sokolov (2014); Bettles et al. (2015); Lee et al. (2016); Schilder et al. (2016); Longo et al. (2016); Jones et al. (2017); Zhu et al. (2016); Sutherland and Robicheaux (2016); Jenkins et al. (2016b); Bettles et al. (2016); Sutherland and Robicheaux (2017). In free space the light-induced resonant dipole-dipole interactions are of finite range, resulting in light-mediated interactions that depend on the positions of the individual atoms in the ensemble. The combination of recurrent scattering and the position-dependent light-mediated interactions can lead to a correlated atom response that dramatically differs from that predicted by mean-field theoretical models of continuous medium electrodynamics in which such effects are ignored Javanainen et al. (2014); Javanainen and Ruostekoski (2016).

In a cavity the optical response can be quite different from that of free space. Even a single atom can experience recurrent scattering by repeatedly absorbing and emitting the same photon. Quantum effects can more clearly manifest themselves owing to the directed light confinement. In contrast to atoms in free space, the sole common electromagnetic mode of a single-mode optical cavity allows long-range scattering of light between atoms without attenuation. As a consequence, collective radiative effects can be important in a cavity even if the interatomic spacing is larger than a wavelength. However, the collective phenomena are simplified due to the absence of attenuation in the scattered light, considerably suppressing the spatial dependence of the light-mediated interactions. Consequently, in some limits, exact solutions may be obtained in cavities, for example when all atoms experience an identical cavity field and trapping potential Tavis and Cummings (1968). However, modern experiments can involve atomic distributions which are broad compared to the cavity wavelength. As a result, spatial variations in coupling strengths, densities and detunings may become important Carmichael (2007), leading to position-dependent correlations between individual atoms.

In this paper, we calculate the optical response of a many-atom cavity system where we include the position-dependent correlations between pointlike atoms that are generated by light-mediated interactions in the cavity. We develop a semiclassical stochastic electrodynamics simulation method that includes recurrent scattering processes between the atoms, atomic saturation, internal level structure, and the atom statistics. Crucially, the description can incorporate the effects of spatial variation of the atom distribution, the cavity mode, and the detuning, so that the atoms no longer are assumed to experience identical fields or potentials. We derive the semiclassical approach from a quantum field-theoretical analysis for the coupled theory of atoms and the cavity field, with an approach analogous to the one we have previously implemented for the interaction of light with an atomic ensemble in free space Lee et al. (2016). In the limit of low light intensity, the stochastic electrodynamics simulation technique exactly reproduce the full quantum field-theoretical description for cold stationary two-level atoms Javanainen et al. (1999); Lee et al. (2016). At higher light intensities where saturation is important, we implement the semiclassical approximation that neglects quantum fluctuations between the levels. For quantum degenerate atomic ensembles the quantum statistical position correlations between the atoms can be incorporated in the simulations of the optical response, provided that one can synthesize the corresponding stochastic ensemble of atomic positions.

Using the stochastic electrodynamics simulations, we discuss the eigenmode structure of the cavity-atom system, with particular emphasis on the collective atomic modes that can be excited. The most commonly studied optical response of the system can be considered to be due to a collective mode with superradiant characteristics - the cavity optical field induces this collective mode to decay much faster than would be the case for an isolated atom, known as cavity-enhanced spontaneous emission Carmichael (2007). We show that, together with this superradiant mode, there exist subradiant collective modes with extremely narrow linewidths, and that these subradiant modes can be coupled to the cavity light by introducing a spatial variation in atomic transition frequency. The subradiant modes can have distinct spatial profiles, allowing targeted excitation by suitably tailored transverse pump beams, and their long decay times also provide an opportunity for cavity light-storage mediated by collective atomic excitations, provided that the decay to other modes than the cavity mode is suppressed, e.g., by a regular arrangement of atoms.

Since our stochastic electrodynamics are semiclassical outside of the low intensity limit, we also compare such simulations against the exact quantum treatment obtained by solving the full coupled equations for the atomic correlation functions governing the optical response Ruostekoski and Javanainen (1997a) for a tractable two-atom system. Such comparisons not only inform about the validity of the numerical approximation, but also reveal quantum features in the two-atom response of the cavity system. The semiclassical stochastic simulations are able to capture the broad behavior of the optical response when it is dominated by the superradiant collective mode, although we find pronounced quantum features in the intensity resonance profile in the full quantum treatment, and these features are beyond the description of the semiclassical treatment. Quantum effects are even more significant when subradiant modes dominate the response, and the stochastic simulations do not give a quantitative agreement outside of the low intensity limit. For frequencies detuned from collective atomic mode resonances, the stochastic simulations agree well with the full quantum results. The existence of quantum features in the collective cavity-atom response indicates the tendency of an optical cavity to enhance the effect of quantum fluctuations via light confinement.

While the general techniques developed here can be used for any cold-atom system in a cavity, the simulations we present in this paper specialize to systems in which the atoms are additionally confined in a lattice potential along the cavity axis. A crucial advantage of our simulation technique is that the quantum statistical position correlations of the atoms are accounted for by means of the stochastic sampling procedure, a necessary feature to consider for cavities containing quantum degenerate atomic gases. We show that the optical response can act as a diagnostic for the atomic many-body phase by considering the different behavior of atoms in Mott insulator (MI) or superfluid states within the lattice potential. We show that both the subradiant spectrum and the resonance shift of the cavity mode are appreciably altered by the different position correlations of those states. By comparing the coherently and incoherently scattered light we also find how the difference between the total scattered light intensity and the coherently scattered light intensity is sensitive to the atom statistics, indicating how the fluctuations of the atomic positions are mapped onto the fluctuations of the scattered light.

We begin in Sec. II by deriving the full quantum field-theoretical formalism for atom-cavity system. Specializing in Sec. III to a system of two-level atoms, we discuss the stochastic method that allows an efficient numerical simulation. In Sec. IV we calculate the eigenmode structure of the cavity-atom system, before studying the two atom case and the comparison between stochastic simulations and the full quantum field-theoretical treatment in Sec. V. Finally, in Sec. VI we show how the different atomic many-body correlations in MI and superfluid states can easily be accounted for in the stochastic simulations, and their effects on the optical response of the system.

II Hamiltonian and equations of motion

II.1 System

We consider atoms confined within an optical cavity of length $L$. We assume tight confinement in the transverse dimensions such that the motion of the atoms is restricted to a single transverse mode – one which is sufficiently narrow that the coupling of atoms to the cavity depends only upon the longitudinal position along the cavity axis $x\in(0,L)$. The system then reduces to a one-dimensional (1D) problem for atomic quantities integrated over the transverse dimensions, such as the atomic polarization $\hat{{\boldsymbol{\mathbf{P}}}}(x)$. We include two pumping mechanisms, both at frequency $\Omega$, a direct cavity axis pump of strength $\eta$ and a beam illuminating the atoms in the cavity from a transverse direction with a strength, shape and polarization incorporated into ${\bf h}(x)$.

Within the dipole approximation, and expressing the Hamiltonian in the length gauge from the Power-Zienau-Woolley transformation Power and Zienau (1959); Woolley (1971); Cohen-Tannaudji et al. (1989), the interaction between the atoms and cavity light modes may be given by the 1D Hamiltonian density

[TABLE]

Here, $\hat{{\boldsymbol{\mathbf{D}}}}(x)$ is the electric displacement, which is the basic dynamical variable for light, and which can be separated into a sum of the positive and negative frequency components $\hat{{\boldsymbol{\mathbf{D}}}}^{+}(x)$ and $\hat{{\boldsymbol{\mathbf{D}}}}^{-}(x)$ respectively, with $\hat{{\boldsymbol{\mathbf{D}}}}^{+}(x)=[\hat{{\boldsymbol{\mathbf{D}}}}^{-}(x)]^{\dagger}$. The positive frequency component can be written in terms of the cavity mode photon annihilation operators $\hat{a}_{q}$

[TABLE]

where $\hat{\boldsymbol{\mathbf{e}}}_{q}$ is the mode polarization, and $f_{q}(x)$ encapsulates the spatial variation of the cavity mode. We specialize now, and in the remainder of this paper, to the case of a cavity with a single mode relevant mode, and for clarity drop the index $q$, however the treatment that we subsequently develop could, in principle, be extended to treat multimode cavities.

Defining ${\bf g}(x)=\zeta f(x)\hat{\boldsymbol{\mathbf{e}}}/\hbar\epsilon_{0}$, such that for a Fabry-Perot cavity ${\bf g}(x)={\bf g}_{0}\sin(qx)$, while for a ring cavity we may have modes ${\bf g}(x)={\bf g}_{0}e^{\pm iqx}$, then the Hamiltonian density of Eq. (1) becomes

[TABLE]

Here we have additionally made the rotating wave approximation. We note that, in principle, counter-rotating terms could be included in a treatment analogous to that below, and may be required to consider very strong coupling.

Factoring out the dominant driving frequency time dependence for all time-dependent operators, the full Hamiltonian for the system in the rotating frame of the pump may be written 111In length gauge of the Power-Zienau-Woolley transformation, there is in addition a term proportional to $\hat{{\boldsymbol{\mathbf{P}}}}(x)\cdot\hat{{\boldsymbol{\mathbf{P}}}}(x)$, however we neglect this term here since for non-overlapping point dipoles this it vanishes.

[TABLE]

where $\Delta_{c}=\Omega-\omega_{c}$ is the detuning of the cavity mode from the pump. Here the second line includes the transverse pumping of the atoms, and the last line describes the cavity mode and its direct pumping with strength $\eta$. In the first line $\mathcal{H}_{g}$ and $\mathcal{H}_{e}$ are Hamiltonian densities for atomic fields governing the ground and excited levels respectively. In many realistic experimental situations the internal sublevel structure may allow multiple transitions to be driven by the cavity mode. We include the sublevel structure by defining quantum field operators $\hat{\psi}_{g\nu}(x)$ and $\hat{\psi}_{e\eta}(x)$ for the sublevels involved in the transition $\left|g,\nu\right\rangle\rightarrow\left|e,\eta\right\rangle$. For the case of the linear Zeeman shift caused by a magnetic field of strength $B(x)$ the effective 1D atomic Hamiltonian densities are then given by

[TABLE]

where $g^{(g)}_{l}$ and $g^{(e)}_{l}$ are the Landé g-factors Bransden and Joachain (2003) for levels $g$ and $e$, and $\omega_{0}$ is the resonance frequency of the $\left|g\right\rangle\leftrightarrow\left|e\right\rangle$ transition in the absence of any magnetic field.

In this notation, the polarization operator $\hat{{\boldsymbol{\mathbf{P}}}}^{+}$ can be represented as a sum over contributions from the different possible transitions

[TABLE]

where $\boldsymbol{\mathbf{d}}_{g\nu e\eta}$ represents the dipole matrix element for the transition $|e,\eta\rangle\rightarrow|g,\nu\rangle$

[TABLE]

Here the sum is over the unit circular polarization vectors $\sigma=\pm 1,0$, and ${\cal C}_{\nu,\eta}^{(\sigma)}$ denote the Clebsch-Gordan coefficients of the corresponding optical transitions. The reduced dipole matrix element is represented by ${\cal D}$ (here chosen to be real) and $\boldsymbol{\mathbf{d}}_{e\eta g\nu}=\boldsymbol{\mathbf{d}}_{g\nu e\eta}^{*}$. The light fields with the polarizations $\hat{\boldsymbol{\mathbf{e}}}_{\pm}$ and $\hat{\boldsymbol{\mathbf{e}}}_{0}$ drive the transitions $|g,\nu\rangle\rightarrow|e,\nu\pm 1\rangle$ and $|g,\nu\rangle\rightarrow|e,\nu\rangle$, respectively, in such a way that only the terms $\sigma=\eta-\nu$ in Eq. (9) are nonvanishing.

The open nature of the cavity system is important to consider, and the effect of transmissive loss of photons from the cavity mode through the end mirrors at a rate $\kappa$ may be included by considering the evolution of the density matrix $\hat{\varrho}$

[TABLE]

The equation of motion for $\hat{\varrho}$ can be used to derive equations of motion for the cavity field $\hat{a}(t)$ also incorporating, e.g., the backaction of quantum measurement Carmichael (2007).

II.2 Equations of motion for quantum field operators

II.2.1 Light

We may now calculate equations of motion for the atomic and cavity mode operators. From the equation of motion for $\hat{a}$, and assuming that $\kappa$ is large compared with the other relevant frequency scales of the system, we adiabatically eliminate the cavity field. In the limit that fluctuations from the emitting atoms dominate, we obtain the dependence of $\hat{a}$ on the atomic polarization

[TABLE]

where

[TABLE]

and $\hat{a}_{F}$ represents the free field that would exist if there were no atoms present.

II.2.2 Matter

Substituting Eq. (11) into the equations of motion for the atomic operators, and assuming that the atomic motional state is unchanged by the scattering of light 222we therefore neglect any change to the atomic motional energy from a scattering event, an assumption equivalent to assuming stationary atoms throughout. See Ref. Lee et al. (2016) for a more detailed discussion of this approximation, leads to

[TABLE]

where repeated indices $\zeta,\tau$ are summed over. The atom-light detuning is denoted by the expressions

[TABLE]

where we have again factored out the dominant pump frequency $\Omega$. We have also defined the cavity propagator

[TABLE]

Due to the ability of the cavity to mediate long-range photon exchange between atoms this propagator depends only upon the positions $x$ and $x^{\prime}$, rather than the distance $|x-x^{\prime}|$ which would appear in free-space Lee et al. (2016). Additionally, we have reordered the free-field contribution $\hat{a}_{F}$ ($\hat{a}_{F}^{\dagger}$) to the right (left) hand side. We will assume the free field to be in a coherent state and this ordering means that the free-field operators will appear as multiplicative classical coherent fields $\langle\hat{a}_{F}\rangle=a_{F}=i\eta/\bar{\kappa}$ after expectation values are taken. The commutator

[TABLE]

which is required for the reordering is readily obtained from Eq. (11) by observing that atomic operators must commute with $\hat{a}$ at equal times.

II.3 Cavity-atom optical response for two-level atoms

II.3.1 The general two-level atom case

Finally, combining the above equations of motion for atom fields, and reordering the free-field contributions, results in the coupled equations for the polarizations and one-body correlation functions which together give the optical response of the system. We give the full expressions in App. A, but for clarity we now simplify to the case of atoms which are well described by a two-level model, coupling to a single mode optical cavity, a case which we use in the remainder of the paper unless otherwise specified. The two-level optical response is governed by the coupled equations

[TABLE]

where we have defined the two-level atom versions of Eqs. (7) and (17) by $\hat{P}^{+}={\cal D}\hat{\psi}^{\dagger}_{g}\hat{\psi}_{e}$ and $G_{c}(x,x^{\prime})=g(x)g^{*}(x^{\prime})$, with $g(x)={\bf g}(x)\cdot\boldsymbol{\mathbf{d}}_{eg}/{\cal D}$ and $h(x)={\bf h}(x)\cdot\boldsymbol{\mathbf{d}}_{eg}/{\cal D}$. The detuning $\bar{\Delta}=\Delta_{e}-\Delta_{g}$ now represents simply the detuning between the pump frequency and the atomic resonance frequency. In the presence of a spatially dependent applied magnetic field strength, $\bar{\Delta}$ can be spatially dependent [see Eqs. (5) and (6)], however for brevity we will not generally indicate the spatial dependence unless relevant.

In this two-level model, the origin of the terms in these equations may be discussed transparently. In Eq. (19), the last line represents the interaction of an atom at point $x$ with the transverse driving field $h(x)$ and the driven field of the cavity $\hat{a}_{F}$. The effect of saturation is accounted for by the factor of $(\hat{\psi}^{\dagger}_{e}\hat{\psi}_{e}-\hat{\psi}^{\dagger}_{g}\hat{\psi}_{g})$. The second and third lines include the effect of a dipole at point $x^{\prime}$ interacting with an atom at point $x$ via the cavity mode, an interaction governed by the propagator $G_{c}(x,x^{\prime})$. Finally, the term proportional to $G_{c}(x,x)$ in the first line arises from the reordering of the free-field contribution. In free-space, the similar procedure resulted in terms corresponding to both the spontaneous emission and the Lamb shift Lee et al. (2016). In the cavity, the term has the form of a self-interaction of a dipole at $x$ via the cavity field. However, the real and imaginary parts could alternatively be expressed as the position dependent cavity-enhanced emission rate ${\cal D}^{2}|g(x)|^{2}\kappa/(\kappa^{2}+\Delta_{c}^{2})$ and Lamb shift ${\cal D}^{2}|g(x)|^{2}\Delta_{c}/(\kappa^{2}+\Delta_{c}^{2})$ Brune et al. (1994), respectively.

II.3.2 The limit of low light intensity

The system simplifies greatly in the limit of low light intensity, where we may work to first order in the excited level operator $\hat{\psi}_{e\nu}$ or driving field operator $\hat{a}_{F}$. In this limit, the excited state density operator vanishes, the ground state density operator is invariant, and it remains only to treat the equation for the atomic polarization operator, which becomes

[TABLE]

Taking expectation values of both sides of this equation gives the equation of motion for the polarization $P(x)\equiv\langle\hat{P}^{+}(x)\rangle$, from which the optical response is readily calculated from Eq. (11). However, the equation of motion for $\langle\hat{P}^{+}(x)\rangle$ depends on the two-body correlation

[TABLE]

The equation of motion for $P_{2}$ may be derived in an analogous manner

[TABLE]

and so depends in turn upon a three-body correlation

[TABLE]

This continues up to the $N$th order – leading to a hierarchy of equations of motion for differing orders of correlations, which we give explicitly in App. B. Similar hierarchies, but with greater complexity, are encountered outside of the low intensity limit or in the case of multilevel atoms Ruostekoski and Javanainen (1997a).

In Eq. (25) the term proportional $P_{2}(x^{\prime};x)$ on the right hand side represents recurrent scattering processes where an excitation between the atoms at the positions $x$ and $x^{\prime}$ is repeatedly exchanged by the scattering of the same photon. In a high-finesse cavity even a single atom can undergo recurrent scattering by repeatedly absorbing and emitting the same photon. In cavity quantum electrodynamics this is referred to as the cooperative regime Carmichael (2007). In the case of interaction of light with atomic ensembles in free space, the cooperative terminology for scattering has traditionally had a notably broader meaning Guerin et al. (2017).

It rapidly becomes prohibitive to tackle the full hierarchy, and instead the problem can be approached by a perturbative expansion, for example in the small parameter $\rho_{3D}/k^{3}$, which is valid at low atom densities Morice et al. (1995); Ruostekoski and Javanainen (1999); van Tiggelen et al. (1990); Ruostekoski (1999). Rather than such a perturbative approach, the field-theoretical model can be mapped onto stochastic electrodynamics that is reminiscent of stochastic Langevin equation approach to a diffusion equation Javanainen et al. (1999); Lee et al. (2016). The formalism of Ref. Lee et al. (2016) was derived for free space light propagation – presenting techniques which exactly solve the free-space counterparts of Eqs. (40) for the case of atoms with a single ground level in the limit of low light intensity, and which provide approximate but numerically practical solutions outside of this limit. Analogous approaches can also be introduced for light propagation in confined 1D waveguides Ruostekoski and Javanainen (2016). Here we derive the stochastic electrodynamics for a cavity system in a semiclassical approximation. We consider explicitly only the case of two-level atoms, however the multi-level generalization follows from Eqs. (40) and the discussions in Ref. Lee et al. (2016).

III Stochastic electrodynamics simulations of two-level atoms

III.1 Stochastic sampling

Taking expectation values of Eqs. (19), and (21) leads to a hierarchy of coupled equations of motion for correlation functions describing the optical response of two-level atoms in a cavity. Computationally efficient solutions, particularly when spatial variation of the atom density, cavity mode and detunings are nontrivial, are provided by stochastic simulations of the optical response: A set of $N$ discrete atomic positions $\{X_{1},\ldots,X_{N}\}$ is sampled from the thermal equilibrium $N$-body joint probability distribution function for the atoms. For each such realization of atomic positions, we solve the optical response for the hypothetical model of $N$ atoms pinned to the fixed positions $\{X_{1},\ldots,X_{N}\}$. Subsequent ensemble averaging over many such independent realizations gives the quantum expectation values for the quantities governing the optical response of the ensemble.

The appropriate $N$-body distribution governing the position sampling is that representing the thermal equilibrium state of the atoms in the absence of light, prior to the introduction of any pump beams into the cavity. Since we have assumed that the atoms are stationary, this distribution function for the discrete position sampling is invariant. The subsequent introduction of light can induce a collective optical response, as governed by the equations of motion (19), and (21). These equations can be integrated, or solved in the steady-state, for each single realization of discrete positions. The optical response of the ensemble, including the effects of quantum statistical density correlations and light-induced atomic correlations, results from the subsequent averaging over many such realizations.

In principle, any density correlations can be included in the distribution function from which positions are drawn, for example Fermi-Dirac statistics can be modeled using a Metropolis algorithm Javanainen et al. (1999); Ruostekoski and Javanainen (2016). However many situations of interest can be sampled more straightforwardly – for an uncorrelated ensemble of classical atoms or an ideal Bose-Einstein condensate (BEC) the sampling reduces to independently sampling the position of atom $i$ from the ground state atomic density distribution of the ensemble in the absence of light. Since we here concentrate on stationary atoms whose center-of-mass position is constant and only the internal electronic state of the atoms evolves as a function of time, the total one-body density $\rho_{1}(x)$ is given by this is initial ground-state atom density at all times.

To carry out this procedure we first need to find the equations of motion for the atoms and light for each stochastic realization of atomic positions, which can be obtained from the expectation values of Eqs. (19), and (21) conditioned upon the set of positions $\{X_{1},\ldots,X_{N}\}$. The conditioned expectation values of one-body operators can be written in terms of one-body internal-level density matrix elements

[TABLE]

where $\rho_{ge}$ is the matrix element $\langle g|\hat{\rho}|e\rangle$ which is related to $\langle\hat{\psi}_{g}^{\dagger}\hat{\psi}_{e}\rangle$.

Eqs. (19), and (21) then lead to the coupled set of $2N$ nonlinear equations for a single realization of atomic positions

[TABLE]

Here we have used a semiclassical factorization approximation for all conditioned expectation values of two-body correlation functions, e.g.,

[TABLE]

Such an approximation means that the stochastic equations are not able to fully reproduce all the light-induced correlations in the system. However, the semiclassical approximation provides the means to include saturation effects in our stochastic electrodynamics while avoiding the need to solve the full hierarchy of equations of motion for quantum correlation functions discussed in the previous section. Furthermore, the semiclassical approximation does not neglect all correlations, since light-induced correlations that depend on the spatial positions of the atoms are included. In Sec. V we compare in detail the effect of the semiclassical approximation to the full quantum field-theoretical solution for the case of two atoms in a cavity, and show that the semiclassical approximation to the stochastic electrodynamics can still describe strong cooperative effects.

Ensemble averaging the single realization results from Eqs. (30) over many stochastic realizations of atom positions allows the calculation of physical observables. We could also extend the stochastic electrodynamics to include quantum effects beyond the semiclassical factorization of Eq. (31). This constitutes of including the correlation functions between the different internal atomic levels, as described in Ref. Lee et al. (2016) for multilevel atoms in free space. In the lowest order such correlation functions include the pair correlations between the internal levels of all atom pairs.

We note that in a cavity, the cavity field itself may be a significant contribution to the confining potential for the atoms, in addition to any external potential. In such a situation a self-consistent solution must be found, such that the cavity amplitude reflects the assumed $N$-body distribution function from which the atomic positions are sampled.

III.2 Limit of low light intensity

A limit which is both useful and important conceptually is that of low light intensity. In this limit, dropping terms of higher than first order in the light field amplitude or excited state operators leads to the simpler linear equation for the atomic ensemble in a single stochastic realization of atomic positions

[TABLE]

To first order in light field amplitude $\langle\hat{\psi}^{\dagger}_{e}\hat{\psi}_{e}\rangle$ vanishes, and $\langle\hat{\psi}^{\dagger}_{g}\hat{\psi}_{g}\rangle$ is invariant. Each stochastic realization in this limit is equivalent to solving the response of a model of linear harmonic oscillators at positions $X_{1},\ldots,X_{N}$ within an optical cavity. The importance of the limit of low light intensity is that, in contrast to the case where saturation is non-negligible, for atoms with a single relevant electronic ground level the stochastic solutions can be shown to fully reproduce the dynamics of the full hierarchy of correlation functions 333This can be shown by arguments that are analogous to those presented in App. B of Ref. Lee et al. (2016), starting from (42)., and so in this limit the accuracy of the stochastic method is limited only by the statistics of the number of realizations that are generated. In the following section we calculate the eigenmodes of the cavity-atom system, analyzing the limit of low light intensity in detail.

The stochastic electrodynamics equations (32) or (30) describe the optical response of an atomic ensemble in a cavity. The initial distribution of the atoms and their quantum statistics before the light enters the system can be sampled from the appropriate joint probability distribution of atomic positions, provided that such a distribution can be synthesized Javanainen et al. (1999). The ensemble-averaging the time evolution is then reminiscent of the stochastic Langevin equation solution to the diffusion equation. Many-atom cavity responses can also be simulated using stochastic phase-space methods based on the Wigner representation of the atomic and light fields Lee and Ruostekoski (2014, 2015). Here the initial quantum distribution of the atoms is also represented by stochastic sampling in the corresponding Langevin equation, obtained from the appropriate phase-space distribution. Each stochastic realization of the atomic wavefunction is then evolved according to the dynamical equations that also include a stochastic Wiener noise increment, resulting from the dynamics that is conditioned on the leaking of photons out of the cavity. However, in the phase-space approach the emphasis is on the center-of-mass motion of the atoms, and the atomic polarization is treated as a continuous field – an approximation that generally neglects the position-dependent light-induced correlations that we investigate in this paper.

IV Eigenmodes and optical response of the cavity system

IV.1 Eigenmodes in the limit of low light intensity

IV.1.1 Structure of the eigenmodes

In the limit of low light intensity, the equation of motion for the polarization density [Eq. (32)] for a single realization of discrete atomic positions can be written in the form

[TABLE]

where $F$ involves only the driving field terms. The matrix ${\cal M}$ therefore accounts for the cavity mediated interparticle interactions, along with the single particle cavity evolution and self-interaction of the atom with the cavity field. The eigenmodes of ${\cal M}$ are then instrumental in understanding the response of the atomic ensemble at a given cavity-pump detuning. In general, we find $N$ eigenmodes, however only one of the $N$ eigenvectors is certain to have a non-zero overlap with the cavity mode $g(x)$ and therefore to couple to the cavity mode. This mode involves all atoms collectively polarized with directions determined by the local sign of $g(x)$, and represents the strongest possible collective coupling of atoms to the cavity mode. It is characterized by a decay rate (linewidth) of the order of $Ng_{0}^{2}/\mbox{Re}[\bar{\kappa}]$ and frequency (line shift) $Ng_{0}^{2}/\mbox{Im}[\bar{\kappa}]$, it therefore decays much faster than would be expected for a single atom, and represents a superradiant mode.

In contrast, the other $N-1$ modes are subradiant, characterized by much narrower linewidths. In fact, in the absence of some non-periodic dependence on $g(x)$ generally the subradiant modes are completely orthogonal to the cavity mode and hence are infinitely long lived (trivially, since we neglect in our model any loss other than to the cavity mode). The exception to this statement occurs when the atomic detuning $\bar{\Delta}(x)$ is spatially varying. Such a detuning is generally introduced by a trapping potential and can be tailored in practice by employing a magnetic field strength gradient. Since ${\cal M}$ includes $\bar{\Delta}(x)$, any change in the detuning perturbs the eigenmodes of the system. If $\bar{\Delta}(x)$ is not uniform or of a form proportional to $g(x)$, then the perturbation of the eigenmodes is not merely a trivial frequency shift. Instead it results in the coupling of subradiant modes to the cavity mode and therefore in non-zero decay rates for the subradiant modes. Consequently, in subradiant mode results presented in this paper we necessarily re-diagonalize ${\cal M}$ for each different $\bar{\Delta}(x)$ considered. Of course, realistic systems will always have some subradiant mode decay even in the absence of coupling to the cavity due to loss to free space modes transverse to the cavity. We do not include any such loss mechanism in our model, which allows easy identification of decay to the cavity mode, and subsequent measurement via the cavity output. The limit of weak transverse decay may be approached, for example, by trapping the atoms in lattices with subwavelength-spacing Sutherland and Robicheaux (2016); Bettles et al. (2016).

In addition to the collective atomic superradiant and subradiant modes, the cavity mode itself has a resonance (at $\Delta_{c}=0$ in the absence of atoms). The full cavity-atom system therefore admits $N+1$ normal modes, whose nature may be mixed. In particular, the superradiant atomic and cavity modes exhibit an avoided crossing near $\bar{\Delta}=\Delta_{c}=0$, leading to the characteristic “vacuum Rabi splitting”Sanchez-Mondragon et al. (1983); Carmichael (2007) behavior illustrated in Fig. 1. In contrast to the superradiant mode, the subradiant modes remain unperturbed by the cavity mode frequency. In the remainder of the paper, when considering collective atomic modes we will generally work in the weak-coupling regime where the cavity mode is far-detuned from resonance such that the atomic modes can be considered in isolation.

IV.1.2 Lattice system

As an example, we consider a small system of $8$ atoms confined in a lattice commensurate with the cavity mode, such that atoms are confined near maxima of $|g(x)|$, and we will assume that the atoms are in a Mott insulator (MI) state Morsch and Oberthaler (2006) of exactly one atom per lattice site. In practice, a harmonic potential is often used in conjunction with an optical lattice and leads to a “wedding cake” MI ground-state of differing occupancies. A system of one atom per site can then be achieved by manipulating the sites with excess occupancy Weitenberg et al. (2011). The quantum phase of the atomic ensemble can be important in determining the optical response, and it plays a central part in our stochastic method since it determines the joint-probability distribution from which single realizations of atomic positions are sampled. For an MI state with one atom per site the stochastic sampling for each site is similar to that for independent atoms Jenkins and Ruostekoski (2012); Lee et al. (2016): We take the atoms to be confined along their one degree of freedom ($x$) by an external lattice potential commensurate with the cavity mode, giving rise to Wannier functions $\phi_{i}(x)=\phi(x-\ell_{i})$ centered at antinodes $\ell_{i}$ of the coupling strength $g(x)=g_{0}\sin(kx)$ in a Fabry-Perot cavity. Here $\phi(x)$ is given by

[TABLE]

where $L_{x}$ governs the confinement of the well, and the linear density of a single atom in site $i$ is $|\phi_{{i}}({x})|^{2}$. A single realization of discrete atomic positions for the stochastic method of Sec. III is then found by sampling a single position from the distribution $|\phi_{{i}}({x})|^{2}$ for each occupied lattice site. Later, in Sec. VI we will contrast these results with those for atoms in an ideal superfluid state.

In order to couple the subradiant atomic modes to the cavity optical response, we apply a simple linear gradient to the detuning, such that $\bar{\Delta}(x)=\bar{\Delta}_{0}+\bar{\Delta}_{1}x/\lambda$. Ensemble averaging the eigenmodes found from many stochastic realizations of discrete atomic positions leads to the distributions of eigenmode decay rates and frequencies shown in Fig. 2. A single superradiant mode appears, whose distribution of decay rates and frequencies mirrors closely the real and imaginary parts, respectively, of the distribution of $\sum_{j}g^{2}(X_{j})/\bar{\kappa}$. In addition, there are also $7$ non-trivial subradiant modes with linewidths three orders of magnitude smaller than the superradiant mode. In this example the subradiant modes are conveniently separable in frequency, but this is not generally the case. The linear detuning gradient used here is responsible for the uniform frequency spacing of the subradiant modes, and increasing the number of atoms in the system will lead to additional subradient modes with the same spacing.

IV.2 Optical response in the low light intensity limit

In the previous section we illustrated the calculation of the many-atom collective eigenmodes in the cavity system, where the atomic detuning is subject to a linear spatial variation that couples the sub-radiant eigenmodes to the cavity mode. Next we show that these modes, when coupled in this manner, can be excited and identified in the optical response to incident light. We then show that these modes can also be preferentially excited by suitable transverse light.

When coupled to the cavity mode, each of the collective modes can be excited by axially pumping the cavity on resonance with the corresponding mode frequency. Figure 3 demonstrates the resonances in the steady-state cavity intensity $|\langle\hat{a}\rangle|^{2}$, directly proportional to the transmitted light intensity, obtained from an ensemble of stochastic realizations of Eq. (32) for each value of the axial pump frequency. A full scan of frequencies for a given atom-cavity detuning shows a broad cavity mode resonance, with a linewidth of approximately $\kappa$, a superradiant mode with a linewidth of $\sim Ng_{0}^{2}/\mbox{Im}[\bar{\kappa}]$, and the very narrow and comparatively weakly excited subradiant modes.

While the subradiant modes only couple weakly to the cavity mode, they can also be excited by a suitably tailored transverse pump $h(x)$. Under an ensemble average of many realizations, the spatial profile of the centermost of the subradiant modes in Fig. 2(d) corresponding to the linearly varying detuning can be determined, and is shown in Fig. 4(a). In contrast to the superradiant mode in which the polarizations of neighboring lattice sites are out of phase, in this case the subradiant mode has the centermost sites in phase (the exact distribution depends on the spatial form of $\bar{\Delta}(x)$ along with the other spatial dependence in the system). A simple tailoring of a transverse beam $h(x)$ which might predominantly excite this mode is suggested by the step function in Fig. 4(a)444We have used a rectangular transverse pump beam for clarity in demonstrating the effect. In practice, a rectangular beam with width $\lambda$ may be difficult to realize experimentally. However, any reasonably narrow beam which predominantly and symmetrically excites the central two lattice sites should suffice.. Using this transverse beam as the only pump mechanism, and calculating the steady-state cavity intensity as a function of pump frequency, we see in Fig. 4(b) that indeed the targeted (centermost) subradiant mode is strongly excited at the appropriate resonance frequency. The neighboring subradiant mode also shows a weak response, but no other subradiant modes are appreciably excited. In order to demonstrate that the subradiant mode has indeed been excited, in preference to the superradiant mode, Fig. 4(c) shows the overlaps of the excited polarization density $P(x)$ with the spatial profile of the superradiant and the targeted subradiant modes. Over the frequency range where the subradiant mode shows a strong response it can clearly be seen that the superradiant mode is only negligibly excited. The subradiant mode excited by this tailored transverse pump beam has a decay rate distribution peaked at $8.5\times 10^{-4}g_{0}^{2}/\mbox{Re}[\bar{\kappa}]$ [see Fig. 2(c)], which is almost four orders of magnitude smaller than that of the superradiant mode.

IV.3 Eigenmodes in the saturated case

Outside of the limit of low light intensity, where it is necessary to include the effects of saturation, the description of the optical response becomes more complicated. We must consider the coupled system of $3N$ nonlinear equations of motion (30) for $\rho_{ge}^{(j)}$, $\rho_{eg}^{(j)}=\rho_{ge}^{(j)*}$ and $\rho_{ee}^{(j)}$. Given a steady-state solution ${\bf v}_{0}$ of Eqs. (30), the response to a small perturbations ${\bf\Delta v}$ about that point are encoded in the Jacobian matrix ${\bf J}({\bf v}_{0})$ evaluated at the steady state

[TABLE]

The eigenmodes of ${\bf J}({\bf v}_{0})$ can be useful to describe the behavior of the system near this steady state. We have $3N$ eigenmodes of the atom system, which must then be coupled to the eigenmode of the cavity. At low light intensities, $\rho_{ee}^{(j)}$ plays a negligible role in the dynamics, and the eigenmodes can be grouped into sets: $N$ conjugate pairs, which are eigenmodes predominantly representing response of $\rho_{ge}^{(j)}$ and its conjugate, which tend in the limit of low intensity to those discussed in the previous section. There are then $N$ eigenmodes which describe the response of the (almost negligible) $\rho_{ee}^{(j)}$. The latter modes have eigenvalues predominantly determined by the values of $g^{2}(X_{j})/\bar{\kappa}$ for a single realization, and so simply represent the self interaction of an excited state atom via the cavity mode at point $X_{j}$.

V Comparison of semiclassical approximation with full

treatment for two atoms

In this section we compare the stochastic results including saturation of Sec. III with the results of the full quantum field-theoretical representation of the optical response for coherent scattering for the simple test case of two atoms.

In Sec. II we derived quantum field-theoretical equations of motion governing the optical response of an ensemble of atoms in an optical cavity. We then discussed how these equations can be simulated by stochastic electrodynamics simulations in Sec. III. We noted that the factorizations similar to (31) in the derivation of electrodynamics equations of motion introduce semiclassical approximations to the calculations that go beyond the limit of low light intensity. Even in the limit of low light intensity, some approximate factorizations are necessary to describe atoms with multiple electronic ground levels Lee et al. (2016). Here we compare the results of the stochastic simulations with the full field-theoretical solution of a coupled set of equations for correlation functions of atomic densities and polarizations that does not rely on the factorization approximation of Eq. (31). The full field-theoretical solution is exact for atoms which are stationary regarding their center-of-mass motion. In particular, we neglect the potential effect of the cavity field on the atomic recoil or on the confining potential experienced by the atoms.

Equations (19) and (21) show that single particle correlation functions, such as the polarization and level populations, in general depend on second-order correlation functions. In turn, each of these depends on higher-order correlation functions, leading to a hierarchy of equations of motion which continue until the $N$th order. Although this is a full quantum field-theoretical treatment for the system, it is not feasible to solve this full hierarchy of coupled equations for large $N$, which motivates the method of stochastic electrodynamics simulations discussed in this paper. However, it is feasible to analyze the field-theoretical treatment for the simple test case of a pair of two-level atoms, as described in App. B. We are therefore able to compare the semiclassical approximation to the stochastic electrodynamics (Sec. III) with the exact solution.

Again, we take the atoms to be confined in an optical lattice, with linear densities in each site governed by the Wannier functions of Eq. (34) for $i={1,2}$. We assume a MI state of the atoms, with a single atom in each site. As explained in the previous section, for the stochastic method the discrete particle positions are then sampled from the corresponding joint probability distribution for the ground state densities in the absence of cavity light

[TABLE]

In the absence of cavity light, this expression must also be equal to $\rho^{\rm tot}_{2}(x,x^{\prime})$ in the atom number conservation condition given by Eq. (48). And since we have assumed that the atoms are stationary with respect to their center of mass motion, this remains the case even once light enters the cavity. With this identification, the full hierarchy of equations of motion for this system can be solved as shown in App. B. Below, we compare the two techniques for the different manifestations of collective optical response.

V.1 Superradiant mode

As described in the previous section, a two-atom system admits two collective atomic modes: the superradiant mode along with a single subradiant mode. Figure 5 shows the spectrum of the steady-state optical cavity response for the two-atom system for different values of the pump strength (or atomic saturation). We display the light intensity inside the cavity close to resonance with the superradiant eigenmode. The cavity is far-detuned from the pump frequency so that the cavity and collective atomic modes are only weakly coupled. In the limit of low light intensity, the resonance peak is Lorentzian, and we confirm that the stochastic electrodynamics converges to the result obtained by solving the correlation functions. At intermediate intensities, when the stochastic electrodynamics is only approximate, the results become more complicated. Firstly, the results from the full treatment of correlation functions show that structure appears in the resonance profile, with two peaks becoming evident. This structure is a signature of quantum effects which are not able to be fully accounted for by the semiclassical approximation to stochastic electrodynamics. The stochastic electrodynamics results do capture the peak at higher $\bar{\Delta}$ reasonably well, however, the approximation to factorizing two-body correlation functions cannot reproduce the structure of the local minima and the peak at lower $\bar{\Delta}$. As the intensity increases still further, the superradiant peak broadens, the structure becomes less significant and the stochastic results once again give good quantitative agreement.

An advantage of the semiclassical approximation to stochastic electrodynamics is the ability to include spatial correlation effects in a tractable manner even in larger systems. Spatial variation in atomic density, cavity coupling strength and detuning can lead to nontrivial effects in the polarization response of the atoms. For example, since the cavity Lamb shift [${\cal D}^{2}|g(x)|^{2}\Delta_{c}/(\kappa^{2}+\Delta_{c}^{2})$] is a spatially varying quantity, for a given pump frequency only certain values of $x$ are exactly in resonance with the local superradiant mode. The polarization response $P(x)$ of the atoms therefore exhibits a nontrivial spatial structure, which changes as a function of detuning from the superradiant mode frequency. Figure 6 shows the polarization density at a single frequency close to the superradiant resonance frequency, with qualitative agreement between the two approaches even at significant saturation.

The spatial nature of the cavity coupling strength also affects the frequency of the eigenmode. One might expect to estimate the mode position by a simplified two-atom model in which spatial variations were averaged. When $\Delta_{c}$ is sufficiently large that $\bar{\kappa}$ may be treated as approximately constant, this estimate of the superradiant resonance frequency is found by integrating over Eq. (32) and solving the resulting eigenvalue problem, resulting in a resonant frequency of $2{\cal D}^{2}\mbox{Im}\left[g^{2}_{\rm av}/\bar{\kappa}\right]$, where $g^{2}_{\rm av}=\int g^{2}(x)\rho_{1}(x)dx/2$ is the average cavity coupling strength experienced by an atom. In fact, this estimate is not particularly good; a much better estimate arises from setting $g^{2}_{\rm av}=g_{0}^{2}$, i.e., setting the coupling parameter to the maximum experienced by an atom, rather than the average.

V.2 Subradiant mode

Excluding the internal structure of the mode spectral response, the semiclassical approximation to stochastic electrodynamics gives reasonable qualitative agreement for the position and width of the superradiant mode. In contrast, however, once saturation becomes important the semiclassical approximation fails to describe the subradiant mode well. Figure 7 shows the spectrum of the steady-state optical cavity response for the two-atom system for different values of the pump strength. Here, in order to provide some coupling between the subradiant collective modes and the cavity mode (see Sec. IV), we have added a differential atomic level shift between the two sites

[TABLE]

where $\theta(x)$ is the Heaviside function and the center of the two wells is at $x=0$. By construction, the semiclassical stochastic approach coincides exactly (within the sampling error) with the full field-theoretical solution in the low intensity limit, displaying a resonance peak over an order of magnitude narrower than that corresponding to the superradiant mode. In contrast to the superradiant mode, increasing the intensity does not lead to any additional internal structure of the resonance profile in the results from the full quantum treatment. However, the subradiant mode does broaden substantially as the saturation increases. The semiclassical approximation clearly overestimates the broadening and resonance shift once saturation is included. At very high intensities, the subradiant mode becomes indistinguishable from the background created by the broadened superradiant mode and in this limit the semiclassical approximation again provides good quantitative agreement.

V.3 Detuned from collective resonances

Detuned from any particular collective resonance, the semiclassical approximation to stochastic electrodynamics works comparatively well at all studied intensities. Figure 8 shows a comparison with the full hierarchy of equations of motion over a range of saturation strengths when the incident light is detuned from any of the collective eigenmode resonances. Agreement is excellent at high values of saturation, although at intermediate intensities the semiclassical approximation does lead to a small discrepancy with the full treatment. It is instructive to compare these with other approximate treatments. For $N$ atoms in a cavity, each of which interact with the cavity mode with identical coupling strengths $g_{0}$, the system reduces to the Tavis-Cummings model Tavis and Cummings (1968), and in the limit of low light intensity this simplifies to give the equivalent response to that of a single atom with coupling strength $g_{0}\sqrt{N}$. Motivated by this, a common simple approximate treatment of spatial variations of coupling strengths $g(x)$ is to solve a similar model of a single atom with coupling strength $g_{0}\sqrt{N_{\mathrm{eff}}}$, where $g_{0}^{2}N_{\mathrm{eff}}=\int g^{2}(x)\rho_{1}(x)\mathrm{d}x$. The results of this approximation are also shown in Fig. 8, and while they tend to the correct limit for low light intensity, they are in general less accurate than the results of the semiclassical electrodynamics simulations.

V.4 Comparison between correlation functions

In the semiclassical model for each stochastic realization we solve the equations of motion for single-particle variables for internal atomic states. Ensemble-averaging over many realizations of atomic positions establishes spatial correlations between the atoms. However, not all the correlations that result from the full quantum dynamics are incorporated. This can be identified by comparing in more detail the equations for internal atomic states for each stochastic realization in the semiclassical approximation with the full quantum dynamics of the internal states for fixed atomic positions.

In Eqs. (45) and (47) we give examples of the full quantum equations of motion for two of the two-body correlation functions $P_{2}(x;x^{\prime})$ and $\langle:\hat{P}^{+}(x)\hat{P}^{-}(x^{\prime}):\rangle$, where $:\cdot:$ represents normal ordering. The semiclassical factorization approximation together with the stochastic sampling technique synthesizes two-body correlation functions from the products of one-body quantities evaluated for each stochastic realization. For a single realisation of atoms at the discrete positions $X_{j}$ and $X_{l}$, one can obtain the implied semiclassical equations of motion for the products $\rho_{gg}^{(j)}\rho_{ge}^{(l)}$ and $\rho_{ge}^{(j)}\rho_{eg}^{(l)}$ [using Eqs. (30)]. On the other hand, full quantum equations of motion for the two-body quantities $\rho_{gg;ge}^{(j,l)}$ and $\rho_{ge;eg}^{(j,l)}$ can be derived from Eqs. (45) and Eqs. (47), respectively, evaluated at the same discrete positions, by substituting

[TABLE]

and

[TABLE]

and similarly for all other two-body correlation functions. Subsequently factorizing all two-body terms in the resultant equations of motion [in the manner of Eq. (31)] leads to terms consisting of products of one body quantities. Comparing these equations with those obtained by making the factorization approximation at the outset gives some insight into the nature of the semiclassical approximation.

For both of these two-body correlation functions, the semiclassical equations of motion reproduce the diagonal and driving terms (those proportional to $a_{F}$ or $h(x)$) from the full quantum treatment. In the case of $d\rho_{ge;eg}^{(j,l)}/dt$ the off-diagonal term proportional to $\rho_{ee}^{(j)}\rho_{ge}^{(l)}$ is included, while the terms missing from the semiclassical factorized version are identical to those included but with a swap of coordinates $X_{j}\leftrightarrow X_{l}$. A crucial aspect of the stochastic treatment is the averaging over many such stochastic realizations. In the low intensity limit it can be shown that, although not present in the factorized version for a single realization, the swapped coordinate term proportional to $P_{2}(x^{\prime};x)$ in Eq. (45) is included by means of averaging over the stochastic samples Lee et al. (2016).

In contrast, the semiclassical approximation to the equation of motion for $\rho_{ge}^{(j)}\rho_{eg}^{(l)}$ does not reproduce any of the off-diagonal couplings present in the full result for $\rho_{ge;eg}^{(j,l)}$. Since in general in the full quantum picture all two-body correlation functions are coupled, the missing terms mean the semiclassical approximation cannot reproduce the quantum correlations, even when averaged over multiple realizations. However, in the low intensity limit the only two-body correlation functions which are important for the cavity optical response are $P_{2}(x;x^{\prime})$ and $P_{2}(x^{\prime};x)$, and hence the same factorization approach is able to fully reproduce the correlations in this limit.

Figure 9 compares two different two-body correlation functions obtained from the full treatment of the hierarchy of correlation functions with factorized semiclassical approximations at two different frequencies near resonant with the superradiant eigenmode and within the range spanned by Fig 5. It can be seen that when the semiclassical approximation to stochastic electrodynamics is least accurate, the full two-body correlation functions exhibit structure that cannot be reproduced by any factorization approximation. In contrast, where the factorized and full correlation functions show similar qualitative features the semiclassical approximation agrees rather well with the full treatment for the optical response of this collective mode.

In summary, while the semiclassical approximation to stochastic electrodynamics has no effect in the limit of low light intensity and the two methods are the same (beyond the sampling noise), at intermediate intensities when saturation effects become important the semiclassical approximation only gives qualitative agreement for the excitation of the superradiant eigenmode and is unable to describe the excitation of the subradiant eigenmode. In the limit of very high intensities, quantitative agreement is restored although the collective features have become substantially power-broadened in this limit. Detuned from the collective modes, the semiclassical approximation works well at all intensities.

The semiclassical approximation to stochastic simulations is designed to capture spatially-dependent correlation effects induced by the scattered light. In cavities such effects are weaker than in free space, owing to the unattenuated long-range dipole-dipole interactions between the atoms. Furthermore, light confinement and directed light emission enhance quantum effects in cavities compared with quantum fluctuations in free space. Consequently, finding differences in the cavity response of a small two-atom system between the semiclassical approximation and the exact quantum result is not entirely surprising, but such deviations are likely to become smaller in larger atomic ensembles and in multimode cavities.

VI Diagnostics of atomic quantum phase

An advantage of the stochastic electrodynamics simulations is that different atom distributions and statistics are incorporated in the joint probability distribution from which the stochastic realizations of discrete atomic positions are sampled. In many situations an independent atom sampling can be employed; for example, for an ideal BEC or for uncorrelated (in the absence of light) classical atoms confined in an optical lattice potential the positions are sampled from a distribution proportional to the total linear density $\rho_{1}(x)=\sum_{i}n_{i}|\phi_{i}(x)|^{2}dx$ where $\phi_{i}(x)$ is the wavefunction for site $i$ with site population $n_{i}$. In contrast, we described how to simulate a Mott Insulator (MI) state in Sec. IV.

The atom statistics can affect the cavity system response when site-to-site spatial variation in the Wannier functions, atomic detunings or cavity coupling strength are present. The optical response can therefore be used in principle as a diagnostic tool to query the quantum phase of the atoms. As another illustration of the stochastic electrodynamics simulations we discuss two examples below where the MI and BEC states differ: firstly, the distribution of subradiant modes, and secondly the case of a cavity transmission spectrum when some sites are masked from the cavity, following the model of Ref. Mekhov et al. (2007). Off-resonance spontaneous scattering in optical lattices has actively studied as a diagnostic tool for quantum and thermal states of atoms; see for instance Refs. Javanainen and Ruostekoski (2003); Ruostekoski et al. (2009); Rist et al. (2010); Łakomy et al. (2009); Douglas and Burnett (2010, 2011); Ye et al. (2011); Sandner et al. (2011); Cordobes Aguilar et al. (2014); Elliott et al. (2015b).

VI.1 Subradiant modes for a BEC

In Fig. 2 we presented the subradiant modes which appear when $\bar{\Delta}(x)$ was not spatially constant but had a weak constant gradient, for a system of $8$ atoms in a MI state of 1 atom per lattice site. Seven subradiant and one superradiant modes are clearly (in this case) resolvable. The subradiant modes are only coupled to the cavity mode when the detuning $\bar{\Delta}(x)$ varies between lattice sites. Since different atomic quantum phases lead to different atom probability distributions, they will in turn lead to different sampled detunings in the possible stochastic realizations of atom positions. The nature of the subradiant modes might therefore be expected to vary with the atomic quantum phase.

Indeed, if we instead consider the same system as Fig. 2 but in a superfluid BEC phase, the site number fluctuations significantly alter both the lifetimes and resonance frequencies of the subradiant modes. Modes are no longer easily separable in frequency, and the distribution of lifetimes acquires a two peak structure, as shown in Fig. 10. In contrast, the superradiant collective mode is not significantly affected by the atomic quantum phase, since it is not greatly affected by any small changes in $\bar{\Delta}(x)$. For simplicity, we have assumed that the lattice potential is unaffected by the spatially dependent $\bar{\Delta}(x)$, so that the Wannier functions for each lattice site are identical.

VI.2 Number fluctuations in the cavity transmission spectrum

A clear example of the importance of the atomic quantum phase was presented by Mekhov, Maschler and Ritsch Mekhov et al. (2007) when they considered the transmission spectra of a cavity containing an optical lattice commensurate with the cavity mode, but such that only $n_{A}\leq n_{L}$ of the total $n_{L}$ sites occupied by atoms interact with the cavity mode (Ref. Mekhov et al. (2007) proposes to achieve this by tilting the axis of the 1D optical lattice with respect to the cavity axis).

In contrast to most of the results in this paper, we concentrate in this section on the dressed mode of the cavity-atom system which is dominated by the bare cavity mode (illustrated by the red dashed line in Fig. 1). We work in the regime where the bare atomic transition is detuned far from resonance with the driving axial cavity pump, but where the cavity mode itself is near resonant with the pump ($|\bar{\Delta}|\gg|\Delta_{c}|$). In the absence of atoms, the empty cavity exhibits a resonance peak in the transmitted photon intensity centered at $\Delta_{c}=0$ with a width dictated by the cavity loss rate $\kappa$. The presence of far-detuned atoms in the cavity shifts the position of this resonance, proportional to the number of atoms in the cavity (with reference to Fig. 1, the shift is that of the dressed state of the cavity-atom system from the bare cavity resonance). This shift, together with the geometry of the system described above, allows the signature of the atomic many-body state to be seen in the cavity transmission spectrum.

For a MI state with exactly one atom per site, $n_{A}$ atoms interact with the cavity mode to shift the cavity resonance frequency. For an ideal BEC superfluid however, site number fluctuations mean that a given experimental realization can involve any number of atoms between [math] and $N$, each realization giving rise to a different resonance shift. The signature of the superfluid state is therefore a comb-like transmission spectra, compared to the single peaked response of the MI state.

Our stochastic simulations are well suited to tackle this system, and the contrasting spectra between the two atomic phases are shown in Fig. 11(a) Here, we have simulated a system of $12$ atoms with in an optical lattice with narrow Wannier functions and with only half of the sites of the optical lattice able to interact with the cavity light. The numerical calculation extends the results of Ref. Mekhov et al. (2007) to situations where the atomic transitions are allowed to be saturated by high pump intensities and the coupling strength is allowed to vary in space. In Fig. 11(b) we show results for a lattice where the spatial widths of Wannier functions are on the order of $\lambda/2$, leading to a range of cavity coupling strengths $g(x)$ experienced by the atoms. In Fig. 11(c), we show the effect of significant saturation. For each case, we show both the coherently scattered light transmission $\propto|\langle\hat{a}\rangle|^{2}$ and the total intensity transmission $\propto\langle|\hat{a}|^{2}\rangle$ that also includes the incoherently scattered light component. The two signals for the MI state are very similar, but they differ more in the case of a BEC due to the fluctuating density, as in the latter case the incoherently scattered light $\propto\langle|\hat{a}|^{2}\rangle-|\langle\hat{a}\rangle|^{2}$ has a stronger effect on the optical response.

In comparison to Fig. 11(a), both widening the widths of the Wannier functions and including saturation can be seen to lead to a loss of resolution in the comb shape for the superfluid response and a broadening for the MI case. The larger width of the Wannier function leads to a broadening of the resonance peaks because the wider atom distribution samples a greater range of cavity coupling strengths. Furthermore, since the geometry of the system leads to a decrease in the total experienced coupling strength, the peak is shifted to smaller $\Delta_{c}$. The distribution of $\sum_{j}g(X_{j})$ sampled is not symmetric, and this is reflected in the appearance of the resonance peak. Higher intensities decrease the resolution through power broadening of the spectra. Nonetheless, although such considerations may mean the distinctive comb-like structure is not always resolvable, the differing widths of the superfluid and MI resonances means that the spectra can still be used as a diagnostic tool for the quantum phase of the atoms.

VII Concluding remarks

We have formulated stochastic electrodynamics for many-atom systems in a cavity. In this work the approach was implemented in a semiclassical approximation that could be extended also to include quantum fluctuations. However, in the limit of low light intensity and for two-level atoms, the presented simulations are limited in accuracy only by the sampling error introduced by the finite number of stochastic samples used. Our formulation of the stochastic simulations has the advantage that atomic position correlations and spatially dependent potentials and couplings are readily included.

We have shown that a system of atoms in a cavity can exhibit a collective optical response with a strongly enhanced superradiant mode, but also with a number of weakly coupled subradiant modes with very narrow linewidths. These subradiant modes can be coupled to the cavity response via spatially dependent detunings, and can exhibit distinct spatial profiles allowing them to be driven by tailoring the shape and phase of a transverse optical pump.

This phenomenon also raises a possible mechanism for the storage of light by many-body atomic excitations in the cavity. Provided that a controllable spatially-dependent level shift can be introduced for the atoms when the cavity mode is driven, the modes that are subradiant in the absence of the spatially-dependent shifts can be directly excited by the driving field. The idea is related to the analysis of Ref. Facchinetti et al. (2016) where a subradiant mode was driven in a planar optical lattice in free space by coupling the different atomic polarization components with the Zeeman level shifts. After a suitable excitation pulse, the level shifts are turned off, such that the finite, nonvanishing resonance linewidth of the excited subradiant mode (of the case with spatially-dependent level shifts) is tuned as close as possible to zero. The excitation may then become trapped in the subradiant mode in a way that a decay via the cavity mode may in principle exhibit an extremely long lifetime. The atoms may still decay via the free space modes perpendicular to the cavity, but such a decay is sensitive to the spatial arrangement of the atoms, and the decay can be weak for regularly spaced array of atoms with a subwavelength lattice constant Sutherland and Robicheaux (2016); Bettles et al. (2016).

The semiclassical stochastic method was compared with the full quantum solution, revealing quantum features in the optical response. These are particularly prominent in the excitation of narrow subradiant resonances at high intensities. We also showed how the atomic quantum phase of the ensemble may be discerned from the the cavity optical response due to the different atomic position correlations in differing phases. The quantum statistics of the atoms inside the cavity is mapped onto the optical signal, and the statistics is reflected in different ways in the coherently and incoherently scattered light.

Appendix A Polarization and correlation function equations of motion for

multi-level atoms

Equations (19) and (21) gave the optical response for a system of two-level atoms. For the case of multilevel atoms the response follows from the general equations of motion (13) and (14), resulting in the coupled equations for the set of all possible polarizations and coherences

[TABLE]

where the repeated indices $\zeta$ and $\tau$ should be implicitly summed over, and $\bar{\Delta}_{a\nu b\eta}=\Delta_{b\eta}-\Delta_{a\nu}$. Here we explicitly indicate only the non-local position dependence of the atomic field operators, and we have also introduced the convenient tensor

[TABLE]

Taking expectation values of Eqs. (40) gives equations of motion for all one-body correlation functions, and knowledge of those correlation functions solves the problem of optical response. However, as in the two-level case, the one-body correlation functions depend in turn upon two-body correlation functions, leading to the hierarchy of equations of motion for multi-level atoms. The number of correlation functions involved in the hierarchy increases rapidly with the number of levels present in the atoms, and the complexity of the multi-level case can therefore be substantially greater than in the two-level case.

Appendix B Hierarchy of equations for the optical response of atoms in a cavity

B.1 Limit of low light intensity

In the limit of low light intensity, Eqs. (13) and (14), along with commutators such as (18), may be used to write the hierarchy of coupled integral equations describing the optical response of an ensemble of two-level atoms in an optical cavity in the compact form

[TABLE]

analogous to the free-space case of Refs. Ruostekoski and Javanainen (1997a) and Lee et al. (2016). Here we have neglected all terms of greater than first order in the light field amplitude of excited state operators. Similar to Ref. Lee et al. (2016), we have defined the $\ell^{\mathrm{th}}$ order one-dimensional correlation functions in the limit of low light intensity as

[TABLE]

and again used the definitions $\hat{P}^{+}=\mathcal{D}\hat{\psi}^{\dagger}_{g}\hat{\psi}_{e}$, $h(x)={\bf h}(x)\cdot\boldsymbol{\mathbf{d}}_{eg}/{\cal D}$, $g(x)={\bf g}(x)\cdot\boldsymbol{\mathbf{d}}_{eg}/{\cal D}$, and $G_{c}(x,x^{\prime})=g(x)g^{*}(x^{\prime})$.

Equation (42) shows that, as in the free space case, the $\ell^{\mathrm{th}}$ order correlation function $P_{\ell}$ depends on the integral over $P_{\ell+1}$, leading to a hierarchy of equations of motion which terminates only at $\ell=N$. Avoiding the need to solve this hierarchy of equations, the stochastic technique presented in Sec. III gives a computationally efficient method to obtain the optical response. For the low light intensity case corresponding to Eq. (42), the stochastic method solves the linear set of equations (32) for stochastic realizations of fixed atom positions. An argument akin to that in App. B of Ref. Lee et al. (2016) shows that subsequent averaging over an ensemble of such realizations reproduces the full dynamics of the correlation functions as dictated by Eq. (42).

B.2 Including saturation for a system of two atoms

Going beyond the limit of low light intensity to include the effects of saturation rapidly increases the complexity of the hierarchy of equations of motion governing the correlation functions, since one must now account for all $4^{n}$ $n$-body correlation functions for every $n<N$. However, for a small system of just two atoms the hierarchy terminates with the $16$ equations of motion for the two-body correlation functions $\langle\hat{\psi}^{\dagger}_{i}(x)\hat{\psi}^{\dagger}_{j}(x^{\prime})\hat{\psi}_{k}(x^{\prime})\hat{\psi}_{m}(x)\rangle$, where $i,j,k,m$ can refer to the ground or excited atomic state. In such a system the direct steady-state solution of the full hierarchy for the optical response is therefore a realistic, if somewhat tedious, approach.

The equations of motion for the two-body correlation functions may be derived using Eqs.  (13) and (14), and terms reordered using commutators similar to Eq. (18). For example, the equations of motion for the correlation functions $P_{2}(x;x^{\prime})=\mathcal{D}\langle\hat{\psi}^{\dagger}_{g}(x)\hat{\psi}^{\dagger}_{g}(x^{\prime})\hat{\psi}_{e}(x^{\prime})\hat{\psi}_{g}(x)\rangle$ and $\langle\hat{\psi}^{\dagger}_{g}(x)\hat{\psi}^{\dagger}_{e}(x^{\prime})\hat{\psi}_{g}(x^{\prime})\hat{\psi}_{e}(x)\rangle$ are governed by

[TABLE]

[TABLE]

Deriving all $16$ such equations of motion, solving for the steady state, and including the atom number conservation relation

[TABLE]

leads to a set of $16$ linear independent equations which can be solved numerically to give the set of all two-body correlation functions. Here $\rho^{\rm tot}_{2}(x,x^{\prime})$ is the two-body density-density correlation function that we assume is known. Consistent with our earlier stochastic treatment, we have assumed that the atoms are stationary, consequently $\rho^{\rm tot}_{2}(x,x^{\prime})$ is invariant and has a form dictated by the trapping potential, and atom statistics. Once the two-body correlation functions are obtained, polarization and excited state densities then follow from the steady state of the expectation values of Eqs. (19) and (21), and the hierarchy of equations has then been solved.