Exact electrodynamics versus standard optics for a slab of cold dense gas

TL;DR

This paper compares standard electrodynamics and atom-by-atom simulations for light propagation in cold dense gases, revealing qualitative differences when density and wavelength conditions lead to strong dipole-dipole interactions.

Contribution

It demonstrates that mean-field electrodynamics may fail in dense gases where atom correlations become significant, highlighting the importance of microscopic simulations.

Findings

Deviations from mean-field theory scale with $ ho k^{-3}$.

Noticeable effects appear at $ ho k^{-3} oughly 10^{-2}$.

Simulations show shifts in resonance lines consistent with Lorentz-Lorenz and cooperative Lamb shifts.

Abstract

We study light propagation through a slab of cold gas using both the standard electrodynamics of polarizable media, and massive atom-by-atom simulations of the electrodynamics. The main finding is that the predictions from the two methods may differ qualitatively when the density of the atomic sample $\rho$ and the wavenumber of resonant light $k$ satisfy $\rho k^{-3}\gtrsim 1$. The reason is that the standard electrodynamics is a mean-field theory, whereas for sufficiently strong light-mediated dipole-dipole interactions the atomic sample becomes correlated. The deviations from mean-field theory appear to scale with the parameter $\rho k^{-3}$, and we demonstrate noticeable effects already at $\rho k^{-3} \simeq 10^{-2}$. In dilute gases and in gases with an added inhomogeneous broadening the simulations show shifts of the resonance lines in qualitative agreement with the predicted Lorentz-Lorenz shift and "cooperative Lamb shift", but the quantitative agreement is unsatisfactory. Our interpretation is that the microscopic basis for the local-field corrections in electrodynamics is not fully understood.

Equations

$$\hat{\mathbf{P}}(\mathbf{r})={\mathbf{d}}_{ge}{\psi }_g^{†}(\mathbf{r}){\psi }_e(\mathbf{r}),$$

$$\hat{\mathbf{E}}(\mathbf{r})={\hat{\mathbf{E}}}_0(\mathbf{r})+\int d^3{r}^{\prime }\mathsf{G}(\mathbf{r},{\mathbf{r}}^{\prime })\hat{\mathbf{P}}({\mathbf{r}}^{\prime }),$$

$${\rho }_1({\mathbf{r}}_1)$$

$${\rho }_2({\mathbf{r}}_1,{\mathbf{r}}_2)$$

$${\mathbf{P}}_1(;{\mathbf{r}}_1)$$

$${\mathbf{P}}_2({\mathbf{r}}_1;{\mathbf{r}}_2)$$

$${\mathbf{P}}_3({\mathbf{r}}_1,{\mathbf{r}}_2;{\mathbf{r}}_3)$$

$$\dot{\mathbf{P}}({\mathbf{r}}_1)$$

$${\dot{\mathbf{P}}}_2({\mathbf{r}}_1;{\mathbf{r}}_2)$$

$${\dot{\mathbf{P}}}_3({\mathbf{r}}_1,{\mathbf{r}}_2;{\mathbf{r}}_3)$$

$$\gamma =\frac{{\mathcal{D}}^2{k}_0^3}{6\pi \hslash {ϵ}_0}.$$

$${\mathbf{P}}_2({\mathbf{r}}_1;{\mathbf{r}}_2)=\rho ({\mathbf{r}}_1)\mathbf{P}({\mathbf{r}}_2).$$

$$\dot{\mathbf{P}}({\mathbf{r}}_1)$$

$$\mathbf{E}(\mathbf{r})={\mathbf{E}}_0(\mathbf{r})+\int d^3{r}^{\prime }\mathsf{G}(\mathbf{r},{\mathbf{r}}^{\prime })\mathbf{P}({\mathbf{r}}^{\prime }).$$

$${\rho }_2({\mathbf{r}}_1,{\mathbf{r}}_2)=⟨\sum _{i\ne j}\delta ({\mathbf{r}}_1-{\mathbf{X}}_i)\delta ({\mathbf{r}}_2-{\mathbf{X}}_j)⟩,$$

$$\dot{\mathbf{d}}({\mathbf{X}}_i)=(i\Delta -\gamma )\mathbf{d}({\mathbf{X}}_i)+i\zeta {\mathbf{E}}_0({\mathbf{X}}_i)+i\zeta \sum _{j\ne i}\mathsf{G}({\mathbf{X}}_i,{\mathbf{X}}_j)\mathbf{d}({\mathbf{X}}_j).$$

$$\mathbf{d}({\mathbf{X}}_i)=\alpha \left({\mathbf{E}}_0({\mathbf{X}}_i)+\sum _{j\ne i}\mathsf{G}({\mathbf{X}}_i,{\mathbf{X}}_j)\mathbf{d}({\mathbf{X}}_j)\right),$$

$$\alpha =-\frac{\zeta }{\Delta +i\gamma }=-\frac{{\mathcal{D}}^2}{\hslash }\frac{1}{\Delta +i\gamma }$$

$$\mathbf{E}({\mathbf{X}}_i)={\mathbf{E}}_0({\mathbf{X}}_i)+\alpha \sum _{j\ne i}\mathsf{G}({\mathbf{X}}_i,{\mathbf{X}}_j)\mathbf{E}({\mathbf{X}}_j).$$

$$\mathbf{E}(\mathbf{r})={\mathbf{E}}_0(\mathbf{r})+\alpha \sum _j\mathsf{G}(\mathbf{r},{\mathbf{X}}_j)\mathbf{E}({\mathbf{X}}_j)$$

$$k=c=\hslash =\frac{1}{4\pi {ϵ}_0}=1$$

$$\mathbf{E}(\mathbf{r})=\mathsf{G}(\mathbf{r},{\mathbf{r}}_0)\mathbf{d},\mathbf{B}(\mathbf{r})=\mathsf{H}(\mathbf{r},{\mathbf{r}}_0)\mathbf{d},$$

$${\mathsf{G}}_{ij}(\mathbf{r},{\mathbf{r}}_0)=$$

$${\hat{\mathbf{e}}}_i\cdot \left\{(\hat{\mathbf{n}}\times {\hat{\mathbf{e}}}_j)\times \hat{\mathbf{n}}+[3\hat{\mathbf{n}}(\hat{\mathbf{n}}\cdot {\hat{\mathbf{e}}}_j)-{\hat{\mathbf{e}}}_j]\left(\frac{1}{r^2}-\frac{i}{r}\right)\right\}\frac{e^{ir}}{r},$$

$${\mathsf{H}}_{ij}(\mathbf{r},{\mathbf{r}}_0)={\hat{\mathbf{e}}}_i\cdot \hat{\mathbf{n}}\times {\hat{\mathbf{e}}}_j\left(1-\frac{1}{ir}\right)\frac{e^{ir}}{r}.$$

$$\mathbf{B}(\mathbf{r})=-i\nabla \times \mathbf{E}(\mathbf{r}),$$

$$=$$

$$\mathbf{S}$$

$$\Delta =\delta \gamma ,$$

$$\mathcal{D}=\sqrt{\frac{3\gamma }{2}},$$

$$\alpha =-\frac{3}{2}\frac{1}{\delta +i}.$$

$$P=\frac{|\alpha {|}^2|{\mathbf{E}}_0{|}^2}{3}=\frac{3|{\mathbf{E}}_0{|}^2}{4(1+{\delta }^2)}.$$

$$I_0=\frac{1}{8\pi }|{\mathbf{E}}_0{|}^2.$$

$$\sigma (\delta )=\frac{8\pi }{3}|\alpha {|}^2=\frac{6\pi }{1+{\delta }^2}.$$

Full text

I Introduction

With laser cooling and trapping, and evaporative cooling, it is now experimentally possible to prepare what is arguably the most elementary medium for light propagation, atoms effectively at standstill. More specifically, a cold enough gas presents what according to the long-standing terminology of laser spectroscopists is a homogeneously broadened medium; the atoms move only a small fraction of the wavelength over the time it takes their internal state to relax to steady state. Each atom in the sample is subject not only to the driving light but also to the light sent by all other atoms. On the microscopic level this makes the problem of light propagation in the medium a major challenge. This subject is obviously very old, but cold atomic samples afford an opportunity for experiments in unprecedented regimes and in unprecedented detail. Correspondingly, new experiments are emerging rapidly Meir et al. (2014); Pellegrino et al. (2014); Kwong et al. (2015); Bromley et al. (2016); Guerin et al. (2016a); Bons et al. (2016); Jenkins et al. (2016a); Jennewein et al. (2016); Roof et al. (2016); Araújo et al. (2016).

On the theoretical side, there is the old idea that one could solve the problem of light propagation in a medium on an atom-by-atom basis directly numerically Javanainen et al. (1999). The growing throughput of computers available to researchers is making such a plan practical. These methods, whether called classical-electrodynamics simulations or coupled-dipole simulations, are now a routine theoretical tool Chomaz et al. (2012); Jenkins and Ruostekoski (2012a); T. Bienaimé and Kaiser (2013); Javanainen et al. (2014); Pellegrino et al. (2014); Javanainen and Ruostekoski (2016); Sutherland and Robicheaux (2016a); Bromley et al. (2016); Bettles et al. (2015); Schilder et al. (2016); Yoo and Paik (2016); Jenkins et al. (2016a); Zhu et al. (2016); Sutherland and Robicheaux (2016b); Bettles et al. (2016); Jennewein et al. (2016). Closely related numerical techniques based on the analysis of the eigenstates of the coupled system of the light and the atoms Rusek et al. (1996); Jenkins and Ruostekoski (2012b, a); Antezza and Castin (2013); Balik et al. (2013); Skipetrov and Sokolov (2014); Facchinetti et al. (2016) or density matrices and quantum trajectories Berhane and Kennedy (2000); Clemens et al. (2003); Jones et al. (2017) are also widely used today. Other ideas drawn from the theory of radiative transfer Chandrasekhar (1960); Ishimaru (1997) and multiple scattering Lagendijk and van Tiggelen (1996); van Rossum and Nieuwenhuizen (1999), amended with numerics, also have potential to make inroads into the questions about light propagation in atomic media Guerin et al. (2016b).

