# Collective coupling of driven multilevel atoms and its effect on four-wave
mixing

P. Yanes-Thomas  Instituto de Física, Universidad Nacional Autónoma de México,
Ciudad de México 04510, México  R. Gutiérrez-Jáuregui  Instituto de Física,
Universidad Nacional Autónoma de México, Ciudad de México 04510, México  P.
Barberis-Blostein  Instituto de Investigaciones en Matemáticas Aplicadas y en
Sistemas, Universidad Nacional Autónoma de México, Ciudad de México 04510,
México  D. Sahagún-Sánchez  Instituto de Física, Universidad Nacional Autónoma
de México, Ciudad de México 04510, México  R. Jáuregui<EMAIL_ADDRESS>Instituto de Física, Universidad Nacional Autónoma de México, Ciudad de México
04510, México  A. Kunold<EMAIL_ADDRESS>Área de Física Teórica y Materia
Condensada, Universidad Autónoma Metropolitana Azcapotzalco, Ciudad de México
02200, México

###### Abstract

Microscopic models based on multilevel atoms are central to controlling non-
linear optical responses and coherent control of light. These models are
traditionally based on single-atom responses that are then parametrically
extrapolated to include collective effects, such as an enhanced response or
propagation within the medium. In this work we present a systematic analysis
of the cooperative effects arising in driven systems composed of multi-level
atoms coupled by a common electromagnetic environment. The analysis is based
on an interplay between dressed states induced by the driving field, photon
exchanges, and the collective decay channels. This theory is applied to the
case of four-wave mixing mediated by a pair of atoms in the diamond
configuration, a nonlinear process commonly used to create correlated photon
pairs. The effects of inter-atomic correlations and collective decay channels
over outgoing photons are then explored. We find that the resulting
correlation function changes from a Lorentzian distribution for independent
atoms to a two-peaked distribution when dipole-dipole interactions are
included. The two-peak structure is related to a joint action between the
dressed states and the collective decay channels. Signatures of these
processes are identified for existing experimental realizations. The intuition
obtained from this connection helps to uncover relevant parameters that could
be exploited for quantum control protocols based on dispersive and dissipative
cooperative effects in multi-level systems.

##  I Introduction

The optical response of an atomic medium is modified by external
electromagnetic fields. Intense laser fields, for example, have been shown to
alter this response directly by creating excitation pathways between
electronic states or indirectly by manipulating the atomic positions. The
excitation pathways affect both matter and its electromagnetic environment.
The resulting dressed states can be controlled in real time to create
transparency windows in otherwise opaque media [1], promote a non-linear
response [2, 3, 4] or tune the interaction between particles [5]. Control over
the atomic positions, by comparison, can be used to promote coherent
scattering between the atomic dipoles. Atoms inside dense ensembles display a
collective response that is manifested as altered scattering rates,
directional scattering or modified line shapes [6], to name a few examples. It
is this ability to fully control the atomic degrees-of-freedom with high
precision that has made atomic systems into a fundamental tool for creating
versatile light matter interfaces.

As experiments scale in both size and complexity, an increasingly rich
interplay between the control of external and internal degrees-of-freedom is
observed. In a regime characterized by strong fields and dense arrays, the
external fields can be used to induce collective pathways opened by coherent
scattering [7]. These pathways can be further exploited to create exotic
responses as non-reciprocal [8, 9, 10, 11] or topological materials [12, 13,
14]. They can also be used to explore the connection between macroscopic
properties of a material with microscopic models. A paradigmatic example of
the latter corresponds to the quantification of properties of the coherent
scattering into compact, macroscopic parameters related to the optical density
[15]. Such an approach extends to modern technologies, where these
possibilities place a premium on identifying the relevant variables for
controlling these pathways in a nonlinear setting.

We present here a systematic approach for describing the interplay between
collective and coherent effects in atomic media. For this, we consider an
ensemble of multi-level atoms driven by arbitrary laser fields. The dynamical
equations describing this system are a generalization of those obtained in
Ref. [16] for incorporating non-linear optical processes. As a specific
example, we consider the case of two four-level atoms driven by pump lasers.
This case provides the basic ingredients needed to describe non-linear and
collective effects in non-dilute atomic clouds, where pair-wise interactions
become prominent. The pump fields are shown to induce both inter- and intra-
atomic coherences. Intra-atomic coherences lead to interference between the
excitation pathways of individual atoms. Inter-atomic coherences, in turn,
extend this effect to collective pathways created from the interaction between
different atoms. We show that these collective pathways manifest as a self-
seeding effect, an apparent line-splitting in the atomic transitions and a
change in frequency of the outgoing photons. This splitting can not be
obtained using standard macroscopic models for multi-level atoms. Our
systematic approach is directed towards identifying the relevant variables to
control these pathways in a more general setting, with potential applications
in quantum control and engineering.

The paper is organized as follows. In Section II we introduce a microscopic
model for an atomic ensemble driven by external laser fields. Dynamical
equations are obtained for the atomic operators with emphasis on the dressed
states induced by the laser fields and how they are affected by dissipation.
Then, in Sec. III, we treat the explicit example of two four-level atoms using
experimental parameters for cold atomic clouds. An analysis of how this
collective response is imprinted on the outgoing photons is presented in Sec.
IV. We identify identify the relevant variables to control the properties of
the outgoing photons working from the dressed states. Section V is left for
summary and conclusions, where we identify signatures of collective processes
in four-wave mixing that can be found with current experimental reazlizations.

##  II Multi-level emitters coupled to a common environment

Consider an ensemble of $n_{a}$ identical, non-overlapping emitters that are
coupled to a common electromagnetic environment. Each emitter is characterized
by its position $\boldsymbol{r}_{\alpha}$, with
$\alpha=\\{1,2,\dots,n_{a}\\}$, and is composed of $n_{l}$ internal levels
$\ket{m}_{\alpha}$ of corresponding energies $\varepsilon_{m}$, with
$m=\\{0,1,\dots,n_{l}-1\\}$. Their common environment is described by the
electric field operator $\mathbf{E}$. When decomposed into a plane-wave
expansion, this operator reads

$\mathbf{E}(\mathbf{r},t)=\sum_{q}\mathcal{E}_{q}\left(\boldsymbol{e}_{q}\mathrm{e}^{i\boldsymbol{k}_{q}\cdot\boldsymbol{r}}a_{q}(t)+\boldsymbol{e}^{*}_{q}\mathrm{e}^{-i\boldsymbol{k}_{q}\cdot\boldsymbol{r}}a_{q}^{\dagger}(t)\right)\,,$
(1)

where $\mathcal{E}_{q}=\sqrt{\hbar\omega_{q}/2\epsilon_{0}\mathcal{V}}$ is the
vacuum field amplitude, $\mathcal{V}$ the spatial volume confining the field,
and creation $a_{q}^{\dagger}$ and annihilation operators $a_{q}$ obey the
commutation relation $[a_{q},a_{q^{\prime}}^{\dagger}]=\delta_{q,q^{\prime}}$.
The plane waves are characterized by a wave vector $\boldsymbol{k}_{q}$,
angular frequency
$\omega_{q}=c\sqrt{\boldsymbol{k}_{q}\cdot\boldsymbol{k}_{q}}$, and
polarization unit vector $\boldsymbol{e}_{q}$. For simplicity, we encompass
all these traits within an index $q$.

The emitters are assumed to move slowly, such that no appreciable changes
occur in their position in the time required to modify their collective
internal states. We then remove this motion from our analysis and describe
this light-matter system by the Hamiltonian

$\mathcal{H}=\mathcal{H}_{A}+\mathcal{H}_{P}+\mathcal{H}_{AP}\,,$ (2)

where

$\displaystyle\mathcal{H}_{A}$ $\displaystyle=$
$\displaystyle\sum_{\alpha=1}^{n_{a}}\sum_{m=1}^{n_{l}}\varepsilon_{m}\sigma_{m,m}^{\alpha}\,,$
(3) $\displaystyle\mathcal{H}_{P}$ $\displaystyle=$
$\displaystyle\sum_{q}\hbar\omega_{q}a_{q}^{\dagger}a_{q}\,,$ (4)

are free terms for emitters and field, respectively. The operators
$\sigma_{mn}^{\alpha}=\ket{m}_{\alpha}\leftindex_{\alpha}{\bra{n}}$ signal the
transitions within internal levels and thus commute for any two distinct
emitters $(\alpha\neq\beta)$:

$\left[\sigma_{mn}^{\alpha},\sigma_{m^{\prime}n^{\prime}}^{\beta}\right]=\left[\sigma_{mn}^{\alpha\dagger},\sigma_{m^{\prime}n^{\prime}}^{\beta}\right]=\left[\sigma_{mn}^{\alpha\dagger},\sigma_{m^{\prime}n^{\prime}}^{\beta\dagger}\right]=0.$

The last term of Eq. (2) describes the coupling between field and emitters.
This coupling is dipolar in character, with each emitter probing its local
electric field via the dipole operator $\mathbf{p}$. The interaction is given
by

$\mathcal{H}_{AP}=-\sum_{\alpha=1}^{n_{a}}\sum_{n,m=0}^{n_{l}-1}\mathbf{E}(\mathbf{r}_{\alpha})\cdot\boldsymbol{s}_{mn}^{\alpha}\,,$
(5)

where both rotating and counter-rotating terms are included via

$\boldsymbol{s}_{mn}^{\alpha}(t)=\boldsymbol{p}_{mn}\sigma_{mn}^{\alpha}+\boldsymbol{p}_{mn}^{*}\sigma_{mn}^{\alpha\dagger}.$
(6)

with $\boldsymbol{p}_{mn}=\left\langle m\left|\mathbf{p}\right|n\right\rangle$
the matrix element connecting two internal states.

###  II.1 Dynamical equations

The field mediates the interaction between the emitters and introduces
boundary conditions via external drives and collective decay channels. We now
include these effects over arbitrary emitter operators.

We begin from the Heisenberg equations for the field operators $a_{q}(t)$.
Under the action of $\mathcal{H}$, these operators evolve as

$\dot{a}_{q}=-i\omega_{q}a_{q}+\frac{i}{\hbar}\mathcal{E}_{q}\sum_{\alpha}\sum_{n,m}\boldsymbol{e}^{*}_{q}\cdot\boldsymbol{s}_{mn}^{\alpha}\mathrm{e}^{-i\boldsymbol{k}_{q}\cdot\boldsymbol{r}_{\alpha}},$
(7)

with formal solution

