# Classical versus quantum intensity-field correlations of scattered light
from extended cold atomic clouds

N. Piovella  Dipartimento di Fisica ”Aldo Pontremoli”, Università degli Studi
di Milano, Via Celoria 16, I-20133 Milano, Italy

###### Abstract

We calculate the intensity-field correlations in the light scattered by $N$
cold atoms driven by a quasi-resonant laser field. Fundamental differences
occur if the atomic state is an entangled single-excitation state or a
coherent factorized state. We provide analytic expressions for the two-time
field and intensity correlation functions for the timed Dicke state and the
quasi-Bloch state. The comparison with multi-atom simulations shows good
agreement between numerical and analytic solutions.

##  I Introduction

Statistical properties of the radiation scattered by $N$ identical two-level
atoms are usually studied by means of the single-time second-order normalized
correlation function $g^{(2)}$ [1]. When the system is in an excited state,
$g^{(2)}$ roughly equals two. This value is typical of the Hanbury Brown and
Twiss effect for thermal or pseudothermal Gaussian fields [2]. When the system
decays from the excited state to the superradiant state [3, 4], the behavior
of $g^{(2)}$ is essentially classical, since it decreases from 2 to nearly 1,
which is the value from a coherent state. However, as the system further
decays to the ground state, the intensity fluctuation function, calculated in
a model where the field is quantized, is quite different from the classical
one [5]. Similar results are expected also for subradiance, the odd-twin of
superradiance, in which the excited atoms stay trapped or decay slower to the
ground state. It has been shown that subradiance can be considered to be a
purely nonlocal, nonclassical phenomenon displayed by quantum sources [6].

More recently, a renewed interest in superradiance and subradiance has
increased in the context of the cooperative light scattering in cold atomic
samples [7, 8]. These studies stem from the seminal work by Scully et al. [9],
which described the single-photon superradiance from $N$ two-level atoms
prepared by the absorption of one photon of wave vector $\mathbf{k}_{0}$ [10,
11, 12]. It has been shown that the photon is spontaneously emitted in the
same direction of the incident photon with a cooperative decay rate
proportional to $N$ and inversely proportional to the size of the atomic
cloud. These studies considered the decay of atoms prepared in the ‘timed
Dicke state’:

$|+\rangle_{\mathbf{k}_{0}}=\frac{1}{\sqrt{N}}\sum_{j=1}^{N}e^{i\mathbf{k}_{0}\cdot\mathbf{r}_{j}}|g_{1},g_{2},\dots,e_{j},\dots,g_{N}\rangle$
(1)

where $|g_{1},g_{2},\dots,e_{j},\dots,g_{N}\rangle$ is a Fock state in which
the $j$th atom is prepared in the excited state $|e_{j}\rangle$ and all the
other atoms are in the ground state $|g_{j}\rangle$, and $\mathbf{r}_{j}$ is
the position of the $j$th-atom. A natural platform for studying this kind of
superradiance is provided by the cooperative scattering, in which the atoms
cooperate to scatter the photons from an incident laser beam with wave number
$\mathbf{k}_{0}$ and frequency $\omega_{0}=ck_{0}$ close to the atomic
resonance, leading to a directional emission. This phenomenon is due to the
synchronization of the atomic dipoles with the laser. By different experiments
and comparing different models [13, 14] it has been shown that in the linear
optics approximation, i.e. with a weak excitation of the atoms, the behavior
is almost classical and can be described by considering the atoms as classical
dipoles. Quantum effects detectable by measuring the scattered intensity or
the cooperative radiation force exerted on the atoms could be observed only
with a strong pump such that saturation effects become important [15], or by
investigating the fluctuations of the system. The objective of this work is to
give evidence of the quantum or classical effects observable in the linear
optics approximation, by evaluating the two-time correlation functions
$g^{(1)}(t,\tau)$ and $g^{(2)}(t,\tau)$ for the photons scattered by cold two-
level atoms in different quantum states, i.e., (a) the entangled single-
excitation state, (b) the timed-Dicke state [9], (c) the product (or coherent)
state, and (d) what we call a quasi-Bloch [16] or Eberly’s [12] state. Whereas
$g^{(1)}(t,\tau)$ is the same for all these states, the intensity-correlations
function $g^{(2)}(t,\tau)$ exhibits fundamental differences related to the
quantum or classical description of the atomic ensemble.

##  II The model

Our system consists of a gas of $N$ two-level atoms (with random, fixed
positions $\mathbf{r}_{j}$, lower and upper states $|g_{j}\rangle$ and
$|e_{j}\rangle$ with $j=1,\dots,N$, and transition frequency $\omega_{a}$ with
linewidth $\Gamma=d^{2}\omega_{a}^{3}/2\pi\hbar\epsilon_{0}c^{3}$, where $d$
is the electric dipole matrix element), driven by a uniform resonant radiation
beam with wave vector $\mathbf{k}_{0}=k_{0}\mathbf{\hat{e}}_{z}$, frequency
$\omega_{0}=\omega_{a}+\Delta_{0}$ and electric field $E_{0}$. The atom-field
interaction Hamiltonian in the rotating-wave approximation (RWA) is

$\hat{H}=\hat{H}_{0}+\hat{H}_{1}$ (2)

where

$\displaystyle\hat{H}_{0}$ $\displaystyle=$
$\displaystyle\hbar\sum_{j=1}^{N}\left\\{-\frac{\Delta_{0}}{2}\hat{\sigma}_{3j}+\frac{\Omega_{0}}{2}\left(\hat{\sigma}_{j}e^{-i\mathbf{k}_{0}\cdot\mathbf{r}_{j}}+{\hat{\sigma}_{j}}^{\dagger}e^{i\mathbf{k}_{0}\cdot\mathbf{r}_{j}}\right)\right\\}$
$\displaystyle\hat{H}_{1}$ $\displaystyle=$
$\displaystyle\hbar\sum_{j=1}^{N}\sum_{\mathbf{k}}g_{k}\left[\hat{a}_{\mathbf{k}}^{\dagger}\hat{\sigma}_{j}e^{i(\omega_{k}-\omega_{0})t-i\mathbf{k}\cdot\mathbf{r}_{j}}+{\hat{\sigma}_{j}}^{\dagger}\hat{a}_{\mathbf{k}}e^{-i(\omega_{k}-\omega_{0})t+i\mathbf{k}\cdot\mathbf{r}_{j}}\right].$
(3)

Here $\Omega_{0}=dE_{0}/\hbar$ is the pump Rabi frequency,
$\hat{a}_{\mathbf{k}}\exp(-i\omega_{k}t)$ is the photon annihilation operator
in the interaction picture, with wavenumber $\mathbf{k}$ and frequency
$\omega_{k}=ck$, $g_{k}=d\sqrt{\omega_{k}/(2\hbar\epsilon_{0}V_{ph})}$,
$V_{ph}$ is the photon volume,
$\hat{\sigma}_{j}=\exp(i\Delta_{0}t)|g_{j}\rangle\langle e_{j}|$ and
$\hat{\sigma}_{3j}=|e_{j}\rangle\langle e_{j}|-|g_{j}\rangle\langle g_{j}|$.
We write the Heisenberg equations of the atomic and field operators as

$\displaystyle\frac{d\hat{\sigma}_{j}}{dt}$ $\displaystyle=$
$\displaystyle\frac{1}{i\hbar}[\hat{\sigma}_{j},\hat{H}]=i\Delta_{0}\hat{\sigma}_{j}+\frac{i\Omega_{0}}{2}\hat{\sigma}_{3j}e^{i\mathbf{k}_{0}\cdot\mathbf{r}_{j}}+i\sum_{\mathbf{k}}g_{k}\hat{\sigma}_{3j}\hat{a}_{\mathbf{k}}e^{-i(\omega_{k}-\omega_{0})t+i\mathbf{k}\cdot\mathbf{r}_{j}}$
(4) $\displaystyle\frac{d\hat{\sigma}_{3j}}{dt}$ $\displaystyle=$
$\displaystyle\frac{1}{i\hbar}[\hat{\sigma}_{3j},\hat{H}]=i\Omega_{0}\hat{\sigma}_{j}e^{-i\mathbf{k}_{0}\cdot\mathbf{r}_{j}}+2i\sum_{\mathbf{k}}g_{k}\hat{a}_{\mathbf{k}}^{\dagger}\sigma_{j}e^{i(\omega_{k}-\omega_{0})t-i\mathbf{k}\cdot\mathbf{r}_{j}}+\textrm{h.c.}$
(5) $\displaystyle\frac{d\hat{a}_{\mathbf{k}}}{dt}$ $\displaystyle=$
$\displaystyle\frac{1}{i\hbar}[\hat{a}_{\mathbf{k}},\hat{H}]=-ig_{k}e^{i(\omega_{k}-\omega_{0})t}\sum_{m=1}^{N}\hat{\sigma}_{m}e^{-i\mathbf{k}\cdot\mathbf{r}_{m}}.$
(6)

We consider the atoms initially in their ground state and we assume weak
excitation ($\Omega_{0}\ll\Gamma$), so that we approximate
$\hat{\sigma}_{3j}(t)\approx-\hat{I}_{j}$, where $\hat{I}_{j}$ is the identity
operator for the $j$th atom. This approximation amounts to neglecting
saturation and multiexcitation, i.e. all the processes generating more than
one photon at the same time (linear regime). Integrating Eq.(6) and
substituting it into Eq.(4), neglecting $a_{k}(0)$ (since the initial field
state is vacuum) we obtain

$\displaystyle\frac{d\hat{\sigma}_{j}}{dt}$ $\displaystyle=$
$\displaystyle\left(i\Delta_{0}-\frac{\Gamma}{2}\right)\hat{\sigma}_{j}-\frac{i\Omega_{0}}{2}\hat{I}_{j}e^{i\mathbf{k}_{0}\cdot\mathbf{r}_{j}}-\sum_{\mathbf{k}}g_{k}^{2}\sum_{m=1}^{N}e^{i\mathbf{k}\cdot(\mathbf{r}_{j}-\mathbf{r}_{m})}\int_{0}^{t}dt^{\prime}\hat{\sigma}_{m}(t-t^{\prime})\,e^{-i(\omega_{k}-\omega_{0})t^{\prime}}.$
(7)