The present work started with our chance observation in numerical light propagation simulations that the density-dependent Lorentz-Lorenz (LL) shift of the atomic resonance Lorentz (2011), a quintessential local-field correction, is absent in cold, dense atomic samples Javanainen et al. (2014). Delving into the problem deeper, we discovered that the standard electrodynamics of polarizable media (EDPM) Jackson (1999); Born and Wolf (1999), and the resulting standard optics, may fail qualitatively in cold, dense atomic samples Javanainen and Ruostekoski (2016). However, adding inhomogeneous broadening that mimics the Doppler shifts of thermal atoms restored the behavior of standard EDPM. Along the way we made a number of additional qualitative and quantitative observations. Among others, in our analysis it emerged that the so-called “cooperative Lamb shift” in a slab of atomic matter Friedberg et al. (1973) can be explained in standard electrodynamics as an etalon effect due to the reflections of light from the faces of the slab of matter. In retrospect the issues with the EDPM are not much of a surprise, as EDPM is an effective-medium mean-field theory (MFT) and is bound to fail when the light-mediated dipole-dipole interactions make the atomic sample strongly correlated. Nonetheless, they beg for questions about the meaning of cooperative light-atom interactions Javanainen and Ruostekoski (2016); Guerin et al. (2016b) and the limits of the predictive power of EDPM.

In this paper we add technical details and new results and discussions that illuminate, support and expand on the observations in Refs. Javanainen et al. (2014); Javanainen and Ruostekoski (2016). We begin in Sec. II by reviewing the theoretical basis of our classical-electrodynamics simulations. In Sec. III we orient the reader to concepts such as coherent and incoherent scattering and cooperative line shifts and linewidths by presenting simple analytically solvable examples. The core of the present paper, however, is about the comparison of the EDPM solutions and numerical simulations of the response of a gas of atoms confined to a slab to light at normal incidence. What exactly is involved here is explained in Sec. IV. The remaining sections V and VI present and discuss the results.

II Background

The purpose of the present Sec. II is threefold. First we briefly summarize our fully quantum mechanical approach to light-matter interactions as in Refs. Ruostekoski and Javanainen (1997a, b), especially as it comes to a hierarchy of equations of motion for the correlation functions that involve polarization of the atoms and densities of the atoms in different points in space. Second, following Ref. Javanainen et al. (1999), we explain how and in what sense we may solve the hierarchy for the correlation functions numerically using classical-electrodynamics simulations. Third, both in our numerical computations and in the discussions of this paper we almost exclusive use certain natural units for microscopic theory of light propagation in dipolar samples. We conclude by introducing these units.

II.1 Quantum theory of light propagation in dipolar medium

Our approach Ruostekoski and Javanainen (1997a) begins with the boson field operators for ground-state and excited-state atoms $\psi_{g}({\bf r})$ and $\psi_{e}({\bf r})$. The labels $g$ and $e$ implicitly include the Zeeman state labels of the angular-momentum degenerate energy levels. We adopt a summation convention whereby repeated indices $g$ and $e$ in a product are summed over. Since we always deal with pairs of atom field operators, we believe that our scheme is also valid for fermionic atoms.

The atoms are coupled to the quantized electrodynamic field via the dipole interaction. We deviate from the dominant practice in that we adopt the Power-Zienau-Woolley viewpoint Power and Zienau (1959); Power (1964); Cohen-Tannaudji et al. (1989), whereby the primary quantized variable having to do with the electromagnetic field is the electric displacement $\hat{\bf D}$ not the electric field $\hat{\bf E}$. The result is a quantum field theory that in appearance closely resembles the usual EDPM. However, here we deliberate phrase our arguments in terms of the electric field.

To begin with, we have the positive frequency part of the polarization operator for the atoms

[TABLE]

where ${\bf d}_{ge}$ are the dipole moment matrix elements. If there is a difference between positive- and negative-frequency parts of the quantity in question, we write down the positive-frequency part without further comment. Correspondingly, when we consider analogous classical quantities, we alway assume a dominant frequency of the driving light $\omega$ in the problem, so that a classical counterpart of a positive-frequency part of a quantity may be written as, say, $\langle{\hat{\bf P}({\bf r},t)}\rangle=e^{-i\omega t}{\bf P}({\bf r},t)$, where ${\bf P}({\bf r},t)$ now assumedly varies little over the time scale $\omega^{-1}$. We then express the classical polarization (and electric field, and dipole moment, and so on) in terms of the slowly varying part ${\bf P}({\bf r},t)$ without further ado. The physical polarization, a real quantity, is {\hbox{\frac{1}{2}}}{\bf P}({\bf r},t)e^{-i\omega t}+{\rm c.c.}. This convention is, of course, deeply ingrained in optical physics and quantum optics.

In analogy to classical electrodynamics, the electric field operator is related to the polarization operator by

[TABLE]

where $\hat{\bf E}_{0}\equiv\hat{\bf D}_{0}/\epsilon_{0}$ is the electric field in the absence of matter, and ${\sf G}({\bf r},{\bf r}^{\prime})$, a $3\times 3$ matrix, is the dipole propagator such that ${\sf G}({\bf r},{\bf r}^{\prime}){\bf d}$ is the electric field at ${\bf r}$ from an oscillating dipole moment ${\bf d}$ at ${\bf r}^{\prime}$ Jackson (1999); Born and Wolf (1999). The dipole propagator should also include a singularity term $-\frac{1}{3}\delta({\bf r}-{\bf r}^{\prime})/\epsilon_{0}$ Jackson (1999) in order for (2) to be the integral representation of the correct Maxwell’s wave equation, although its presence in the equations of motion for matter is a subtle matter Ruostekoski and Javanainen (1997b) to which we briefly return later. Integrals involving ${\sf G}$ are typically not absolutely convergent either at small or large $|{\bf r}-{\bf r}^{\prime}|$. The values of such integrals depend on how they are done. This type of ambiguities are widespread in the theory of the electrodynamics of dipolar media, and are often difficult to resolve. This is one of the reasons why we think that one should be suspicious of any and all “physical” approximations in this field.

One can have the electric field radiated by matter fall back on the matter and change the atomic dipole moments, hence polarization. The self-field giving rise to radiation reaction and transition linewidths can be handled with the Markov and Born approximations of quantum optics as usual, but otherwise the ensuing operator equations are (most likely) impossible to solve directly. Instead we go to expectation values.

Here we proceed under the limit of low light intensity, only keeping the leading nontrivial contribution in the strength of the incoming field $\hat{\bf E}_{0}$, and specialize to the case when the angular momenta of the levels are $J_{g}=0$ and $J_{e}=1$, c.f. Ref. Lee et al. (2016). Specifically, introduce normally ordered correlation functions for ground state density and correlations between polarization and ground state density as

[TABLE]

A rigorous quantum mechanical analysis Ruostekoski and Javanainen (1997a) finds a hierarchy of equation of motions for these expectation values beginning with

[TABLE]

[TABLE]

[TABLE]

and continuing along these lines all the way up to the order equal to the number of the atoms $N$. Here the detuning $\Delta=\omega-\omega_{0}$ is the difference of the frequency of the driving light $\omega$ from the atomic resonance frequency $\omega_{0}$ and $\gamma$ is the HWHM linewidth of the transition. Further, we have $\zeta={\cal D}^{2}/\hbar$, where ${\cal D}$ is the reduced dipole matrix element related to the linewidth and the wave number of resonant light $k_{0}=\omega_{0}/c$ by

[TABLE]

In the usual way we assume that the incoming quantum light is in a coherent state, and replace it with the classical electric field ${\bf E}_{0}$.

Let us factor in Eq. (5) the second-order correlation function ${\bf P}_{2}$ as in

[TABLE]

The resulting approximate equation,

[TABLE]

tells us that the polarization at ${\bf r}_{1}$, basically the dipole moment of an atom at ${\bf r}_{1}$, evolves under the joint influence of the incoming field and the electric field radiated from the polarization of the atoms, as if the atoms were smeared out continuously in space. This is an effective-medium MFT. Moreover, it is easy to see in explicit examples that this MFT is the same as the standard EDPM for the atoms.

The equation of motion of the second-order correlation function ${\bf P}_{2}({\bf r}_{1};{\bf r}_{2})$ has a similar structure except for one crucial point, the term on the right-hand side $\propto{\bf P}_{2}({\bf r}_{2};{\bf r}_{1})$. This is the effect of the dipolar field of the atom at the position ${\bf r}_{1}$ on the atom at ${\bf r}_{2}$. The equation for ${\bf P}_{2}({\bf r}_{2};{\bf r}_{1})$ has an analogous term, effect of the dipolar field of the atom at the position ${\bf r}_{2}$ on the atom at ${\bf r}_{1}$. This is our first glimpse of recurrent scattering, repeated photon exchange between a group of atoms; in this case, two atoms. The equation for the correlation function ${\bf P}_{3}({\bf r}_{1},{\bf r}_{2};{\bf r}_{3})$ similarly exhibits recurrent scattering between atoms at ${\bf r}_{1}$, ${\bf r}_{2}$ and ${\bf r}_{3}$, and so forth.

We also argued in Ref. Ruostekoski and Javanainen (1997b) essentially as follows: First, assume that the density equals a constant $\rho$, and all position correlations factor as $\rho_{n}=\rho^{n}$. Second, state an ansatz for all position-polarization correlations ${\bf P}_{k}({\bf r}_{1},\ldots,{\bf r}_{k-1};{\bf r}_{k})=\rho^{k-1}{\bf P}({\bf r}_{k})$. Third, ignore all recurrent-scattering cross terms such as the direct coupling of ${\bf P}_{2}({\bf r}_{2};{\bf r}_{1})$ to ${\bf P}_{2}({\bf r}_{1};{\bf r}_{2})$. Then ${\bf P}({\bf r})$ from Eq. (10) gives the exact solution to the entire hierarchy of the equations of the correlation functions. While this was not claimed in Ref. Ruostekoski and Javanainen (1997b), we believe that the same argument is valid even if the density is not constant as long as all density correlation functions factorize to a product of one-particle densities. Thus, there would be two possible reasons for beyond-MFT effect: Repeated exchange of photons between the atoms (recurrent scattering), and pre-existing correlations in the positions between the atoms.