$\displaystyle a_{q}(t$ $\displaystyle)=a_{q}(0)\mathrm{e}^{-i\omega_{q}t}$
$\displaystyle+\frac{i}{\hbar}\mathcal{E}_{q}\sum_{\alpha,n,m}\int_{0}^{t}dt^{\prime}\mathrm{e}^{i\omega_{q}(t^{\prime}-t)-i\boldsymbol{k}_{q}\cdot\boldsymbol{r}_{\alpha}}\boldsymbol{e}^{*}_{q}\cdot\boldsymbol{s}_{mn}^{\alpha}(t^{\prime}).$
(8)

When substituted into Eq. (1), this solution naturally divides the electric
field into free and scattered terms
$\mathbf{E}=\mathbf{E}_{\text{f}}+\mathbf{E}_{\text{s}}$. Their positive
frequency components are

$\displaystyle\mathbf{E}^{(+)}_{\text{f}}(\mathbf{r},t)=\sum_{q}\mathcal{E}_{q}\boldsymbol{e}_{q}\mathrm{e}^{i\left(\boldsymbol{k}_{q}\cdot\boldsymbol{r}-\omega_{q}t\right)}a_{q}(0)\,,$
(9)
$\displaystyle\mathbf{E}^{(+)}_{\text{s}}(\mathbf{r},t)=\frac{i}{\hbar}\sum_{q}\mathcal{E}_{q}^{2}\sum_{\alpha,n,m}\mathrm{e}^{i\boldsymbol{k}_{q}\cdot(\boldsymbol{r}-\boldsymbol{r}_{\alpha})}\boldsymbol{e}_{q}\,$
$\displaystyle\times\int_{0}^{t}dt^{\prime}\mathrm{e}^{i\omega_{q}(t^{\prime}-t)}\boldsymbol{e}^{*}_{q}\cdot\boldsymbol{s}_{mn}^{\alpha}(t^{\prime}).$
(10)

The scattered field describes the field radiated by the emitter ensemble. In
free-space, where the ensemble couples to a large number of modes, the sum
over $q$ turns into an integral over all possible frequencies $d\omega$ and
propagation angles $d\Omega_{\mathbf{k}}$, and a sum over polarization vectors
$e_{q;\mathbf{k}}$. If the coupling between emitters and field is weak enough,
the coupling strength and mode densities can be evaluated at the resonance
frequency

$\Delta_{mn}=(\varepsilon_{n}-\varepsilon_{m})/\hbar\,.$ (11)

for each transition. Under this first Markov approximation [17], the integrals
are readily evaluated to give

$\displaystyle\mathbf{E}^{(+)}_{s}(\mathbf{r}_{\beta},t)=\sum_{\alpha}\sum_{\\{n,m\\}}$
$\displaystyle\left(\frac{\Delta_{mn}^{3}}{4\pi\epsilon_{0}c^{3}}\right)\overleftrightarrow{\mathbf{F}}^{(+)}(\mathbf{r}_{\alpha\beta},\Delta_{mn})\cdot\mathbf{p}_{n,m}$
$\displaystyle\times$
$\displaystyle\sigma_{mn}^{\alpha}(t-t_{\alpha\beta})\Uptheta(t-t_{\alpha\beta})\,,$
(12)

where the sum $\\{n,m\\}$ is taken over positive frequencies $\Delta_{mn}$
only. Here, $ct_{\alpha\beta}=|\mathbf{r}_{\alpha\beta}|$ describes the time
delay between the emission of a photon at $\mathbf{r}_{\alpha}$ and its
detection at $\mathbf{r}_{\beta}$, with
$\mathbf{r}_{\alpha\beta}=\mathbf{r}_{\alpha}-\mathbf{r}_{\beta}$, and
$\Uptheta$ the Heaviside function. The tensor

$\overleftrightarrow{\mathbf{F}}^{(+)}(\mathbf{r},\omega)=\frac{1}{\xi^{3}}\left[(\xi^{2}+i\xi-1)\mathbbm{1}-(\xi^{2}+3i\xi-3)\frac{\mathbf{r}\otimes\mathbf{r}}{r^{2}}\right]$
(13)

accounts for the propagation of the scattered field with
$\xi=\omega|\mathbf{r}|/c$.

By evaluating source and field operators at different times, Eq. (II.1)
describes the evolution of an emitter as its output field propagates towards a
target. The spatial separation between source and target is effectively
removed by transforming into an interaction picture where the sources are
retarded. For a small ensemble, where the only changes in the emitter as the
photon propagates from one end of the ensemble to the other are given by its
free evolution, this retardation only results in a local phase [18]. The
scattered field under this assumption is

$\displaystyle\mathbf{E}^{(+)}_{s}(\mathbf{r}_{\beta},t)=\sum_{\alpha}\sum_{\\{n,m\\}}$
$\displaystyle\left(\frac{\Delta_{mn}^{3}}{4\pi\epsilon_{0}c^{3}}\right)\overleftrightarrow{\mathbf{F}}^{(+)}(\mathbf{r}_{\alpha\beta},\Delta_{mn})\cdot\mathbf{p}_{nm}$
$\displaystyle\times$
$\displaystyle\mathrm{e}^{i\kappa_{nm}{r}_{\alpha\beta}}\sigma_{mn}^{\alpha}(t)\Uptheta(t)\,,$
(14)

where $\kappa_{mn}=\Delta_{mn}/c$ is the transition wave-number.

We can now describe the evolution of an arbitrary emitter operator $Q(t)$
coupled to this common environment. Under the action of $\mathcal{H}$ and
within the rotating wave approximation these operators evolve as

$\displaystyle\dot{Q}$
$\displaystyle=\frac{i}{\hbar}\left[\mathcal{H}_{A},Q\right]$
$\displaystyle+\frac{i}{\hbar}\sum_{\alpha}\sum_{\\{n,m\\}}\left[\sigma_{mn}^{\alpha},Q\right]\mathbf{p}_{mn}\cdot\mathbf{E}^{(+)}_{\text{f}}(r_{\alpha})$
$\displaystyle\qquad-\mathbf{p}_{mn}^{*}\cdot\mathbf{E}^{(-)}_{\text{f}}(r_{\alpha})\left[Q,\sigma_{mn}^{\alpha\dagger}\right]$
$\displaystyle+\frac{i}{\hbar}\sum_{\alpha}\sum_{\\{n,m\\}}\left[\sigma_{mn}^{\alpha},Q\right]\mathbf{p}_{mn}\cdot\mathbf{E}^{(+)}_{\text{s}}(r_{\alpha})$
$\displaystyle\qquad-\mathbf{p}_{mn}^{*}\cdot\mathbf{E}^{(-)}_{\text{s}}(r_{\alpha})\left[Q,\sigma_{mn}^{\alpha\dagger}\right]$
(15)

where all operators are evaluated at time $t$. For simplicity, we consider
that the free internal energies $\epsilon_{m}$ are already Lamb shifted [19].

By inserting the fields from Eqs. (9-10) into Eq. (II.1), the dynamical
equations are obtained following the standard approach [20, 17]. Here, terms
that involve operators $\sigma_{mn}^{\alpha}$ and
$\sigma_{m^{\prime}n^{\prime}}^{\beta}$ with
$\\{m,n\\}\neq\\{m^{\prime},n^{\prime}\\}$ are removed unless the separation
between their transition frequencies,
$||\Delta_{mn}|-|\Delta_{m^{\prime}n^{\prime}}||$, is smaller than the natural
linewidths of the relevant internal levels. After removing these fast
oscillating terms and performing the summations explicitly, Eq. (II.1) becomes

$\dot{Q}=-(i\hbar)^{-1}\left[\mathcal{H}_{A}+\mathcal{H}_{E}+\mathcal{H}_{dd},Q\right]+\mathcal{L}\left[Q(t)\right],$
(16)

where free-field effects and dipole-dipole interactions occur via

$\displaystyle\mathcal{H}_{E}$
$\displaystyle=-\sum_{\alpha,n,m}{\sigma}_{mn}^{\alpha}(t)\mathbf{p}_{mn}\cdot\mathbf{E}^{(+)}_{\text{f}}(\mathbf{r}_{\alpha},t)+\text{H.c.}\,,$
(17) $\displaystyle\mathcal{H}_{dd}$
$\displaystyle=\sum_{\alpha\neq\beta}\sum_{m,n}{\sum_{m^{\prime},n^{\prime}}}\hbar\Omega^{\alpha\beta}(m,n;m^{\prime},n^{\prime})\sigma_{mn}^{\alpha\dagger}(t)\sigma_{m^{\prime}n^{\prime}}^{\beta}(t),$
(18)

while damping is included via the superoperator

$\mathcal{L}\left[\bullet\right]=\tfrac{1}{2}\sum_{\alpha\beta}{\sum_{mn}}{\sum_{m^{\prime}n^{\prime}}}\gamma^{\alpha\beta}(m,n;m^{\prime},n^{\prime})\mathfrak{L}_{mn;m^{\prime}n^{\prime}}^{\alpha\beta}\left[\bullet\right]$
(19)

with

$\mathfrak{L}_{mn;m^{\prime}n^{\prime}}^{\alpha\beta}\left[\bullet\right]=2\sigma_{mn}^{\alpha\dagger}\bullet\sigma_{m^{\prime},n^{\prime}}^{\beta}-\left\\{\sigma_{mn}^{\alpha\dagger}\sigma_{m^{\prime}n^{\prime}}^{\beta}\bullet\right\\}\,.$
(20)

The dipole-dipole coupling $\Omega^{\alpha\beta}$ and decay rates
$\gamma^{\alpha\beta}$ are related to the imaginary and real parts of
$\mathbf{F}^{(\pm)}$ by

$\displaystyle\Omega^{\alpha\beta}(m,n;m^{\prime},n^{\prime})+i\tfrac{1}{2}\gamma^{\alpha,\beta}(m,n;m^{\prime},n^{\prime})=$
$\displaystyle\tfrac{3}{4}\Gamma_{mn;m^{\prime}n^{\prime}}e^{i\kappa_{mn}r_{\alpha\beta}}\hat{\mathbf{p}}_{mn}^{*}\cdot\overleftrightarrow{\mathbf{F}}(\mathbf{r}_{\alpha\beta},\Delta_{m^{\prime}n^{\prime}})\cdot\hat{\mathbf{p}}_{m^{\prime}n^{\prime}}\,,$
(21)

with

$\Gamma_{mn;m^{\prime}n^{\prime}}=\left(\frac{1}{4\pi\epsilon_{0}}\right)\frac{4\Delta^{3}_{mn}\left|\mathbf{p}_{mn}^{*}\right|\left|\mathbf{p}_{m^{\prime}n^{\prime}}\right|}{3\hbar
c^{3}}$ (22)