The last term in Eq.(7) describes the effect of the spontaneously emitted
photons on the atoms. In the Markov approximation (i.e. when the photon
transit time through the atomic sample is much shorter than the excitation
decay time), we may approximate under the integral
$\hat{\sigma}_{m}(t-t^{\prime})\approx\hat{\sigma}_{m}(t)$. Then, the
remaining time integral yields a real part [with a term proportional to
$\delta(k-k_{0})$] and an imaginary part (corresponding to the principal part
of the integral). We transform the sum over the modes $\mathbf{k}$ into an
integral, $\sum_{\mathbf{k}}\rightarrow(V_{ph}/8\pi^{3})\int d\mathbf{k}$. The
real and imaginary parts of the double integral over $t$ and $\mathbf{k}$
yield the cooperative decay and frequency shift (collective Lamb shift),
respectively. The proper expression of the cooperative frequency shift has
been obtained adding to the Hamiltonian (3) the not-RWA contributions
associated to virtual photons exchanged between different atoms. It results in
the following relation [17]:

$\sum_{\mathbf{k}}g_{k}^{2}e^{i\mathbf{k}\cdot\mathbf{R}}\int_{0}^{\infty}dt^{\prime}e^{-ic(k-k_{0})t^{\prime}}\longrightarrow\frac{\Gamma}{2ik_{0}|\mathbf{R}|}e^{ik_{0}|\mathbf{R}|}$
(8)

where $\Gamma=V_{ph}g_{k_{0}}^{2}k_{0}^{2}/(\pi c)$. Using Eq.(8) in Eq.(7) we
obtain [18],

$\displaystyle\frac{d\hat{\sigma}_{j}(t)}{dt}$ $\displaystyle=$
$\displaystyle\left(i\Delta_{0}-\frac{\Gamma}{2}\right)\hat{\sigma}_{j}(t)-\frac{i\Omega_{0}}{2}\hat{I}_{j}e^{i\mathbf{k}_{0}\cdot\mathbf{r}_{j}}-\frac{\Gamma}{2}\sum_{m\neq
j}^{N}\gamma_{jm}\hat{\sigma}_{m}(t).$ (9)

where $\gamma_{jm}=\exp(ik_{0}r_{jm})/(ik_{0}r_{jm})$ and
$r_{jm}=|\mathbf{r}_{j}-\mathbf{r}_{m}|$. Equation (9) describes the time
evolution of the atomic operators $\hat{\sigma}_{j}$ of $N$ weakly excited
atoms. The real part of $\gamma_{jm}$ describes the spontaneous emission decay
and the imaginary part of $\gamma_{jm}$ describes the energy shift due to
resonant dipole-dipole interactions. Note that we use a scalar model for the
field, neglecting thus any polarization and near field dependence. Detailed
calculations for small and large samples of various geometries show that near-
field and far-field contributions as well as resonant and antiresonant terms
need to be taken properly into account for quantitative predictions [19, 20],
and the present model thus needs to be considered with care illustrating only
a part of the dipole-dipole coupling for real systems.

The positive-frequency part of the electric field is defined as

$\hat{E}^{(+)}(\mathbf{r},t)=i\sum_{\mathbf{k}}{\cal
E}_{k}\hat{a}_{\mathbf{k}}(t)e^{-i\omega_{k}t+i\mathbf{k}\cdot\mathbf{r}}$
(10)

where ${\cal E}_{k}=\sqrt{\hbar\omega_{k}/2\epsilon_{0}V_{ph}}$ is the single-
photon electric field. By integrating Eq.(6) and inserting it in Eq.(10) we
obtain

$\hat{E}^{(+)}(\mathbf{r},t)=\sum_{\mathbf{k}}{\cal
E}_{k}g_{k}\sum_{m=1}^{N}e^{i\mathbf{k}\cdot(\mathbf{r}-\mathbf{r}_{m})-i\omega_{0}t}\int_{0}^{t}dt^{\prime}e^{-i(\omega_{k}-\omega_{0})t^{\prime}}\hat{\sigma}_{m}(t-t^{\prime})$
(11)

Using Eq.(8), the Markov approximation leads to

$\hat{E}^{(+)}(\mathbf{r},t)\approx-i\frac{dk_{0}^{2}}{4\pi\epsilon_{0}}\sum_{j=1}^{N}\frac{e^{-i\omega_{0}(t-|\mathbf{r}-\mathbf{r}_{j}|/c)}}{|\mathbf{r}-\mathbf{r}_{j}|}\hat{\sigma}_{j}(t)$
(12)

which has a transparent interpretation as the sum of wavelets scattered by $N$
dipoles of position $\mathbf{r}_{j}$ and detected at distance $\mathbf{r}$ and
time $t$. In the far field limit, $|\mathbf{r}-\mathbf{r}_{j}|\approx
r-(\mathbf{r}\cdot\mathbf{r}_{j})/r$ and [21]

$\hat{E}^{(+)}(\mathbf{r},t)\approx-i\frac{dk_{0}^{2}}{4\pi\epsilon_{0}r}e^{-i\omega_{0}(t-r/c)}\sum_{j=1}^{N}e^{-i\mathbf{k}\cdot\mathbf{r}_{j}}\hat{\sigma}_{j}(t)$
(13)

where $\mathbf{k}=k_{0}(\mathbf{r}/r)$. The intensity of scattered radiation
by $N$ atoms measured at distance $\mathbf{r}$ and time $t$ is

$I_{N}(\mathbf{r},t,\mathbf{r},t)=\frac{c\epsilon_{0}}{2}\langle\hat{E}^{(-)}(\mathbf{r},t)\hat{E}^{(+)}(\mathbf{r},t)\rangle=\frac{cd^{2}k_{0}^{4}}{32\pi^{2}\epsilon_{0}r^{2}}\sin^{2}\theta_{\mathbf{k}}\cdot
I(t)$ (14)

where $\theta_{\mathbf{k}}$ is the detection angle and

$I(t)=\overline{\sum_{j,m}e^{-i\mathbf{k}\cdot(\mathbf{r}_{j}-\mathbf{r}_{m})}\langle\hat{\sigma}_{m}^{\dagger}(t)\hat{\sigma}_{j}(t)\rangle}$
(15)

is the dimensionless intensity, where the bar calls for an average over
positions, assuming that the locations of the $N$ atoms are not controlled. In
that case an average over a series of experiments with many otherwise
identical samples will be needed to lead to a stable observed intensity.

##  III Correlation functions

The first-order and second-order coherence of light can be described by the
normalized two-time and equal position correlation functions introduced by
Glauber [22, 23]:

$\displaystyle g^{(1)}(t,\tau)$ $\displaystyle=$ $\displaystyle\frac{\langle
E^{(-)}(\mathbf{r},t)E^{(+)}(\mathbf{r},t+\tau)\rangle}{\langle
E^{(-)}(\mathbf{r},t)E^{(+)}(\mathbf{r},t)\rangle}$ (16) $\displaystyle
g^{(2)}(t,\tau)$ $\displaystyle=$ $\displaystyle\frac{\langle
E^{(-)}(\mathbf{r},t)E^{(-)}(\mathbf{r},t+\tau)E^{(+)}(\mathbf{r},t+\tau)E^{(+)}(\mathbf{r},t)\rangle}{\langle
E^{(-)}(\mathbf{r},t)E^{(+)}(\mathbf{r},t)\rangle^{2}}$ (17)

(we drop the hats over the operators) where $\tau$ is the time difference
between two-photon detection events within a two-photon coincidence count.
Two-photon bunching is defined as $g^{(2)}(0)>g^{(2)}(\tau)$ ($\tau\neq 0$),
whereas antibunching is defined as $g^{(2)}(0)<g^{(2)}(\tau)$ ($\tau\neq 0$),
which is usually regarded as a nonclassical effect [24]. Using (13) we obtain

$\displaystyle g^{(1)}(t,\tau)$ $\displaystyle=$
$\displaystyle\frac{1}{I(t)}\overline{\sum_{j,m}e^{-i\mathbf{k}\cdot(\mathbf{r}_{j}-\mathbf{r}_{m})}\langle\sigma_{m}^{\dagger}(t)\sigma_{j}(t+\tau)\rangle}$
(18) $\displaystyle g^{(2)}(t,\tau)$ $\displaystyle=$
$\displaystyle\frac{1}{I^{2}(t)}\overline{\sum_{j,m,p,q}e^{-i\mathbf{k}\cdot(\mathbf{r}_{j}-\mathbf{r}_{m}+\mathbf{r}_{q}-\mathbf{r}_{p})}\langle\sigma_{p}^{\dagger}(t)\sigma_{m}^{\dagger}(t+\tau)\sigma_{j}(t+\tau)\sigma_{q}(t)\rangle}.$
(19)

The determination of the correlation functions $g^{(1)}(t,\tau)$ and
$g^{(2)}(t,\tau)$ requires the evaluation of the two-time quantum averages
$\langle\sigma_{m}^{\dagger}(t)\sigma_{j}(t+\tau)\rangle$ and
$\langle\sigma_{p}^{\dagger}(t)\sigma_{m}^{\dagger}(t+\tau)\sigma_{j}(t+\tau)\sigma_{q}(t)\rangle$,
which can be related to that of averages evaluated at a single time using the
quantum regression theorem [25, 26].

##  IV Single-excitation state

Let us restrict the Hilbert space of the $N$ atoms to the subspace spanned by
the ground state $|g\rangle=|g_{1},\dots,g_{N}\rangle$ and the single-excited-
atom states $|j\rangle=|g_{1},\dots,e_{j},\dots,g_{N}\rangle$ with
$j=1,\dots,N$:

$|\Psi(t)\rangle=\alpha(t)|g\rangle+e^{-i\Delta_{0}t}\sum_{j=1}^{N}\beta_{j}(t)|j\rangle$
(20)

where we approximate $\alpha\approx 1$, so that
$\langle\hat{\sigma}_{j}\rangle=\alpha^{*}\beta_{j}\approx\beta_{j}$ where,
from Eq.(9),

$\displaystyle\frac{d\beta_{j}(t)}{dt}$ $\displaystyle=$
$\displaystyle\left(i\Delta_{0}-\frac{\Gamma}{2}\right)\beta_{j}(t)-\frac{i\Omega_{0}}{2}e^{i\mathbf{k}_{0}\cdot\mathbf{r}_{j}}-\frac{\Gamma}{2}\sum_{m\neq
j}\gamma_{jm}\beta_{m}(t)$ (21)

with initial conditions $\beta_{j}(0)=0$. Since
$\langle\sigma_{m}^{\dagger}(t)\sigma_{j}(t)\rangle=\beta_{m}^{*}(t)\beta_{j}(t)$,
the dimensionless average intensity at the time $t$ is

$\displaystyle I(t)$ $\displaystyle=$
$\displaystyle\overline{\left|\sum_{j}e^{-i\mathbf{k}\cdot\mathbf{r}_{j}}\beta_{j}(t)\right|^{2}}.$
(22)

Using the quantum regression theorem, it is possible to show that (see
Appendix A)

$\displaystyle g^{(1)}(t,\tau)$ $\displaystyle=$
$\displaystyle\frac{1}{I(t)}\overline{\sum_{j,m}e^{-i\mathbf{k}\cdot(\mathbf{r}_{j}-\mathbf{r}_{m})}\beta_{m}^{*}(t)\beta_{j}(t+\tau)}e^{-i\omega_{0}\tau}.$
(23)

and

$\displaystyle g^{(2)}(t,\tau)$ $\displaystyle=$
$\displaystyle\frac{1}{I^{2}(t)}\overline{\sum_{j,m,p,q}e^{-i\mathbf{k}\cdot(\mathbf{r}_{j}-\mathbf{r}_{m}+\mathbf{r}_{q}-\mathbf{r}_{p})}\beta_{p}^{*}(t)H_{mj}(\tau)\beta_{q}(t)}$
(24)

where $H_{mj}(\tau)$ is the solution of the equation

$\displaystyle\frac{d}{d\tau}H_{mj}(\tau)$ $\displaystyle=$
$\displaystyle-\Gamma
H_{mj}(\tau)-\frac{i\Omega_{0}}{2}\left[e^{i\mathbf{k}_{0}\cdot\mathbf{r}_{j}}\beta^{*}_{m}(\tau)-e^{-i\mathbf{k}_{0}\cdot\mathbf{r}_{m}}\beta_{j}(\tau)\right]-\frac{\Gamma}{2}\left[\sum_{k\neq
j}\gamma_{jk}H_{mk}(\tau)+\sum_{k\neq m}\gamma^{*}_{mk}H_{kj}(\tau)\right].$
(25)

with $H_{mj}(0)=0$. We observe that, since $g^{(1)}(t,0)=1$ and
$\lim_{t\rightarrow\infty}|g^{(1)}(t,\tau)|=1$, the system has full first-
order coherence. We can interpret this result also saying that the scattering
is elastic, as expected by a system with a linear response to a cw excitation.
In this case the randomness of the system has no effect on the field
correlation function. Instead, $g^{(2)}(t,0)=0$ since the state (20) has a
single excitation and it is not possible to detect two photons at the same
time. The numerical solution of Eq.(24) will be discussed in Sec. VIII.

##  V Timed Dicke state

Let’s now consider the timed Dicke state of Eq. (1), assuming
$\beta_{j}(t)=\beta(t)e^{i\mathbf{k}_{0}\cdot\mathbf{r}_{j}}$ in Eq. (20).
Then,

$I(t)=N^{2}\overline{|S_{N}(\mathbf{k}-\mathbf{k}_{0})|^{2}}|\beta(t)|^{2}$
(26)

where

$S_{N}(\mathbf{k}-\mathbf{k}_{0})=\frac{1}{N}\sum_{j=1}^{N}e^{-i(\mathbf{k}-\mathbf{k}_{0})\cdot\mathbf{r}_{j}}$
(27)

is the structure factor. Furthermore, since the dynamics is the same for all
the atoms, with $\langle\sigma_{j}(t)\rangle=\langle\sigma(t)\rangle
e^{i\mathbf{k}_{0}\cdot\mathbf{r}_{j}}$,

$\displaystyle g^{(1)}(t,\tau)$ $\displaystyle=$
$\displaystyle\frac{\langle\sigma^{\dagger}(t)\sigma(t+\tau)\rangle}{\langle\sigma^{\dagger}(t)\sigma(t)\rangle}e^{-i\omega_{0}\tau}$
(28) $\displaystyle g^{(2)}(t,\tau)$ $\displaystyle=$ $\displaystyle
R\frac{\langle\sigma^{\dagger}(t)\sigma^{\dagger}(t+\tau)\sigma(t+\tau)\sigma(t)\rangle}{\langle\sigma^{\dagger}(t)\sigma(t)\rangle^{2}}$
(29)

where the factor $R$ is (assuming $N\gg 1$)

$R=\frac{\overline{|S_{N}|^{4}}}{\left[\overline{|S_{N}|^{2}}\right]^{2}}\approx\frac{2+4N|S_{\infty}(\mathbf{k}-\mathbf{k}_{0})|^{2}+N^{2}|S_{\infty}(\mathbf{k}-\mathbf{k}_{0})|^{4}}{(1+N|S_{\infty}(\mathbf{k}-\mathbf{k}_{0})|^{2})^{2}}$
(30)

and we have approximated the structure factor by a continuous distribution
with number density $n(\mathbf{r})$,

$S_{\infty}(\mathbf{k}-\mathbf{k}_{0})=\int
n(\mathbf{r})e^{-i(\mathbf{k}-\mathbf{k}_{0})\cdot\mathbf{r}_{j}}d\mathbf{r}.$
(31)

The term $R$ describes the classical intensity correlation function at zero
delay, due to the spatial distribution of the scatterers. It ranges from $R=1$
to $R=2$ for coherent and chaotic light, respectively. For instance, for a
Gaussian spherical distribution with rms size $\sigma_{r}$,
$S_{\infty}(\mathbf{k}-\mathbf{k}_{0})=\exp[-2(k_{0}\sigma_{r})^{2}\sin^{2}(\theta/2)]$.
Hence $R=1$ for $k_{0}\sigma_{r}|\sin(\theta/2)|\ll 1$ and $R=2$ for
$k_{0}\sigma_{r}|\sin(\theta/2)|\gg 1$ (see Fig. 1).

Figure 1: $R$ vs $\sigma=k_{0}\sigma_{r}$ for $N=10^{3}$ and
$\theta=90^{\circ}$, averaged over $10^{4}$ iterations; solid black line:
numerical simulation; dashed red line: analytic solution.

For a timed Dicke state,
$\langle\sigma^{\dagger}(t)\sigma(t)\rangle=|\beta(t)|^{2}$, with

$\beta(t)=\beta^{st}\left\\{1-e^{(i\Delta-\Gamma_{N}/2)t}\right\\},$ (32)
$\beta^{st}=\frac{\Omega_{0}}{2\Delta+i\Gamma_{N}},$ (33)

$\Gamma_{N}=\Gamma+(N\Gamma/4\pi)\int_{0}^{2\pi}d\phi\int_{0}^{\pi}d\theta\sin\theta|S_{\infty}(k_{0},\theta,\phi)|^{2}$
is the collective decay rate, $\Delta=\Delta_{0}-\Delta_{N}$ and
$\Delta_{N}=(N\Gamma/8\pi^{2})\mathsf{P}\int_{0}^{\infty}d\kappa\kappa^{3}/(\kappa-1)\int_{0}^{2\pi}d\phi\int_{0}^{\pi}d\theta\sin\theta|S_{\infty}(k_{0}\kappa,\theta,\phi)|^{2}$
is the collective Lamb shift [18]. Since $\langle\sigma(t)\rangle=\beta(t)$,
then

$\langle\sigma(t+\tau)\rangle=\beta(t+\tau)=\beta^{st}+[\beta(t)-\beta^{st}]e^{(i\Delta-\Gamma_{N}/2)\tau}=\langle\sigma(\tau)\rangle+e^{(i\Delta-\Gamma_{N}/2)\tau}\langle\sigma(t)\rangle.$
(34)

From the quantum regression theorem,

$\langle\sigma^{\dagger}(t)\sigma(t+\tau)\rangle=\langle\sigma^{\dagger}(t)\rangle\langle\sigma(\tau)\rangle+e^{(i\Delta-\Gamma_{N}/2)\tau}\langle\sigma^{\dagger}(t)\sigma(t)\rangle$
(35)

and

$g^{(1)}(t,\tau)=\frac{\beta^{*}(t)\beta(\tau)+e^{(i\Delta-\Gamma_{N}/2)\tau}|\beta(t)|^{2}}{|\beta(t)|^{2}}e^{-i\omega_{0}\tau}=\frac{1-e^{(i\Delta-\Gamma_{N}/2)(t+\tau)}}{1-e^{(i\Delta-\Gamma_{N}/2)t}}e^{-i\omega_{0}\tau}.$
(36)

As expected, $\lim_{t\rightarrow\infty}|g^{(1)}(t,\tau)|=1$ and the scattered
field has full first-order coherence. In order to evaluate $g^{(2)}(t,\tau)$,
we note that