Another mathematical point we made in Ref. Ruostekoski and Javanainen (1997b) about the cross terms runs as follows: Since the dipolar kernel ${\sf G}({\bf r}_{1},{\bf r}_{2})$ diverges for ${\bf r}_{1}={\bf r}_{2}$, it follows from (the steady-state version of) Eq. (6) that $P_{2}({\bf r}_{1};{\bf r}_{2})$ tends to zero when ${\bf r}_{1}$ and ${\bf r}_{2}$ tend to the same value, and this observation is independent of whether the delta function divergence is present in ${\sf G}({\bf r}_{1},{\bf r}_{2})$. Therefore the delta function divergence in ${\sf G}$, if any, should have no effect on ${\bf P}({\bf r}_{1})$ solved from Eq. (5). By an analogous argument, the delta function divergence has no effect on ${\bf P}_{2}({\bf r}_{1};{\bf r}_{2})$, and so forth. In short, any delta function divergence in ${\sf G}$ should have no effect on the polarization of the gas. This observation is, coincidentally, compatible with the physical notion that it is impossible to overlay two atomic dipoles exactly, so the delta function in ${\sf G}$ should never fire. Put mathematically, the actual position correlations $\rho_{k}$ for $k\geq 2$ should tend to zero when any two of the positions get close. One then surmises from the hierarchy that all correlation functions ${\bf P}_{k}$ for $k\geq 2$ have the same property.

Given the solution to the hierarchy of the equations of the correlation functions, one may obtain the expectation value of the electric field ${\bf E}({\bf r})=\langle{\hat{\bf E}({\bf r})}\rangle$ from Eq. (2),

[TABLE]

Notably absent from the formulation is the motion of the atoms. At the quantum level one may treat the center-of-mass motion of the atoms as a quantized degree of freedom, as one often does in the theory of the mechanical effects of light. We have not done so, however, as the ensuing theory would be cumbersome and opaque. The point to remember here is that the motion of the atoms, forces of light, and effects of photon recoil are all ignored in our present analysis.

II.2 Classical-electrodynamics solution for light propagation

Obviously it is in general impossible to solve the hierarchy for polarization correlation functions starting with Eqs. (5)-(7), and so on, directly numerically. The present subsection is mostly a brief summary of Refs. Javanainen et al. (1999); Lee et al. (2016) that present a workaround. The idea is that, just as one might solve a Fokker-Planck equation (diffusion equation) numerically using stochastic Langevin equations for individual particles Gardiner (2009), one may solve for the polarization correlations using stochastic classical-electrodynamics simulations.

Consider a gas of $N$ atoms. The simulations start with generation of random positions ${\bf X}_{1}$, …, ${\bf X}_{N}$ for the atoms in such a way that the positions are drawn from the probability distribution that gives the prescribed density correlations $\rho_{k}$. Incidentally, the normally ordered quantum correlations correspond to the properly defined classical density correlation functions for point-like particles from which the singularities corresponding to counting of the same particle more than once have been eliminated. For instance, the classical two-particle correlation function would be the stochastic average

[TABLE]

where the restriction $i\neq j$ removes the singularity.

Four cases of atomic probability distribution have been relevant in our work. First, we deal with classical atoms that are assumedly distributed completely independently of one another with the given density $\rho({\bf r})$. This strictly speaking cannot be true, as, for instance, atoms attract or repel each other at short distances, and the positions of the atoms may also be correlated on the length scale of the thermal de Broglie wavelength. However, we assume that the characteristic distance between the atoms is much larger than this type of correlations lengths. Second, maybe paradoxically, the same independent-atom model applies as the leading approximation also to the Bose-Einstein condensate. Third, atoms confined in optical lattices provide structured arrays where the positions of atoms in the Mott-insulator states can be sampled Jenkins and Ruostekoski (2012a); Yoo and Paik (2016). Fourth, we have also done one-dimensional simulations of a noninteracting one-component Fermi gas at zero temperature Javanainen et al. (1999); Ruostekoski and Javanainen (2016). The joint probability density for the positions of the atoms is then given by the absolute square of the many-body wavefunction, a Slater determinant. The characteristic feature of the Fermi-Dirac statistics is that the atoms tend to avoid each other, and end up with more evenly-spaced positions than classical atoms. Atoms with this position distribution may be sampled using the Metropolis algorithm.

Given the positions of the atoms, corresponding to the hierarchy for the correlation functions we next have the equations of motion for the dipole moments of the atoms. Denoting the positions of the dipoles explicitly, we have

[TABLE]

The dipole at ${\bf X}_{i}$ is driven by the external field ${\bf E}_{0}$, and by the dipolar fields from all other atoms. This may be seen even more graphically from the steady-state version of Eqs. (13),

[TABLE]

where

[TABLE]

is the polarizability of an atom; in fact, the well-known polarizability of the proverbial two-level atom at low light intensity. In the present case the steady-state dipole moment aligns with the net field at the position of the atom. This isotropy is because of our underlying assumption of the $J_{g}=0\rightarrow J_{e}=1$ transition. However, for other types of transitions Jenkins et al. (2016b) the appropriate polarization tensor $\alpha_{ij}$ could be defined such that the relation between the vector components of the the dipole moment and the electric field reads $d_{i}=\sum_{j}\alpha_{ij}E_{j}$.

In the present paper we deal solely with the steady-state version of the theory, Eqs. (14). These are a closed inhomogeneous set of linear equations for the dipole moments, or, thinking about it in another way, for the electric fields at the positions of the dipoles ${\bf E}({\bf X}_{i})$:

[TABLE]

Given a sample of the positions of the atoms, we solve these equations numerically for ${\bf E}({\bf X}_{i})$. Regarding the analogy to solving diffusion equations using particle simulations, this solution would be the counterpart of a stochastic trajectory obtained from the Langevin equation.

In the end we are interested in the total electric (and possibly also magnetic) field everywhere in space. It is formally given by

[TABLE]

everywhere except at the exact positions of the atoms, where we have a divergence in the dipolar kernel. In an indirect way, even this divergence has been taken into account: The dipolar field acting back on the atom that sends the field is formally infinite, but the action of this self-field is already included in the damping rate $\gamma$, and the associated level shift (Lamb shift) is incorporated into the energies of the levels. We also caution that the electric field as in Eq. (17) is not necessarily the most practical quantity to calculate. We return to this point below in Sec. IV.2.

The remaining step in the simulations is to repeat the process for a large number of samples of the atomic positions, and average the results. In the limit of an infinite number of samples, the stochastic average of, say, the electric field converges to the corresponding quantum mechanical average that would be obtained by solving the entire quantum hierarchy for the correlation functions.

Several analyses assuming that there is at most one photon present at any time  Prasad and Glauber (2000); Svidzinsky et al. (2010); Balik et al. (2013); T. Bienaimé and Kaiser (2013) also in effect show that in the limit of low light intensity quantum theory of light-matter interactions reduces to classical electrodynamics. There has been an argument along these lines that found some deviations from the standard EDPM Svidzinsky et al. (2015); specifically, a result that was traditionally thought to apply for the displacement was derived for the electric field. We emphasize, though, that by strictly following the Power-Wolley-Zienau procedure Power and Zienau (1959); Power (1964); Cohen-Tannaudji et al. (1989), with the inclusion of the polarization self-energy, we got results that were in a complete agreement with the structure of the usual EDPM Ruostekoski and Javanainen (1997a); compare Eq. (13) of Ref. Ruostekoski and Javanainen (1997a) and Eq. (14) of Ref. Svidzinsky et al. (2015) with Eq. (31) of Ref. Svidzinsky et al. (2015). There are also caveats to the quantum-classical agreement, cases such as 1D nanofibers or photonic crystals Thompson et al. (2013); Douglas et al. (2015) in which the atom-field coupling can be so strong that one photon may saturate an atom. We will not analyze such situations any further in the present paper.

From this point on we leave quantum mechanics behind, and pretend that the electrodynamics of the dipolar medium is, in fact, entirely classical. For instance, for any given atomic sample we find the electric field which we may average over many samples to get the averaged field. Likewise, given the electric field, we square it, average, and obtain the average of the square of the electric field, which is an intensity-like quantity. This is not a trivial point in quantum mechanics: If we were to compute the quantum average of $\hat{\bf E}^{\dagger}({\bf r})\cdot\hat{\bf E}({\bf r})$ exactly from quantum mechanics, we would in principle have to do something like develop a hierarchy of equations for some other correlation functions than those we have dealt with so far, and, if possible, develop a corresponding simulation. The difference between the quantum mechanical average of the square of the quantum field and the classical average of the square of the classical field is in the quantum fluctuations.

Nonetheless, we do not expect significant quantum fluctuations in the kind of situations we consider here. In fact, for a model atom such as ours, in the limit of low light intensity, there are no quantum fluctuations in the scattered light. We may have an issue for instance if the light intensity is increased, or if more than one electronic Zeeman ground level is involved Lee et al. (2016); Hammerer et al. (2010), although observing such quantum fluctuations typically require delicate experimental setups.

II.3 Units and conventions

All numerical computations described here were done in units such that the numerical values

[TABLE]

apply. Here $k=\omega/c$ is the wave number of the driving light. For parameters typical in laser spectroscopy the difference between the wave number of the driving light and of resonant light $k_{0}=\omega_{0}/c$ is negligible, and we henceforth ignore it. The unit of length $k^{-1}$ is related to the wavelength of the driving light by $k^{-1}=\lambdabar=\lambda/2\pi$, and the unit of quantities such as area and density follow accordingly. Below all discussions are in these units, unless explicitly stated otherwise.

For a dipole ${\bf d}$ at ${\bf r}_{0}$, the electric and magnetic fields of dipole radiation at position $\bf r$ are

[TABLE]

where ${\sf G}$ is again the dipolar field propagator, and $\sf H$ gives the magnetic field from a dipole. Expressed in Cartesian coordinates, these are matrices with the components

[TABLE]

Here $r$ and $\hat{\bf n}$ are the distance from the source point to the field point and the unit vector directed from the source point to the field point, and $\hat{\bf e}_{i}$ are the cartesian unit vectors. We have dropped the contact term in (20) as we here observe the light outside the sample, and in the interactions between the atomic dipoles it is inconsequential Ruostekoski and Javanainen (1997b). The relation between the positive frequency parts of the electric and magnetic fields reads

[TABLE]

and the energy density and Poynting vector at the given field position are

[TABLE]

Another convention here is that, by default, we express the detuning in units of the linewidth of the transition,

[TABLE]

with the dimensionless detuning $\delta$. Be virtue of Eq. (8) and our conventions, the relation between dipole moment matrix element and linewidth reads

[TABLE]

and the polarizability (15) is

[TABLE]