When the electric dipole moments coincide, $\Gamma_{mn;mn}=:\Gamma_{mn}$, this
rate corresponds to the spontaneous decay rate between levels $\\{m,n\\}$ of
an individual atom.

Equation (16) describes the evolution of an arbitrary operator of this
ensemble of multi-level emitters. The superoperator $\mathcal{L}$ introduces
damping into the system through collective decay channels formed from the
interference of individual decay paths and obtained by diagonalizing the
$\gamma$ matrix of Eq. (19). Destructive interference gives way to subradiant
channels whose eigenvalue is smaller than the corresponding individual decay
rate. Constructive interference, in turn, gives way to superradiant channels
whose eigenvalues are larger than the individual decay rate. The multi-level
structure is made apparent at transitions of similar energies. By allowing for
these transitions to interfere, we open the door for coherences among
different levels to arise from spontaneous emission [21].

The relative phases that give rise to these sub- and superradiant channels are
usually washed away by the dipole-dipole couplings and free fields of Eqs.
(17,18). As each atom probes a different local environment, this dephasing
usually mixes the subradiant and superradiant channels of the system. In the
absence of free fields, this mixing can be understood from the matrices
$\gamma$ and $\Omega$ not commuting with one another. The underlying physical
problem can then be better understood in terms of the basis that diagonalizes
each one of these matrices.

###  II.2 Dressed states and collective decays

We have placed emphasis on the collective frequency shifts
$\Omega^{\alpha\beta}$ and decay channels $\gamma^{\alpha\beta}$ introduced by
the electromagnetic environment. The environment, however, also carries
external fields via the free modes of Eq. (17). These input fields— along with
the photons that are coherently interchanged by the emitters—can dress the
emitter-photons system [22, 23] or alter their decays rates [24, 25] depending
on their nature. We finish this section by introducing these coherent fields.

When driven by a coherent source, such as a laser, the operators $\hat{a}_{q}$
in Eq. (9) can be replaced by a coherent state amplitude $\alpha_{q}$ [26].
Their effect on the emitters is accounted for by Eq. (17), which now reads

$\mathcal{H}_{E}(t)=-\hbar\sum_{\alpha}\sum_{mn}\sigma_{mn}^{\alpha}\Lambda_{mn}\mathrm{e}^{i\boldsymbol{k}\cdot\boldsymbol{r}_{\alpha}}\mathrm{e}^{-i\omega_{\ell}t}+\text{H.c.}$
(23)

with $\omega_{\ell}$ the laser frequency and
$\Lambda_{mn}=\mathcal{E}_{\ell}\alpha_{\ell}\mathbf{e}_{\ell}\cdot\mathbf{p}_{mn}/\hbar$
the Rabi frequencies characterizing the coupling. The driving fields couple
predominantly to states whose transition frequency is near the laser frequency
$\omega_{\ell}\simeq\Delta_{mn}$, thus allowing for the sum in Eq. (23) to be
restricted to these quasi-resonant transitions.

The periodic oscillation of the driving field in a multilevel system can be
removed working in the interaction picture with respect to the Hamiltonian

$\displaystyle\mathcal{H}_{0}$ $\displaystyle=$
$\displaystyle\sum_{\alpha}\sum_{mn}\tfrac{1}{2}\hbar\varphi_{nm}\left(\sigma_{m,m}^{\alpha}-\sigma_{n,n}^{\alpha}\right).$
(24)

where the frequencies $\varphi_{mn}$ are chosen to cancel out the laser
frequencies. The explicit time dependence of the total Hamiltonian and the
decay super-operators is eliminated in this picture, and the dynamical Eqs.
(16) take the form

$\dot{\tilde{Q}}(t)=\frac{i}{\hbar}\left[\mathcal{H}_{sys},\tilde{Q}(t)\right]+\mathcal{L}\left[\tilde{Q}(t)\right],$
(25)

where $\tilde{Q}=\mathcal{R}(t)Q\mathcal{R}^{\dagger}(t)$ with
$\mathcal{R}(t)=\text{exp}[i\mathcal{H}_{0}t/\hbar]$ are the emitter operators
in the interaction picture and

$\mathcal{H}_{sys}=\tilde{\mathcal{H}}_{A}-\tilde{\mathcal{H}}_{0}+\tilde{\mathcal{H}}_{dd}+\tilde{\mathcal{H}}_{E}\,.$

includes both the effect of driving fields and dipole-dipole couplings.

The states that diagonalize $\mathcal{H}_{sys}$ are the collective dressed
states. They are formed from the superposition of bare states, containing
phases imprinted onto the emitters and their distribution across the ensemble.
For intense driving fields, these dressed states provide a natural picture for
describing this collective system that complements the insight acquired from
the collective decay channels. The latter are obtained by diagonalizing the
decay tensor $\gamma^{\alpha\beta}$ that defines the strength of the
contribution of each super-operator
$\mathcal{L}^{\alpha\beta}_{mn;m^{\prime}n^{\prime}}$ to the super-operator
$\mathcal{L}$, Eq. (19). The collective decay channels are the main focus in
spontaneous decay experiments.

The following Sections are concerned with a particular example where the
interplay between dressed states and collective decay channels alters the
optical response of an atomic system. We explore this response by solving Eq.
(16) using a general algorithm described in the Appendix. The physical
picture, however, is enriched by exploring the dressed states and their
collective decays.

##  III Two four level atoms

We now explore the interplay between dressed states and collective decay for
two atoms, $\alpha\in\\{1,2\\}$, with four-levels each, $m\in\\{0,1,2,3\\}$.
Their internal levels are distributed in the diamond configuration shown in
Fig. 1, where states $|0\rangle$ and $|2\rangle$ are connected by two paths: a
two-photon transition along $|1\rangle$ and a cascaded decay through
$|3\rangle$. The two-photon transition is induced by a pair of co-propagating
beams. These beams have orthogonal polarizations set to match the dipole
moment of their respective atomic transitions. Their intensity leads to the
Rabi frequencies $\Lambda_{01}$ and $\Lambda_{12}$ and are detuned for the
single- and two-photon transitions by $\delta_{1}$ and $\delta_{2}$,
respectively. Each atom interacts with its neighbor via the dipole-dipole
coupling. We have selected the case of two atoms since it presents the
simplest case where collective effects are already made apparent. This pair
interaction can then be used as a building block for larger order effects.

Figure 1: Schematics and energy level diagram for two driven, four-level atoms
($\alpha$ and $\beta$) in the diamond configuration. Two pump beams orthogonal
to the atomic plane couple states $|0\rangle$, $|1\rangle$ and $|2\rangle$ via
a two photon transition (straight arrows) characterized by Rabi frequencies
$\Lambda_{n,m}$ and detuned from the single- and two-photon transitions by
$\delta_{1}$ and $\delta_{2}$, respectively. De-excitation occurs via
spontaneous emission (curled arrows) with individual decay rate
$\gamma_{n,m}$. The atoms interact with one another via the environment, which
induces a dipolar coupling and collective decay at respective rates
$\Omega^{\alpha\beta}_{i,j}$ and $\gamma^{\alpha\beta}_{i,j}$ (pink arrows).
Blue and red lines denote two different polarizations.

To be more specific, we have in mind four states of atomic 87Rb:
$\\{5S_{1/2},5P_{3/2},5D_{3/2},5P_{1/2}\\}$, which have been used to implement
quantum control protocols and create squeezed light in alkali atomic gases
[27, 28, 29]. These states correspond, respectively, to
$\\{\ket{0}_{\alpha},\ket{1}_{\alpha},\ket{2}_{\alpha},\ket{3}_{\alpha}\\}$ in
our notation. They are separated by the transition wavelengths
$\\{\lambda_{01},\lambda_{12},\lambda_{23},\lambda_{30}\\}=\\{780,776,762,795\\}\,\text{nm}$
and are connected radiatively by the individual decay rates
$\\{\Gamma_{0,1},\Gamma_{1,2},\Gamma_{2,3},\Gamma_{0,3}\\}=\\{36.2,0.641,2.86,17.2\\}$
MHz. Since the transition energies differ significantly compared to the decay
rates, cooperative effects are dominated by transitions between the same
atomic levels, as denoted by pink arrows in Fig. 1. The parameters for the
pumping beams are taken from recent experiments with atomic clouds [30, 31,
32]. The atoms are then driven by two pump lasers (blue and red arrows in the
Figure) of powers $P_{1}=0.1\,$mW and $P_{2}=4\,$mW with laser spot diameter
of $1.1\,$mm and variable detunings. When joined by the dipole matrix elements
[33]
$\\{p_{0,1},p_{1,2},p_{2,3},p_{0,3}\\}=\\{5.956,0.787,1.616,4.11\\}\,ea_{0}$
these beams give way to the Rabi frequencies $\Lambda_{01}=7.5$ MHz and
$\Lambda_{12}=6.3$ MHz. The dipole moments are oriented such that
$\hat{p}_{0,1}=\hat{p}_{2,3}$ and $\hat{p}_{3,0}=\hat{p}_{1,2}$ while
$\hat{p}_{0,1}^{*}\cdot\hat{p}_{1,2}=0$ (to guarantee a closed loop transition
we take $\hat{p}_{m,n}$ to be circularly polarized). In writing the Rabi
frequencies we have further considered that the two atoms are located in a
plane perpendicular to the laser propagation vector, thus allowing for their
frequencies to be in phase.

We consider important to point out two key differences between our description
below and the standard approach used to describe the light generated by this
system [34]. First, we maintain all states throughout the description. The
state $|1\rangle_{\alpha}$ is usually adiabatically eliminated for the large
detunings $\delta_{1}$ considered here [35, 36]. We show that, once dressed by
the pump beams and dipole-dipole couplings, this state leads to a coupled
state whose frequency is brought close to resonance to other relevant states,
thus affecting the optical response of the system [Sec. III.2]. Second, we
make explicit that there is no external beam connecting $|3\rangle$ to the
other states. It has been shown that the population inversion between levels
$|2\rangle$ and $|3\rangle$ induced by the two-photon transition leads to
lasing scheme in a many-atom system [37]. The generated field can then be used
to generate an efficient frequency up-conversion when two infrared low power
laser fields are converted spontaneously into infrared and blue radiation in a
rubidium gas [27]. Usually, the transition is closed by a weak seeding-field
[28, 29]. We are interested in exploring how atom-atom correlations give rise
to a self-seeding field, i.e., a four wave mixing effect without an external,
weak pump field coupling the $|3\rangle$ and $|0\rangle$ levels to enhance the
cascade transition $|2\rangle\rightarrow|3\rangle\rightarrow|0\rangle$.