$\displaystyle\langle\sigma^{\dagger}(t+\tau)\sigma(t+\tau)\rangle$
$\displaystyle=$ $\displaystyle|\beta(t+\tau)|^{2}$ (37) $\displaystyle=$
$\displaystyle\left|\beta^{st}+[\beta(t)-\beta^{st}]e^{(i\Delta-\Gamma_{N}/2)\tau}\right|^{2}$
$\displaystyle=$
$\displaystyle|\beta^{st}|^{2}+|\beta(t)-\beta^{st}|^{2}e^{-\Gamma_{N}\tau}+\left\\{(\beta(t)-\beta^{st}){\beta^{st}}^{*}e^{i\Delta\tau}+\mathrm{c.c.}\right\\}e^{-\Gamma_{N}\tau/2}.$

where $\mathrm{c.c.}$ stands for the complex conjugate. It can be rewritten in
the form

$\displaystyle\langle\sigma^{\dagger}(t+\tau)\sigma(t+\tau)\rangle$
$\displaystyle=$ $\displaystyle
A(\tau)+B(\tau)\langle\sigma(t)\rangle+B^{*}(\tau)\langle\sigma^{\dagger}(t)\rangle+C(\tau)\langle\sigma^{\dagger}(t)\sigma(t)\rangle$
(38)

where

$\displaystyle A(\tau)$ $\displaystyle=$
$\displaystyle|\beta^{st}|^{2}\left\\{1+e^{-\Gamma_{N}\tau}-2e^{-\Gamma_{N}\tau/2}\cos(\Delta\tau)\right\\}$
(39) $\displaystyle B(\tau)$ $\displaystyle=$ $\displaystyle
e^{-\Gamma_{N}\tau/2}{\beta^{st}}^{*}\left\\{e^{i\Delta\tau}-e^{-\Gamma_{N}\tau/2}\right\\}$
(40) $\displaystyle C(\tau)$ $\displaystyle=$ $\displaystyle
e^{-\Gamma_{N}\tau}.$ (41)

From the quantum regression theorem,

$\displaystyle\langle\sigma^{\dagger}(t)\sigma^{\dagger}(t+\tau)\sigma(t+\tau)\sigma(t)\rangle$
$\displaystyle=$ $\displaystyle
A(\tau)\langle\sigma^{\dagger}(t)\sigma(t)\rangle+B(\tau)\langle\sigma^{\dagger}(t)\sigma(t)\sigma(t)\rangle+B^{*}(\tau)\langle\sigma^{\dagger}(t)\sigma^{\dagger}(t)\sigma(t)\rangle$
(42) $\displaystyle+$ $\displaystyle
C(\tau)\langle\sigma^{\dagger}(t)\sigma^{\dagger}(t)\sigma(t)\sigma(t)\rangle.$

Since $\sigma(t)\sigma(t)=0$ and $\sigma^{\dagger}(t)\sigma^{\dagger}(t)=0$,

$g^{(2)}(t,\tau)=R\frac{A(\tau)}{|\beta(t)|^{2}}=R\frac{1+e^{-\Gamma_{N}\tau}-2e^{-\Gamma_{N}\tau/2}\cos(\Delta\tau)}{1+e^{-\Gamma_{N}t}-2e^{-\Gamma_{N}t/2}\cos(\Delta
t)}$ (43)

Hence, $\lim_{t\rightarrow\infty}g^{(2)}(t,\tau)=g^{(2)}(\tau)$ where

$g^{(2)}(\tau)=R\left|1-e^{(i\Delta-\Gamma_{N}/2)\tau}\right|^{2}=R\left\\{1+e^{-\Gamma_{N}\tau}-2e^{-\Gamma_{N}\tau/2}\cos(\Delta\tau)\right\\}.$
(44)

In conclusion, the timed Dicke state behaves as a single driven atom, but with
a collective decay rate $\Gamma_{N}$ and Lamb shift $\Delta_{N}$. As expected,
$g^{(2)}(0)=0$ and $\lim_{\tau\rightarrow\infty}g^{(2)}(\tau)=R$. Figure 2
shows Eq.(44) for $R=2$ and two different values of detuning,
$\Delta=5\Gamma_{N}$ (solid black line) and $\Delta=0$ (red dashed line).
Since $g^{(2)}(\tau)>g^{(2)}(0)$, the system exhibits antibunching. Notice
that, contrary to $g^{(1)}$, $g^{(2)}$ is proportional to the factor $R$,
equal to 2 when the photons are emitted isotropically and randomly out of the
diffraction cone, with aperture $\Delta\theta\sim\lambda_{0}/\sigma_{r}$, and
equal to 1 in the opposite case.

Figure 2: $g^{(2)}(\tau)$ vs $\Gamma_{N}\tau$ for $R=2$ and two different
values of detuning: $\Delta=5\Gamma_{N}$ (solid black line) and $\Delta=0$
(red dashed line).

In the case of the timed Dicke state, $H_{mj}(\tau)=H(\tau)$ and
$g^{(2)}(\tau)=H(\tau)/|\beta^{st}|^{2}$, where $H(\tau)$ is the solution of
the equation

$\displaystyle\frac{dH(\tau)}{d\tau}$ $\displaystyle=$
$\displaystyle\frac{i\Omega_{0}}{2}[\beta(\tau)-\beta^{*}(\tau)]-\Gamma_{N}H(\tau),$
(45)

yielding the solution (44).

##  VI Product state

It has been noted that the results for the intensity $I(t)$ emitted by $N$
weakly-excited two-level atoms, obtained for a symmetric timed Dicke state as
described in Sec.V, could be obtained assuming a product state of $N$ two-
level atoms [12] (named also ‘pure Bloch state’ by some authors [16]). More
specifically, pure Bloch states are product states in which every one of $N$
two-level atoms is in the same superposition of ground and excited state. Such
states are easily produced experimentally. As it happens for the timed Dicke
state, the driving field imposes a coherence in the photons emitted
spontaneously by each atom, so that superradiance arises because the state is
symmetric under exchange of particles. However, it is expected that the
quantum statistic of the timed Dicke state will be quite different from that
of the pure Bloch state. The aim of this section is to obtain the stationary
correlation functions $g^{(1)}(\tau)$ and $g^{(2)}(\tau)$ for such product
state, in the presence of a cw driving field. To be more general, we assume
first that the excitation probability amplitudes $\beta_{j}$ are not the same
for every atom,

$|\Psi(t)\rangle=\prod_{j=1}^{N}\left\\{\alpha(t)|g_{j}\rangle+\beta_{j}(t)e^{-i\Delta_{0}t}|e_{j}\rangle\right\\}$
(46)

with $|\alpha(t)|^{2}+|\beta_{j}(t)|^{2}=1$ for every $j$. Assuming
$\alpha\sim 1$, we obtain $\langle\sigma_{j}\rangle=\beta_{j}$,
$\langle\sigma_{m}^{\dagger}\sigma_{j}\rangle=\beta_{m}^{*}\beta_{j}$ and