Assuming a single atom at the origin and an incoming field with the vector amplitude ${\bf E}_{0}$, one may straightforwardly obtain the dipole moment, the electric and magnetic fields, the Poynting vector, and finally the radiated power. The result, a useful reference, is

[TABLE]

The intensity of the incoming plane wave is

[TABLE]

Writing the radiated power in terms of the intensity and the scattering cross section $\sigma$ as $P=\sigma I_{0}$, we find

[TABLE]

The on-resonance light scattering cross section therefore is $\sigma(0)=6\pi$.

III Simple examples

In this section we discuss independent-atom response and cooperative response to light in simple analytically solvable cases. Among other things we demonstrate the difference between coherent and incoherent scattering, and show how even the seemingly starkly contrasting concepts of cooperativity from Dicke states and radiation from independent atoms may be difficult to tell apart. Our examples serve as a general reference for more comprehensive numerical-simulation studies that follow. We also introduce several pieces of physics that are absent from the present simplest models, but might well figure in real experiments and complicate the interpretation of simulation results.

III.1 Radiation from a Gaussian clouds of atoms

III.1.1 Continuous medium

For the problem of $N$-atom gases, let us start with a hypothetical model with a continuous spatial distribution of atoms. In terms of macroscopic electromagnetism, there is a monochromatic polarization of the sample ${\bf P}({\bf r})=\rho({\bf r})\,{\bf d}({\bf r})$, where $\rho({\bf r})$ is the density of the sample and $\bf d({\bf r})$ is the electric dipole of an atom at the position $\bf r$. Taking the atoms to reside around the origin of the coordinates, in the far field at the distance $r\gg 1$ and with $r$ much larger than the size of the sample, the terms $\propto 1/r^{2}$ and $\propto 1/r^{3}$ in Eq. (20) are negligible and the field radiated (“scattered”) by this polarization is

[TABLE]

$\hat{\bf r}={\bf r}/r$ is the unit vector that points from the source at $\simeq 0$ toward the field point at the distance $r$.

For easy analysis, we model the density with a Gaussian,

[TABLE]

where $N$ is the atom number and $R$ is the size scale of the sample. The parametrization is chosen in such a way that the rms value of $|{\bf r}|$ equals $R$. In the limit $R\gg 1$ there will be a narrow cone of radiation around the direction of the incoming beam; let us denote the angle from the incident beam by $\theta$.

For a tangible example we take a $\sigma_{+}$ circularly polarized plane wave propagating in the $z$ direction, writing

[TABLE]

The assumption is that the incoming light dominates even inside the sample, i.e., that each atom responds to the incoming light only. Accordingly, we write the dipole moment of an atom at ${\bf r}$ as

[TABLE]

The radiated field from Eq. (31) is then

[TABLE]

In the far field the light locally makes a plane wave, the Poynting vector points radially outwards and has the magnitude

[TABLE]

and the total power in the radiation is readily obtained as

[TABLE]

The intensity of the scattered light scales with the square of the atom number, $N^{2}$. This is similar to an important characteristic of superradiance. However, it is early for conclusions yet. Instead, we will next inspect a more realistic problem with discrete atoms.

III.1.2 Independent discrete radiators

Take a collection of $N$ identical dipoles sitting at the positions ${\bf r}_{i}$ in the incoming field, and assume that each of these dipoles radiates a field ${\bf E}_{i}({\bf r})$ independently. In other words, we again assume that only the incoming field ${\bf E}_{0}({\bf r})$ drives each dipole. In terms of scattering theory, one might say that a photon scatters from an atom at most once, so this model is occasionally called single-scattering approximation.

The total dipolar field at the point ${\bf r}$ is

[TABLE]

with

[TABLE]

In the far field the radial component of the Poynting vector is

[TABLE]

We take the position of each atom to be a random variable independent of the positions of the other atoms, governed by the probability density function $f({\bf r})$. Then the average outward energy flux (average over many samples of the gas) is determined from

[TABLE]

The first term represents coherent scattering, as if the atom was spread out to a continuous dielectric material with the spatial shape specified by $f({\bf r})$. It arises from adding the fields of different radiators, and is essentially proportional to $N^{2}$. The second term $\propto N$ is for incoherent scattering, the sum of the intensities radiated by the individual atoms. It is present because the gas is not a continuous (nonfluctuating) dielectric medium, but consists of discrete scatterers.

Incidentally, while the above argument might not be as widely known as it deserves to be, the basic message is far from novel. If air were a continuous dielectric medium, it would not scatter sunlight sideways and the sky would be black. The blue sky comes from incoherent scattering that results because air consists of discrete molecules. This observation goes back to (at least) Lord Rayleigh Lord Rayleigh (1899).

For a comparison, we apply the same incoming light as in Eq. (33), and the position distribution for each atom is taken to be the same Gaussian,

[TABLE]

Given the usual polarizability $\alpha$, the far field (the $1/r$ part of dipole radiation) averaged over the positions of an atom, the absolute square of the former, and the absolute square of the field averaged over the positions give

[TABLE]

and the total radiated power becomes

[TABLE]

The part $\propto N(N-1)$ in Eq. (46) is the same as it would be for the continuous atom density in Eq. (32), except for the factor of $N(N-1)$ instead of $N^{2}$. In the forward direction $\theta=0$, the intensities of coherently and incoherently scattered components of light add up to exactly $N^{2}$ times the intensity from a single atom. The bigger is the sample, the narrower is the cone in the direction $\theta\simeq 0$ for coherent scattering. This aspect is demonstrated in Fig. 1. The sample acts as an antenna that directs the radiation in the forward direction. The total power of scattered light decreases with increasing size of the cloud.

Conversely, in the limit $R\rightarrow 0$ the intensity in all directions is enhanced by the factor $N^{2}$, and of course so is the total power. Given the well-known Dicke cooperative regime, a reader might erroneously interpret such an enhancement as a cooperative phenomenon. It cannot be, since in this example we have simply added the fields from independent radiators.

III.2 Radiation from two atoms

For two atoms (or ions, as things might be) the radiation field can be solved explicitly, and there have even been experiments already a while ago Eichmann et al. (1993); DeVoe and Brewer (1996). We discuss as an example the special case when a plane wave polarized in the $x$ direction and propagating in the $z$ direction strikes two atoms sitting on the $x$ axis separated by the distance $\ell$. Specifically, we have the incoming field and the two positions for the dipoles

[TABLE]

In this case the fields at the positions of the dipoles as solved from Eqs. (16) are

[TABLE]

Since the dipolar field diverges with decreasing distance from the dipole, one might expect that the fields at the positions of the dipoles should diverge when the dipoles approach one another. However, the exact opposite holds true: For a fixed detuning and hence fixed polarizability $\alpha$, ${\bf E}({\bf r}_{\pm})$ actually tend to zero as $\ell^{3}$ when the distance $\ell$ between the dipoles tends to zero. Maybe counterintuitively, when the detuning is kept constant and the atoms approach each other, they decouple from the light altogether Jenkins et al. (2016b). That is why we are not overly concerned about some atoms being close to one another in the steady-state numerical simulations.

Fixed detuning, however, may not be the most useful way of viewing the result. Instead, we insert the explicit expression of the polarization. In this subsection III.2 only, we find it expedient not to scale the detuning to the linewidth $\gamma$, and write

[TABLE]

This shows our first instance of cooperative shift and broadening of the resonance of the atoms as a result of the radiation from one dipole falling on the other. The sines and cosines originate from retardation, propagation delay of light between the atoms. In the limit $\ell\rightarrow 0$ we have the expansions, keeping the leading terms,

[TABLE]

The shift $\Delta(\ell)-\Delta$ diverges as $\ell^{-3}$, which clearly reflects the dipole-dipole interactions between the atoms. It is this shift that leads to the decoupling of two closely spaced atoms from the light. On the other hand, the linewidth doubles.

Moving on to the energy flux in the far field and to the radiated power, we find after some tedious mathematics the expressions

[TABLE]

The angular distribution of the radiation, normalized in such a way that $\int_{0}^{\pi}d\theta\,\sin\theta\,P(\theta)=1$, reads

[TABLE]

This shows the dipole radiation pattern modulated by the interference of the radiation from the two dipoles; see the demonstration in Fig. 2. For $\ell=0$ we have the usual dipole radiation. With increasing $\ell$ the interference first concentrates the radiation more to the $\theta=\pi$ plane perpendicular to the dipoles. With increasing $\ell$ the side lobes grow numerous, and the overall angular distribution pattern rounds out.

As to the radiated power, the difference of the laser frequency from the atomic resonance shifted by the cooperative effects is the true gauge of the detuning. We momentarily assume that this shift of reference point is implicit in the expression of the power in the two-atom radiation $P_{S}$, and simply replace $\Delta(\ell)\rightarrow\Delta$. On the other hand, if there were no cooperative effects or interference between the radiations from the two dipoles, the total power would be twice the power radiated by one dipole under the same driving field. The ratio of the actual two-atom power and the power from two independent dipoles

[TABLE]

is a quantitative measure for the effects of the presence of two dipoles. In the limit $\ell\rightarrow\infty$, $\gamma(\ell)\rightarrow\gamma$, and we have $C\rightarrow 1$, as expected; the radiated powers from the two dipoles simply add.

The situation is more intriguing in the opposite case $\ell\rightarrow 0$ with $\gamma(\ell)\rightarrow 2\gamma$, and the nature of the result depends in an interesting manner on the detuning. With $|\Delta|\gg\gamma$, we have $C\simeq 2$, and the two-atom sample radiates twice as much power as two separate atoms would. If $|\Delta|\ll\gamma$, we have C\simeq{\hbox{\frac{1}{2}}}. This means that on resonance the two atoms together emit the same power as one atom would. In fact, in the limit $\ell\rightarrow 0$ the radiated power can be expressed at all detunings as

[TABLE]

as if we had a single dipole with the dipole moment matrix element equal to $\sqrt{2}$ times the original dipole moment matrix element, hence the linewidth $2\gamma$. Now, if the two dipoles were at the same place but completely independent, the total induced dipole moment would be twice the dipole moment induced on one atom. But the induced dipole moment is proportional to the square of the dipole moment matrix element, hence the multiplier $\sqrt{2}$ in the dipole moment matrix element is consistent with the notion of independently radiating atoms. There is a similar $\sqrt{2}$ in the quantum mechanics of the Dicke states. From the present angle, this factor is classical physics in disguise.