Figure 2: Dressed-state energies for a single atom in the diamond
configuration as as a function of the two photon detuning $\delta_{2}$ (solid
lines). The existence of several levels causes the bare energies (dashed
lines) to be shifted from the dressed states while far from resonance. Here,
$\delta_{1}=-70$ MHz with the remaining parameters given in the main text.
Figure 3: Dressed states energies for a pair of four-level atoms as a function
of the inverse inter-atomic distance. At large distances the states approach
the symmetric (solid lines) and antisymmetric (dashed) states formed from the
individual dressed states of Fig. 1 [see Eqs. (28)]. At short interatomic
distances, dipole-dipole interactions couple these asymptotic states and bring
several lines close to resonance as shown near $50$MHz. The parameters follow
from Fig. 2 at two photon resonance $\delta_{2}=0$.

###  III.1 Individual dressed states

We begin by reviewing the dressed states of an individual atom and their
dependence on the external parameters in the absence of decay. To do so, we
move into the interaction picture and carry out a complementary transformation
to a rotating frame according to Eqs. (24-25) to obtain a representation of
the Hamiltonian that is time independent. This transformation involves
rotation phases determined by the laser detunings and given by

$\displaystyle\varphi_{01}=\frac{2}{3}(\delta_{1}+\delta_{2})\,$ ,
$\displaystyle\varphi_{12}=\frac{2}{3}(2\delta_{2}-\delta_{1}).$ (26)

The resulting time-independent Hamiltonian is

$\mathcal{H}_{sys}=\left(\begin{array}[]{cccc}-\frac{1}{3}\left(\delta_{1}+\delta_{2}\right)&\Lambda_{0,1}&0&0\\\
\Lambda_{0,1}&\frac{1}{3}\left(2\delta_{1}-\delta_{2}\right)&\Lambda_{1,2}&0\\\
0&\Lambda_{1,2}&\frac{1}{3}\left(-\delta_{1}+2\delta_{2}\right)&0\\\
0&0&0&0\\\ \end{array}\right).$ (27)

Figure 2 shows the dressed states $|a\rangle,|b\rangle,|c\rangle,|d\rangle$
that diagonalize $\mathcal{H}_{sys}$ as function of the two photon detuning
$\delta_{2}$. The results suggest that the bare states $|0\rangle$ and
$|2\rangle$ (dashed lines) couple via a two-photon transition to create states
$|b\rangle$ and $|c\rangle$. These dressed states display an anticrossing when
the two-photon transition is brought to resonance $(\delta_{2}=0)$. The
anticrossing is related to the Autler-Townes doublet [38, 39, 23] measured in
experimental settings. By comparison, state $|d\rangle$ remains mostly
decoupled from these two states due to the large single-phoon detuning
$\delta_{1}=-70$ MHz. The state $|a\rangle$ corresponds to the $|3\rangle$ and
is completely decoupled from the other states in the absence of spontaneous
decay. We use the energy of this state to define a zero of energy.

The dressed states can be obtained in closed form nearby this two-photon
resonance. When written within the bare basis, they read

$\displaystyle\left|a\right\rangle=$ $\displaystyle\left|3\right\rangle,$
(28a) $\displaystyle\left|b\right\rangle=$
$\displaystyle\mathcal{N}_{b}^{-1}\Big{[}\Lambda_{1,2}\ket{0}-\Lambda_{0,1}\ket{2}\Big{]},$
(28b) $\displaystyle\left|c\right\rangle=$
$\displaystyle\mathcal{N}_{c}^{-1}\Big{[}{\Lambda_{0,1}}\ket{0}+\tfrac{1}{4}({\delta_{1}+\zeta})\ket{1}+{\Lambda_{1,2}}\ket{2}\Big{]},$
(28c) $\displaystyle\left|d\right\rangle=$
$\displaystyle\mathcal{N}_{d}^{-1}\Big{[}{\Lambda_{0,1}}\ket{0}-\tfrac{1}{4}({\delta_{1}+\zeta})\ket{1}+{\Lambda_{1,2}}\ket{2}\Big{]}],$
(28d)

where $\mathcal{N}_{b}$, $\mathcal{N}_{c}$ and $\mathcal{N}_{d}$ are
normalization constants. The corresponding eigenvalues are

$\displaystyle\varepsilon_{a}=0,$
$\displaystyle\varepsilon_{b}=-\delta_{1}/3,$ (29a)
$\displaystyle\varepsilon_{c}=(\delta_{1}+3\zeta)/6,$
$\displaystyle\varepsilon_{d}=(\delta_{1}-3\zeta)/6,$ (29b)

where $\zeta=\sqrt{\delta_{1}^{2}+16\Lambda_{0,1}^{2}+16\Lambda_{1,2}^{2}}$.
From Eqs. (28) it is shown that the state $|b\rangle$ remains dark to an
experiment probing the $|1\rangle$ state. And, when
$\Lambda_{0,1}\approx\Lambda_{1,2}$ and
$\left|\delta_{1}\right|\gg\Lambda_{0,1}$ the state $\ket{1}$ can be
adiabatically eliminated to give

$\displaystyle\left|b\right\rangle$ $\displaystyle\approx$
$\displaystyle\frac{1}{\sqrt{2}}\Big{[}\ket{0}-\left|2\right\rangle\Big{]},$
(30a) $\displaystyle\left|c\right\rangle$ $\displaystyle\approx$
$\displaystyle\frac{1}{\sqrt{2}}\Big{[}\ket{0}+\left|2\right\rangle\Big{]},$
(30b) $\displaystyle\left|d\right\rangle$ $\displaystyle\approx$
$\displaystyle\left|1\right\rangle.$ (30c)

Figure 4: Energies of the collective dressed states of two atoms as a function
of the inverse of the inter-atomic distance for $\delta_{1}=-70$ MHz and
$\delta_{2}=0$ MHz. Traces of the oscillatory behavior of $\Omega^{1,2}_{mn}$
inherited from the spherical Bessel functions composing the tensor
$\mathbf{F}^{(+)}(\mathbf{r},\omega)$, are observed. In this particular case
all involved wavelengths are very similar, around 780 nm.

###  III.2 Collective dressed states

As more atoms are included, these dressed states begin to mix. The mixing is
induced by dipole-dipole interactions and, as such, depends on the inter-
atomic distance and the underlying symmetries of the system. At large
distances ($\kappa_{mn}r_{1,2}\gg 1$) the interaction consists mainly of
photon exchanges associated with the $1/\kappa_{mn}r_{1,2}$ term that
dominates the behavior of $\Omega^{12}$ in the radiation zone [see Eq. (21)].
At short inter-atomic distances this tensor acquires a much richer structure
that strongly modifies the dressed state energies. The effect is greater as
the distance is reduced until the dipole approximation begins to break down as
$r_{1,2}$ approaches the atomic size [40].

Figure 3 shows the dressed-state energies obtained for two atoms as a function
of their interatomic distance. The results are obtained without spontaneous
decay. There are 16 different collective states accounting for the possible
combinations of individual dressed states. We then label these states
according to the state they reach asymptotically at large distances. These
asymptotic states separate into symmetric and antisymmetric combinations
(drawn as solid and dashed lines) that take the form

$\ket{i,j}_{\pm}=\frac{1}{\sqrt{2}}\bigg{[}\ket{i}_{1}\otimes\ket{j}_{2}\pm\ket{j}_{1}\otimes\ket{i}_{2}\bigg{]},\,$
(31)

with $i,j=\\{a,b,c,d\\}$. Each pair shares the same energy when the atoms are
far apart, which is given by the sum of the individual dressed-state energies
of Eq. (29). As the atoms are brought together, this degeneracy is lifted. The
division into symmetric and antisymmetric states follows from the invariance
of $\mathcal{H}_{sys}$ under atomic exchange and, as such, is mantained at
small distances. This division leads to anticrossings between states of the
same symmetry and crossing among states of opposite symmetries.

Crossings and anticrossings speak of the mixture between dressed states and
can be used to understand the collective response of the atomic pair. The
states $|{a},{a}\rangle_{+}$ and $|{b},{b}\rangle_{+}$, for example, do not
couple to other dressed states through the dipole-dipole interaction. Their
energies remain constant throughout. By contrast, the energies of the
asymptotic states $|{c},{c}\rangle_{+}$ and $|{d},{d}\rangle_{+}$ are modified
due to the presence of the bare state $\ket{1}$ in their wave function. This
bare state mediates dipole transitions between $|0\rangle$ and $|2\rangle$.
Had it been adiabatically eliminated, the displaced lines would have been
lost.

Symmetric and antisymmetric states are also affected by the possible dipole
transitions. For $r_{1,2}>350$ nm ($1/r_{12}<2.85\mu\text{m}^{-1}$), the
states that asymptotically correspond to $\ket{c,i}_{\pm}$ and
$\ket{b,i}_{\pm}$, ($i=\\{a,d\\}$), show essentially the same oscillations on
their energy shifts but with opposite sign. This similarity is attributed to
the structure of $|b\rangle$ and $|c\rangle$, which is almost equal at the
large detuning $\delta_{1}$ considered here [Eqs. (30)]. As the distance
continues to decrease the similarity between the lines is broken. The
existence of additional states at positive frequencies leads to strong
interactions denoted by additional anticrossings. These anticrossings couple
asymptotic states and lead to a relative change between $|b,i\rangle_{\pm}$
and $|c,i\rangle_{\pm}$ lines.

For example, the energy difference between the dressed states
$|c,d\rangle_{+}$ and $|b,b\rangle$ which was about 70 MHz at
$r_{12}\rightarrow\infty$ becomes 2.1 MHz at $r_{12}=120$ nm causing these
levels to interact strongly. Of relevance for the optical properties shown
below (Sec. IV) is the second avoided crossing observed between the asymptotic
$|c,d\rangle_{+}$ and $|c,c\rangle$ around $r_{12}\simeq 200$ nm
($1/r_{12}\simeq 5$ $\mu\text{m}^{-1}$).

The reason why these dipoles transition start to strongly couple the dressed
states at these short distances is that they become of the order of the Rabi
frequencies for distances such that $\kappa_{mn}r_{1,2}>1$ and then dominate.
The coupling depends on the relative wavelength and the population of
different bare states within the asymptotic states. By numerically switching
each of the couplings on and off, we now monitor the changes in each dressed
state to identify the most relevant parameters.