$\langle\sigma_{m}\sigma_{j}\rangle=\left\\{\begin{array}[]{ccc}\displaystyle\beta_{m}\beta_{j}&\mathrm{if}&j\neq
m\\\ &&\\\ 0&\mathrm{if}&j=m\end{array}\right..$ (47)

Hence, for the product state all the expectation values of the operators
factorize and the dynamics is determined solely by $\beta_{j}$, solution of
Eq.(21). As a consequence, following the same reasoning adopted in Sec.IV, we
arrive at the same expression (23) for $g^{(1)}(t,\tau)$ obtained assuming the
entangled state (20). Differences between the product state and the single-
excitation state appear when higher-order expectation values are observed. We
report the details of the demonstration in Appendix B, leading to

$\displaystyle g^{(2)}(t,\tau)$ $\displaystyle=$
$\displaystyle\frac{1}{I^{2}(t)}\overline{\sum_{j,m,p,q}e^{-i\mathbf{k}\cdot(\mathbf{r}_{j}-\mathbf{r}_{m}+\mathbf{r}_{q}-\mathbf{r}_{p})}\beta_{p}^{*}(t)H_{mj}(t,\tau)\beta_{q}(t)}$
(48)

where $H_{mj}(\tau)$ is the solution of the equation

$\displaystyle\frac{d}{d\tau}H_{mj}(t,\tau)$ $\displaystyle=$
$\displaystyle-\Gamma
H_{mj}(t,\tau)-\frac{i\Omega_{0}}{2}\left[e^{i\mathbf{k}_{0}\cdot\mathbf{r}_{j}}F^{*}_{m}(t,\tau)-e^{-i\mathbf{k}_{0}\cdot\mathbf{r}_{m}}F_{j}(t,\tau)\right]-\frac{\Gamma}{2}\left[\sum_{k\neq
j}\gamma_{jk}H_{mk}(t,\tau)+\sum_{k\neq
m}\gamma^{*}_{mk}H_{kj}(t,\tau)\right]$

with initial condition

$H_{mj}(t,0)=\left\\{\begin{array}[]{ccc}\displaystyle\beta_{m}^{*}(t)\beta_{j}(t)&&\mathrm{if}\quad
p\neq m\quad\mathrm{or}\quad j\neq q\\\ &&\\\
0&&\mathrm{otherwise}\end{array}\right.$ (50)

and $F_{j}(t,\tau)$ is the solution of

$\displaystyle\frac{d}{d\tau}F_{j}(t,\tau)$ $\displaystyle=$
$\displaystyle\left(i\Delta_{0}-\frac{\Gamma}{2}\right)F_{j}(t,\tau)-\frac{i\Omega_{0}}{2}e^{i\mathbf{k}_{0}\cdot\mathbf{r}_{j}}-\frac{\Gamma}{2}\sum_{k\neq
j}\gamma_{jk}F_{k}(t,\tau).$ (51)

with

$F_{j}(t,0)=\left\\{\begin{array}[]{ccc}\displaystyle\beta_{j}(t)&\mathrm{if}&j\neq
q\\\ &&\\\ 0&&\mathrm{otherwise}\end{array}\right..$ (52)

The numerical solution of Eq. (48) will be discussed in Sec. VIII.

##  VII The Eberly’s state

As done for the timed Dicke state, we assume now that the excitation
probability is the same for all the atoms [12],

$|\Psi(t)\rangle_{E}=\prod_{j=1}^{N}\left\\{\alpha(t)|g_{j}\rangle+\beta(t)e^{-i\Delta_{0}t+i\mathbf{k}_{0}\cdot\mathbf{r}_{j}}|e_{j}\rangle\right\\}$
(53)

with $|\alpha(t)|^{2}+|\beta(t)|^{2}=1$. In the following we will assume the
linear approximation, so that $\alpha(t)\sim 1$ and $\beta(t)$ is given by
Eq.(32). Then, defining

$E(t)=\frac{1}{N}\sum_{j=1}^{N}e^{-i\mathbf{k}\cdot\mathbf{r}_{j}}\sigma_{j}(t)$
(54)

we have

$\displaystyle\langle E(t)\rangle$ $\displaystyle=$ $\displaystyle
S_{N}\beta(t)$ (55) $\displaystyle\langle E(t+\tau)\rangle$ $\displaystyle=$
$\displaystyle
S_{N}\beta(t+\tau)=S_{N}\beta(\tau)+e^{(i\Delta-\Gamma_{N}/2)\tau}\langle
E(t)\rangle$ (56) $\displaystyle\langle E^{\dagger}(t)E(t)\rangle$
$\displaystyle=$ $\displaystyle|S_{N}|^{2}|\beta(t)|^{2}$ (57)
$\displaystyle\langle E^{\dagger}(t+\tau)E(t+\tau)\rangle$ $\displaystyle=$
$\displaystyle A(\tau)|S_{N}|^{2}+B(\tau)S_{N}^{*}\langle
E(t)\rangle+B^{*}(\tau)S_{N}\langle E^{\dagger}(t)\rangle+C(\tau)\langle
E^{\dagger}(t)E(t)\rangle$ (58)

where $S_{N}$ has been defined in Eq.(27) and $A(\tau)$, $B(\tau)$ and
$C(\tau)$ have been defined in Eqs.(39)-(41). From the quantum linear
regression theorem,

$\displaystyle\langle E^{\dagger}(t)E(t+\tau)\rangle$ $\displaystyle=$
$\displaystyle S_{N}\beta(\tau)\langle
E^{\dagger}(t)\rangle+e^{(i\Delta-\Gamma_{N}/2)\tau}\langle
E^{\dagger}(t)E(t)\rangle$

and

$\displaystyle g^{(1)}(t,\tau)$ $\displaystyle=$
$\displaystyle\frac{\beta(\tau)\beta^{*}(t)+e^{(i\Delta-\Gamma_{N}/2)\tau}|\beta(t)|^{2}}{|\beta(t)|^{2}}e^{-i\omega_{0}\tau}$
(60)

which coincides with Eq.(36). From Eq.(58) we obtain

$\displaystyle\langle E^{\dagger}(t)E^{\dagger}(t+\tau)E(t+\tau)E(t)\rangle$
$\displaystyle=$ $\displaystyle A(\tau)|S_{N}|^{2}\langle
E^{\dagger}(t)E(t)\rangle$ (61) $\displaystyle+$
$\displaystyle\left[B(\tau)S_{N}^{*}\langle
E^{\dagger}(t)E(t)E(t)\rangle+\textrm{c.c.}\right]$ $\displaystyle+$
$\displaystyle C(\tau)\langle E^{\dagger}(t)E^{\dagger}(t)E(t)E(t)\rangle$

The last two terms make the difference from the timed Dicke state, since in
general they are not zero. It is possible to demonstrate that (see Appendix C)

$\displaystyle\langle E^{\dagger}(t)E(t)E(t)\rangle$ $\displaystyle=$
$\displaystyle\beta(t)|\beta(t)|^{2}\left\\{\frac{2(N-1)}{N^{2}}+K_{2}\right\\}$
(62) $\displaystyle\langle E^{\dagger}(t)E^{\dagger}(t)E(t)E(t)\rangle$
$\displaystyle=$
$\displaystyle|\beta(t)|^{4}\left\\{\frac{2(N-1)}{N^{3}}+\frac{4(N-2)}{N^{2}}K_{2}+K_{4}\right\\}$
(63)

where

$\displaystyle K_{2}$ $\displaystyle=$
$\displaystyle\frac{1}{N^{2}}\sum_{j}\sum_{p\neq
j}e^{i(\mathbf{k}_{0}-\mathbf{k})\cdot(\mathbf{r}_{p}-\mathbf{r}_{j})}\approx|S_{\infty}|^{2}$
(64) $\displaystyle K_{4}$ $\displaystyle=$
$\displaystyle\frac{1}{N^{4}}\sum_{m}\sum_{j\neq m}\sum_{p\neq j,m}\sum_{q\neq
p,j,m}e^{i(\mathbf{k}_{0}-\mathbf{k})\cdot(\mathbf{r}_{p}+\mathbf{r}_{q}-\mathbf{r}_{j}-\mathbf{r}_{m})}\approx|S_{\infty}|^{4}$
(65)

For $N\gg 1$ and after the statistical average,

$\displaystyle\langle E^{\dagger}(t)E^{\dagger}(t+\tau)E(t+\tau)E(t)\rangle$
$\displaystyle=$
$\displaystyle\left(\frac{2}{N^{2}}+\frac{4}{N}|S_{\infty}|^{2}+|S_{\infty}|^{4}\right)\left(A(\tau)|\beta(t)|^{2}+C(\tau)|\beta(t)|^{4}\right)$
(66) $\displaystyle+$
$\displaystyle\left(\frac{1}{N}+|S_{\infty}|^{2}\right)\left(\frac{2}{N}+|S_{\infty}|^{2}\right)\left[\beta(t)B(\tau)+\textrm{c.c.}\right].$

Using the definitions of $A(\tau)$, $B(\tau)$ and $C(\tau)$ and taking the
limit $t\rightarrow\infty$, we obtain

$\displaystyle g^{(2)}(\tau)$ $\displaystyle=$ $\displaystyle
R+2Q\left[e^{-\Gamma_{N}\tau}-e^{-\Gamma_{N}\tau/2}\cos(\Delta\tau)\right],$
(67)

where

$Q=R-\frac{2+N|S_{\infty}|^{2}}{1+N|S_{\infty}|^{2}}$ (68)

We observe that in the ”chaotic” limit $N|S_{\infty}|^{2}\ll 1$, $R\approx 2$,
$Q\ll 1$, and $g^{(2)}(\tau)\approx 2$. Conversely, in the ”coherent” limit
$N|S_{\infty}|^{2}\gg 1$, $R\approx 1$, $Q\ll 1$, and $g^{(2)}(\tau)\approx
1$. More generally, $g^{(2)}(\tau)$ depends on $N$, $\sigma=k_{0}\sigma_{r}$,
and the detection angle $\theta$. The parameter $Q$ takes its maximum value
$Q_{max}=1/2$ for $N|S_{\infty}|^{2}=1$, with $R=7/4$. As an example, Fig.3
shows $g^{(2)}(\tau)$ vs $\Gamma_{N}\tau$ for a spherical Gaussian
distribution with $N=10^{6}$ and $k_{0}\sigma_{r}=20$, a laser beam with
$\Delta=5\Gamma_{N}$, and detection angles
$\theta=9^{\circ},10^{\circ},11^{\circ},12^{\circ}$. Within a few degrees, the
value of $g^{(2)}(\tau)$ changes from $1$ to $2$, with damped oscillations as
a function of $\tau$.

Figure 3: $g^{(2)}(\tau)$ vs $\Gamma_{N}\tau$ from Eq.(67) for
$\Delta=5\Gamma_{N}$, $N=10^{6}$, $k_{0}\sigma_{r}=20$ and angles
$\theta=12^{\circ}$ (dashed black line), $\theta=11^{\circ}$ (continuous red
line), $\theta=10^{\circ}$ (dashed-dotted blue line) and $\theta=9^{\circ}$
(dotted green line).

##  VIII Numerical results

In this section we compare the exact solution of $g^{(1)}(t,\tau)$ and
$g^{(2)}(t,\tau)$, calculated numerically from Eqs. (23), (24), and (48) with
the analytic expressions of Eqs. (36), (43), and (67). The expression of
$g^{(2)}(t,\tau)$ assuming $N$ classical dipoles is

$\displaystyle g^{(2)}(t,\tau)$ $\displaystyle=$
$\displaystyle\frac{1}{I^{2}(t)}\overline{\sum_{j,m,p,q}e^{-i\mathbf{k}\cdot(\mathbf{r}_{j}-\mathbf{r}_{m}+\mathbf{r}_{q}-\mathbf{r}_{p})}\beta_{p}^{*}(t)\beta_{m}^{*}(t+\tau)\beta_{j}(t+\tau)\beta_{q}(t)}.$
(69)

where $I(t)$ is given by Eq.(22). In this case
$g^{(2)}(\tau)=\lim_{t\rightarrow\infty}g^{(2)}(t,\tau)=R$, i.e. is
independent on $\tau$. Hence the oscillations observed in Eq.(67) for Eberly’s
state seem to have a quantum nature (see Fig.3). Figure 4 shows the result of
a numerical simulation calculating $g^{(2)}(t,\tau)$ from Eq.(24) for a
spherical Gaussian distribution with $N=100$, $k_{0}\sigma_{r}=5$, detuning
$\Delta=5\Gamma$, detection angle $\theta=90^{\circ}$, and $\Gamma t=5$,
averaged over 20 iterations (solid blue line). The numerical result is
compared with the analytic expression of Eq.(44) obtained assuming the timed
Dicke state (red dashed line). The numerical value $R\approx 1.7$ obtained for
$\Gamma\tau\rightarrow\infty$ is less than the chaotic value 2 because of the
small number of particles and iterations. The exact result differs from the
approximated timed Dicke solution because of the spread of values of
$\beta_{j}$ around the average value, causing a decoherence which reduces the
oscillation amplitudes.

Figure 4: $g^{(2)}(\tau)$ vs $\Gamma\tau$, calculated from Eq. (24), for
$N=100$, $k_{0}\sigma_{r}=5$, $\Delta=5\Gamma$, $\theta=90^{\circ}$ and
average after 20 iterations (continuous blue line). The dashed red line is the
analytic expression of Eq. (44), with $R=1.7$.

Figure 5 shows the result of a numerical simulation calculating
$g^{(2)}(t,\tau)$ from Eq.(48) for a spherical Gaussian distribution with
$N=100$, $k_{0}\sigma_{r}=5$, detuning $\Delta=5\Gamma$, detection angle
$\theta=16.26^{\circ}$, and $\Gamma t=5$, averaged over 20 iterations (solid
blue line). The numerical result is compared with the analytic expression of
Eq. (67) obtained assuming Eberly’s state with $R\approx 1.025$,
$Q\approx-0.042$, and $\Gamma_{N}\approx 2\Gamma$ (dashed red line). These
parameters have been calculated numerically for the spatial distribution. Also
in this case, the ideal case of Eberly’s state shows oscillations with
amplitudes larger than the ones of the exact solution. Nevertheless, the
transient oscillations are clearly visible and detectable. We outline again
that these oscillations disappear in the classical limit, when the atoms are
described as classical dipoles.

Figure 5: $g^{(2)}(\tau)$ vs $\Gamma\tau$, calculated from Eq. (44), for
$N=100$, $k_{0}\sigma_{r}=5$, $\Delta=5\Gamma$, $\theta=16.26^{\circ}$ and
average after 20 iterations (continuous blue line). The dashed red line is the
analytic expression of Eq. (67), with R=1.025,$,Q=-0.042$, and
$\Gamma_{N}=2\Gamma$.

##  IX Conclusions

We have calculated the two-time field and intensity correlation functions of
the light scattered by an ensemble of two-level atoms driven by a cw resonant
laser beam, in the linear optics regime. The atoms have fixed random
positions. We have calculated the correlation functions for two different
quantum atomic states, i.e., the single-excitation state and the product
state, obtaining exact expressions to be evaluated numerically. Furthermore,
we have obtained analytic expressions of $g^{(1)}(t,\tau)$ and
$g^{(2)}(t,\tau)$ in the case of uniform excitation. This approximation leads
to the well-known timed Dicke state [9] for the single-excitation state and to
the pure Bloch state [16] for the product state (named here Eberly’s state, in
honor of Eberly [12] who first discussed the differences between these two
states). Our conclusions are that both these states lead, as expected, to
$g^{(1)}(\tau)=\lim_{t\rightarrow\infty}g^{(1)}(t,\tau)=1$, as it results in
also describing the atoms as classical dipoles. This result is independent on
the statistical properties of the atomic distribution, and corresponds to
elastic scattering at the frequency of the incident driving beam. Differences
between the single-excitation state and the product state appear when the
stationary intensity correlation function is considered,
$g^{(2)}(\tau)=\lim_{t\rightarrow\infty}g^{(2)}(t,\tau)$. In the classical
limit (i.e., atoms as classical dipoles), $g^{(2)}(\tau)=R$, where $R$ depends
on the randomness and spatial distribution of the atomic sample, varying from
$R=1$ when the emission is coherent (i.e., within the diffraction cone for an
extended cloud) to $R=2$ for a chaotic random distribution. For the timed
Dicke state, $g^{(2)}(0)=0$ and $\lim_{\tau\rightarrow\infty}g^{(2)}(\tau)=R$,
with a transient time of the order of $1/\Gamma_{N}$, where $\Gamma_{N}$ is
the cooperative decay rate for $N$ atoms. In this case we observe
antibunching, since $g^{(2)}(0)<g^{(2)}(\tau)$. For Eberly’s state,
$g^{(2)}(0)=R$ and $\lim_{\tau\rightarrow\infty}g^{(2)}(\tau)=R$. However, it
is possible to observe transient oscillations as a function of $\tau$ in an
intermediate regime with $1<R<2$ and in the detuned case, where the coherent
and chaotic emission are competing. These oscillations are purely quantum and
are not visible when the atoms are treated as classical dipoles.

We outline again that the aim of this work is to propose a method to
distinguish between the possible states generated in the cooperative
scattering, by measuring the second-order correlation function
$g^{(2)}(\tau)$. In fact, it gives different results for the single-excitation
entangled state or the factorized coherent state. Generally, these states
should necessitate a different preparation, as discussed in the original
papers by Scully and co-workers [9, 27] and more recently studied
experimentally by Felinto and co-workers [28]. Intuitively, it is likely that
an ensemble of $N$ two-level atoms driven by a classical field will be
described by the product state. However, this has not been proved yet, and an
experiment measuring $g^{(2)}(\tau)$ can do it.

This study has assumed a cw driving beam and neglected any atomic motion,
either due to temperature or recoil [29]. It would be interesting in the
future to extend it to include the temperature inducing a decay of the
correlations (see, for instance, Ref. [30]). Also, the statistical properties
of subradiance [8] should deserve attention, which may be investigated with
the present formalism just switching off the driving laser and observing the
fluorescence light emitted by the excited atoms at sufficiently long time,
such that only the subradiance component survives. All these points will be
the object of a future publication.

###### Acknowledgements.

We thank R. Bachelard for useful discussions and R. Gaiba for the analytic
calculations of Appendix C done during his stage in Milan. This work was
performed within the framework of the European Training Network ColOpt, which
is funded by the European Union (EU) Horizon 2020 program under the Marie
Sklodowska-Curie actions, Grant Agreement No. 721465 R.A.

##  Appendix A Evaluation of $g^{(1)}(\tau)$ and $g^{(2)}(\tau)$ for the
single-excitation state (20)

Using Eqs. (9) with the quantum regression theorem,

$\displaystyle\frac{d\langle\sigma_{m}^{\dagger}(t)\sigma_{j}(t+\tau)\rangle}{d\tau}$
$\displaystyle=$
$\displaystyle\left(i\Delta_{0}-\frac{\Gamma}{2}\right)\langle\sigma_{m}^{\dagger}(t)\sigma_{j}(t+\tau)\rangle-\frac{i\Omega_{0}}{2}e^{i\mathbf{k}_{0}\cdot\mathbf{r}_{j}}\beta_{m}^{*}(t)-\frac{\Gamma}{2}\sum_{k\neq
j}\gamma_{jk}\langle\sigma_{m}^{\dagger}(t)\sigma_{k}(t+\tau)\rangle,$ (70)

where $\beta_{m}^{*}(t)$ is the solution of Eq. (21). Equation (70) is
integrated with the initial condition, at $\tau=0$,
$\langle\sigma_{m}^{\dagger}(t)\sigma_{j}(t)\rangle=\beta_{m}^{*}(t)\beta_{j}(t)$.
Let us note that, setting
$\langle\sigma_{m}^{\dagger}(t)\sigma_{j}(t+\tau)\rangle=\beta_{m}^{*}(t)G_{mj}(t,\tau)$,
Eq. (70) becomes

$\displaystyle\frac{dG_{mj}(t,\tau)}{d\tau}$ $\displaystyle=$
$\displaystyle\left(i\Delta_{0}-\frac{\Gamma}{2}\right)G_{mj}(t,\tau)-\frac{i\Omega_{0}}{2}e^{i\mathbf{k}_{0}\cdot\mathbf{r}_{j}}-\frac{\Gamma}{2}\sum_{k\neq
j}\gamma_{jk}G_{mk}(t,\tau).$ (71)

with $G_{mj}(t,0)=\beta_{j}(t)$. It is clear from (71) that $G_{mj}$ is
independent on $m$ and, comparing Eq.(71) with Eq.(21),
$G_{mj}(t,\tau)=\beta_{j}(t+\tau)$, so that

$\displaystyle g^{(1)}(t,\tau)$ $\displaystyle=$
$\displaystyle\frac{1}{I(t)}\overline{\sum_{j,m}e^{-i\mathbf{k}\cdot(\mathbf{r}_{j}-\mathbf{r}_{m})}\beta_{m}^{*}(t)\beta_{j}(t+\tau)}e^{-i\omega_{0}\tau}.$
(72)

In order to obtain $g^{(2)}(t,\tau)$, we need to evaluate
$g_{pmjq}(t,\tau)=\langle\sigma_{p}^{\dagger}(t)\sigma_{m}^{\dagger}(t+\tau)\sigma_{j}(t+\tau)\sigma_{q}(t)\rangle$.
First, we obtain the equation

$\displaystyle\frac{d}{d\tau}(\sigma^{\dagger}_{m}\sigma_{j})$
$\displaystyle=$
$\displaystyle-\Gamma\sigma^{\dagger}_{m}\sigma_{j}-\frac{i\Omega_{0}}{2}\left[e^{i\mathbf{k}_{0}\cdot\mathbf{r}_{j}}\sigma^{\dagger}_{m}-e^{-i\mathbf{k}_{0}\cdot\mathbf{r}_{m}}\sigma_{j}\right]-\frac{\Gamma}{2}\left(\sum_{k\neq
j}\gamma_{jk}\sigma^{\dagger}_{m}\sigma_{k}+\sum_{k\neq
m}\gamma^{*}_{mk}\sigma^{\dagger}_{k}\sigma_{j}\right).$ (73)

and then, from the quantum regression theorem,

$\displaystyle\frac{d}{d\tau}g_{pmjq}(t,\tau)$ $\displaystyle=$
$\displaystyle-\Gamma
g_{pmjq}(t,\tau)-\frac{i\Omega_{0}}{2}\left[e^{i\mathbf{k}_{0}\cdot\mathbf{r}_{j}}f^{*}_{qmp}(t,\tau)-e^{-i\mathbf{k}_{0}\cdot\mathbf{r}_{m}}f_{pjq}(t,\tau)\right]$
(74) $\displaystyle-$ $\displaystyle\frac{\Gamma}{2}\left[\sum_{k\neq
j}\gamma_{jk}g_{pmkq}(t,\tau)+\sum_{k\neq
m}\gamma^{*}_{mk}g_{pkjq}(t,\tau)\right],$

where
$f_{pjq}(t,\tau)=\langle\sigma_{p}^{\dagger}(t)\sigma_{j}(t+\tau)\sigma_{q}(t)\rangle$
are solutions of the equations:

$\displaystyle\frac{d}{d\tau}f_{pjq}(t,\tau)$ $\displaystyle=$
$\displaystyle\left(i\Delta_{0}-\frac{\Gamma}{2}\right)f_{pjq}(t,\tau)-\frac{i\Omega_{0}}{2}e^{i\mathbf{k}_{0}\cdot\mathbf{r}_{j}}\beta^{*}_{p}(t)\beta_{q}(t)-\frac{\Gamma}{2}\sum_{k\neq
j}\gamma_{jk}f_{pkq}(t,\tau).$ (75)

The initial conditions for Eqs.(74) and (75) are $g_{pmjq}(t,0)=0$ and
$f_{pjq}(t,0)=0$. By setting
$f_{pjq}(t,\tau)=\beta_{p}^{*}(t)\beta_{q}(t)F_{pjq}(t,\tau)$, Eq.(75) yields

$\displaystyle\frac{d}{d\tau}F_{pjq}(t,\tau)$ $\displaystyle=$
$\displaystyle\left(i\Delta_{0}-\frac{\Gamma}{2}\right)F_{pjq}(t,\tau)-\frac{i\Omega_{0}}{2}e^{i\mathbf{k}_{0}\cdot\mathbf{r}_{j}}-\frac{\Gamma}{2}\sum_{k\neq
j}\gamma_{jk}F_{pkq}(t,\tau).$ (76)

with $F_{pjq}(t,0)=0$. Hence $F_{pjq}(t,\tau)$ is independent on $p,q,t$ and,
comparing Eq. (76) with Eq. (21), $F_{pjq}(t,\tau)=\beta_{j}(\tau)$ and
$f_{pjq}(t,\tau)=\beta_{p}^{*}(t)\beta_{j}(\tau)\beta_{q}(t)$, so that Eq.
(74) becomes

$\displaystyle\frac{d}{d\tau}g_{pmjq}(t,\tau)$ $\displaystyle=$
$\displaystyle-\Gamma
g_{pmjq}(t,\tau)-\frac{i\Omega_{0}}{2}\beta_{p}^{*}(t)\left[e^{i\mathbf{k}_{0}\cdot\mathbf{r}_{j}}\beta^{*}_{m}(\tau)-e^{-i\mathbf{k}_{0}\cdot\mathbf{r}_{m}}\beta_{j}(\tau)\right]\beta_{q}(t)$
(77) $\displaystyle-$ $\displaystyle\frac{\Gamma}{2}\left[\sum_{k\neq
j}\gamma_{jk}g_{pmkq}(t,\tau)+\sum_{k\neq
m}\gamma^{*}_{mk}g_{pkjq}(t,\tau)\right].$

Setting $g_{pmjq}(t,\tau)=\beta_{p}^{*}(t)H_{pmjq}(t,\tau)\beta_{q}(t)$, Eq.
(77) becomes

$\displaystyle\frac{d}{d\tau}H_{pmjq}(t,\tau)$ $\displaystyle=$
$\displaystyle-\Gamma
H_{pmjq}(t,\tau)-\frac{i\Omega_{0}}{2}\left[e^{i\mathbf{k}_{0}\cdot\mathbf{r}_{j}}\beta^{*}_{m}(\tau)-e^{-i\mathbf{k}_{0}\cdot\mathbf{r}_{m}}\beta_{j}(\tau)\right]$
(78) $\displaystyle-$ $\displaystyle\frac{\Gamma}{2}\left[\sum_{k\neq
j}\gamma_{jk}H_{pmkq}(t,\tau)+\sum_{k\neq
m}\gamma^{*}_{mk}H_{pkjq}(t,\tau)\right]$

with $H_{pmjq}(t,0)=0$. It is evident from (78) that $H_{pmjq}(t,\tau)$ does
not depend on $p$, $q$ or $t$, i.e., $H_{pmjq}(t,\tau)=H_{mj}(\tau)$. Hence
the second-order correlation function is

$\displaystyle g^{(2)}(t,\tau)$ $\displaystyle=$
$\displaystyle\frac{1}{I^{2}(t)}\overline{\sum_{j,m,p,q}e^{-i\mathbf{k}\cdot(\mathbf{r}_{j}-\mathbf{r}_{m}+\mathbf{r}_{q}-\mathbf{r}_{p})}\beta_{p}^{*}(t)H_{mj}(\tau)\beta_{q}(t)}$
(79)

where $H_{mj}(\tau)$ is the solution of the equation

$\displaystyle\frac{d}{d\tau}H_{mj}(\tau)$ $\displaystyle=$
$\displaystyle-\Gamma
H_{mj}(\tau)-\frac{i\Omega_{0}}{2}\left[e^{i\mathbf{k}_{0}\cdot\mathbf{r}_{j}}\beta^{*}_{m}(\tau)-e^{-i\mathbf{k}_{0}\cdot\mathbf{r}_{m}}\beta_{j}(\tau)\right]-\frac{\Gamma}{2}\left[\sum_{k\neq
j}\gamma_{jk}H_{mk}(\tau)+\sum_{k\neq m}\gamma^{*}_{mk}H_{kj}(\tau)\right].$
(80)

with $H_{mj}(0)=0$.

##  Appendix B Evaluation of $g^{(2)}(\tau)$ for the product state (46)

In the chain of the derivation of the expectation values
$g_{pmjq}(t,\tau)=\langle\sigma_{p}^{\dagger}(t)\sigma_{m}^{\dagger}(t+\tau)\sigma_{j}(t+\tau)\sigma_{q}(t)\rangle$
with Eq.(74) and
$f_{pjq}(t,\tau)=\langle\sigma_{p}^{\dagger}(t)\sigma_{j}(t+\tau)\sigma_{q}(t)\rangle$
with Eq.(75), their initial conditions are

$g_{pmjq}(t,0)=\left\\{\begin{array}[]{ccc}\displaystyle\beta^{*}_{p}(t)\beta_{m}^{*}(t)\beta_{j}(t)\beta_{q}(t)&&\mathrm{if}\quad
p\neq m\quad\mathrm{or}\quad j\neq q\\\ &&\\\
0&&\mathrm{otherwise}\end{array}\right..$ (81)

and

$f_{pjq}(t,0)=\left\\{\begin{array}[]{ccc}\displaystyle\beta_{p}^{*}(t)\beta_{j}(t)\beta_{q}(t)&\mathrm{if}&j\neq
q\\\ &&\\\ 0&&\mathrm{otherwise}\end{array}\right..$ (82)

By setting $f_{pjq}(t,\tau)=\beta_{p}^{*}(t)\beta_{q}(t)F_{pjq}(t,\tau)$,
Eq.(75) yields

$\displaystyle\frac{d}{d\tau}F_{pjq}(t,\tau)$ $\displaystyle=$
$\displaystyle\left(i\Delta_{0}-\frac{\Gamma}{2}\right)F_{pjq}(t,\tau)-\frac{i\Omega_{0}}{2}e^{i\mathbf{k}_{0}\cdot\mathbf{r}_{j}}-\frac{\Gamma}{2}\sum_{k\neq
j}\gamma_{jk}F_{pkq}(t,\tau).$ (83)

with

$F_{pjq}(t,0)=\left\\{\begin{array}[]{ccc}\displaystyle\beta_{j}(t)&\mathrm{if}&j\neq
q\\\ &&\\\ 0&&\mathrm{otherwise}\end{array}\right..$ (84)

Hence $F_{pjq}(t,\tau)$ is independent on $p,q$,
$F_{pjq}(t,\tau)=F_{j}(t,\tau)$. Comparing Eq. (83) with Eq. (21), we deduce
that $F_{j}(t,\tau)=\beta_{j}(t+\tau)$ if $j\neq q$ and
$F_{j}(\tau)=\beta_{j}(\tau)$ if $j=q$. The same reasoning will be true for
$g_{pmjq}(t,\tau)$: from Eq.(74),

$\displaystyle\frac{d}{d\tau}g_{pmjq}(t,\tau)$ $\displaystyle=$
$\displaystyle-\Gamma
g_{pmjq}(t,\tau)-\frac{i\Omega_{0}}{2}\beta_{p}^{*}(t)\left[e^{i\mathbf{k}_{0}\cdot\mathbf{r}_{j}}F^{*}_{m}(t,\tau)-e^{-i\mathbf{k}_{0}\cdot\mathbf{r}_{m}}F_{j}(t,\tau)\right]\beta_{q}(t)$
(85) $\displaystyle-$ $\displaystyle\frac{\Gamma}{2}\left[\sum_{k\neq
j}\gamma_{jk}g_{pmkq}(t,\tau)+\sum_{k\neq
m}\gamma^{*}_{mk}g_{pkjq}(t,\tau)\right].$

Setting $g_{pmjq}(t,\tau)=\beta_{p}^{*}(t)H_{pmjq}(t,\tau)\beta_{q}(t)$, Eq.
(85) becomes

$\displaystyle\frac{d}{d\tau}H_{pmjq}(t,\tau)$ $\displaystyle=$
$\displaystyle-\Gamma
H_{pmjq}(t,\tau)-\frac{i\Omega_{0}}{2}\left[e^{i\mathbf{k}_{0}\cdot\mathbf{r}_{j}}F^{*}_{m}(t,\tau)-e^{-i\mathbf{k}_{0}\cdot\mathbf{r}_{m}}F_{j}(t,\tau)\right]$
(86) $\displaystyle-$ $\displaystyle\frac{\Gamma}{2}\left[\sum_{k\neq
j}\gamma_{jk}H_{pmkq}(t,\tau)+\sum_{k\neq
m}\gamma^{*}_{mk}H_{pkjq}(t,\tau)\right]$

with

$H_{mj}(t,0)=\left\\{\begin{array}[]{ccc}\displaystyle\beta_{m}^{*}(t)\beta_{j}(t)&&\mathrm{if}\quad
p\neq m\quad\mathrm{or}\quad j\neq q\\\ &&\\\
0&&\mathrm{otherwise}\end{array}\right..$ (87)

It is evident from (86) that $H_{pmjq}(t,\tau)$ depends on $p,q$ only for the
initial condition at $\tau=0$, so that $H_{pmjq}(t,\tau)=H_{mj}(t,\tau)$.
Hence, the second-order correlation function is

$\displaystyle g^{(2)}(t,\tau)$ $\displaystyle=$
$\displaystyle\frac{1}{I^{2}(t)}\overline{\sum_{j,m,p,q}e^{-i\mathbf{k}\cdot(\mathbf{r}_{j}-\mathbf{r}_{m}+\mathbf{r}_{q}-\mathbf{r}_{p})}\beta_{p}^{*}(t)H_{mj}(t,\tau)\beta_{q}(t)}$
(88)

where $H_{mj}(\tau)$ is the solution of the equation

$\displaystyle\frac{d}{d\tau}H_{mj}(t,\tau)$ $\displaystyle=$
$\displaystyle-\Gamma
H_{mj}(t,\tau)-\frac{i\Omega_{0}}{2}\left[e^{i\mathbf{k}_{0}\cdot\mathbf{r}_{j}}F^{*}_{m}(t,\tau)-e^{-i\mathbf{k}_{0}\cdot\mathbf{r}_{m}}F_{j}(t,\tau)\right]$
(89) $\displaystyle-$ $\displaystyle\frac{\Gamma}{2}\left[\sum_{k\neq
j}\gamma_{jk}H_{mk}(t,\tau)+\sum_{k\neq
m}\gamma^{*}_{mk}H_{kj}(t,\tau)\right].$

##  Appendix C Evaluation of $\langle E^{\dagger}(t)E(t)E(t)\rangle$ and
$\langle E^{\dagger}(t)E^{\dagger}(t)E(t)E(t)\rangle$ for the Eberly’s state
(53)

From the definition, we have

$\displaystyle\langle E^{\dagger}(t)E(t)E(t)\rangle$ $\displaystyle=$
$\displaystyle\frac{1}{N^{3}}\sum_{jmp}e^{i(\mathbf{k}_{0}-\mathbf{k})\cdot(\mathbf{r}_{m}+\mathbf{r}_{p}-\mathbf{r}_{j})}\langle\sigma^{\dagger}_{j}\sigma_{m}\sigma_{p}\rangle$
(90) $\displaystyle\langle E^{\dagger}(t)E^{\dagger}(t)E(t)E(t)\rangle$
$\displaystyle=$
$\displaystyle\frac{1}{N^{4}}\sum_{jmpq}e^{i(\mathbf{k}_{0}-\mathbf{k})\cdot(\mathbf{r}_{p}+\mathbf{r}_{q}-\mathbf{r}_{j}-\mathbf{r}_{m})}\langle\sigma^{\dagger}_{j}\sigma^{\dagger}_{m}\sigma_{p}\sigma_{q}\rangle$
(91)

and the not vanishing expectation values are

$\langle\sigma^{\dagger}_{j}\sigma_{m}\sigma_{p}\rangle=\beta|\beta|^{2}\left\\{\begin{array}[]{cc}\displaystyle\mathrm{if}&m=j\neq
p\\\ &\\\ \mathrm{if}&p=j\neq m\\\ &\\\ \mathrm{if}&p\neq j\neq
m\end{array}\right.$ (92)

and

$\langle\sigma^{\dagger}_{j}\sigma^{\dagger}_{m}\sigma_{p}\sigma_{q}\rangle=|\beta|^{4}\left\\{\begin{array}[]{cc}\displaystyle\mathrm{if}&j=p\neq
m=q\\\ &\\\ \mathrm{if}&j=q\neq m=p\\\ &\\\ \mathrm{if}&j=p,j\neq m,j\neq
q,m\neq q,\\\ &\\\ \mathrm{if}&j=q,j\neq m,j\neq p,p\neq m,\\\ &\\\
\mathrm{if}&m=p,m\neq q,m\neq j,q\neq j,\\\ &\\\ \mathrm{if}&m=q,p\neq m,p\neq
j,m\neq j,\\\ &\\\ \mathrm{if}&p\neq j\neq m\neq q\end{array}\right.$ (93)

From them we obtain:

$\displaystyle\langle E^{\dagger}(t)E(t)E(t)\rangle$ $\displaystyle=$
$\displaystyle\beta|\beta|^{2}\left\\{\frac{2(N-1)}{N^{2}}+K_{2}\right\\}$
(94) $\displaystyle\langle E^{\dagger}(t)E^{\dagger}(t)E(t)E(t)\rangle$
$\displaystyle=$
$\displaystyle|\beta|^{4}\left\\{\frac{2(N-1)}{N^{3}}+\frac{4(N-2)}{N^{2}}K_{2}+K_{4}\right\\}$
(95)

where

$\displaystyle K_{2}$ $\displaystyle=$
$\displaystyle\frac{1}{N^{2}}\sum_{j}\sum_{p\neq
j}e^{i(\mathbf{k}_{0}-\mathbf{k})\cdot(\mathbf{r}_{p}-\mathbf{r}_{j})}\approx|S_{\infty}|^{2}$
(96) $\displaystyle K_{4}$ $\displaystyle=$
$\displaystyle\frac{1}{N^{4}}\sum_{m}\sum_{j\neq m}\sum_{p\neq j,m}\sum_{q\neq
p,j,m}e^{i(\mathbf{k}_{0}-\mathbf{k})\cdot(\mathbf{r}_{p}+\mathbf{r}_{q}-\mathbf{r}_{j}-\mathbf{r}_{m})}\approx|S_{\infty}|^{4}.$
(97)

## References

  * [1] C. Cohen-Tannoudji, J. Dupont-Roc, G. Grynberg, Atom-Photon Interactions: Basic Processes and Applications, Wiley Ed., 1998.
  * [2] R. Hanbury Brown and R. Q. Twiss, Nature 178,1046 (1956).
  * [3] R.H. Dicke, Phys. Rev. 93, 99 (1954).
  * [4] J.H. Eberly and N.E. Rehler, Phys. Rev. A 2, 1607 (1970).
  * [5] E. Ressayre and A. Tallet, Phys. Rev. A5, 1572 (1972).
  * [6] D. Bhatti, R. Schneider, S. Oppel, and J. von Zanthier, Phys. Rev. Lett. 120, 113603 (2018).
  * [7] T. Bienaimé, S. Bux, E. Lucioni, Ph.W. Courteille, N. Piovella, and R. Kaiser, Phys. Rev. Lett. 104, 183602 (2010).
  * [8] W. Guerin, M. O. Araújo, and R. Kaiser, Phys. Rev. Lett. 116, 083601 (2016).
  * [9] M.O. Scully, E.S. Fry, C.H.R. Ooi, K. Wodd́iewicz, Phys. Rev. Lett. 2006, 96, 010501.
  * [10] A.A. Svidzinsky, J.T. Chang and M.O. Scully, Phys. Rev. Lett. 100, 160504 (2008).
  * [11] M.O. Scully, A.A. Svidzinsky, Science 325, 1510 (2009).
  * [12] J.H. Eberly, J. Phys. B: At. Mol. Opt. Phys. 39, S599 (2006).
  * [13] J. Chabé, M.T. Rouabah, L. Bellando, T. Bienaimé, N. Piovella, R. Bachelard, R. Kaiser, Physical Review A 89, 043833 (2014).
  * [14] R. Bachelard, N. Piovella, W. Guerin, R. Kaiser, Physical Review A 94, 033836 (2016).
  * [15] L. Pucci, A. Roy, T. Santiago do Espirito Santo, R. Kaiser, M. Kastner, and R. Bachelard Phys. Rev. A 95, 053625 (2017).
  * [16] R. Friedberg and J.T. Manassah, Laser Phys. Lett. 4 900 (2007).
  * [17] A.A. Svidzinsky, J.T. Chang, M.O. Scully, Phys. Rev.A 81, 053821 (2010).
  * [18] T. Bienaimé, M. Petruzzo, D. Bigerni, N. Piovella and R. Kaiser, Journal of Modern Optics 58, 1942 (2011).
  * [19] R. Friedberg, S.R. Hartmann, J.T. Manassah, Phys. Rep. 7, 101 (1973).
  * [20] R. Friedberg, J.T. Manassah, Phys. Lett. A 374, 1648 (2010).
  * [21] N.E. Rehler and J.H. Eberly, Phys. Rev. A 3, 1735 (1971).
  * [22] R. J. Glauber, Phys. Rev. 130, 2529 (1963).
  * [23] R. J. Glauber, Phys. Rev. 131, 2766 (1963).
  * [24] L. Mandel and E. Wolf, Optical Coherence and Quantum Optics (Cambridge University Press, New York, 1995), p. 714.
  * [25] R. Loudon, The Quantum Theory of Light, Oxford Univ. Press (2000), p. 333.
  * [26] S.M. Barnett, and P.M. Radmore, Methods in Theoretical Quantum Optics, Oxford Univ. Press (2002), p. 158.
  * [27] M.O. Scully, Laser Physics 17, 635 (2007).
  * [28] R.A. de Oliveira et al., Phys. Rev. A 90, 023848 (2014).
  * [29] F. Robicheaux, and Deepak A. Suresh, Phys. Rev. A 104, 023702 (2021).
  * [30] A. Eloy, Z. Yao, R. Bachelard, W. Guerin, M. Fouché, and R. Kaiser, Phys. Rev. A 97, 013810 (2018).