Sufficiently far off resonance the dipoles, even if close to one another, are independent, and each radiates the same amplitude as one dipole would. This means twice the amplitude and four times the power, which is the result we already noted. This is again an interference effect, and has nothing to do with cooperativity. Cooperativity is clearly responsible for the shift of the resonance.

Regarding the modifications of the linewidth, there is some ambiguity. As is well known, on resonance the power radiated from a two-level atom is independent of the dipole matrix element. This may be thought of as a consequence of energy conservation: With increasing dipole moment the atom tends to radiate more, but at the same time the increased radiation damps the resonant response more and these effects exactly balance. Viewed in this way, on resonance even two independent atoms should radiate the same power as one atom, and we might call this an interference effect; four times as much radiation, but also four times as much damping. Nonetheless, we may regard the resonance behavior as cooperative as well. Light from both atoms falls back on both of them, and we have a cooperative radiation reaction that determines the altered damping rate and linewidth. For one thing, the variation of the linewidth $\gamma(\ell)$ with $\ell$ shows that the propagation of light from one atom to the other is involved. The limit $\gamma(\ell)=2\gamma$ for $\ell\rightarrow 0$ appears to allow one to think of the resonance linewidth both as an independent-atom phenomenon and as a cooperative phenomenon. This is a remarkable coincidence, if a coincidence it is.

III.3 What’s missing?

There are several obvious pieces of physics missing from our picture that may figure in the interpretation of the experiments and simulations alike. We mention a few most notable items here, and amplify as we go along. An extended and somewhat complementary account is given in Ref. Guerin et al. (2016b).

First, there is the interference of the scattered light with the incoming light. This is behind the “absorption” of light. Atomic samples that scatter light elastically return all of the light energy back to the light field, and there is no genuine absorption. Instead, the light from the incident driving field, say, a laser beam, and the forward-scattered light interfere destructively. The energy that gets removed from the incident beam is directed elsewhere.

Second, suppose we actually did have a continuous and nonfluctuating distribution of polarization ${\bf P}({\bf r})$. A typical microscopic model would state that for a continuous density $\rho({\bf r})$ the polarization is ${\bf P}({\bf r})=\rho({\bf r}){\bf d}({\bf r})$ if the dipole moment of an atom at ${\bf r}$ were ${\bf d}({\bf r})$. The standard method to analyze this situation is to use the EDPM. It is not an independent-atom theory, but takes into account the effects of the radiation from the atoms on each other in some averaged way. In fact, EDPM is a MFT. Whether one can solve it accurately is another matter, but one can apply general intuition. For instance, a Gaussian cloud might act like a (poor-quality) converging or diverging lens depending on the sign of the dielectric constant, which in turn depends on the sign of the detuning. EDPM is not an exact and possibly not even a quantitatively useful description of the response of an atomic sample to light, but optics-like effects should be expected to be present in the results of both experiments and simulations.

Ordinarily, when one thinks of light propagation through a sample as a standard optics problem, the initial and scattered field are dealt with at the same time. There is the remarkable Ewald-Oseen extinction theorem Born and Wolf (1999), which roughly says that inside a dielectric medium the electric field has a component that cancels the incoming field. That is why the light inside a dielectric medium has the wavelength $\lambda/n$ appropriate for the dielectric constant of the medium $n$, even if from a microscopic standpoint one also concludes that the field inside is the sum of the incoming field and the field scattered from the atoms.

Finally, suppose one solves the light propagation problem numerically using classical-electrodynamics simulations. Since the positions of the atoms are random, there are fluctuations in the scattered radiation, and perforce, in the interference of the scattered and incoming radiation. Smooth radiation patterns as in Fig. 1 are averages over a large number of atomic distributions, but at least over time scales such that the atoms may be regarded as being at standstill there is no such averaging. An individual sample of atomic positions may give a radiation pattern that looks quite ragged. Eventually one has to confront the possibilities of spatial fluctuations in the scattered light, incoherent scattering, and generalizations of incoherent scattering beyond the single-scattering framework.

IV The slab

In this paper we study mostly a slab of matter, with the light coming in to a face of the slab at normal incidence. In this case EDPM may be solved exactly in what amounts to a student exercise. The idea is to compare these exact solutions of the MFT with ab-initio numerical simulations.

IV.1 Elementary optics

Figure 3 illustrates the standard-optics problem. The light is coming in from vacuum with the refractive index 1. We denote the refractive index of the medium by $n$, so the wave number inside is $K=n$ ($nk$ in SI units). The incoming light with the reference electric field amplitude at the entrance $E_{0}$ gets either reflected or transmitted at the entrance face, with the corresponding amplitude reflection and transmission coefficients being $(1-n)/(n+1)$ and $2/(n+1)$ Jackson (1999); Born and Wolf (1999). Inside the medium we have two amplitudes, $E_{r}$ corresponding to the right-going wave with the propagation factor $e^{iKz}$, and the left-going amplitude $E_{l}$. By matching the incoming, reflected and transmitted waves at the front face we have

[TABLE]

Similar matching can be made at the exit face, which leads to the relation between the incoming and transmitted amplitudes

[TABLE]

To complete the exercise we note that, according to the local-field corrections the effective electric field inside the sample is $E_{e}=E+\frac{4\pi}{3}P$ Jackson (1999); Born and Wolf (1999), the polarization is $P=4\pi\chi E$, where $\chi$ is the susceptibility, and also $P=\rho\alpha E_{e}$, where $\rho$ is the atom density and $\alpha$ the polarizability (27) of an atom. We then have

[TABLE]

where $\delta_{L}=-2\pi\rho$ is the Lorentz-Lorenz (LL) shift of the resonance. Simple algebra gives the power transmission coefficient, optical thickness (depth, density), and the conventional absorption coefficient defined as

[TABLE]

Standard scattering theory says that if light propagates in a medium with density $\rho$ for a distance $h$, for the scattering cross section $\sigma$ the fraction of light energy that makes it through and the corresponding optical thickness are

[TABLE]

This is Beer’s law, and the reason why we usually state our results in terms of optical thickness: If Beer’s law were valid, the line shape of the optical thickness, e.g., its variation with the tuning of the driving light, would be independent of the thickness of the sample. The physical thickness would simply be a multiplicative factor.

Clearly, the EDPM solution (58)-(60) cannot agree with Beer’s law exactly. There are interesting lessons to be learned from this discrepancy. Suppose we have an electric field propagating in the $z$ direction, of the form $E(z)={\cal E}(z)e^{iz}$, where ${\cal E}(z)$ varies little with $z$ over the scale length 1; basically, over the scale of the wavelength. The so-called slowly varying envelope approximation would then say

[TABLE]

By assuming a dominant propagation direction and a slowly varying electric field amplitude, the wave-equation for the electric field is converted to a first-order differential equation, and Beer’s law follows. What gives in a slab is that the dielectric medium is taken to have an abrupt face, so the electric field cannot be expected to vary slowly over the length scale of a wavelength everywhere. Indeed, there is a counterpropagating reflected wave inside the slab, and not just a single direction of propagation. In physical terms, the problems with Beer’s law can be attributed to etalon effects, reflections of light from the faces of the dielectric slab.

When we deal with samples that have relevant features of the size of a wavelength like the slab Keaveney et al. (2012), or that are of a size comparable to a wavelength Jenkins et al. (2016a); Jennewein et al. (2016), the usual approximations of optics such as slowly-varying envelope approximation and paraxial approximation tend to break down, not to mention ray optics. One then has to solve the full Maxwell equations, an onerous requirement even numerically Davidson (2011); Wriedt (2012); Kahnert (2016); Jennewein et al. (2016). For the slab the standard optics gives an exact solution to Maxwell’s equations. The standard-optics result for a slab may also be derived straightforwardly Ruostekoski and Javanainen (1997b) from Eq. (10).

It should be noted that we have inserted the local-field correction Jackson (1999); Born and Wolf (1999) by hand. Given the specific form of the polarizability, the result is then exactly the LL redshift of the resonance as a function of the density of the sample, $\delta_{L}=-2\pi\rho$. In the standard SI units it would read $\Delta_{L}=-2\pi(\rho/k^{3})\gamma$.

What is the shift of the resonance is a difficult question operationally. Beer’s law (61) says that the line shape of optical thickness is a Lorentzian and the position of the resonance can be easily determined. However, in general we find a line shape that is not Lorentzian, and worse, not symmetric about any particular tuning of the driving light. Keeping this in mind, we next discuss the line shift for the EDPM solution of the slab (58)-(60).

Let us assume that the line shape is of the form $D=\rho/(K_{0}+K_{1}\rho)$, where $K_{0}$ and $K_{1}$ are independent of the density $\rho$. We expand this form in density $\rho$, also expand the EDPM optical thickness from Eqs. (58)-(60) in $\rho$, and choose the coefficients $K_{0}$ and $K_{1}$ in such a way that up to second order in $\rho$ we have the same expansions. We find that, up to this order in $\rho$, the line shape is Lorentzian, and is shifted from the one-atom resonance by

[TABLE]

This is the “cooperative Lamb shift” of Friedberg, Hartmann and Manassah Friedberg et al. (1973). Here the first term is the LL shift that we put in by hand, and the second, oscillatory, part comes from the etalon effects. The theoretical result (63) was recently tested experimentally, albeit in a hot gas with moving atoms Keaveney et al. (2012), and found to work quite well, apart from a shift between the theory and the observations proportional to the density of the gas. In a cold and dense trapped cloud of atoms an analogous expression was shown to fail Jennewein et al. (2016), but in a very dilute limit analogous expressions, derived from the standard optics, are expected to provide qualitative estimates for the shift even for cold atoms Roof et al. (2016).

For an arbitrary density we use the maximum of the resonance line as a proxy of the position of the resonance, hence as the shift from the one-atom resonance. The results obtained numerically from Eqs. (58)-(60) are plotted in Fig. 4. We show the shift in units of the (absolute value) of the LL shift as a function of sample thickness for various densities as solid red lines, and the low-density limit (63) as the dashed black line. At low atom densities the oscillatory etalon effect of the “cooperative Lamb shift” is clearly observable. However, as the density increases, the sample becomes optically thicker and the fraction of the light that propagates from one face to the other decreases. Etalon effects are reduced, and the MFT resonance shift tends to the LL shift.

IV.2 Numerical simulations