Figure 4 shows the energies resulting from the predominant couplings
$\Omega_{mn}^{1,2}$ for different states as a function of the inter-atomic
distance compared with the intesnsity of the dominant dipole-dipole coupling.
In Figure 4(a), $\Omega_{0,1}^{1,2}$ is kept on and all other coupling are
artificially turned off. We observe the lifting of the degeneracy for the
states labeled $\ket{b,d}_{\pm}$ and $\ket{c,d}_{\pm}$. Similarly, in Figs.
4(b) and 4(c) only $\Omega_{0,3}^{1,2}$ is kept on and the degeneracy is
lifted from the states labeled $\ket{a,c}_{\pm}$, $\ket{a,d}_{\pm}$ and
$\ket{a,b}_{\pm}$. The modified energies follows the same shape as the
dominant $\Omega_{mn}^{1,2}$, it is also affected by the avoided crossings
between states of the same symmetry.

###  III.3 Collective decay rates

Up to now we have focused on the dressed states and neglected the effect of
spontaneous emission, where an excitiation stored within the atoms is set to
escape in the form of free photons. The emission is accounted for by the
superoperator $\mathcal{L}$ of Eq. (25) where the collective decay channels
are obtained by diagonalizing the decay tensor $\gamma^{\alpha\beta}$.

For the two atom case worked above, there are two decay channels for each
$(m,n)$-transition. A symmetric, superradiant channel with a decay rate
$\Gamma_{mn}^{+}=\gamma_{mn}^{1,1}+\gamma_{mn}^{1,2}$; and a subradiant,
antisymmetric channel with decay rate
$\Gamma_{mn}^{-}=\gamma_{mn}^{1,1}-\gamma_{mn}^{1,2}$. These channels are
constructed from the interference of individual decay paths of each atom, as
such, their values tend to the ones of isolated atoms $\gamma_{mn}^{1,1}$ at
large distances. Table 1 shows relevant values for this work.

Decay via these channels is accounted for by the operators

$\Pi_{mn}^{\pm}=\frac{1}{\sqrt{2}}\big{(}\sigma_{mn}^{1}\otimes\mathbbm{1}_{2}\pm\mathbbm{1}_{1}\otimes\sigma_{mn}^{2}\big{)}$
(32)

For the case of two atoms studied here, the symmetric operators drive
transitions between states having the same symmetry under atom-permutation,
while antisymmetric ones connect states with the opposite symmetry only.

$r_{1,2}$ nm | $\Gamma_{0,1}^{+}$ | $\Gamma_{1,2}^{+}$ | $\Gamma_{3,2}^{+}$ | $\Gamma_{0,3}^{+}$ | $\Gamma_{0,1}^{-}$ | $\Gamma_{1,2}^{-}$ | $\Gamma_{3,2}^{-}$ | $\Gamma_{0,3}^{-}$
---|---|---|---|---|---|---|---|---
$60$ nm | $59.7$ | $1.04$ | $4.71$ | $28.1$ | $12.6$ | $0.234$ | $0.998$ | $6.23$
$120$ nm | $58.1$ | $0.991$ | $4.58$ | $26.6$ | $14.2$ | $0.290$ | $1.13$ | $7.68$
$180$ nm | $55.6$ | $0.907$ | $4.37$ | $24.4$ | $16.7$ | $0.376$ | $1.34$ | $9.87$
$240$ nm | $52.4$ | $0.805$ | $4.11$ | $21.8$ | $19.9$ | $0.477$ | $1.60$ | $12.5$
$\infty$ nm | $36.2$ | $0.641$ | $2.85$ | $17.2$ | $36.2$ | $0.641$ | $2.85$ | $17.2$
Table 1: Sample values of the collective decay rates in MHz as a function of
the inter-atomic distance.

##  IV Collective effects imprinted on photon pairs

The collective effects that are present in the dipole shifts and spontaneous
decay are also reflected on the outgoing photons. In this last section we
explore these effects over photon pairs through the second order correlation
function

$\mathcal{G}^{(2)}(x_{1},x_{2})=\expectationvalue{\boldsymbol{E}^{\dagger}(x_{1})\boldsymbol{E}^{\dagger}(x_{2})\boldsymbol{E}(x_{2})\boldsymbol{E}(x_{1})}\,,$
(33)

which is proportional to the probability of detecting a photon at point
$x_{1}$ along with a photon at point $x_{2}$. Here $x=(\boldsymbol{r},t)$ are
spatio-temporal coordinates. The following results remains grounded on
experiments where correlations between photon pairs are controlled externally
by the pump beams. We are particularly interested in the effect of dipole-
dipole interactions over this control.

To explore these effects we consider an array of photon detectors that are
used to characterize the generated light. The detectors are placed at
positions $\mathbf{R}=(R,\theta,\phi)$ in the far field to cover a surface
area ${\mathcal{A}}=R^{2}\Delta\Omega$ of solid angle $\Delta\Omega$ each. The
positive frequency component of the electric field reaching a detectors is

$\boldsymbol{E}^{(+)}(\boldsymbol{R},t_{R})=\mathbf{E}_{\text{f}}^{(+)}(\boldsymbol{R},t_{R})+\mathbf{E}_{\text{s}}^{(+)}(\boldsymbol{R},t_{R})$
(34)

where

$\boldsymbol{E}_{\text{s}}^{(+)}(\boldsymbol{R},t_{R})=\frac{\hbar}{\epsilon_{0}c^{2}}\sum_{\alpha}\sum_{mn}\frac{\boldsymbol{p}_{mn}-\hat{\boldsymbol{R}}_{\alpha}(\hat{\boldsymbol{R}}_{\alpha}\cdot\boldsymbol{p}_{mn})}{R_{\alpha}}\\\
\times\Delta_{mn}^{2}\exp\left[i\kappa_{mn}\left(R_{\alpha}-R\right)\right]\sigma_{mn}^{\alpha}(t)\,,$
(35)

as obtained from Eqs. (9-14) with $t_{R}=t+R/c$ the retarded time and
$\boldsymbol{R}_{\alpha}=\boldsymbol{R}-\boldsymbol{r}_{\alpha}$ the relative
position between each atoms and the detector.

Figure 5: Correlation function $\text{G}^{(2)}_{2,3;3,0}$ as a function of the
two-photon detuning $\delta_{2}$ for increasing laser power $P_{2}$ [panels
(a) and (b)] and interatomic distance [panel (c)]. (a) Two isolated atoms
($\Omega^{1,2}_{mn}=0$) display a Lorentzian response which is power-broadened
and shifted for increasing powers. (b) When interactions are included, a
double-peaked correlation emerges. The peaks are separated for increasing
powers. In both cases $P_{2}=2-10$ mW increasing in steps of 2 mW (lower to
higher curve) and $r_{1,2}=120$ nm. (c) For $P_{2}=4.0$ mW the double peak is
maintained as the inter-atomic distance grows. Here $r_{12}=80-160$ nm
increasing in steps of 20nm (lower to higher curve). In all cases
$\delta_{1}=-70$ MHz.

From equations (33) and (35) the correlation function of the emitted photons
is traced back to particular atomic transitions. We are interested in the
correlations between photons produced by the cascaded decay
$|2\rangle\rightarrow|3\rangle\rightarrow|0\rangle$. For simplicity, we focus
on cases where the photons are detected in the direction of the laser
propagation ($\hat{R}\cdot\vec{p}_{23}=\hat{R}\cdot\vec{p}_{30}=0$) within a
small time-window ($x_{1}=x_{2}=(R,t)$). These events are very likely due to
phase-matching conditions of the system. The probability to detect two-photon
coincidences within a small time window is proportional to

$\mathcal{G}^{(2)}_{2,3;3,0}{\mathcal{A}}^{2}=\mathcal{K}_{c}\Big{(}\expectationvalue{\sigma^{1}_{22}+\sigma^{1}_{22}\sigma^{2}_{33}+\sigma^{1}_{33}\sigma^{2}_{22}+\sigma^{2}_{22}}\\\
+2\expectationvalue{\sigma^{1}_{23}\sigma^{2}_{32}e^{i\kappa_{23}r_{1,2}}+\sigma^{1}_{32}\sigma^{2}_{23}e^{-i\kappa_{23}r_{1,2}}}\\\
+\expectationvalue{\sigma^{1}_{20}\sigma^{2}_{02}e^{i\kappa_{20}r_{1,2}}+\sigma^{1}_{02}\sigma^{2}_{20}e^{-i\kappa_{20}r_{1,2}}}\Big{)}\\\
:=\mathcal{K}_{c}\text{G}^{(2)}_{2,3;3,0}\\\ $ (36)

with $\kappa_{20}=\kappa_{23}+\kappa_{30}$ and

$\mathcal{K}_{c}=\Big{[}\frac{\hbar^{4}\omega^{4}_{23}\omega_{30}^{4}}{\epsilon_{0}^{2}c^{4}}\Big{]}|{\mathbf{p}}_{2,3}|^{2}|{\mathbf{p}}_{3,0}|^{2}\Delta\Omega.$

The contribution $\text{G}^{(2)}_{2,3;3,0}$ on the right-hand side depends on
the expectation values of strings of atomic operators only. We will use Eq.
(LABEL:eq:G2(0)) to describe the correlations between photonic pairs below.
The results are obtained numerically by solving Eq. (25) for atoms originally
prepared in the ground state until a steady-state is reached (at $t\sim 800$
ns). The computational details are described thoroughly in Appendix A.

Figure 5 shows the correlation $\text{G}^{(2)}_{2,3;3,0}$ as a function of the
two-photon detuning $\delta_{2}$ for increasing laser power $P_{2}$. Panels
(a) and (b) show, respectively, the case where the atoms are independent
($\Omega^{1,2}_{mn}=0$ and $\gamma^{1,2}_{mn}=0$) and that where dipole-dipole
interactions are present. For the independent atoms shown in 5(a) the
correlation displays a Lorentzian shape. The Lorentzian is power broadened and
its center is shifted as the intensity of the laser is increased. By
comparison, when all dipole-dipole interactions are present in Fig. 5(b), we
find a significant reduction in the overall value of $\text{G}_{2,3;3,0}$ and
an emerging two-peak structure. The reduction in the overall value of
$\text{G}^{(2)}_{2,3;3,0}$ is attributed to the photon exchange between nearby
atoms that prevents light from reaching at least one detector. The two peaks,
in turn, arise from the collective channels described above. We observe a
shift of the local maxima of each peak from one side of the two photon
resonance to the other as $P_{2}$ increases. In order to parameterize this
response, we artificially turned on and off different dipole interactions and
found that the coupling $\Omega^{1,2}_{0,1}$ dominates this response.