A slab extending to infinity in the transverse directions would correspond to an infinite number of atoms, an impossibility for numerical analysis. In our simulations we attempt to do the next best thing: We study an atomic sample confined to a circular disk of radius $R$ and thickness $h$, and make the radius as large as practicable. The area of the disk is then $A=\pi R^{2}$. We put some given number $N$ atoms evenly distributed inside the disk, which gives the number density $\rho=N/(hA)$. We again assume a plane wave of light propagating along the axis of the disk, and denote the direction of propagation by $z$. In keeping with the symmetry of the disk, we take the incoming light to be circularly polarized, so that it again reads ${\bf E}_{0}({\bf r})=E_{0}\,e^{iz}\,\hat{\bf e}_{+}$.

Simulating the scattered field for an individual sample of atoms and averaging over the samples is straightforward per se. In contrast to the standard-optics solution of Sec. IV.1, we do not put in any ad hoc local field corrections or LL shifts. As far as the microscopic model of the dipolar medium is concerned, they simply do not belong there.

Here we mostly discuss the transmission of light through the sample. Unfortunately, for the reasons we already touched upon in Sec. III.3, the general simulation scheme we have described previously would not work satisfactorily. The problem is graphically illustrated in Fig. 5. Here we have a disk with thickness $h=1$ and area $A=1024$, and with $N=2048$ atoms inside in random positions. The figure shows the ratio of the transmitted intensity $I=\frac{1}{8\pi}|{\bf E}_{T}|^{2}$ to the incoming intensity $I_{0}$ in a plane parallel to the disk and at the distance $10\pi$ (five wavelengths) downstream from the center of the disk. One can see the depression at the center, the shadow cast by the disk, but also diffraction rings and fluctuations of the intensity as a function of position. Except for very low atom densities, we cannot handle numerically large enough disks to materially eliminate the diffraction, which would have a large effect on the computed transmission. Besides, the spatial fluctuations of intensity will be present regardless.

There is, however, a shortcut Chomaz et al. (2012) that appears to expedite the approach to the limit of large area of the disk enormously. We compute the transmitted intensity outside of the disk as if the atoms only scattered light in the forward direction, and hence call this the forward-scattering approximation. To begin with, let us take a dipole in the $xy$ plane in a disk with radius $R$ centered at the origin, and an observation point in the far field with $x=0$, $y=0$, $z=\xi$ and $\xi\geq 1$. We also assume that $R\gg\xi$. In the argument we take the dipole to be polarized in the $x$ direction, ${\bf d}=d\,\hat{\bf e}_{x}$, and have it reside somewhere in the disk with the coordinates $\{\varrho\cos\phi,\varrho\sin\phi,0\}$.

In fact, we take the position of the dipole to be random and evenly distributed inside the disk, so that the average of the field at the observation point is

[TABLE]

Here ${\sf G}_{F}$ retains only the far-field contribution $\propto 1/r$ to the dipole field propagator ${\sf G}$ of Eq. (20). We first carry out the angular integral, and in the remaining integral over $\varrho$ make the substitution $x=\sqrt{\xi^{2}+\varrho^{2}}$. This results in the electric field

[TABLE]

This boils down to two integrals. We first have

[TABLE]

This estimate says that for $R\gg\xi$ the value of the second term oscillates rapidly with a large $R$, and we use its average value 0. The other integral we approximate as

[TABLE]

where the second form is found numerically for the limit $\xi\gg 1$. The same argument could just as well be made for the $y$ polarization of the dipole, so for a dipole $\bf d$ in the $xy$ plane we simply have

[TABLE]

Here $e^{i\xi}$ is a phase-matching factor as dictated by the driving plane wave of light.

Average over the circle with the radius $R$ eliminates the longitudinal component of the dipolar field, but it may be present in the field of a dipole that is not on the axis of the circle. In that case we remove the longitudinal component by hand. Given the incoming plane wave of light ${\bf E}_{0}e^{iz}$, we therefore write the total transmitted light as a sum over the dipoles at their positions ${\bf r}_{k}$ as

[TABLE]

This is the same prescription as given in Chomaz et al. (2012), albeit in our system of units.

We have compared the transmission coefficient calculated from Eq. (69) with the analytically known result for one atom in the disk obtained from scattering theory. Suppose the light is on resonance so that the scattering cross section is $6\pi$, then to the leading order the analytical approximation of the absorption coefficient is ${\cal A}=1-6\pi/A$. For a disk with the area $A=256$, the difference from this prediction and the forward-scattering approximation is about $2\%$, and the difference decreases inversely proportionally to the area of the disk. In the limit of a dilute sample the forward-scattering approximation also reproduces the MFT results, the main difference being in the shift of the resonance. We see this even if the disk is so thick that most of the light is absorbed. Given Fig. 5 and the ugly approximations in the derivation, Eq. (69) reproduces the MFT results amazingly well when MFT is expected to be valid.

We have interspersed test cases among our simulations where we have increased the area of the disk and looked for convergence of the transmission coefficient. To verify the convergence has proven to be exceedingly expensive in computer time, but we obtain order-of-magnitude estimates of the simulation error due to the finite size of the disk. We occasionally quote them with our results. These truncation errors are usually the largest known numerical errors in our computations, surpassing the statistical fluctuations that result from the necessarily finite number of samples used in the averaging over the atomic positions.

The forward-scattering approximation together with the increasing size of the disk can evidently be used to mitigate the complications due to the optics of the finite-size disk. However, the forward-scattering approximation also removes the spatial fluctuations from the transmitted light. To quantify the fluctuations we momentarily discuss the reradiated dipolar field only, the sum on the right-hand side of Eq. (17). We imagine placing a probe disk of the same radius $R$ as the simulation sample, at the distance $R^{2}/2$ downstream from the atomic sample. This distance is analogous to the Rayleigh range, where the light radiated by the atomic sample starts to transition from the near-field form of a beam of light to the far-field form of a cone with a constant opening angle, and as such gives a natural place where to observe the scattered field. We integrate the component of the electric field with the same polarization $\hat{\bf e}_{+}$ as the incoming beam over the probe disk, and study fluctuations of the integrated field over the atomic samples. For a disk with $h=1,\,A=1024,\,N=1024$, and on resonance, the fluctuations are about 6%.

For an infinitely large disk absorption would correspond to the interference of the incoming and scattered light. In the case when the optics of the disk has a significant effect, the interference between the scattered field and a field with the wavefront matched to the diffraction pattern of a disk-shape aperture might lead to a more meaningful measure of absorption, but the diffraction pattern is difficult to calculate accurately as this would require solving full vectorial Maxwell’s equations. Instead, our studies of fluctuations give us an indirect estimate of the limitations of the calculations of the absorption coefficient as we have done them in this paper: No matter what wave front, in our particular example a random residual field of about $6/100$ of the original field amplitude must remain after the incoming and scattered fields have canceled each other to the maximum extent allowed by optics. A fraction of the intensity of about $(6/100)^{2}$ invariably gets through as a result of the fluctuations. Such fluctuations are not present in EDPM (here we always assume that EDPM refers to a static, continuous medium). We might conceivably think of the transmission of the random component of the field as diffusion of light through a sample of randomly spaced scatterers Chandrasekhar (1960); Ishimaru (1997); van Rossum and Nieuwenhuizen (1999). In our example a meaningful comparison with the MFT would only be possibly at optical thicknesses of less than $D\sim-\ln(6/100)^{2}\sim 6$. From this kind of arguments we surmise that diffusion of light does not materially affect the conclusions of the present paper.

V Results for the slab

Our main qualitative result is displayed in Fig. 6 that compares the optical thickness from the MFT and from numerical simulations for a dense ($\rho=1$) slab. Even if one discounts the LL shift that was put in by hand to the MFT anyway, there is a large difference between the curves. Our interpretation is that the standard EDPM fails. This is not hard to fathom: EPMD is a MFT for the light-mediated interactions between the dipoles, and as garden variety MFTs do, it goes bad when the interactions between the atoms increase; here, as a result of density.

Our next question is, what kind of densities are needed for notable discrepancies between the MFT and the simulations? For a demonstration, we have developed the following scheme. We first fix the thickness of the slab at $h=\pi$. For various values of density $\rho$ we then find from the MFT exactly, numerically, the corresponding positive detuning $\delta(\rho)$ for which the absorption coefficient is ${\cal A}_{M}=0.01$. There is a corresponding negative detuning such that the absorption is the same, not precisely the negative of the positive detuning because the line shape is not symmetric about $\delta=0$, but we arbitrarily pick the positive detunings. We also compute the absorption coefficients ${\cal A}_{C}$ from the simulations for the same densities $\rho$ and detunings $\delta(\rho)$. In Fig. 7 we plot the ratio ${\cal A}_{C}/{\cal A}_{M}$, the cooperative enhancement of absorption, for a number of sample densities. The enhancement exceeds two already for a sample as dilute as $\rho=0.1$, with $\delta(\rho)\,=\,24.2\,$. The message here is that, depending on what the experiment might be, significant deviations from MFT may be found at unexpectedly low densities even for the off-resonance case.

The results of Fig. 7 also tally with the simulation results in Fig. 6. The detunings are large in Fig. 7, about $\delta\sim 80$ for the case of $\rho=1$, and it is evident from Fig. 6 that this far in the wings the MFT understates optical thickness and absorption.

There are aspect in Fig. 7 that also bear on the interpretation of our simulations. The overall absorption in the figure varies between 0.01 and 0.1 depending on density, so that between 90% and 99% of the light gets through. The sample is optically thin, the thickness being at most 1/10 of the “mean free path” of a photon. The deviations from MFT evidently cannot be attributed to diffuse scattering or radiation trapping Chandrasekhar (1960); Ishimaru (1997); van Rossum and Nieuwenhuizen (1999); Bourgain et al. (2013); Guerin et al. (2016b). The area of the disk was also large, approximately $A=2000$, and we looked at the convergence of the results with the disk area extensively. Even at its worst, at the highest density $\rho=1$ in the figure, the relative error in the results should be no more than on the order of 5%. The finite size of the disk should not be a major contributor to the deviations of the simulations from MFT either.

The sample in Fig. 7 is optically thin because of the large detuning. For the MFT, Eqs. (58)-(60), an expansion in $1/\delta$ appropriate to the limit of a large detuning $\delta$ gives

[TABLE]