Figure 5(c) shows how the two-peaked structure changes with the interatomic
distance $r_{1,2}$. The numerical results show that the two peak structure
remains up to $\sim 180$ nm. As expected, when the distance $r_{1,2}$
increases way past this point the shape of the Lorentzian shape found for
independent atoms is recovered.

###  IV.1 Correlations from a dressed-states perspective

To further characterize this spectral signature, we consider the dressed-
states developed in Sec. III.2 above. To do so, we expand the atomic operators
of $\text{G}^{(2)}$ within the dressed basis

$\sigma_{nm}^{1}\sigma^{2}_{m^{\prime},n^{\prime}}=\sum_{\nu,\mu}\sum_{s=\pm}\zeta^{nmm^{\prime}n^{\prime}}_{\nu,\mu,s}|\nu,\mu\rangle_{s}\leftindex_{s}{\langle}\nu,\mu|$
(37)

with $\nu,\mu$ labeling the dressed states. Given the invariance of the system
under atomic exchange there is no coupling between symmetric and antisymmetric
states.

For small atomic distances raging between $120$ and $200$ nm [corresponding to
5-8 $\mu$m-1 in Fig. 3], we find that the main contributions to the
correlation function involve transitions between the asymptotic states
$|a\rangle$, $|b\rangle$, and $|c\rangle$

$\displaystyle\sigma_{22}^{1}\sigma_{33}^{2}+\sigma_{33}^{1}\sigma_{22}^{2}$
$\displaystyle\approx$ $\displaystyle\zeta_{ac+}^{2233}|ac\rangle_{+}\langle
ac|+\zeta_{ab-}^{2233}|ab\rangle_{-}\langle ab|\,$ (38)
$\displaystyle\sigma_{23}^{1}\sigma_{32}^{2}+\sigma_{32}^{1}\sigma_{23}^{2}$
$\displaystyle\approx$ $\displaystyle\zeta_{ac+}^{2332}|ac\rangle_{+}\langle
ac|+\zeta_{ab-}^{2332}|ab\rangle_{-}\langle ab|\,.$ (39)

where the first line refers to contributions by atomic populations and the
second to contributions by photon-exchange. The specific contributions to the
other relevant operators
$\sigma_{20}^{1}\sigma_{02}^{2}+\sigma_{02}^{1}\sigma_{20}^{2}$ and
$\sigma_{22}^{1}+\sigma_{22}^{2}$ are more complicated and are plotted in Fig.
6. In all cases, the coefficients $\zeta$ related to the asymptotic state
$|a\rangle$ remain virtually unaffected by the two-photon detuning
$\delta_{2}$. Operators that are independent of $|a\rangle$, however, change
profoundly with the detuning. This change is attributed to the dressed-states
energies plotted in Fig. 3. Transitions within $|b\rangle$ and $|c\rangle$
states are brought close to resonance for the interatomic distances
considering here and become more likely to interact.

Figure 6: Contributions $\zeta$ to the correlation function $\text{G}^{2}$
arising from transitions among different dressed states [see Eqs.
(LABEL:eq:G2(0)) and (37)]. (a) Dominant coefficients $\zeta_{\nu\mu}^{2222}$
in the expansion of the operator $(\sigma_{22}^{1}+\sigma_{22}^{2})$. (b)
Dominant coefficients $\zeta_{\nu\mu}^{2002}$ in the expansion of
$(\sigma_{20}^{1}\sigma_{02}^{2}+\sigma_{02}^{1}\sigma_{20}^{2})$.

From Eqs. (38)-(39) and the data shown in Fig. 6, it can be concluded that the
structure of the correlation function $\text{G}^{(2)}_{2,3;3,0}$ is dictated
by the $|b,c\rangle_{-}$, $|b,c\rangle_{+}$ and $|b,b\rangle_{+}$ states. An
increased contribution of the state $|b,b\rangle_{+}$ for negative detunings
leads to a diminishing value of $\text{G}^{(2)}_{2,3;3,0}$, while an increased
contribution of $|b,c\rangle_{+}$ and $|c,d\rangle_{+}$, near the positive
detuning at $2.1$ MHz, leads to a maximum value of $\text{G}^{(2)}_{2,3;3,0}$
at that detuning. This behavior underlies the competition between excitation
$|2\rangle\rightarrow|1\rangle\rightarrow|0\rangle$ and de-excitation
$|2\rangle\rightarrow|3\rangle\rightarrow|0\rangle$ paths that control the
correlations between outgoing photon pairs.

The question regarding the double peaked structure of
$\text{G}^{(2)}_{2,3;3,0}$ in Fig. 5 remains. The contribution of each dressed
states is dependent on both the dressed-state matrix element and the relevant
transition probabilities caused by decay. A careful analysis reveals that
these conditions are primarily met by two different sets of transitions, each
one associated with a particular peak of $\text{G}^{(2)}_{2,3;3,0}$. The peak
located at $\delta_{2}\sim 0$ corresponds to the symmetric decay channels
$\Pi_{0,1}^{+}$ and $\Pi_{1,2}^{+}$ of Eq. 32 that mediate transitions with
matching symmetry. These superradiant decay channels are the only two channels
that allow zero energy transitions between the same collective dressed states.
By contrast, the second peak, located at around $-3$ MHz is found to be
produced by the antisymmetric decay channels $\Pi_{0,1}^{-}$ and
$\Pi_{1,2}^{-}$. These subradiant channels couple states with opposite
symmetry. In fact, it is found that the most relevant collective dressed-
states for the structure of $\text{G}^{(2)}_{2,3;3,0}$ are $\ket{b,c}_{+}$ and
$\ket{b,c}_{-}$. Figure 5(c) suggests that the transition energies between
these dressed states remains stable over a range of inter-atomic distances
between $80$ nm and $160$ nm. This is congruent with the results shown in
Figure 3 where the transition energies of the relevant states are seen to be
in the range of 5 MHz.

The relevance of the states $\ket{b,c}_{+}$ and $\ket{b,c}_{-}$ is
corroborated by the following argument. The second order correlation function
$\text{G}^{(2)}_{2,3;3,0}$, Eq. (LABEL:eq:G2(0)), can be split into two parts.
The first, which we denote by $\text{G}_{\text{pp}}^{(2)}(0)$, contains all
the photon interchange processes that do not alter the atomic population of a
given atomic level,

$\text{G}_{\text{pp}}^{(2)}(0)=\expectationvalue{\sigma^{1}_{22}+\sigma^{1}_{22}\sigma^{2}_{33}+\sigma^{1}_{33}\sigma^{2}_{22}+\sigma^{2}_{22}}.$
(40)

The rest of the terms in Eq. (LABEL:eq:G2(0)), correspond to terms that in a
density matrix approach would be directly related terms with an inter-atomic
photon exchange

$\displaystyle\text{G}_{\text{pe}}^{(2)}(0)=$ $\displaystyle
2\expectationvalue{\sigma^{1}_{23}\sigma^{2}_{32}e^{i\kappa_{23}r_{1,2}}+\sigma^{1}_{32}\sigma^{2}_{23}e^{-i\kappa_{23}r_{1,2}}}$
$\displaystyle+$
$\displaystyle\expectationvalue{\sigma^{1}_{20}\sigma^{2}_{02}e^{i\kappa_{20}r_{1,2}}+\sigma^{1}_{02}\sigma^{2}_{20}e^{-i\kappa_{20}r_{1,2}}}.$
(41)

The functions $\text{G}_{\text{pp}}^{(2)}(0)$ and
$\text{G}_{\text{pe}}^{(2)}(0)$ are plotted, respectively, in Figs. 7(a) and
7(c) as a function of $\delta_{2}$ and $r_{1,2}$. Curved dashed lines
corresponding to the energy difference between the dressed states
$\ket{b,c}_{+}$ and $\ket{b,c}_{-}$ and a straight dashed white line at
$\delta_{2}=0$ are also plotted as guide. We observe that
$\text{G}_{\text{pp}}^{(2)}(0)>0$ while $\text{G}_{\text{pe}}^{(2)}(0)<0$ and
its absolute value is smaller than $\text{G}_{\text{pp}}^{(2)}(0)$ so that
$\text{G}^{(2)}(0)>0$, as expected. The negativity of
$\text{G}_{\text{pe}}^{(2)}(0)$ induces the nontrivial structure of the
correlation function and is a consequence of photon mediated transitions
$\ket{0}_{\alpha}\leftrightarrow\ket{1}_{\beta}\leftrightarrow\ket{2}_{\gamma}$,
$\alpha,\beta,\gamma=\\{1,2\\}$ that do not include level $|3\rangle$.

To further emphasize the role of the dressed states $\ket{b,c}_{+}$ and
$\ket{b,c}_{-}$ we also plotted the outcome of normalizing
$\text{G}_{\text{pp}}^{(2)}(0)$ and $\text{G}_{\text{pe}}^{(2)}(0)$ to their
maximum absolute value for each inter-atomic distance $r_{12}$. This is shown
in Figs. 7(c) and 7(d). A strong correlation between the local maximum at a
given inter-atomic distance and the line showing the difference in energies of
the two-atom dressed states $\ket{b,c}_{+}$ and $\ket{b,c}_{-}$ appears. The
varying width of each peak also becomes evident. With a larger distinction
between superradiant and subradiant decays appearing at lower distances [see
Table 1].

Figure 7: (a) $\text{G}_{\text{pp}}^{(2)}(0)$ as a function of $\delta_{2}$
and $r_{1,2}$; (b) $|\text{G}_{\text{pp}}^{(2)}(0)|$ normalized to its maximum
absolute value for each inter-atomic distance as a function of $\delta_{2}$
and $r_{1,2}$; (c) $\text{G}_{\text{pe}}^{(2)}(0)$ as a function of
$\delta_{2}$ and $r_{1,2}$; (d) $|\text{G}_{\text{pp}}^{(2)}(0)|$ normalized
to its maximum absolute value for each inter-atomic distance as a function of
$\delta_{2}$ and $r_{1,2}$. The curbed white line correspond to the energy
difference $\mathcal{E}_{b,c_{+}}-\mathcal{E}_{b,c_{-}}$ between the
corresponding dressed states. In this calculation the phases
$\kappa_{23}r_{1,2}$ and $\kappa_{20}r_{1,2}$ lead to a quasi periodic
modulation in the ordinate axis in Fig. (a).

##  V Conclusions

In this work we have studied the cooperative effects of multi-level atoms
coupled by the radiation field. To accomplish this, the known formalism for
the dipole-dipole coupling of two-level atomic systems developed by Lehmberg
was thoroughly extended to many-level systems.

As an illustrative example we considered a pair of four level atoms in the
diamond configuration. The parameters of the calculations were chosen to be
compatible with experimental realizations. Two external lasers were set to
couple only one side of the diamond configuration. As a consequence we were
able to study the coherent phenomena generated by the inter-atomic coupling in
two distinct and coupled branches, one defined by the levels linked by the
external lasers and another by those connected only via spontaneous decays.

We found that two atoms are enough to generate non trivial correlations
between both branches at short inter-atomic distances. These correlations lead
to a structured coincidence function $\text{G}^{(2)}(0)$ as a function of the
two-photon detuning $\delta_{2}$, the inter-atomic distance and the input
power of the external lasers, Fig. 5. In order to fully understand this
structure, we performed a dressed state description of the system as well as
we identified collective decay channels. We have also analyzed the role of
$\text{G}_{\text{pp}}^{(2)}(0)$ and $\text{G}_{\text{pe}}^{(2)}(0)$.

This careful analysis allowed to understand the presence of a minimum of
$\text{G}^{(2)}_{2,3;3,0}$ in terms of several factors. One of which is the
prominence of the two-atom dark state $|b,b\rangle_{+}$ in the terms
$\sigma_{2,0}^{1}\sigma_{0,2}^{2}+\sigma_{0,2}^{1}\sigma_{2,0}^{1}$ and
$\sigma_{2,2}^{1}+\sigma_{2,2}^{2}$ for $\delta_{2}<0$ [see Figs. 6(a) and 6
(b)]. Another factor behind the presence of such a minimum is the behavior of
the eigenenergies of the two atoms dressed states $\mathcal{E}_{b,c_{+}}$ and
$\mathcal{E}_{b,c_{-}}$. Finally we recognize that all of these factors show
that the dynamics of the branch
$|0\rangle_{\alpha}\leftrightarrow|1\rangle_{\beta}\leftrightarrow|2\rangle_{\gamma}$
controls the structure of the branch
$|2\rangle\alpha\leftrightarrow|3\rangle_{\beta}\leftrightarrow|0\rangle_{\gamma}$,
that for isolated atoms would be a spontaneous cascade decays branch.

It is worthwhile to mention that the influence of the atomic state $\ket{1}$
is enhanced by the coherent interaction between the atoms. As a consequence,
this atomic state should not be adiabatically eliminated when studying
cooperative effects, even though the first pump laser is highly detuned from
the $|0\rangle\rightarrow|1\rangle$ transition.

We must mention that in the spectroscopy of atomic gases the presence of
structured $\text{G}_{2}(0)$ coincidence functions has been experimentally
observed [30].These experiments involve general parameters compatible with our
calculations even though the atoms do not have a constant inter-atomic
separation. Nevertheless, our results are a proof of principle that self-
seeding is a cooperative effect that can be induced by a single pair of atoms
at an adequate distance. More detailed calculations involving statistical
distributions of inter-atomic distances are required to test quantitatively
that the self-seeding evolves in an atom cloud to the macroscopic effect
observed in the microscopic QED description. The feasibility of observing a
macroscopic manifestation of this effect is supported by the stability of the
relevant energy differences between the associated dressed states as functions
of the inter-atomic distance for $\kappa_{nm}r_{12}<1$.

The most important qualitative consequence of our calculations is that if the
distance between the atoms can be controlled, so can the creation of inter-
atomic correlations via dipole-dipole interactions with clear consequences for
the nonlinear optical observables.

## ACKNOWLEDGEMENTS

RG-J, DS and RJ gratefully acknowledge support by the PAPIIT-UNAM grants
IA103024, IN106821, IN112624, and IN104523. AK is indebted to IFUNAM for their
hospitality and was financially supported by Departamento de Ciencias Básicas
UAM-A grant number 2232218. PY-T was financially supported by CONAHCYT
Estancias Posdoctorales por México 2022 (3). PB-B gratefully aknowledges
support from the PAPIIT-UNAM grant IG101324.

##  Appendix A Solution of the dynamical equations

In this Section we describe a general approach to solve Eq. (16), written in
short-hand notation as

$\dot{Q}(t)=L\left(Q(t),t\right)\,.$ (42)

Any atomic operator $Q(t)$ can be written as a linear combination of the
elements of the set of matrices
$\mathcal{Q}_{N}=\\{Q_{1},Q_{2},\dots,Q_{N}\\}$ that spans the whole set of
atomic operators. The operator $\dot{Q}(t)$ as well as the operator
$L\left(Q(t),t\right)$ in Eq. (42) can also be written in terms of such a
basis. Since the total number of states for a system of emitters is
$n_{t}=n_{l}^{n_{a}}$, the set $\mathcal{Q}_{N}$ must be formed of
$N=n_{t}^{2}=n_{l}^{2n_{a}}$ matrices.

It can be shown that any such group of matrices forms an internal vector space
and can be arranged so that they are orthogonal under the trace scalar product

$\mathrm{Tr}\left[Q_{i}^{\dagger}Q_{j}\right]=\mathrm{Tr}\left[Q_{i}^{\dagger}Q_{i}\right]\delta_{i,j}.$
(43)

Thus, any atomic operator $A$ can be expanded in terms of the elements of
$\mathcal{Q}_{N}$ as

$A=\sum_{i}c_{i}Q_{i},$ (44)

where $c_{i}=\mathrm{Tr}[Q_{i}^{\dagger}A]/\mathrm{Tr}[Q_{i}^{\dagger}Q_{i}]$,
according to the orthogonality condition (43).

There is a wide variety of sets that can be used as the basis
$\mathcal{Q}_{N}$. To construct a complete set of matrices of dimension
$n_{t}$ we first build a complete set of matrices of dimension $n_{l}$ that
span the space of a single emitter. Altogether this set will have $n_{l}^{2}$
elements. The elements of this set belong to one of three categories of
matrices. The first $n_{l}$ matrices $q_{1}$, $q_{2}$, $\dots$, $q_{n_{l}}$
are diagonal matrices with only one diagonal element equal to $1$ and all
other elements equal to $0$. These are given explicitly by

$\left(q_{k}\right)_{mn}=\delta_{m,k}\delta_{n,k},$ (45)

where $i,m,n=1,2,\dots,n_{l}$. Next is a set of matrices that are
generalizations of the $\sigma_{x}$ Pauli matrix. These are give by

$\displaystyle(q_{k})_{mn}$ $\displaystyle=$
$\displaystyle\delta_{m,j}\delta_{n,l}+\delta_{m,l}\delta_{n,j},$ (46)

where $j=1,2,\dots,n_{l}-1$, $l=m+1,m+2,\dots,n_{l}$ and
$k=n_{l}+1,n_{l}+2,\dots,n_{l}(n_{l}+1)/2$. The final matrices in the third
set are generalizations of the $\sigma_{y}$ Pauli matrix. Similarly to the
previous set, these are give by

$\displaystyle(q_{k})_{mn}$ $\displaystyle=$ $\displaystyle
i\left(\delta_{m,j}\delta_{n,l}-\delta_{m,l}\delta_{n,j}\right),$ (47)

where $j=1,2,\dots,n_{l}-1$, $l=m+1,m+2,\dots,n_{l}$ and
$k=n_{l}(n_{l}+1)/2+1,n_{l}(n_{l}+1)/2+2,\dots,n_{l}^{2}$. We can now build
the complete base as a Kronecker product of the $q_{k}$ matrices

$Q_{k}=q_{k_{n_{a}}}\otimes q_{k_{n_{a}-1}}\otimes\dots\otimes q_{k_{1}}$ (48)

where $k_{1},k_{2},\dots,k_{n_{a}}=1,2,\dots,n_{l}^{2}$ and
$k=\sum_{i=0}^{n_{a}-1}k_{i}n_{l}^{i}$. This base has the advantage of being
composed of Hermitian matrices therefore any Hermitian matrix expanded in this
base will have exclusively real coefficients.

We can write $N$ differential equations, one for each one of the elements of
$\mathcal{Q}_{N}$ as

$\frac{d}{dt}Q_{i}(t)=L\left(Q_{i}(t),t\right).$ (49)

The formal solution to this equation is

$Q_{i}(t)=U(t)Q_{i}U^{-1}(t),$ (50)

where $U(t)$ is, in general, non-unitary due to the decoherence terms in the
generalized Lehmberg equation.

The set of all the $Q_{i}(t)$’s is a complete set of matrices since the
orthogonality relation is preserved under the transformation $U(t)$:

$\mathrm{Tr}\left[Q_{i}^{\dagger}(t)Q_{j}(t)\right]=\mathrm{Tr}\left[U(t)Q_{i}^{\dagger}Q_{j}U^{-1}(t)\right]\\\
=\mathrm{Tr}\left[Q_{i}^{\dagger}(t)Q_{i}(t)\right]\delta_{i,j}.$ (51)

This allows us to expand the right-hand side of (49) as

$\frac{d}{dt}Q_{i}(t)=\sum_{j}\frac{\mathrm{Tr}\left[Q_{j}^{\dagger}(t)L\left(Q_{i}(t),t\right)\right]}{\mathrm{Tr}\left[Q_{j}^{\dagger}(t)Q_{j}(t)\right]}Q_{j}(t)\\\
=\sum_{j}\frac{\mathrm{Tr}\left[Q_{j}^{\dagger}L\left(Q_{i},t\right)\right]}{\mathrm{Tr}\left[Q_{j}^{\dagger}Q_{j}\right]}Q_{j}(t),$
(52)

since due to the properties of the trace,

$\displaystyle\mathrm{Tr}\left[Q_{j}^{\dagger}(t)L\left(Q_{i}(t),t\right)\right]$
$\displaystyle=$
$\displaystyle\mathrm{Tr}\left[Q_{j}^{\dagger}L\left(Q_{i},t\right)\right],$
(53) $\displaystyle\mathrm{Tr}\left[Q_{j}^{\dagger}(t)Q_{i}(t)\right]$
$\displaystyle=$ $\displaystyle\mathrm{Tr}\left[Q_{j}^{\dagger}Q_{i}\right].$
(54)

Now lets assume that initially the ensemble of emitters is in the state

$\left|\Psi_{0}\right\rangle=\left|\psi_{n_{a}}\right\rangle\otimes\left|\psi_{n_{a}-1}\right\rangle\otimes\dots\otimes\left|\psi_{1}\right\rangle,$
(55)

where $\left|\psi_{i}\right\rangle$ is the initial state of the $i-$th
emitter. Taking the expectation value of both sides of (52) and defining
$w_{i}(t)=\left\langle\Psi_{0}\left|Q_{i}(t)\right|\Psi_{0}\right\rangle$ we
obtain