The factor in front is the absorption coefficient of independent atoms far-off resonance when the cross section for photon scattering $\simeq{6\pi}/\delta^{2}$ is asymptotically small. However, even in MFT and far-off resonance, the sample does not behave like a collection of independent radiators. The reason is etalon effects, reflections of light from the faces of the slab: For a large detuning the refractive index $n$ is close to one and deviations from unity scale like $\rho/\delta$. A fraction of light $\simeq 1-n\propto\rho/\delta$ that made it through the slab, almost all of it for a large $\delta$, is reflected from the back face, almost all of the reflected light propagates to the front face, gets reflected again with a reflection coefficient $\propto\rho/\delta$, and finally interferes with the light that goes straight through. That is the reason for the two terms in Eq. (70) proportional to $(\rho h/\delta^{2})$ (light straight through) and $(\rho/\delta)^{2}$ (interference).

Aside from pointing out that run-of-the-mill optics can pop up in quite unexpected places, we use Fig. 7 and Eq. (70) as a springboard for dimensional analysis. Our hypothesis is that there are two in principle independent dimensionless density parameters in the light-propagation problem, the on-resonance optical thickness $6\pi\rho h$ ($6\pi\rho hk^{-2}$ in terms of full dimensional quantities) that is a MFT parameter, and the plain $\rho$ ($\rho k^{-3}$) that governs the role of dipole-dipole interactions beyond MFT. In several recent experiments Kwong et al. (2015); Guerin et al. (2016a); Bromley et al. (2016); Araújo et al. (2016) the on-resonance optical thickness has proven to be the dimensionless parameter that governs the density dependence of the results. In contrast, our interpretation is that the nontrivial results in Fig. 7 are attributed to the beyond-MFT parameter $\rho$.

Figure 7 does not separate the two dimensionless quantities cleanly as we kept the MFT absorption coefficient constant, not the on-resonance optical thickness. It also seems that optics, in this case due to the reflection of light from the surfaces of the slab, inevitably imposes some ambiguity in the interpretation of the results. Here we attempted to minimize the effects of the optics by choosing the sample thickness $h=\pi$, whereupon the etalon-effect term in Eq. (70) vanishes.

We next describe a numerical experiment in which we literally keep the on-resonance optical thickness constant and vary the density. The basic idea is to keep the area of the disk ($A=4096$) and atom number ($N=512$) constant while varying the thickness $h$. The resulting optical thickness as a function of detuning from the independent-atom scattering theory (30) would be

[TABLE]

$D=2.4$ on resonance with $\delta=0$, which means less than 10% power transmission. In practice we had to deviate from this ideal to keep the error due to the finite area of the disk somewhat under control, so that at the high end of the thickness range we finished with $N=2048$, $A=16384$ and $h=16$, but of course without altering the preset optical thickness. The resulting absorption line shapes, optical thickness vs. detuning, are shown in Fig. 8 for a range of thicknesses varying from $h=0.25$ to $h=16$ in multiples of two (solid red lines from bottom to top), corresponding to the densities ranging from $\rho=0.25$ to $\rho=0.0078125$ decreasing by factors of $0.5$. Also shown is the corresponding independent-atom prediction (71) (dashed black line), which basically differs from the lowest-density simulation graph by a very small shift of the resonance frequency.

As before, there are etalon effects that affect the results, although the optical thickness should depress them some as the light reflected from the back face of the slab should complete a back-and-forth trip before interfering with the light that gets through on the first try. This caveat notwithstanding, we see substantial effects of the varying density on the absorption line shapes. The dipole-dipole interactions make a noticeable difference already at the density $\rho=0.015625$, the second red curve from the top.

There are two other observations to be made here. First, at higher densities the resonance shifts as large as in Fig. 4 should be plain visible to the eye in Fig. 8, and the resonances should move to the red. There are visible density shifts of the resonance alright, but much smaller than one would surmise from Fig. 4 and with the opposite sign. The line shifts large enough that one can plainly see them are not compatible with the MFT, especially not so when the LL shift is included in the MFT. Second, one should note that the graphs for the thicknesses $h=0.25$, $h=0.5$ and $h=1$ are close to one another; the first two are hard to distinguish at all. This holds true even though the density increases by a factor of two for each curve. Our interpretation is that at decreasing thicknesses, here evidently below $h=1$, the physics must eventually become two-dimensional and then it is the area density $N/A$ instead of the volume density $\rho=N/(Ah)$ that governs the scaling. The area density, of course, is the same for all curves in this figure. We have yet another complication to take into account when interpreting the numerical simulations Chomaz et al. (2012).

Even though the “cooperative Lamb shift” (63) was seen in experiments with dense hot gases Keaveney et al. (2012), our simulations for dense cold gases have produced no comparable shifts, and no sign of the LL shift either Javanainen et al. (2014). The line shifts from our simulations simply do not agree with the standard expectations. We next discuss two of our quantitative studies of the line shift from Refs. Javanainen et al. (2014) and Javanainen and Ruostekoski (2016) in added detail.

The “cooperative Lamb shift” (63), in our view, is an etalon effect calculated analytically for the limit of asymptotically low atom density. We therefore study the line shifts as a function of sample thickness $h$ at low atom densities, $\rho=0.01$ and $\rho=0.005$. For the latter the LL shift is 3% of the natural linewidth, so the position of the resonance has to be found accurately. Fortunately, for such low densities the resonance lines $D(\delta)$ are well approximated by Lorentzians, so we can fit them with a Lorentzian with an adjustable center (and width). We did the fits over the detuning range $\delta\in[-2,2]$, and in a few explicit tests found that the resonance positions for data from different simulation runs were reproducible on a few-percent level. The results are shown in Fig. 9. Also shown is the prediction (63), albeit shifted up by the absolute value of the LL shift. The maximum optical thickness in these simulations was $D\simeq 1$.

The oscillatory dependence of the line position on the thickness of the sample is evident, but with two caveats. First, the numerical results track the oscillations in the theory fairly well, but only at thicknesses comparable to or larger than $h=1$. Below that, the approximately constant difference between numerical results and the shifted analytical results shrinks to zero, as the numerically computed shifts clearly tend to 0 for $h\rightarrow 0$. We attribute the thin-sample behavior to transition from 3D physics to 2D physics, as already discussed in conjunction with Fig. 8. For a fixed density $\rho$, the area density $h\rho$ tends to zero as $h\rightarrow 0$, and since in 2D physics it is the area density that counts, the density shift of the line tends to zero as well. Second, there is a large additive constant in the shifts compared to the MFT prediction; recall that we have already removed the LL shift from the theory curve in Fig. 9. It appears that even in the limit of low atom density there is little quantitative validity to the LL shift.

On the other hand, the experiments that found a variation of the line shift in accordance with the “cooperative Lamb shift” and/or the LL shift Maki et al. (1991); Keaveney et al. (2012) were carried out in hot atomic vapors in which the atoms move at thermal speeds, and also collide. We have coded classical-electrodynamics simulations of moving atoms, but the convergence of the results is not yet adequately under control. We therefore adopt a shortcut. Namely, to the lowest order of approximation the motion of the atoms causes Doppler shifts, and as a result the resonance frequencies of the atoms appear to have a corresponding random distribution. We simply add such inhomogeneous broadening to our simulations: While generating a random position for each atom, we also add a random shift to the resonance frequency drawn from a Gaussian distribution with zero average and the rms width $\Omega=100$. This in fact is a reasonable estimate for the D lines in near-room temperature alkali vapors.

The result is a spectrum $D(\delta)$ that under a casual inspection looks like a Gaussian with the rms width $\Omega$. The assignment is to find the center of the resonance line. We resort to a standard method in experimental spectroscopy: We define what is known as the Voigt profile, convolution of a Lorentzian (width $\Gamma$, unit height) and Gaussian with a width $\Omega$,

[TABLE]

where erfc is the complement of the error function as defined, say, in Mathematica, and fit the observed line shape $D(\delta)$ to a Voigt profile $H\,V(\delta-s,\Gamma,\Omega)$. Here we regard the overall height $H$, the shift $s$, and the widths of the Lorentzian $\Gamma$ and of the Gaussian $\Omega$ all as adjustable parameters.

Figure 10 presents the shift of resonance line $s$ as a function of the sample thickness for a disk with the density $\rho=1$, given inhomogeneous broadening with the rms value $\Omega=100$. The dots are from numerical simulations, the sizes being comparable to our estimate of the statistical errors. The solid line is the “cooperative Lamb shift” as from Eq. (63), and the dashed line is a vertically translated version of Eq. (63) that gives the best fit to the numerical data points with $h\geq 1$. The fits of the simulated line shapes to the Voigt profile turn out to be excellent, and in our examples we obtain reproducible results for the shifts $s$ that are on the order 1% of the width of the Gaussian $\Omega$. However, practical constraints forced us to a rather small disk area of $A=256$, which contributes an unknown truncation error. Again, the maximum optical thickness in the samples used to prepare the figure was on the order of $D\simeq 1$.

We once more attribute the exceptional behavior of thin samples, $h\lesssim 1$, to the transition from 3D to 2D physics. Other than this, the simulation results are quite close to the prediction Eq. (63), the deviation being about $0.4\,|\delta_{L}|$. There was a similar “collision shift” in the recent experiments Keaveney et al. (2012) that verified the prediction (63), so the agreements of our numerical experiments and real laboratory experiments with the theory (63) are on a similar footing.

The effect of inhomogeneous broadening is to modify the optical response by emphasizing the mean-field phenomenology via the suppression of light-induced correlations between the atoms. The basic principle is simple to understand: with increasing inhomogeneous broadening the atoms are farther away from resonance with the light that mediates the interactions between the atoms. We can illustrate the interplay between the inhomogeneous broadening and light-mediated interactions by a simple two-atom example Javanainen et al. (2014). The atoms 1 and 2 are assumed to have different resonance frequencies, hence different polarizabilities $\alpha_{1}$ and $\alpha_{2}$. The field amplitude at the atom 2 is then the sum of the incident field amplitude and the field scattered by the atom 1. Formally, we can write it as

[TABLE]

The operator expression in the denominator in the first line is expanded in a power series, as illustrated in Fig. 11. The first term is the free field on atom 2; in the second term the free field excites atom 1, which sends its dipolar field back on atom 2; in the third term the free field excites atom 2, which sends a dipolar field to excite atom 1, which sends a dipolar field back on atom 2. Further terms in the expansion come out the same way reflecting repeated photon exchanges between the atoms. The last term shown is also the first example of a recurrent scattering process in which a light wave interacts more than once with the same atom.