$\frac{d}{dt}w_{i}(t)=\sum_{j}\mathrm{Tr}\left[Q_{j}^{\dagger}L\left(Q_{i},t\right)\right]w_{j}(t).$
(56)

Defining the vector of expectation values
$\boldsymbol{w}(t)=(w_{1}(t),w_{2}(t),\dots,w_{N}(t))$ we note that this
procedure yields quite a simple linear system of differential equations of the
form

$\displaystyle\frac{d\boldsymbol{w}(t)}{dt}=\Lambda(t)\boldsymbol{w}(t),$ (57)

with initial conditions
$w_{i}(0)=\left\langle\Psi_{0}\left|Q_{i}\right|\Psi_{0}\right\rangle$ and
matrix elements given by
$\Lambda_{i,j}(t)=\mathrm{Tr}\left[Q_{j}^{\dagger}L\left(Q_{i},t\right)\right]$.
Furthermore, the matrix elements of $\Lambda$ are all real numbers, because,
as we mentioned above, $\mathcal{Q}_{N}$ is constituted of Hermitian matrices.

Once the $w_{i}(t)$ are known, the expected value of any operator $A$ can be
expressed by using expansion (44) as

$\left\langle\Psi_{0}\left|A(t)\right|\Psi_{0}\right\rangle=\sum_{i}\frac{\mathrm{Tr}\left[Q_{i}^{\dagger}A\right]}{\mathrm{Tr}\left[Q_{i}^{\dagger}Q_{i}\right]}w_{i}(t).$
(58)

In this manner one need only solve Eq. (16) once for a given set of parameters
and can employ that solution to calculate the expectation value of any
arbitrary atomic operator.

## References

  * Finkelstein  _et al._ [2023] R. Finkelstein, S. Bali, O. Firstenberg, and I. Novikova, A practical guide to electromagnetically induced transparency in atomic vapor, New Journal of Physics 25, 035001 (2023).
  * Ardelt  _et al._ [2016] P.-L. Ardelt, M. Koller, T. Simmet, L. Hanschke, A. Bechtold, A. Regler, J. Wierzbowski, H. Riedl, J. J. Finley, and K. Müller, Optical control of nonlinearly dressed states in an individual quantum dot, Phys. Rev. B 93, 165305 (2016).
  * Bounouar  _et al._ [2017] S. Bounouar, M. Strauß, A. Carmele, P. Schnauber, A. Thoma, M. Gschrey, J.-H. Schulze, A. Strittmatter, S. Rodt, A. Knorr, and S. Reitzenstein, Path-controlled time reordering of paired photons in a dressed three-level cascade, Phys. Rev. Lett. 118, 233601 (2017).
  * Alexanian and Bose [2006] M. Alexanian and S. K. Bose, Two-photon resonance fluorescence, Phys. Rev. A 74, 063418 (2006).
  * Yan  _et al._ [2020] Z. Z. Yan, J. W. Park, Y. Ni, H. Loh, S. Will, T. Karman, and M. Zwierlein, Resonant dipolar collisions of ultracold molecules induced by microwave dressing, Phys. Rev. Lett. 125, 063401 (2020).
  * Bettles  _et al._ [2016] R. J. Bettles, S. A. Gardiner, and C. S. Adams, Enhanced optical cross section via collective coupling of atomic dipoles in a 2d array, Phys. Rev. Lett. 116, 103602 (2016).
  * Jenkins and Ruostekoski [2013] S. D. Jenkins and J. Ruostekoski, Metamaterial transparency induced by cooperative electromagnetic interactions, Phys. Rev. Lett. 111, 147401 (2013).
  * Hadad and Steinberg [2010] Y. Hadad and B. Z. Steinberg, Magnetized spiral chains of plasmonic ellipsoids for one-way optical waveguides, Phys. Rev. Lett. 105, 233904 (2010).
  * Hadad and Steinberg [2013] Y. Hadad and B. Z. Steinberg, One way optical waveguides for matched non-reciprocal nanoantennas with dynamic beam scanning functionality, Opt. Express 21 (2013).
  * Gutiérrez-Jáuregui and Asenjo-Garcia [2022] R. Gutiérrez-Jáuregui and A. Asenjo-Garcia, Directional transport along an atomic chain, Phys. Rev. A 105, 043703 (2022).
  * Joshi  _et al._ [2023] C. Joshi, F. Yang, and M. Mirhosseini, Resonance fluorescence of a chiral artificial atom, Phys. Rev. X 13, 021039 (2023).
  * de Léséleuc al. [2019] S. de Léséleuc al., Observation of a symmetry-protected topological phase of interacting bosons with rydberg atoms, Science 365, 775 (2019).
  * al. [2021] G. S. al., Probing topological spin liquids on a programmable quantum simulator, Science 374, 1242 (2021).
  * Kiczynski  _et al._ [2022] M. Kiczynski, S. K. Gorman, H. Geng, _et al._ , Engineering topological states in atom-based semiconductor quantum dots, Nature 606, 694 (2022).
  * Eberly and Rehler [1969] J. Eberly and N. Rehler, Dynamics of superradiant emission, Physics Letters A 29, 142 (1969).
  * Lehmberg [1970a] R. H. Lehmberg, Radiation from an $n$-atom system. i. general formalism, Phys. Rev. A 2, 883 (1970a).
  * Gardiner and Collett [1985] C. W. Gardiner and M. J. Collett, Input and output in damped quantum systems: Quantum stochastic differential equations and the master equation, Phys. Rev. A 31, 3761 (1985).
  * Lehmberg [1970b] R. H. Lehmberg, Radiation from an n-atom system. i. general ormalism, Phys. Rev. A 2, 883 (1970b).
  * Agarwal [1973] G. S. Agarwal, Rotating-wave approximation and spontaneous emission, Phys. Rev. A 7, 1195 (1973).
  * Gross and Haroche [1982] M. Gross and S. Haroche, Superradiance: An essay on the theory of collective spontaneous emission, Physics Reports 93, 301 (1982).
  * Javanainen [1992] J. Javanainen, Effect of state superpositions created by spontaneous emission on laser-driven transitions, Europhys. Lett. 17, 407 (1992).
  * Juzeliūnas and Carmichael [2002] G. Juzeliunas and H. J. Carmichael, Systematic formulation of slow polaritons in atomic gases, Phys. Rev. A 65, 021601(R) (2002).
  * Fleischhauer  _et al._ [2005] M. Fleischhauer, A. Imamoglu, and J. P. Marangos, Electromagnetically induced transparency: Optics in coherent media, Rev.ModPhys.633 77 (2005).
  * Gardiner [1986] C. W. Gardiner, Inhibition of atomic phase decays by squeezed light: A direct effect of squeezing, Phys. Rev. Lett. 56, 1917 (1986).
  * Carmichael  _et al._ [1987] H. J. Carmichael, A. S. Lane, and D. F. Walls, Resonance fluorescence from an atom in a squeezed vacuum, Phys. Rev. Lett. 58, 2539 (1987).
  * Carmichael [1999] H. J. Carmichael,  _Statistical Methods in Quantum Optics 1: Master Equations and Fokker–Planck Equations_ (Springer, Berlin Chapter VII, 1999).
  * Zibrov  _et al._ [2002] A. S. Zibrov, M. D. Lukin, L. Hollberg, and M. O. Scully, Efficient frequency up-conversion in resonant coherent media, Phys. Rev. A 65, 051801(R) (2002).
  * Becerra  _et al._ [2008] F. E. Becerra, R. T. Willis, S. L. Rolston, and L. A. Orozco, Nondegenerate four-wave mixing in rubidium vapor: The diamond configuration, Phys. Rev. A 78, 013834 (2008).
  * Willis  _et al._ [2010] R. T. Willis, F. E. Becerra, L. A. Orozco, and S. L. Rolston, Correlated photon pairs generated from a warm atomic ensemble, Phys. Rev. A 82, 053842 (2010).
  * Cerè  _et al._ [2018] A. Cerè, B. Srivathsan, G. K. Gulati, B. Chng, and C. Kurtsiefer, Characterization of a photon-pair source based on a cold atomic ensemble using a cascade-level scheme, Phys. Rev. A 98, 023835 (2018).
  * Arias-Tellez  _et al._ [2022] N. Arias-Tellez, I. F. Angeles-Aguillón, D. Martínez-Cara, A. Martínez-Vallejo, L. Y. Villegas Aguilar, L. A. Mendoza-López, Y. M. Torres, R. A. Gutiérrez-Arenas, R. Jáuregui, I. Pérez-Castillo, A. Ceré, and D. Sahagún Sánchez, An experimental setup to generate narrowband bi-photons via four-wave mixing in cold atoms, Revista Mexicana de Física 68, 031303 1– (2022).
  * Martínez-Cara  _et al._ [2023] D. Martínez-Cara, I. F. Ángeles-Aguillón, N. Arias-Téllez, R. Gutiérrez-Jáuregui, R. Jáuregui, and D. S. Sánchez, Quantum coherence and collective processes during four-wave mixing in cold atoms, Archive (2023).
  * Safronova  _et al._ [2004] M. S. Safronova, C. J. Williams, and C. W. Clark, Relativistic many-body calculations of electric-dipole matrix elements, lifetimes, and polarizabilities in rubidium, Phys. Rev. A 69, 022509 (2004).
  * Reitz  _et al._ [2022] M. Reitz, C. Sommer, and C. Genes, Cooperative quantum phenomena in light-matter platforms, PRX Quantum 3, 010201 (2022).
  * Brion  _et al._ [2007] E. Brion, L. H. Pedersen, and K. Mølmer, Adiabatic elimination in a lambda system, Journal of Physics A: Mathematical and Theoretical 40, 1033 (2007).
  * Breuer and Petruccione [2007] H.-P. Breuer and F. Petruccione,  _The Theory of Open Quantum Systems_ (Oxford University Press, 2007).
  * Malcuit  _et al._ [1985] M. S. Malcuit, D. J. Gauthier, and R. W. Boyd, Suppression of amplified spontaneous emission by the four-wave mixing process, Phys. Rev. Lett. 55, 1086 (1985).
  * Autler and Townes [1955] S. H. Autler and C. H. Townes, Stark effect in rapidly varying fields, Phys. Rev. 100, 703 (1955).
  * Cohen-Tannoudji [1996] C. N. Cohen-Tannoudji, The autler-townes effect revisited, in  _Amazing Light: A Volume Dedicated To Charles Hard Townes On His 80th Birthday_ (Springer, 1996) pp. 109–123.
  * Anderson and Raithel [2013] S. E. Anderson and G. Raithel, Ionization of rydberg atoms by standing-wave light fields, Nature Communications 4, 2967 (2013).