Let us now regard atom 2 as the spectator and imagine averaging over the position of atom 1. This operation faces major mathematical obstacles because of the divergence of ${\sf G}({\bf r}_{1},{\bf r}_{2})$, but we do not attempt to sort them out because these problems are evidently similar for homogeneously and inhomogeneously broadened samples. Next add the inhomogeneous broadening $\omega_{D}$. To the order of magnitude, averaging over the resonant frequencies suppresses the polarizability by a factor of $\gamma/\omega_{D}$. Thus, the first nontrivial term in the expansion corresponding the mean-field polarization gets suppressed by this small factor, and the higher terms by higher powers of the small quantity $\gamma/\omega_{D}$. Qualitatively, repeated photon exchanges are de-emphasized because in such processes both the emitter and the absorber are off resonance.

Our numerical simulations confirm analogous behavior in many-atom ensembles. In fact, when the inhomogeneous broadening $\omega_{D}$ exceeds the resonance linewidth $\gamma$ of the atoms, the results begin to approach the mean-field phenomenology of standard optics, indicating that macroscopic EDPM is an emergent theory, resulting from the suppression of light-induced correlations. The same effect was demonstrated experimentally in the case of fluorescence where the resonance shifts of a cold, dense gas of atoms substantially differed from those predicted for thermal atomic ensembles Jenkins et al. (2016a). Both experimental observations and numerical simulations revealed the absence of any notable shift in cold trapped atomic ensemble. However, introducing inhomogeneous broadening in the simulations restored a large value for the shift.

The suppression of light-mediated interactions by inhomogeneous broadening is a generic effect in coupled resonant emitter systems. For instance, electromagnetic interactions between solid-state radiators, such as plasmonic circuit resonators, may be described by analogous coupled-dipole model simulations Jenkins and Ruostekoski (2012b). Inhomogeneous broadening in such a system can result, e.g., from fabrication imperfections, and have been shown to notably suppress strong radiative interactions between the resonators Jenkins and Ruostekoski (2012c).

We conclude with an after-the-fact test that reinforces our interpretations. We take simulation data for the same parameters we used to demonstrate the qualitative failure of the MFT in Fig. 6, optical thickness $D$ for the sample density $\rho=1$ and thickness $h=1$ as a function of the tuning of the driving light $\delta$, and plot on the same figure also the optical thickness divided by two for a disk that is twice as thick. The result is shown in Fig. 12. The curve with $h=2$ is beset with visible numerical noise since with these atom numbers, $N=2048$, the runs are getting expensive and we have used a reduced number of samples for the atomic positions. Nevertheless, the obvious conclusion is that doubling the thickness from $h=1$ to $h=2$ to a good approximation doubles the optical thickness.

This firstly means that two slabs of thickness $h=1$ back to back would basically behave like one slab with thickness $h=2$. In other words, $h=1$ already represents bulk, 3D, behavior, as we have concluded three times already under different conditions.

Second, within our shortcut to compute the transmission of light, for these parameters the transmission still decreases exponentially with sample thickness. Now, in such a sample the excitation of the dipoles obviously decrease approximately exponentially downstream in the sample as well. If there were a transition to diffuse optics or if the optics of the finite-size sample ($A=1024$) drastically changed with the increasing thickness, there would be changes in the functional dependence of the dipole moments on the distance downstream, which even our approximate way of calculating the transmission would presumably have picked up. MFT fails when diffuse optics sets in, and the forward-scattering approximation does not fully include the optics of the finite-size disk, but neither of these complications apparently is substantial even for the thicker sample with the maximum optical thickness of $D\simeq 5$. We again surmise that our simulations are in the regime when our comparisons with MFT are meaningful.

VI Concluding remarks

Our basic observation is that textbook EDPM and the ensuing usual optics may fail qualitatively as dipole-dipole interactions between the atoms get stronger with increasing density of the atoms Javanainen and Ruostekoski (2016). An effective-medium MFT that spreads out the neighboring atoms into a continuous polarization no longer suffices to describe the influence of the other atoms on an each “spectator” atom. Instead, the effect of the other atoms depends on where exactly they are.

In several examples we have studied the question of when the MFT starts showing strain. The scaling of the whole problem we have employed throughout this paper suggests that the relevant scale is on the order $\rho\simeq 1$, or $\rho\simeq k^{3}$ in dimensional units. This is the kind of a density one would see in experiments with Bose-Einstein condensates, or with tightly trapped cold atoms. However, the scaling-away of dimensional quantities in itself does not give any particular numerical criterion. The observation from our simulations is that MFT may be off by quite a lot already at $\rho\sim 0.01$.

Here we would like to force the issue of (on-resonance) optical thickness versus density Javanainen and Ruostekoski (2016), $6\pi\rho h$ versus $\rho$ in our examples. Optical thickness is a characteristic dimensionless parameter of MFT, and while MFT remains valid, optical thickness may be expected to be the dimensionless parameter. Numerous theoretical analyses are phrased in terms of optical thickness, and a scaling with optical thickness has been demonstrated in recent experiments, e.g. Kwong et al. (2015); Guerin et al. (2016a); Bromley et al. (2016); Araújo et al. (2016). It is not a surprise that one can observe superradiance even in standard optics, and that it scales with optical thickness; optical thickness makes optical resonances broader, whereupon the conventional wisdom about Fourier transformations automatically predicts shortertening time scales. However, from our perspective the more interesting case would be when MFT fails, whereupon, we hypothesize, the density becomes an independent parameter governing the deviations from the MFT. We demonstrate such behavior in Fig. 8 obtained from our simulations, but at present there apparently are no real experiments showing this type of $\rho$ scaling. On the contrary, the scaling of subradiance with optical thickness in a dilute sample as observed experimentally Guerin et al. (2016a) severely challenges our picture, as it is unclear if subradiance can exist in MFT in the first place. At the moment we have no resolution to this issue.

There are phenomena for which EDPM and standard optics with their continuous polarization field do not apply as a matter of principle. Incoherent scattering sideways, as in the two bands at the base of the angular distribution of forward scattering on the right panel of Fig. 1, is an example. In this case, though, we could amend standard optics and still make predictions for sideways scattering: In the single-scattering approximation we would simply add the intensities (not amplitudes) of the light scattered from different atoms. On the other hand, if one studies resonance fluorescence from a few ions, EDPM is a meaningless as a starting point. One can easily imagine intermediate scenarios. What are the predictions from standard optics may also be very difficult to determine per se: If the atomic sample is comparable to the wavelength in size, EDPM boils down to solving the full Maxwell’s equations, which remains a challenge even numerically. All of these caveats notwithstanding, we propose the criterion that a phenomenon should not be called cooperative if standard optics cannot reasonably be excluded as the cause. To give an example, we would object to the notion that the functioning of eyeglasses reflects cooperative response to light of the molecules that make the lenses.

As we have already noted, our recent interest in this research area was triggered by our observation that we did not see the predicted LL shift in numerical simulations of disks of dense, cold gas. The absence of the LL shift has since been demonstrated in light scattering experiments from a small and dense trapped cloud of atoms Pellegrino et al. (2014); Jenkins et al. (2016a, b). These experiments were about sideways scattering, however, which does not directly belong to the MFT framework. From our present viewpoint it is particularly relevant that experiments have also been carried out with forward scattered light under similar conditions that could be directly compared with optics solved numerically from Maxwell’s equations Jennewein et al. (2016). The general result was that at higher atom numbers ($\sim 180$) ab-initio simulation analogous to the ones we have described here came closer to the experimental results than the predictions from optics. However, “[t]he remaining difference with the microscopic model shows that a quantitative understanding of the light-induced interactions even in a relatively simple situation is still a challenge” Jennewein et al. (2016).

The line shifts still present a puzzle. We found the oscillatory dependence of the line shift in accordance with the etalon effects in our simulations of both dilute and inhomogeneously broadened samples, but the LL shift is a more delicate affair. Dimensional analysis and the experience in spectroscopy suggest that at asymptotically low densities there should be a line shift proportional to sample density $\rho$ ($\propto\rho k^{-3}$ in terms of full dimensional quantities). The LL shifts amounts to a specific prediction for the numerical factor that cannot be deduced from dimensional analysis alone. In dilute homogeneously broadened samples we found a LL type shift that is on the order of $\rho$, but even has the opposite sign than the LL shift. Now, if we expand susceptibility of the gas as a power series in density, the LL shift produces a term proportional to $\rho^{2}$. In the usual way of MFTs, EDPM apparently is not a systematic expansion in density Ruostekoski and Javanainen (1997a). However, there are indications that going beyond MFT in an ensemble of randomly distributed atoms could produce corrections proportional to $\rho^{2}$ in quantities such as susceptibility Morice et al. (1995); Ruostekoski and Javanainen (1997a, 1999) as well, and corresponding density-dependent line shifts. Our tentative conclusion is that, if there is any validity to the usual concept of local-field corrections in homogeneously broadened samples to begin with, beyond-MFT effects probably overwhelm them.

The case of inhomogeneously broadened samples is also intriguing. The low-density phenomenology persisted in our examples at least up to $\rho=1$, which in and of itself is not a surprise as inhomogeneous broadening reduces the dipole-dipole interactions. Qualitatively, only a fraction on the order of the ratio of the homogeneously and inhomogeneously broadened linewidths of the atoms has a frequency that can be on resonance with the light propagating in the sample. This reduces the dipole-dipole interactions and extends the range of validity of the MFT Javanainen et al. (2014); Jenkins et al. (2016a). On the other hand, when we fitted the Voigt profile to inhomogeneously broadened absorption lines, we found that the etalon-effect oscillations reside on top of a base lineshift that is about 60% of the LL shift.

The local-field corrections, of which the LL shift is a particular example, have been an enormously successful concept in the physics of electricity and magnetism for well over a century. The best we could do was to get to within 60% of the LL shift. If we posit that the notion of local-field corrections is quantitatively sound, the question is, why did we never do better than 60%? We think that there is a significant piece of physics missing here, but so far it has eluded us.

Acknowledgements.

We acknowledge support from NSF, Grant Nos. PHY-0967644 and PHY-1401151, EPSRC, and 2016 Hongik University Research Fund. Most of the computations were done on Open Science Grid, VO Gluex, and on the University of Southampton Iridis 4 computer cluster.