# Driven-dissipative phase separation in free-space atomic ensembles

D. Goncalves  ICFO-Institut de Ciencies Fotoniques, The Barcelona Institute of
Science and Technology, 08860 Castelldefels (Barcelona), Spain.  L. Bombieri
ICFO-Institut de Ciencies Fotoniques, The Barcelona Institute of Science and
Technology, 08860 Castelldefels (Barcelona), Spain.  Institute for Theoretical
Physics, University of Innsbruck, 6020 Innsbruck, Austria  Institute for
Quantum Optics and Quantum Information of the Austrian Academy of Sciences,
6020 Innsbruck, Austria  G. Ferioli  Universite Paris-Saclay, Institut
d’Optique Graduate School, CNRS, Laboratoire Charles Fabry, 91127, Palaiseau,
France.  S. Pancaldi  Universite Paris-Saclay, Institut d’Optique Graduate
School, CNRS, Laboratoire Charles Fabry, 91127, Palaiseau, France.  I.
Ferrier-Barbut  Universite Paris-Saclay, Institut d’Optique Graduate School,
CNRS, Laboratoire Charles Fabry, 91127, Palaiseau, France.  A. Browaeys
Universite Paris-Saclay, Institut d’Optique Graduate School, CNRS, Laboratoire
Charles Fabry, 91127, Palaiseau, France.  E. Shahmoon  Department of Chemical
& Biological Physics, Weizmann Institute of Science, Rehovot 7610001, Israel.
D. E. Chang  ICFO-Institut de Ciencies Fotoniques, The Barcelona Institute of
Science and Technology, 08860 Castelldefels (Barcelona), Spain.  ICREA-
Institució Catalana de Recerca i Estudis Avançats, 08015 Barcelona, Spain.

###### Abstract

The driven Dicke model, wherein an ensemble of atoms is driven by an external
field and undergoes collective spontaneous emission due to coupling to a leaky
cavity mode, is a paradigmatic example of a system exhibiting a driven-
dissipative phase transition as a function of driving strength. Recently, a
similar phenomenon was experimentally observed, not in a cavity setting, but
rather in a free-space atomic ensemble. The reason why similar behavior should
emerge in free space is not obvious, as the system interacts with a continuum
of optical modes, which encodes light propagation effects. Here, we present
and solve a simple model to explain the behavior of the free-space system,
based on the one-dimensional Maxwell-Bloch equations. On one hand, we show
that a free-space ensemble at a low optical depth can exhibit similar behavior
as the cavity system, as spatial propagation effects are negligible. On the
other hand, in the thermodynamic limit of large atom number, we show that
certain observables such as the transmittance or the atomic excited population
exhibit non-analytic behavior as a function of the driving intensity,
reminiscent of a phase transition. However, a closer analysis reveals that the
atomic properties are highly inhomogeneous in space, and based on this we
argue that the free-space system does not undergo a phase transition but
rather a “phase separation,” roughly speaking, between saturated and
unsaturated regions.

##  I Introduction

The driven-dissipative Dicke model stands as a paradigmatic example of a
driven, open quantum system. It represents a collection of coherently driven,
two-level systems uniformly coupled to a single, lossy electromagnetic field
mode. Despite its apparent simplicity, it approximately describes an atomic
ensemble coupled to a cavity Hepp and Lieb (1973); Agarwal  _et al._ (1977);
Ritsch  _et al._ (2013); Norcia  _et al._ (2016), like in Fig. 1a, and
displays several rich features that have stimulated sustained research. A
notable aspect of this model is the existence of a driven-dissipative phase
transition Carmichael and Walls (1977); Carmichael (1980); Schneider and
Milburn (2002); González-Tudela and Porras (2013); Kirton and Keeling (2017);
Hannukainen and Larson (2018); Kirton  _et al._ (2019); Barberena  _et al._
(2019); Somech and Shahmoon (2022), governed by the parameter
$\beta=2\Omega/N\Gamma_{\rm 1D}$, with $\Omega$ and $\Gamma_{\rm 1D}$ being
the coherent driving strength and atomic emission rate into the cavity,
respectively. The system exhibits a magnetized ($\beta<1$) and non-magnetized
($\beta>1$) phase, with the latter manifesting non-ergodic dynamics Sánchez
Muñoz  _et al._ (2019); Barberena  _et al._ (2019) and so-called time
crystalline behavior Drummond and Carmichael (1978); Iemini  _et al._ (2018);
Zhu  _et al._ (2019). The magnetized phase exhibits spin squeezing Puri and
Lawande (1979); González-Tudela and Porras (2013), which is known to be a
valuable resource in quantum metrology. In realistic cavity setups, atoms can
also decay independently into free space at a rate $\Gamma$, leading to first-
order phase transition with modified properties, like a shifted critical value
of $\beta$ Drummond and Carmichael (1978); Carmichael and Walls (1977);
Carmichael (1980); Walls  _et al._ (1978); Gegg  _et al._ (2018).

Figure 1: (a) Sketch of the driven-Dicke model, where an ensemble of $N$ two-
level atoms couples identically to a single cavity mode $\hat{a}$. The
ensemble is driven by a coherent field with strength $\Omega$ and decays
collectively into the lossy cavity at a rate $\Gamma_{\rm 1D}$. In realistic
setups, atoms can also decay independently outside of the cavity, i.e. into
the $4\pi$ continuum of modes, at a rate $\Gamma$. Here, the system is
permutational invariant and spatially homogeneous. (b) Sketch of the Maxwell-
Bloch model representing an ensemble in free space, where atoms couple to the
driving mode at rate $\Gamma_{\rm 1D}$ and to the continuum of modes at a rate
$\Gamma$. Here, light propagation effects prevent the field $\hat{a}(z)$ and
atomic operators $\hat{\sigma}^{\alpha}(z)$ from being spatially uniform in
general.

Surprisingly, a recent experiment Ferioli  _et al._ (2023) has shown that an
elongated ensemble of atoms in free space, driven along its main axis,
exhibits features that look remarkably similar to the driven-dissipative Dicke
model, even though there are no cavity mirrors. For example, the magnetization
was observed to be suppressed as a function of driving strength, in a way
resembling the phase transition in the Dicke model. The atomic emission also
appeared to be superradiant, scaling quadratically with system size for strong
driving. However, formally, the cavity and free space systems appear quite
different. The Dicke model involves a discrete cavity mode and thus has no
spatial field degrees of freedom, resulting in permutational invariance of
atomic observables and conservation of total angular momentum. Contrarily, a
free space atomic ensemble couples to a continuum of field modes, and indeed
spatiotemporal propagation effects are behind various well-known phenomena in
such systems Phillips  _et al._ (2001); Liu  _et al._ (2001); Fleischhauer
_et al._ (2005); Firstenberg  _et al._ (2013).

Here, we examine the behavior of a driven elongated atomic ensemble in free
space to elucidate the possible existence of a phase transition, and its
connection to that of the iconic “zero-dimensional” driven Dicke model. To
this end, we present the model illustrated in Fig. 1b, which is based on the
well-established 1D Maxwell-Bloch (M-B) equations for light propagation in
dilute atomic clouds McCall and Hahn (1967); Arecchi and Bonifacio (1965);
Haake  _et al._ (1979); Polder  _et al._ (1979); Gross and Haroche (1982);
Castin and Mølmer (1995); Hammerer  _et al._ (2010); Gorshkov  _et al._
(2011). First, taking a discrete version of this theory, we show the atomic
master equation is identical to the cavity problem in terms of driving and
dissipation, but features an extra coherent (dipole-dipole) interaction
encoding light propagation. This term captures the key difference between
cavity and free-space scenarios and suggests their behavior converges when the
atomic system is approximately spatially uniform (e.g., for low optical depths
or strong saturation). Second, solving the M-B model in the thermodynamic
limit within the mean-field approximation, we find a non-analytic behavior in
the light transmission and magnetization, and identify the parameter governing
an apparent “transition”. A closer inspection, however, reveals that the
system exhibits what we term a “phase separation”, where different spatial
regions of the ensemble acquire different properties (saturated vs.
magnetized) due to the spatial degree of freedom. Finally, we compare our M-B
model with the experimental results of Ref. Ferioli  _et al._ (2023),
obtaining a good agreement.

We highlight the complementary work from Ref. Argawal  _et al._ (2024), which
directly models a three-dimensional atomic ensemble, and arrives at similar
conclusions to ours. Together, these works offer significant insight into the
physics of driven-dissipative free-space atomic ensembles and the
manifestation of cavity-like behavior.

##  II Driven-Dissipative Dicke model

To facilitate the comparison with the free space results discussed later, we
briefly recall the driven-dissipative Dicke model of $N$ two-level atoms
interacting identically with a single-mode cavity, and summarize its main
features. Assuming that the cavity mode is lossy, and starting from the full
atom-light master equation, one can trace out the cavity field within a Born-
Markov approximation to obtain a master equation for the atomic density matrix
$\hat{\rho}$ that reads Agarwal  _et al._ (1977); Drummond and Carmichael
(1978); Carmichael and Walls (1977); Carmichael (1980); Schneider and Milburn
(2002)

$\dot{\hat{\rho}}=-\frac{i}{\hbar}\left[\hat{H}_{\rm
drive},\hat{\rho}\right]+\mathcal{L}_{\Downarrow}[\hat{\rho}].$ (1)

The Hamiltonian term

$\hat{H}_{\rm drive}=\frac{\hbar\Omega}{2}(\hat{S}_{+}+\hat{S}_{-}),$ (2)

describes a resonant coherent field driving the ensemble with Rabi amplitude
$\Omega$. Here, $\hat{S}_{\pm}=\sum_{n=1}^{N}\hat{\sigma}_{n}^{\pm}$ are
collective spin raising and lowering operators, which can be written in terms
of the raising/lowering operators of the individual two-level atoms
$\hat{\sigma}_{n}^{\pm}$. The Liouvillian

$\mathcal{L}_{\Downarrow}[\rho]=\frac{\Gamma_{\rm
1D}}{2}\left(2\hat{S}_{-}\hat{\rho}\hat{S}_{+}-\\{\hat{\rho},\hat{S}_{+}\hat{S}_{-}\\}\right)$
(3)

captures the collective atomic emission at a rate $\Gamma_{\rm 1D}$, which
stems from the coupling to the lossy cavity.

The master equation (1) is invariant under the permutation of two atoms and
conserves total angular momentum $\langle\hat{S}^{2}\rangle$. These symmetries
allow $\hat{\rho}$ to be expressed in terms of a limited number of collective
angular momentum states $|j,m\rangle$, which satisfy

$\hat{S}^{2}|j,m\rangle=j(j+1)\ |j,m\rangle,\vspace{-1em}$ (4a)
$\hat{S}_{z}|j,m\rangle=m\ |j,m\rangle,$ (4b)

with $\hat{S}_{\alpha}=\sum_{n=1}^{N}\hat{\sigma}_{n}^{\alpha}/2$,
$j\in[0,N/2]$ and $m\in[-j,j]$. The conservation of angular momentum breaks
the Hilbert space into different $j$ manifolds, which makes dynamics efficient
to solve.

Figure 2: (a) Normalized steady-state magnetization
$\langle\hat{S}_{z}\rangle/N$ as a function of $\beta=2\Omega/N\Gamma_{\rm
1D}$ for $N=30$ atoms. Top: Solid line is a numerical simulation of the Dicke
model [Eq. (1)], while the dashed line is the mean-field prediction of Eq.
(6). Bottom: Solid lines are a simulation of the Cavity model [Eq. (7)] for
values of local dissipation strengths of $\Gamma=[0.5,1,5,10]\Gamma_{\rm 1D}$
(dark to lighter blue), while the dashed line is Eq. (12). (b) Schematic
ordering of the $|j,m\rangle$ states by their quantum numbers $m$ and $j$.
Without local dissipation, $j$ is conserved and any dynamics starting from the
ground state $|g\rangle^{\otimes N}=|j=N/2,m=-N/2\rangle$ is confined to the
purple box. Instead, local dissipation allows for transitions
$\Gamma^{\rightarrow}_{j,m}$ and $\Gamma^{\leftarrow}_{j,m}$ that modify $j$.
(c) Ratio between the angular momentum-changing decay rates
$\Gamma^{\rightarrow}_{j,m}/\Gamma^{\leftarrow}_{j,m}$ for different $j$ and
$m$, which balance at the dashed line given by Eq. (10). The yellow points are
the results of numerical simulations of Eqs. (7) for different driving
strengths, where we convert the expectation values of angular momentum
operators to $j,m$ using the formulas $m\approx\langle\hat{S}_{z}\rangle$ and
$j(j+1)\approx\langle\hat{S}^{2}\rangle=\langle\hat{S}_{x}^{2}+\hat{S}_{y}^{2}+\hat{S}_{z}^{2}\rangle$.
Here, $N=30$ and $\Gamma=[0.1,1,10]\Gamma_{\rm 1D}$ (squares, points and
triangles).

###  II.1 Driven-dissipative phase transition

As discussed in the introduction, the competition between the coherent driving
$\Omega$ and the collective dissipation $\Gamma_{\rm 1D}$ in Eq. (1) yields a
driven-dissipative phase transition Drummond and Carmichael (1978); Carmichael
and Walls (1977); Carmichael (1980); Schneider and Milburn (2002); González-
Tudela and Porras (2013); Kirton and Keeling (2017); Hannukainen and Larson
(2018); Kirton  _et al._ (2019); Barberena  _et al._ (2019); Somech and
Shahmoon (2022). As a function of the ratio $\beta=2\Omega/N\Gamma_{\rm 1D}$,
the steady-state solution of the density matrix features a magnetization that
goes from non-zero ($\beta<1$) to zero ($\beta>1$), a transition that can
already be seen at the mean-field level. In the following, we focus on the
$j=N/2$ manifold, which contains the state $|j=N/2,m=-N/2\rangle$ with all
atoms unexcited (the typical initial state). Keeping only the terms of order
$\mathcal{O}(N^{2})$, the mean-field Heisenberg equations for the collective
spin operators in the steady-state read

$\Omega\langle\hat{S}_{z}\rangle=\Gamma_{\rm
1D}\langle\hat{S}_{z}\rangle\langle\hat{S}_{y}\rangle,\vspace{-1em}$ (5a)
$\Omega\langle\hat{S}_{y}\rangle=\Gamma_{\rm
1D}(\langle\hat{S}^{2}\rangle-\langle\hat{S}_{z}\rangle^{2}),$ (5b)

where $\langle\hat{S}_{x}\rangle=0$. To solve Eqs. (5), one usually exploits
the conservation of angular momentum to substitute
$\langle\hat{S}^{2}\rangle=N(N/2+1)/2\approx N^{2}/4$, yielding

$\langle\hat{S}_{z}\rangle=-\frac{N}{2}\sqrt{1-\beta^{2}},\qquad\beta<1.$ (6)

Figure 2a shows the magnetization of the system $\langle\hat{S}_{z}\rangle$ as
a function of $\beta$. The mean-field prediction in Eq. (6) is consistent with
a full master equation simulation of the Dicke model (1). The non-analytic
behavior of $\langle\hat{S}_{z}\rangle$ at the critical point $\beta_{c}=1$ is
a hallmark of the phase transition. We note this mean-field result is
supported by other methods, such as numerics Hannukainen and Larson (2018);
Barberena  _et al._ (2019); González-Tudela and Porras (2013), exact solution
Drummond and Carmichael (1978); Puri and Lawande (1979), or an analysis based
on approximate eigenstates Somech and Shahmoon (2022). Additional hallmarks of
the phase transition are the steady-state density matrix going from a pure to
completely mixed state González-Tudela and Porras (2013); Puri and Lawande
(1979) and from squeezed to unsqueezed González-Tudela and Porras (2013);
Dalla Torre  _et al._ (2013) at the same critical point.

The unmagnetized phase $(\beta>1)$ also features superradiant atomic emission,
with an emission rate scaling quadratically with atom number. Indeed, the
atomic steady-state for $\beta>1$ is completely mixed within the $j=N/2$
manifold, which implies $\Gamma_{\rm
1D}\langle\hat{S}^{+}\hat{S}^{-}\rangle=\Gamma_{\rm
1D}\text{Tr}(\hat{S}^{+}\hat{S}^{-})/(N+1)\propto N^{2}$. Notably, this
superradiance is preserved at arbitrarily large values of $\beta$ thanks to
angular momentum conservation.

Although not central to the discussion, we point out that even if the steady-
state solution for $\beta>1$ is a mixed density matrix, a large system can
undergo non-ergodic dynamics with persistent oscillations Sánchez Muñoz  _et
al._ (2019), an effect which has been termed as a boundary time crystal Iemini
_et al._ (2018); Drummond and Carmichael (1978); Zhu  _et al._ (2019).

###  II.2 Realistic Cavity model

In realistic cavity QED setups, the cavity is not closed and atoms can also
emit into the continuum of $4\pi$ modes. The typical assumption in quantum
optics consists in treating these other modes as independent emission, which
modifies the original Dicke model by a local dissipation term
$\mathcal{L}_{\downarrow}[\hat{\rho}]$, such that

$\dot{\hat{\rho}}=-\frac{i}{\hbar}\left[\hat{H}_{\rm
drive},\hat{\rho}\right]+\mathcal{L}_{\Downarrow}[\hat{\rho}]+\mathcal{L}_{\downarrow}[\hat{\rho}],\vspace{-1em}$
(7a)
$\mathcal{L}_{\downarrow}[\hat{\rho}]=\frac{\Gamma}{2}\sum_{n=1}^{N}\left(2\hat{\sigma}^{-}_{n}\hat{\rho}\hat{\sigma}^{+}_{n}-\\{\hat{\rho},\hat{\sigma}^{+}_{n}\hat{\sigma}^{-}_{n}\\}\right).$
(7b)

Equations (7) will be referred to as the Cavity model, where atoms decay
collectively into the cavity mode at a rate $\Gamma_{\rm 1D}$, and
independently into free space at a rate $\Gamma$.

The local dissipation term $\mathcal{L}_{\downarrow}[\hat{\rho}]$ breaks the
conservation of angular momentum and allows the system to explore multiple $j$
manifolds during its evolution. Even so, Eqs. (7) remain invariant under the
exchange of two atoms. This permutational symmetry can be exploited to aid
numerical simulations, as it significantly reduces the Hilbert space size,
only requiring $\mathcal{O}(N^{3})$ elements instead of $\mathcal{O}(2^{N})$
Sarkar and Satchell (1987); Shammah  _et al._ (2018); Xu  _et al._ (2013);
Chase and Geremia (2008).

###  II.3 Dynamics in the Cavity model

The non-conservation of angular momentum makes the system depart from the
initial $j=N/2$ manifold. To understand why, it is convenient to rewrite the
cavity master Eq. (7) in terms of the angular momentum states by projecting
the density matrix $\hat{\rho}$ into the $|j,m\rangle$ basis,

$\langle
j^{\prime},m-1|\dot{\hat{\rho}}|j^{\prime},m-1\rangle=\sum_{j,m}\Gamma^{\alpha}_{j,m}\langle
j,m|\hat{\rho}|j,m\rangle+...$ (8)

This reveals the rates $\Gamma_{j,m}^{\alpha}$ at which population is
transferred from the state $|j,m\rangle$ into the states
$|j^{\prime},m-1\rangle$ via spontaneous emission, where
$\alpha=\\{\leftarrow,\rightarrow,\downarrow\\}$ indicates whether the angular
momentum $j^{\prime}$ of the final state has increased or decreased by one
unit or stayed the same, as illustrated in Fig. 2b. In the following, we focus
exclusively on the angular momentum-changing decay rates
$\Gamma^{\leftarrow}_{j,m}$ and $\Gamma^{\rightarrow}_{j,m}$, which change the
value of $j$ by $+1$ and $-1$, respectively, and whose explicit form can be
found in the literature to be Shammah  _et al._ (2018)

$\Gamma^{\leftarrow}_{j,m}=\Gamma\frac{(j-m+1)(j-m+2)(\frac{N}{2}-j)}{2(j+1)(2j+1)},\vspace{-0.5em}$
(9a)
$\Gamma^{\rightarrow}_{j,m}=\Gamma\frac{(j+m-1)(j+m)(j+1+\frac{N}{2})}{2j(2j+1)}.$
(9b)

Figure 2c illustrates how the ratio between these decay rates depends on $j$
and $m$. In the red region,
$\Gamma^{\rightarrow}_{j,m}\gg\Gamma^{\leftarrow}_{j,m}$ and a given state
evolves toward lower $j$ manifolds. The opposite happens in the blue region,
where the angular momentum increases.

Intuitively, in the steady state, one expects the system to lie on the white
curve where these decay rates balance. Let us treat $j$ and $m$ as continuous
variables and consider the limit $N\gg 1$. Imposing
$\Gamma^{\rightarrow}_{j,m}=\Gamma^{\leftarrow}_{j,m}$, one finds a universal
relation between angular momentum and magnetization in the steady state

$j=\frac{N}{2}\sqrt{1-(1+2m/N)^{2}},$ (10)

which interestingly is independent of $\Gamma$ (provided it is not exactly
zero). We supplement this intuitive argument with full master equation
simulations. In particular, in Fig. 2c, we plot with points the calculated
steady-state values of $\langle\hat{S}_{z}\rangle$ and
$\langle\hat{S}^{2}\rangle$ for $N=30$ atoms and various ratios of
$\Gamma/\Gamma_{\rm 1D}$. We find that, after associating the expectation
values of the observables with $m\approx\langle\hat{S}_{z}\rangle$ and
$j(j+1)\approx\langle\hat{S}^{2}\rangle$, the steady-state values indeed lie
on the predicted universal curve of Eq. (10).

As an aside, one can use similar arguments to derive conditions under which
the steady state in the Cavity model becomes mixed. Indeed, since Eq. (7)
never creates coherences between different $j$ manifolds, one expects a pure
state only if the population never leaves the initial $j$ manifold, or
equivalently, if $\Gamma^{\rightarrow}_{j,m}<\Gamma^{\leftarrow}_{j,m}$. Let
us focus on states close to the $j=N/2$ manifold, such that
$j=\frac{N}{2}-\epsilon_{j}$ with $\epsilon_{j}\sim\mathcal{O}(1)$. Using
again Eqs. (9) in the $N\gg 1$ limit, this inequality relation between the
decay rates yields

$n\lesssim\frac{1}{\sqrt{N}},$ (11)

where $n=m/N+0.5$ is the fraction of excited population. Interestingly, the
condition (11) cannot be satisfied in the thermodynamic limit
$(N\rightarrow\infty)$ for a non-zero $n$. This is in stark contrast with the
Dicke model, where the system is always pure in the magnetized phase
$(\beta<1)$.

###  II.4 Phase transition in the Cavity model

A phase transition still persists in the presence of local dissipation, albeit
with modified properties. Let us consider a small amount of local dissipation
$\Gamma<N\Gamma_{\rm 1D}$ in Eq. (7). Since independent emission only
introduces corrections of order $\mathcal{O}(N)$ in the Heisenberg equations,
Eqs. (5) remain valid. However, now $\langle\hat{S}^{2}\rangle\not=N^{2}/4$
because angular momentum is not conserved. Instead, we substitute the value of
angular momentum using Eq. (10), where we again approximate the values of $j$
and $m$ by the expectation values of $\langle\hat{S}^{2}\rangle$ and
$\langle\hat{S}_{z}\rangle$, respectively. This self-consistent mean-field
theory yields three solutions, from which only $\langle\hat{S}^{z}\rangle=0$
and

$\langle\hat{S}_{z}\rangle=-\frac{N}{4}-\frac{N}{4}\sqrt{1-2\beta^{2}},\qquad\beta<\frac{1}{\sqrt{2}},$
(12)

are stable. Equation (12) indicates a first-order phase transition at a
shifted critical point $\beta_{c}^{\prime}=1/\sqrt{2}$ compared to the
original Dicke model ($\beta_{c}=1$). Figure 2a shows how Eq. (12) agrees with
a master equation simulation of the Cavity model (7) for a fixed $N$ and
sufficiently small $\Gamma$. On the other hand, for increasing $\Gamma$ (but
fixed $N$), the transition from magnetized to unmagnetized is displaced to
larger values of $\beta$, far away from the $\beta^{\prime}_{c}$ prediction in
Eq. (12). Again, these findings are supported by other methods such as a
factorization ansatz Walls  _et al._ (1978); Carmichael (1980); Drummond and
Carmichael (1978) and numerics Sarkar and Satchell (1987); Gegg  _et al._
(2018).

##  III Ensemble in free space

We now introduce our minimal model for the experiment reported in Ref. Ferioli
_et al._ (2023), where an elongated ensemble in free space is driven quasi-
uniformly by a focused laser beam (approximately a Gaussian transverse mode).
Our goal is to describe the dynamics of both the atomic operators and the
light propagating in this quasi-one-dimensional mode. To this end, we utilize
the 1D Maxwell-Bloch (M-B) equations, which constitute the predominant theory
to treat these systems McCall and Hahn (1967); Arecchi and Bonifacio (1965);
Haake  _et al._ (1979); Polder  _et al._ (1979); Gross and Haroche (1982);
Castin and Mølmer (1995); Hammerer  _et al._ (2010); Gorshkov  _et al._
(2011). We first introduce these equations in a somewhat different form to
what appears in textbooks – in particular taking the atoms as discrete point-
like scatterers – before showing their equivalence to the standard form of the
M-B equations, where atoms appear as a continuous density. This discrete
formulation has the advantage that it highlights the mathematical difference
with the Cavity model, once the field is integrated out.

###  III.1 Maxwell-Bloch model

We start from the atom-light Hamiltonian of $N$ atoms coupled to a one-
dimensional continuum of propagating photons in free space, described by the
field operator $\hat{a}(z)$. This field represents the right-propagating mode
by which atoms are driven and probed. The system Hamiltonian, including both
atoms and photons, reads Calajó  _et al._ (2019); Shen and Fan (2009, 2005);
Calajó and Chang (2022); Mahmoodian  _et al._ (2020)

$\hat{H}=-i\hbar
c\int\text{d}z\left[\hat{a}^{\dagger}(z)\partial_{z}\hat{a}(z)\right]+\\\
+\hbar\sqrt{\frac{c\Gamma_{\rm
1D}}{2}}\sum_{n=1}^{N}\left[\hat{a}^{\dagger}(z_{n})\hat{\sigma}^{-}_{n}+h.c.\right].$
(13)

The first term describes the energy of the free-propagating fields, where
$\hat{a}^{\dagger}(z)$ and $\hat{a}(z)$ create and annihilate a right-
propagating photon at position $z$. They satisfy standard bosonic commutation
relations $[\hat{a}(z),\hat{a}^{\dagger}(z^{\prime})]=\delta(z-z^{\prime})$.
The second term in Eq. (13) accounts for the interaction between atoms at
discrete positions $z_{n}$ and photons, where we have assumed an identical,
frequency-independent coupling rate $\Gamma_{\rm 1D}$.

Next, we use Eq. (13) to calculate the Heisenberg equations of motion for the
field and atomic operators

$\left[\frac{\partial}{\partial z}+\frac{1}{c}\frac{\partial}{\partial
t}\right]\hat{a}(z)=-i\sqrt{\frac{\Gamma_{\rm
1D}}{2c}}\sum_{n=1}^{N}\hat{\sigma}^{-}_{n}\delta(z-z_{n}),$ (14a)
$\frac{\partial}{\partial t}\hat{\sigma}^{-}_{n}=i\sqrt{\frac{c\Gamma_{\rm
1D}}{2}}\hat{a}(z_{n})\hat{\sigma}_{n}^{z},$ (14b) $\frac{\partial}{\partial
t}\hat{\sigma}^{z}_{n}=i\sqrt{2c\Gamma_{\rm
1D}}\left[\hat{a}^{\dagger}(z_{n})\hat{\sigma}_{n}^{-}-h.c.\right].$ (14c)

The system above constitutes the M-B model, which we later use to investigate
the possibility of a phase transition in free space. However, before doing so,
we show the result of integrating out the field to arrive at an atomic master
equation, as in the Cavity model (7).

Starting from Eq. (14a), we assume the photon propagation occurs in a time-
scale much shorter than the atomic dynamics, as $\Gamma_{\rm 1D}\gg L/c$,
where $L$ is the system length. This allows us to disregard the time
derivative in Eq. (14a), which is equivalent to a Markov approximation Shi
_et al._ (2015); Caneva  _et al._ (2015); Gross and Haroche (1982).
Integrating over space yields the input-output relation

$\hat{a}(z)=a_{\rm in}-i\sqrt{\frac{\Gamma_{\rm
1D}}{2c}}\sum_{n=1}^{N}\hat{\sigma}^{-}_{n}\theta(z-z_{n}),$ (15)

where $\theta(z-z_{n})$ is the Heaviside step function with $\theta(0)=1/2$.
In essence, Eq. (15) shows the field at position $z$ is the coherent sum of
the input light $a_{\rm in}$ (taken to be a $c$-number corresponding to a
coherent state) and the atomic emission. Consistently, only the atoms at the
left of the position $z$ contribute to the total field, as $\hat{a}(z)$ is
right propagating.

The fact that the total field driving the atoms in Eq. (14) can be expressed
purely in terms of the input field and other atoms suggests that an atom-only
description is possible (see, e.g., Refs. Asenjo-Garcia  _et al._ (2017);
Meystre and Sargent (2007) for more details). In addition to the propagating
mode of interest, we assume that the atomic interaction with the remaining
free-space modes can be modeled as independent emission $\Gamma$, as in the
Cavity model of Sec. II.2. The resulting master equation for the atomic
density matrix $\hat{\rho}$ reads Calajó and Chang (2022)

$\dot{\hat{\rho}}=\frac{-i}{\hbar}\left[(\hat{H}_{\rm drive}+\hat{H}_{\rm
dd}),\hat{\rho}\right]+\mathcal{L}_{\Downarrow}^{\prime}[\hat{\rho}]+\mathcal{L}_{\downarrow}[\hat{\rho}],$
(16a) $\hat{H}_{\rm dd}=-\frac{i\hbar\Gamma_{\rm
1D}}{4}\sum_{n>m}\left[\hat{\sigma}^{+}_{n}\hat{\sigma}^{-}_{m}-\hat{\sigma}^{+}_{m}\hat{\sigma}^{-}_{n}\right],$
(16b) $\mathcal{L}_{\Downarrow}^{\prime}[\hat{\rho}]=\frac{\Gamma_{\rm
1D}}{4}\left[2\hat{S}^{-}\hat{\rho}\hat{S}^{+}-\left\\{\hat{\rho},\hat{S}^{+}\hat{S}^{-}\right\\}\right],$
(16c)

where $\hat{H}_{\rm drive}=\hbar\sqrt{c\Gamma_{\rm 1D}/2}\sum_{n}(a_{\rm
in}^{*}(z_{n})\hat{\sigma}^{-}_{n}+h.c.)$ describes the driving by the input
field and $\hat{H}_{\rm dd}$ describes coherent dipole-dipole interactions
mediated by photons in the quasi-1D mode.

The master equation (16) reveals that the only difference between our free
space model and the Cavity model (7) is the coherent dipole-dipole interaction
term $\hat{H}_{\rm dd}$. Note that permutationally invariant dipole-dipole
interactions can also appear in the Cavity model, when the atomic and cavity
resonance frequencies do not match Ritsch  _et al._ (2013). However, the
$\hat{H}_{\rm dd}$ in Eq. (16b) is different in nature, as it depends on the
relative atomic positions and cannot be made to vanish by some choice of
atomic resonance frequency. Intuitively, this term is necessary to capture
pulse propagation within the ensemble, as is allowed in free space, and thus
breaks the permutational symmetry of the Cavity (or Dicke) model.

We emphasize that the M-B model shown in Eqs. (14) is simply a minimal (but
oft-used) model for light propagation in atomic ensembles. As the emission
into the modes not of explicit interest is considered to be independent, it
neglects effects associated with wave interference and multiple scattering of
light in these modes. For example, this model cannot capture the “speckled”
pattern of the scattered far-field intensity when the atomic positions are
disordered Lagendijk and van Tiggelen (1996); Schilder  _et al._ (2017), nor
can it capture Bragg diffraction Birkl  _et al._ (1995) or the perfect
reflection of light from ordered arrays Bettles  _et al._ (2015); Shahmoon
_et al._ (2017); Rui  _et al._ (2020), which has recently been realized to be
a valuable resource for applications Manzoni  _et al._ (2018); Bekenstein  _et
al._ (2020); Moreno-Cardoner  _et al._ (2021). In dense ensembles, the dipole-
dipole interactions that result from atoms experiencing the scattered fields
of other atoms can result in significant shifts in the collective atomic
energies and decay rates Asenjo-Garcia  _et al._ (2017); Glicenstein  _et al._
(2020); Ferioli  _et al._ (2021a); Jennewein  _et al._ (2018, 2016),
fundamental dephasing mechanisms Ji  _et al._ (2023), or limitations to the
ensemble refractive index Andreoli  _et al._ (2021).

While the M-B model certainly misses the physics above, these effects are
generally believed to be negligible in dilute atomic clouds. Our focus remains
to describe the propagation of the input field and its interaction with the
ensemble, which justifies our quasi-one-dimensional model where the remaining
$4\pi$ of modes is treated as independent emission. We also again point to the
complementary work of Ref. Argawal  _et al._ (2024), where the three-
dimensional nature of the problem is treated in more detail, and which arrives
at similar conclusions.

###  III.2 Mean-field steady-state solution

The M-B equations generally have no exact solution beyond the limit of weak
driving (linear response), so we resort to a mean-field treatment to look for
a phase transition in the thermodynamic limit. We include the photonic degrees
of freedom back in the discussion and go back to our Eqs. (14) with the
addition of the independent emission term
$\mathcal{L}_{\downarrow}[\hat{\rho}]$. To further simplify the equations, we
replace the granular ensemble with a macroscopically smooth medium. This is
done by substituting the atomic operators with continuous quantum fields
satisfying
$\sum_{n=1}^{N}\hat{\sigma}^{-}_{n}\delta(z-z_{n})\approx\frac{N}{L}\hat{\sigma}^{-}(z)$.
The resulting equations of motion have the usual M-B form Scully and Zubairy
(1997)

$\partial_{z}\hat{E}(z)=-i\sqrt{\frac{\Gamma_{\rm
1D}}{2}}\frac{N}{L}\hat{\sigma}^{-}(z),$ (17a)
$\dot{\hat{\sigma}}^{-}(z)=-\frac{\Gamma}{2}\hat{\sigma}^{-}(z)+i\sqrt{\frac{\Gamma_{\rm
1D}}{2}}\hat{E}(z)\hat{\sigma}^{z}(z),$ (17b)
$\dot{\hat{\sigma}}^{z}(z)=-\Gamma(\hat{\sigma}^{z}(z)+1)+i\sqrt{2\Gamma_{\rm
1D}}(\hat{E}^{\dagger}(z)\hat{\sigma}^{-}(z)-h.c.),$ (17c)

where the spatiotemporal dependence of the fields is implicit and we have
included the effects of local dissipation $\Gamma$. Here, we have also
normalized the field $\hat{E}=\sqrt{c}\,\hat{a}$ such that
$\langle\hat{E}^{\dagger}\hat{E}\rangle$ represents the number of photons
passing through the plane $z$ per unit of time. In addition, we have neglected
the quantum Langevin noise operators that accompany the dissipation in Eqs.
(17b) and (17c), as they play no role at the mean-field level.

The mean-field solution is obtained by replacing any operator $\hat{O}$ by the
sum of their expectation value and fluctuation $\hat{O}=O+\hat{\delta}_{O}$,
and then neglecting $\hat{\delta}_{O}$ to achieve a solution to lowest order
Gross and Haroche (1982); Kusmierek  _et al._ (2023). Taking
$\Gamma\gg\Gamma_{\rm 1D}$ as in realistic experimental scenarios, the steady-
state solutions of Eqs. (17b) and (17c) read

$\sigma^{-}(z)=-\frac{i\sqrt{2\eta\Gamma}E(z)}{\Gamma+4\eta
E^{2}(z)},\quad\sigma^{z}(z)=\frac{-\Gamma}{\Gamma+4\eta E^{2}(z)},$ (18)

where $\eta=\Gamma_{\rm 1D}/\Gamma$ is the single-atom cooperativity, which
physically represents the ratio of the emission rates into the Gaussian mode
and the remaining free space modes. In terms of physical parameters Tey  _et
al._ (2009); Goncalves  _et al._ (2021),

$\eta=\frac{3}{8\pi}\frac{\lambda^{2}}{A},$ (19)

where $A$ is the mode area of the driving beam and $\lambda$ is the resonant
wavelength of the atomic transition. This guess will be refined later in Sec.
IV to account for the specific could geometry and natural emission mode of the
ensemble, while we keep a more general description here.

Substituting Eq. (18) into Eq. (17a) yields a first-order differential
equation with an analytic solution. Imposing the boundary condition
$E(0)=E_{\rm in}$, the amplitude of the field within a 1D ensemble spanning
$z\in[0,L]$ obeys

$E(z)=E_{\rm in}\ \text{exp}\Bigg{[}\frac{D}{2}\left(s-\frac{z}{L}\right)-\\\
-\frac{1}{2}W\left(\text{log}(D\
s)+D\left(s-\frac{z}{L}\right)\right)\bigg{]}.$ (20)

Here, $W(x)$ is the Wright Omega function Corless and Jeffrey (2002);
Kusmierek  _et al._ (2023), whose specific form will not be relevant. The only
important property is that $W(x)$ exhibits linear behavior for large positive
values of $x\gg 1$ and exponential behavior for negative ones $x\ll-1$. The
solution in Eq. (20) is mathematically equivalent to the problem of light
propagation through a classical saturable absorber Veyron  _et al._ (2022);
Kim and Kwak (2020).

Equation (20) suggests that the mean-field behavior of our system is entirely
determined by two experimentally relevant parameters: the optical depth
$D=2\eta N$ and the saturation parameter $s=2|E_{\rm in}|^{2}/N\Gamma$. The
former characterizes the degree of exponential attenuation of a weak field in
transmission, while the latter is the ratio between the input photon flux and
the photon flux radiated by a completely saturated ensemble of $N$ atoms into
$4\pi$. (Note that a very strongly driven atom becomes completely mixed in
steady-state, and thus $(N/2)\Gamma$ is the maximum rate of photons that could
be scattered by the ensemble.)

For $s\ll 1$, the system is weakly driven and the input field is exponentially
attenuated across the ensemble. One can explicitly see this from the
transmittance in the M-B model, which reads

$T(s,D)=\frac{E^{2}(L)}{E^{2}(0)}.$ (21)

For small $s$, Eqs. (21) and (20) lead to $T\approx e^{-D}$, consistently with
the definition of optical depth. In the opposite regime, where $s\gg 1$, the
system is driven strongly and the atomic response saturates
($\sigma^{-},\sigma^{z}\approx 0$). Physically, the input field greatly
exceeds the atomic emission and the field $E(z)\approx E_{\rm in}$ remains
roughly constant throughout the ensemble, yielding $T\approx 1$.

###  III.3 Existence of phase transition?

Phase transitions are characterized by the non-analytic behavior of
observables (or their derivatives) in the thermodynamic limit upon the change
of a parameter. For example, in the Dicke model, the derivative of $S_{z}$ in
Eq. (6) becomes non-analytic at $\beta=1$. In the M-B model, the thermodynamic
limit corresponds to $D\rightarrow\infty$ while maintaining constant $N/AL$
(fixed density). From Eq. (21), we verify numerically that the transmission
for large optical depth satisfies

$T_{\infty}(s)\approx\begin{cases}0&\text{for $s<1$,}\\\
1-\frac{1}{s}&\text{for $s\geq 1$},\\\ \end{cases}$ (22)

to leading order in $1/s$. We clearly observe two distinct regimes. In the
first regime ($s<1$), the system is magnetized, and the infinite optical depth
leads to zero transmission. In other words, the input photons are scattered
into free space before they reach the end of the ensemble. In the second
regime $(s>1)$, the influx of photons surpasses the maximum atomic emission
rate, resulting in atomic saturation and a non-zero transmission. We compare
the transmission calculated with Eq. (21) and the prediction from Eq. (22) in
Fig. 3a as a function of $s$, where we observe a convergence for increasing
system sizes. Here, we plot the observables in terms of $\sqrt{s}\propto
E_{\rm in}$ to compare later with the Dicke model.

We can also calculate the averaged magnetization
$\bar{s}^{z}=\frac{1}{2L}\int\sigma^{z}(z)\text{d}z$ in the ensemble by
substituting Eq. (20) into Eq. (18). In the thermodynamic limit, one has

$\bar{s}^{z}_{\infty}(s)\approx\begin{cases}\frac{s-1}{2}&\text{for $s<1$,}\\\
0&\text{for $s\geq 1$}.\\\ \end{cases}$ (23)

We again observe two distinct regimes as a function of $s$, where the average
magnetization is non-zero for $s<1$ and zero for $s>1$. This is illustrated in
Fig. 3b, where we notice a striking similarity to the Dicke phase transition
in Fig. 2a.

Figure 3: (a,b) Transmission $T$ and averaged magnetization $\bar{s}^{z}$ in
the M-B model as a function of the square root of the saturation parameter
$\sqrt{s}$ for different optical depths $D=\\{1,5,10,20\\}$ (light to darker
colors). We also plot the respective thermodynamic limits from Eqs. (22) and
(23) in dashed black. Both observables exhibit a sharp transition at $s=1$.
(c) Local magnetization $s^{z}(z)$ as a function of position $z$ and
saturation parameter $s$ for $D=20$. In the magnetized phase ($s<1$), the
system splits into saturated and unsaturated regions at $z\approx z_{c}$.
Here, $\Delta z$ is the length of the boundary between the two regimes. In the
unmagnetized phase ($s>1$), the system is completely saturated.

Upon closer inspection, one notices significant differences with the results
of the Cavity model. In the magnetized phase, the spatial distribution of the
magnetization $s^{z}(z)=\sigma^{z}(z)/2$ is extremely inhomogeneous. Figure 3c
provides a visual representation of a system with large optical depth $D=20$,
showing that for $s<1$ the side of the ensemble where light enters is
saturated and has zero magnetization. However, these atoms scatter all of the
incident light into free space, leaving the atoms on the output side un-
excited.

Indeed, the behavior in Eq. (23) can be easily reproduced by splitting the
ensemble into two regions: fully saturated $(\sigma^{z}(z)=0)$ and not
illuminated $(\sigma^{z}(z)=-1)$. Defining the “critical” position $z_{c}$
between these two regions as the position where the variable in the function
$W(x)$ of Eq. (20) changes sign, i.e. $z_{c}\equiv sL$, the average
magnetization reads

$\bar{s}_{\infty}^{z}=\frac{1}{2L}\int_{0}^{z_{c}}\sigma^{z}(z)\text{d}z+\frac{1}{2L}\int_{z_{c}}^{L}\sigma^{z}(z)\text{d}z=\frac{s-1}{2},$
(24)

which is exactly Eq. (23). This is in stark contrast with the magnetized phase
existing in the Cavity and Dicke models, where each individual atom had the
same magnetization.

We note that although the simple model above predicts well the average
magnetization, the position-dependent magnetization $s^{z}(z)$ is actually a
smooth function, as illustrated in Fig. 3c. The transition between the
magnetized and saturated regions of the system occurs not at a single point
$z_{c}$ but over a length scale $\Delta z\propto L/D$. However, in the
thermodynamic limit, this length scale becomes infinitesimally small compared
to the total length $\Delta z/L\propto 1/D\rightarrow 0$, and leads to the
non-analytic behavior of the average magnetization.

Furthermore, even in the global observables like $\hat{S}_{z}$ and $T$ that
exhibit non-analytic behavior, the “critical value” $s=2E_{\rm
in}^{2}/N\Gamma=1$ depends on field intensity, rather than field amplitude as
in the Cavity model. From these arguments, it thus appears that the physics of
the “transition” is fundamentally different from the ones of the Dicke and
Cavity models. For lack of better terminology, one might say that the free-
space ensemble exhibits a phase separation, with distinct regions of the
system exhibiting different behaviors, rather than a phase transition.

###  III.4 Mapping to Cavity QED

Above, we have argued that the behavior of a free-space ensemble in the
thermodynamic limit (high $D$) is fundamentally different than in the case of
a cavity. This is due to propagation effects in the former system, which
allows for spatial inhomogeneity in local observables. On the other hand, for
low $D$, the atoms cannot have a significant effect on the input field. This
causes the atoms to respond nearly homogeneously, and in this regime, the
system behavior can resemble cavity QED, as we now argue in more detail.

Let us start from the M-B model in Eqs. (17). We first integrate Eq. (17a)
with boundary conditions $E(0)=E_{\rm in}$ and substitute the result in Eq.
(17b). Then, we use the relation $\Omega=\sqrt{2\Gamma_{\rm 1D}}E_{\rm in}$ to
rewrite the input field in the units of the Cavity model, and perform a mean-
field approximation to decorrelate the operators. Finally, we assume that at
low $D$ the ensemble is spatially homogeneous, which allows us to approximate
$\langle\hat{\sigma}^{\alpha}_{i}\rangle\approx\langle\hat{\sigma}^{\alpha}\rangle$.
The resulting mean-field Heisenberg equation for $\hat{S}^{-}$ looks
remarkably similar to the one of the Cavity model for the same assumptions,
where

$\langle\dot{S}_{-}^{\rm Cav}\rangle\approx-\frac{\Gamma
N}{2}\langle\hat{\sigma}^{-}\rangle+\frac{i\Omega
N}{2}\langle\hat{\sigma}^{z}\rangle+\frac{\Gamma_{\rm
1D}N^{2}}{2}\langle\hat{\sigma}^{z}\rangle\langle\hat{\sigma}^{-}\rangle,$
(25a) $\langle\dot{S}_{-}^{\rm MB}\rangle\approx-\frac{\Gamma
N}{2}\langle\hat{\sigma}^{-}\rangle+\frac{i\Omega
N}{2}\langle\hat{\sigma}^{z}\rangle+\frac{\Gamma_{\rm
1D}N^{2}}{4}\langle\hat{\sigma}^{z}\rangle\langle\hat{\sigma}^{-}\rangle.$
(25b)

The only difference is a factor $1/2$, which arises from the fact that within
the M-B model, the radiation of a given atom into the quasi-1D mode can only
affect atoms to the right ($\sim N^{2}/2$ combinations instead of $\sim
N^{2}$), a so-called “chiral” interaction.

To further support this mapping, we use Eqs. (18) and (20) to calculate the
steady-state magnetization $\langle\hat{S}_{z}\rangle$ and total angular
momentum $\langle\hat{S}^{2}\rangle$ in the M-B model for varying optical
depths $D$ and driving strengths. As we illustrate in Fig. 4, the results from
the M-B model obey the same universal relation (10) from the Cavity model at
low $D$. To calculate the expectation values of observables, such as
$\hat{S}^{2}$, we used a product state ansatz, which we discuss in more detail
in the next section.

The fact that an ensemble becomes homogeneous at low $D$ has important
implications for the dipole-dipole interactions. Starting from the effective
Hamiltonian $\hat{H}_{\rm dd}$ in Eq. (16b), one can see that the combination
of mean-field and spatial homogeneity renders the coherent part of the dipole-
dipole interactions negligible. Explicitly,

$\langle\hat{H}_{\rm dd}\rangle\approx\frac{\Gamma_{\rm
1D}N^{2}}{8}(\langle\hat{\sigma}^{+}\rangle\langle\hat{\sigma}^{-}\rangle-\langle\hat{\sigma}^{+}\rangle\langle\hat{\sigma}^{-}\rangle)=0.$
(26)

This again illustrates the role of the dipole-dipole interactions in encoding
field propagation (and thus spatially inhomogeneous response) in the ensemble.

In addition to the low $D$ limit, the Cavity and M-B models converge for
sufficiently large driving strengths. In this limit, the atoms become
saturated, the atomic coherences $\langle\hat{\sigma}^{\pm}_{n}\rangle$ vanish
and the chiral dipole-dipole interactions become negligible
$\langle\hat{H}_{\rm dd}\rangle\approx 0$. This convergence can also be
appreciated in Fig. 4, where the Cavity and M-B curves coincide for large
enough driving (small enough $m$) regardless of the optical depth. In Appendix
B, we provide a complementary analysis based on how $\hat{H}_{\rm dd}$ acts on
the angular momentum states in the strong driving limit.

Figure 4: Steady-state relationship between magnetization and total angular
momentum in the M-B model. The points correspond to the estimated quantum
numbers $m=\langle\hat{S}_{z}\rangle$ and $j(j+1)=\langle\hat{S}^{2}\rangle$
from a numerical simulation of the M-B model for different optical depths $D$
(light to dark blue). For a fixed optical depth, the different points are
obtained by varying the driving strength. We also plot the steady-state
universal relation (10) from the Cavity model with black dashed lines. The
models agree when the optical depth is low, and/or the driving strength is
sufficiently large.

##  IV Modeling of experimental data

We now reproduce some key experimental results from Ref. Ferioli  _et al._
(2023) using the M-B model. Instead of the natural parameters $D$ and $s$ of
the M-B model, the experiment had access to atom numbers $N$ and Rabi
frequencies $\Omega$ Ferioli  _et al._ (2023). Thus, to relate these
quantities, we need to establish first a value for the cooperativity $\eta$ of
the experiment.

To do this, we need to point out a subtlety that arises when modeling a 3D
ensemble by the 1D M-B model. Within the 1D M-B model, the cooperativity
$\eta$ actually has two physical meanings. First, from Eq. (19), we see that
$\eta$ defines the emission rate of a single atom into the 1D mode of the M-B
model (the Gaussian mode of the driving beam). However, from Eq. (16), $\eta$
also characterizes the ratio of independent to collective dissipation
strengths of a collective atomic spin wave excitation. Within the 1D M-B
model, where there is no explicit concept of transverse mode shapes, these two
meanings naturally coincide.

However, in a 3D geometry such as an elongated ensemble, a collective spin
wave excitation will have an emission pattern or transverse mode structure
that is not necessarily the same as the input mode Gross and Haroche (1982).
If the collective atomic excitation is constructively emitting into a
significantly larger solid angle (as is the case of the experiment), then the
$\eta$ as calculated from Eq. (19), where only the contribution into the solid
angle of the driving mode is considered, will significantly underestimate the
collective dissipation rate.

In order to arrive at a more realistic model for $\eta_{\rm eff}$ that
accounts for this effect, we consider the total power $P$ emitted by a cloud
of $N$ atoms into a solid angle $\Delta\Theta$ in free space, Allen and Eberly
(1987)

$P=\int_{\Delta\Theta}\text{d}\Omega_{\textbf{k}}\ I_{N}(\textbf{k}),$ (27)

where $I_{N}(\textbf{k})$ is the intensity emitted in the direction of the
wavevector k. For a single atom, and a full solid angle of
$\Delta\Theta=4\pi$, Eq. (27) can be written in terms of the free-space
emission rate $\Gamma$ as

$P=\int_{4\pi}\text{d}\Omega_{\textbf{k}}I_{1}(\textbf{k})\equiv\Gamma\langle\hat{\sigma}^{+}\hat{\sigma}^{-}\rangle,$
(28)

where $I_{1}(\textbf{k})$ is the far-field radiation pattern of a single atom
in free space. This normalization is equivalent to that of the M-B model,
where $P$ represents the number of photons emitted per unit of time.

We take a simple model for the total power, where we assume the ensemble has a
spatially homogeneous response
$\langle\hat{\sigma}_{n}^{-}\rangle\approx\langle\hat{\sigma}^{-}\rangle
e^{ik_{\rm in}z_{n}}$, where $k_{\rm in}$ is the input field wavevector. For
$N\gg 1$, the intensity $I_{N}(\textbf{k})$ emitted in the direction k can be
separated into coherent and incoherent components as Ferioli  _et al._ (2023);
Tanji-Suzuki  _et al._ (2011); Sutherland and Robicheaux (2017)

$I_{\rm
coh}(\textbf{k})=I_{1}(\textbf{k})N^{2}|\langle\hat{\sigma}^{-}\rangle|^{2}F(\textbf{k},\textbf{k}_{\rm
in}),$ (29a) $I_{\rm
incoh}(\textbf{k})=I_{1}(\textbf{k})N\langle\hat{\sigma}^{+}\hat{\sigma}^{-}\rangle,$
(29b)

with only the former having a well-defined phase relation with the input
field. Here, $F(\textbf{k},\textbf{k}_{\rm in})$ is a form factor that depends
on the angle $\theta$ between k and $\textbf{k}_{\rm in}$ and the distribution
of atomic positions. Assuming the atomic positions $z_{n}$ lie on the optical
axis, it reads

$F(\textbf{k},\textbf{k}_{\rm in})=\frac{1}{N^{2}}\sum_{n\not=m}e^{ik_{\rm
in}(1-\cos\theta)(z_{n}-z_{m})}.$ (30)

Physically, $F(\textbf{k},\textbf{k}_{\rm in})$ captures the angle-dependent
interference between atoms in the coherent atomic emission. The key
observation is that for random $z_{n}$, Eq. (30) predicts constructive
interference only within a small solid angle centered around $\textbf{k}_{\rm
in}$. Indeed, assuming $z_{n}$ follow a Gaussian distribution
$\rho(z;\sigma_{L})$ with variance $\sigma_{L}^{2}$ and zero mean, Eq. (30) in
the continuous limit becomes Ferioli  _et al._ (2023)

$F(\theta)\approx\left|\int_{-\infty}^{\infty}\text{d}z\ \rho(z;\sigma_{L})\
e^{ik_{\rm in}(1-\cos\theta)z}\right|^{2}=\\\
=\text{exp}[-\sigma_{L}^{2}k_{\rm in}^{2}(1-\cos\theta)^{2}].$ (31)

Figure (5) shows Eq. (31) as a function of the angle $\theta$, where we see
that $F(\theta)$ vanishes for angles larger than $\theta\lesssim 0.2$ rad.
Here, we have considered a standard deviation of $\sigma_{L}\approx
12.5\lambda$. We note that estimating $\sigma_{L}$ experimentally is quite
challenging, and the value chosen here is roughly half the measured axial size
$l_{\rm ax}\sim 20-25\lambda$ in Ref. Ferioli  _et al._ (2023), consistently
with the discussion surrounding the Appendix C of Ref. Ferioli  _et al._
(2023).

Figure 5: Form factor $F(\textbf{k},\textbf{k}_{\rm in})$ as a function of the
angle $\theta$ between k and the input field $\textbf{k}_{\rm in}$, and
calculated for typical ensemble sizes of the experiment of Ref. Ferioli  _et
al._ (2023). The ensemble exhibits coherent emission for small angles
$\theta\lesssim 0.2$ rad, while emission is mainly incoherent for larger
$\theta$. Inset: Sketch of the incoherent (blue) and coherent (red) components
of the light emitted by an elongated ensemble.

We then estimate the effective single-atom cooperativity $\eta_{\rm eff}$ as
the fraction of coherently-scattered light in the forward direction
$(k_{z}>0)$, in the weak driving limit
$\langle\hat{\sigma}^{+}\hat{\sigma}^{-}\rangle\approx|\langle\hat{\sigma}^{-}\rangle|^{2}$,

$\eta_{\rm
eff}=\frac{\int_{k_{z}>0}\text{d}\Omega_{\textbf{k}}F(\textbf{k},\textbf{k}_{\rm
in})I_{1}(\textbf{k})}{\int_{k_{z}>0}\text{d}\Omega_{\textbf{k}}I_{1}(\textbf{k})}.$
(32)

Intuitively, $\eta_{\rm eff}=1$ when atoms constructively interfere everywhere
in the forward direction, like in the traditional Dicke scenario where all
atoms are localized at the same point ($\sigma_{L}=0$) Dicke (1954). On the
other hand, $\eta_{\rm eff}=0$ when atoms can have arbitrary relative phases
and do not constructively interfere. In the case of the experiment Ferioli
_et al._ (2023), where atoms were circularly polarized
$I_{1}(\textbf{k})=3\Gamma(1+\sin^{2}\theta\cos^{2}\phi)/16\pi$ and
$\sigma_{L}\approx 12.5\lambda$, Eq. (32) yields an effective cooperativity of
$\eta_{\rm eff}\approx 0.01$. We note this is significantly larger than the
value of $\eta\sim 0.002$ if one uses Eq. (19) for the driving field in Ref.
Ferioli  _et al._ (2023), describing the emission efficiency of a single atom
into the input mode.

Finally, let us detail how we calculate the parameters of the M-B model.
First, the optical depth $D$ is obtained from the atom numbers $N$ measured in
Ref. Ferioli  _et al._ (2023), using the relation $D=2\eta_{\rm eff}N$. For
the saturation parameter $s$, we use the relation

$\beta=\frac{2\Omega}{N\eta_{\rm eff}\Gamma}\equiv\sqrt{\frac{8s}{D}},$ (33)

where we substitute the calculated value of $D$ along with the measured values
of $\Omega$ and $N$. This expression also relates the parameter $\beta$ that
dictates the Dicke phase transition with the M-B model.

###  IV.1 Transmission and magnetization

We start by calculating the average magnetization
$\bar{s}^{z}=\langle\hat{S}_{z}\rangle/N$ for various optical depths $D$. In
Fig. 6a, we plot the prediction from Eqs. (18) and (20) as a function of
$\sqrt{s}$, in solid lines. The results agree with both the thermodynamic
limit from Eq. (23), in dashed-black, and the measurements (points in the
figure) in Ref. Ferioli  _et al._ (2023). Again, we plot the average
magnetization in terms of $\sqrt{s}$ because it is proportional to the field
amplitude, like $\beta$ in the Cavity model.

In the inset of Fig. 6a, we compare the average steady-state magnetization in
both the M-B and Dicke models (solid and dashed curves, respectively) as a
function of $\beta$ for the same number of atoms. We see a remarkable
similarity between the predictions of both models in the parameter regime
studied by Ref. Ferioli  _et al._ (2023). We also note that, for larger $D$,
the transition point in the M-B model approaches $\beta_{c}\sim 0$, instead of
converging to $\beta_{c}=1$ as in the Dicke model.

We now calculate the transmittance $T$ using Eq. (21). In Fig. 6b, we plot $T$
as a function of driving amplitude $\Omega$ and for an optical depth of
$D=36$. Importantly, $T$ was not measured in Ref. Ferioli  _et al._ (2023),
and thus we cannot directly compare theory and experiment.

Nevertheless, one can infer an approximated value for $T$ from the
measurements of $\Omega_{\rm eff}$ in Ref. Ferioli  _et al._ (2023). There,
$\Omega_{\rm eff}$ is the oscillation frequency of the averaged magnetization
$\bar{s}^{z}(t)$ at early times ($t\lesssim 5\Gamma$). Assuming the ensemble
is homogeneous enough, it is given by the total field inside the ensemble
$\Omega_{\rm eff}\approx\Omega-i\Gamma_{\rm 1D}\langle\hat{S}^{z}\rangle$,
which is our Eq. (15). In Fig. 6b, we plot this frequency normalized by the
input $\Omega$, yielding an approximate transmission $T\approx(\Omega_{\rm
eff}/\Omega)^{2}$.

Comparing the transmittance from the M-B model and the one inferred from the
measurements of $\Omega_{\rm eff}$, we observe an overall good agreement,
including the sharp transition around $\Omega\approx 5\Gamma$ (see Eq. (22)).
We also observe a good agreement between the M-B model prediction and the
experiment when studying the transmittance as a function of $D\propto N$ and
for a fixed driving strength $\Omega$, as shown in Fig. 6c.

Figure 6: (a) Averaged steady-state magnetization $\bar{s}^{z}$ as a function
of $\sqrt{s}\propto\Omega$ for $D=[4,12,32]$ (light to dark purple). The
points are experimental data from Ref. Ferioli  _et al._ (2023), with
$D=[4,12,32]$ corresponding to the squares, circles, and diamonds,
respectively. The black-dashed curve is the large $D$ limit from Eq. (23).
Inset: Same averaged magnetization, but now as a function of the $\beta$ as
defined in Ref. Ferioli  _et al._ (2023). Solid lines are obtained with the
M-B model, while we used dashed lines for the Dicke model. (b) Steady-state
transmittance $T$ as a function of driving strength $\Omega$ and fixed optical
depth $D$. The red squares are experimental values of $T\approx(\Omega_{\rm
eff}/\Omega)^{2}$ from Ref. Ferioli  _et al._ (2023), and the solid line is
the theoretical prediction of Eq. (21), with the optical depth of $D=36$
inferred from experimental parameters. Phase separation appears at
$\Omega\approx 5\Gamma$. (c) Same as (b) but for fixed $\Omega=4.5\Gamma$ and
varying $D$.

###  IV.2 Atomic emission

We now discuss the intensity of the light emitted by the atomic cloud in the
forward (transmitted) direction, $\gamma_{\rm at}$. In the M-B model, this is
given by the atomic component of the input-output Eq. (15) at the end of the
ensemble ($z>L$). After normalizing $\hat{E}=\sqrt{c}\hat{a}$ so that the
emitted light has units of photon number per unit time, one has

$\gamma_{\rm at}=\frac{\Gamma_{\rm
1D}}{2}\langle\hat{S}^{+}\hat{S}^{-}\rangle.$ (34)

Our theoretical calculation cannot be directly compared with the experimental
measurement. This is due to a technicality of the setup used to measure this
quantity, which is illustrated in Figure 7a Ferioli  _et al._ (2023). Since
the total forward field consists of both input light and atomic emission, a
spatial filter (SF) is inserted to block the input Gaussian component from
detection. Inevitably, this filter also blocks the atomic emission within the
same solid angle covered by the Gaussian mode. Thus, the measured
$\gamma_{at}$ corresponds to the atomic radiation emitted slightly out of the
optical axis, yet still falling within the numerical aperture of the lens. We
will later discuss how to approximately account for this filter in our M-B
theory. For completeness, though, we first present the predictions for
$\gamma_{\rm at}$ within the ideal M-B model.

####  IV.2.1 Product state ansatz

To calculate quantities like atomic emission $\gamma_{\rm at}$ or second-order
photon correlations $g^{(2)}(0)$ within the M-B model, it is convenient to
undo the continuous limit transformation
$\hat{\sigma}_{n}\leftrightarrow\frac{N}{L}\hat{\sigma}(z_{n})$ on the
operators in Eq. (17) and return to a discrete atom picture. The mean-field
approximation used in the previous sections is equivalent to assuming that the
atomic density matrix $\hat{\rho}$ is in the product state Kusmierek  _et al._
(2023)

$\hat{\rho}=\bigotimes_{n=1}^{N}\hat{\rho}_{n},\qquad\hat{\rho}_{n}=\frac{1}{2}\left(\begin{array}[]{cc}1-\langle\hat{\sigma}^{z}_{n}\rangle&2\langle\hat{\sigma}_{n}^{-}\rangle\\\
2\langle\hat{\sigma}_{n}^{-}\rangle^{*}&1+\langle\hat{\sigma}^{z}_{n}\rangle\end{array}\right).$
(35)

Here, $\hat{\rho}_{n}$ is a local density matrix constructed from the mean-
field result of the M-B model at position $z_{n}=nL/N$ with integer
$n\in[1,N]$. As an aside, applying the mean-field approximation
$\hat{E}\rightarrow\langle\hat{E}\rangle$ in Eq. (17) also assumes that the
atoms only interact with the coherent part of the total field.

Substituting the product state ansatz from Eq. (35), the atomic emission in
Eq. (34) becomes

$\gamma_{\rm at}=\frac{\Gamma_{\rm
1D}}{2}\sum_{n=1}^{N}\frac{\langle\hat{\sigma}^{z}_{n}\rangle+1}{2}+\frac{\Gamma_{\rm
1D}}{2}\sum_{n\not=m}\langle\hat{\sigma}^{+}_{n}\rangle\langle\hat{\sigma}^{-}_{m}\rangle.$
(36)

Equation (36) represents the total atomic emission into the 1D mode. The first
and second terms are also referred to as incoherent and coherent components of
atomic emission, as characterized by the absence (presence) of a well-defined
phase relationship with the input field, analogous to our definitions in Eqs.
(29).

Figure 7: Sketch of the experimental setup in Ref. Ferioli  _et al._ (2023)
for the detection of the atomic field intensity $\gamma_{\rm at}$. A spatial
filter (SF) blocks the input field component (dark red) by eliminating modes
with $\textbf{k}\approx\textbf{k}_{\rm in}$. The measured atomic emission
comprises modes slightly out of the optical axis
$\textbf{k}\not\approx\textbf{k}_{\rm in}$ that still fall within the
numerical aperture of the lens, defining the solid angle $\Delta\Theta$. The
measured $\gamma_{\rm at}$ thus contains both coherent (light red) and
incoherent (lighter red) atomic emission.

####  IV.2.2 Steady-state atomic emission

In the following, we describe the properties of the steady-state atomic
emission $\gamma_{\rm at}^{\rm st}$ in Eq. (36), while we discuss its early-
time dynamics in Appendix A. Figure 8a shows $\gamma_{\rm at}^{\rm st}$ as a
function of the saturation parameter $s$, calculated with Eqs. (18), (20) and
(36) for different optical depths. In the magnetized regime ($s\ll 1$), the
phase separation takes place, and only the region with illuminated atoms can
radiate, such that $\gamma_{\rm at}^{\rm st}\propto\Gamma
Nz_{c}/L\propto\Gamma Ds\propto|E_{\rm in}|^{2}$ independently of system size.
Conversely, deep in the unmagnetized regime ($s\gg 1$), atomic coherence
vanishes and $\gamma_{\rm at}^{\rm st}=\Gamma_{\rm 1D}N/4$ like saturated
independent atoms, scaling linearly with the number of atoms.

Between these two extremes, the steady-state atomic emission features a large
maximum around $s\sim 1$. As shown in the inset of Fig. 8a, these maxima seem
to scale quadratically with system size for small optical depth, a feature
reminiscent of the $N^{2}$ superradiant scaling of the Dicke model. However,
two key differences arise. First, it only appears for $s\sim 1$. Unlike the
Dicke model, where the conservation of angular momentum preserves
superradiance at arbitrarily strong drivings, the free-space ensemble lacks
any symmetry and becomes completely saturated for sufficiently large $s\gg 1$,
which in this limit ultimately yields $\gamma_{\rm at}^{\rm st}\propto N$.
Second, while superradiance in the Dicke model prevails for arbitrarily large
systems, the phase separation for large $D$ splits the free-space ensemble
into saturated and unsaturated regions (neither of which being superradiant),
and destroys the $\sim N^{2}$ scaling.

Let us use the Cavity model to elucidate the origin of the $N^{2}$ scaling at
small optical depths, given independent emission. Similarly to Eq. (9), we
calculate the total decay rate $\Gamma^{\rm out}_{m,j}$ from $|j,m\rangle$ to
any other state $|j^{\prime},m-1\rangle$ Shammah  _et al._ (2018), and rewrite
the result in terms of the fraction of excited population $n=m/N+0.5$, using
Eq. (10) to relate $j$ to $n$. We obtain

$\Gamma^{\rm out}_{j,m}(n)\approx\Gamma_{\rm
1D}N^{2}\left(n-2n^{2}\right)+\Gamma Nn,\quad n\leq\frac{1}{2},$ (37)

where we assumed $N\gg 1$ and a non-zero $n$. The first term represents the
collective emission into the cavity mode, while the second captures the atomic
emission into free space. Equation (37) exhibits superradiant scaling $\sim
N^{2}$ at its maximum in $n=1/4$, while it becomes linear with $N$ at
saturation ($n=1/2$).

From our discussion in Sec. III.4, we map Eq. (37) to the M-B model by
replacing $\Gamma$ by $\Gamma_{\rm 1D}$ in the second term to project the
free-space emission into the 1D mode, and by adding the factor $1/2$ to
account for chiral dipole-dipole interactions. This yields the guess

$\gamma_{\rm at}^{\rm st}(\bar{n})\approx\frac{\Gamma_{\rm
1D}}{2}\left(\frac{N^{2}}{2}\left(\bar{n}-2\bar{n}^{2}\right)+N\bar{n}\right),$
(38)

where we recognize the coherent and incoherent components by their scaling
with $N$. Here, $\bar{n}=\sum_{i=1}^{N}n_{i}/N$ is the average excited
population in the ensemble. Despite the simplicity of this calculation, Eq.
(38) is consistent with the M-B results shown in Fig. 8a. Of course, the
prediction from Eq. (38) breaks down at larger optical depths, where the
mapping between the Cavity and M-B models is no longer valid.

Figure 8: (a) Steady-state atomic emission $\gamma_{\rm at}^{\rm st}$ as a
function of the saturation parameter $s$. Solid lines represent the M-B model
[Eq. (36)] with $D=\\{4,8,12,24\\}$ (light to dark red). Dashed lines are
obtained form the mapping to the Cavity model [Eq. (38)] for the same optical
depths. Inset: Scaling of the steady-state maxima $\max[\gamma^{\rm st}_{\rm
at}]$ with $D$. The purple curve is a guide to the eye indicating quadratic
scaling $\propto D^{2}$. (b) Second-order correlations $g^{(2)}(0)$ as a
function of $s$ for the optical depths $D=\\{4,12,24,40\\}$ (light to dark
blue). Atomic emission is bunched for $s>1$ and follows a universal curve
towards fully chaotic light at very large $s$. The dashed black line is the
prediction made from the mapping to the Cavity model [Eq. (41)]. The inset
shows the details in the correlations around $s\sim 1$.

####  IV.2.3 Second-order photon correlations

We now discuss the second-order photon correlations in the forward-emitted
atomic radiation, which are characterized by

$g^{(2)}(0)=\frac{\langle\hat{S}^{+}\hat{S}^{+}\hat{S}^{-}\hat{S}^{-}\rangle}{|\langle\hat{S}^{+}\hat{S}^{-}\rangle|^{2}}.$
(39)

In the Dicke model, and for maximal angular momentum, $g^{(2)}(0)=1$ in the
magnetized phase because the steady state is a coherently radiating spin state
Somech and Shahmoon (2022). Meanwhile, $g^{(2)}(0)=6/5$ in the unmagnetized
phase, which can be calculated by assuming the system is completely mixed
within the $j=N/2$ manifold. Similarly, in the M-B model, we expect
$g^{(2)}(0)=1$ in the magnetized regime ($s<1$), where the coherent component
dominates the atomic emission. However, we instead expect purely chaotic light
$g^{(2)}(0)=2$ deep in the saturated regime ($s\gg 1$), since the angular
momentum is not conserved and the ensemble becomes a set of independent
emitters decaying at random times.

Figure 8b shows the $g^{(2)}(0)$ from Eq. (39) calculated using the product
state ansatz from Eq. (35) on Eqs. (18) and (20). We observe the predicted
behavior above, where correlations go from coherent to chaotic as a function
of $s$. Interestingly, the transition between these two regimes sharpens
around $s=1$, similar to the magnetization and the transmission. This is
consistent with the maxima in Fig. 8a, which marks the point where the
fraction of incoherent light starts to become relevant.

We note that the calculated correlations are not exactly $g^{(2)}(0)=1$ in the
magnetized regime ($s<1$). This subtlety is due to the incoherent component of
the emission. In our mean-field approach, the amount of transmitted incoherent
light is over-estimated, because it is assumed that atoms are only driven by
the coherent component of the field. Thus, incoherently scattered photons are
“free” to transmit through the system without any additional scattering. We
note that in the limit of very weak driving (at most 2 photons sent through
the system), $g^{(2)}(0)$ can be calculated by other means, including the
effects of re-scattering Mahmoodian  _et al._ (2018).

Figure 9: (a) Coherent $I_{\rm coh}(\theta)$ and incoherent $I_{\rm
incoh}(\theta)$ intensity components of the atomic emission as a function of
the angle $\theta$ between the detected mode k and the optical axis
$\textbf{k}_{\rm in}$ defined by the input field. The grey region
$\theta\leq\theta_{\rm min}$ represents the range of angles blocked by the
spatial filter SF, while the one for $\theta\geq\theta_{\rm max}$ represents
the range falling outside the NA of the lens. Here, we take values of $s=1$,
$\sigma_{L}=12.5\lambda$ and $D=20$. (b) Steady-state atomic emission as a
function of $\sqrt{s}$ for $D=\\{8,24,36\\}$ (light to dark red). Dashed lines
correspond to the modified M-B model and points are the measurements from Ref.
Ferioli  _et al._ (2023), with the global rescaling discussed in the text. (c)
Second-order photon correlations in the atomic emission as a function of
$\sqrt{s}$ and for $D=28$. The points are the data from Ref. Ferioli  _et al._
(2023), the light blue line is the original (on-axis) M-B model and the dark
blue line is the modified prediction when accounting for light measured out of
the optical axis.

To conclude, we again go back to the Cavity model to derive an approximation
for the correlations in the M-B model. We start from the definition of
$g^{(2)}(0)$ in Eq. (39) and substitute the product state ansatz from Eq.
(35). The result can be expressed in terms of sums over the expectation values
$\langle\hat{\sigma}^{z}_{n}\rangle$ and
$\langle\hat{\sigma}^{+}_{m}\rangle\langle\hat{\sigma}^{+}_{n}\rangle$ with
$m\not=n$, which we connect to the coherent $P_{\rm coh}$ and incoherent
$P_{\rm incoh}$ components of the atomic emission using Eq. (36). In the limit
$N\gg 1$, the $g^{(2)}(0)$ approximates to

$g^{(2)}(0)\approx\frac{2P_{\rm incoh}^{2}+4P_{\rm incoh}P_{\rm coh}+P_{\rm
coh}^{2}}{(P_{\rm coh}+P_{\rm incoh})^{2}},$ (40)

to leading order in $N$. From here, we use Eq. (38) to replace each component
of the atomic emission with their respective guess from the Cavity model, such
that $P_{\rm incoh}\approx\Gamma_{\rm 1D}N\bar{n}/2$ and $P_{\rm
coh}\approx\Gamma_{\rm 1D}N^{2}(\bar{n}-2\bar{n}^{2})/4$. The result is an
expression of $g^{(2)}(0)$ that only depends on the fraction of excited
population

$g^{(2)}(0;\bar{n})\approx\frac{\left[N\bar{n}+f(\bar{n})\right]^{2}-2N^{2}\bar{n}^{2}}{\left[f(\bar{n})\right]^{2}},$
(41)

where the function $f(\bar{n})=2\gamma_{\rm at}^{\rm st}(\bar{n})/\Gamma_{\rm
1D}$ using Eq. (38). It is straightforward to check that Eq. (41) yields
coherent correlations for small non-zero $\bar{n}\ll 1$, and chaotic light for
$\bar{n}=1/2$. In practice, the specific $\bar{n}$ is calculated using the M-B
model Eq. (18) for a given $s$ and $D$.

As we show in Fig. 8b, Eq. (41) reproduces remarkably well the M-B simulations
for large $s$. This is because the system becomes spatially homogeneous at
large driving intensities ($s\gg 1$), and one can represent the ensemble by
averaged quantities such as $\bar{n}$. Naturally, this approximation breaks
down at lower driving intensities $s\sim 1$, where the ensemble lacks spatial
homogeneity and undergoes phase separation, in agreement with the breakdown of
Eq. (38) observed in Fig. 8a for large optical depths.

###  IV.3 Deviations from the M-B model

The measurements Ferioli  _et al._ (2023) of the steady-state atomic emission
$\gamma_{\rm at}^{\rm st}$ and its photon correlations $g^{(2)}(0)$ deviate
significantly from what the M-B model predicts. Indeed, as we will see, the
coherent maxima of $\gamma_{\rm at}^{\rm st}$ around $s\sim 1$ in Fig. 8a are
mostly absent, and photon bunching is much stronger than what Fig. 8b
illustrates. Overall, these discrepancies suggest the coherent component of
the atomic emission in the experiment was much weaker than what the M-B model
predicts. The fact that the experimental setup fails to detect part of the
coherent atomic emission due to filtering (Fig. 7) provides a natural
explanation. Here, we present adjusted predictions of the M-B model to
qualitatively account for the filtering, such that we can compare with
experimental data from Ref. Ferioli  _et al._ (2023).

####  IV.3.1 Atomic emission out of the optical axis

Our starting point is the experimental setup from Fig. 7, where a detector
captures the atomic emission $\gamma_{\rm at}$ within a solid angle
$\Delta\Theta=2\pi\Delta\theta$. The range of $\theta_{\rm
min}<\theta<\theta_{\rm max}$ is limited from below by the spatial filter and
from above by the numerical aperture of the lens. For an aspheric lens with
NA$=0.5$ Ferioli  _et al._ (2023), the maximal captured angle is approximately
$\theta_{\rm max}\approx 0.35$ rad. Conversely, for a Gaussian beam with
$w_{0}=6.4\lambda$, the angular spread in the paraxial approximation is given
by $\theta_{\rm min}\approx\lambda/\pi w_{0}=0.05$ rad Ferioli  _et al._
(2023). We note that $\theta_{\rm min}$ was most likely larger than this
value, to ensure that no input light reached the detector. While the specific
$\theta_{\rm min}$ is unknown, we estimate that it could have been as much as
$\theta_{\rm min}\lesssim 0.2\text{ rad}$, and assume such a value going
forward.

Our goal is to calculate the measured atomic emission by integrating Eqs.
(29), which require specifying
$\langle\hat{\sigma}^{+}\hat{\sigma}^{-}\rangle$ and
$\langle\hat{\sigma}^{-}\rangle$. To streamline the discussion, we substitute
the one-dimensional M-B prediction in Eqs. (18) for the specific $D$ and $s$,
averaged over system length. We note that this is a rough approximation to
avoid dealing explicitly with the overly complex diffraction mode in three
dimensions. From the discussion surrounding Fig. 8a, where the cavity model
qualitatively agrees with the original M-B model, we expect that treating the
ensemble as a uniform spin wave only introduces quantitative differences for
the optical depths we discuss.

Figure 9a illustrates the resulting coherent and incoherent atomic intensities
as a function of angle $\theta$, for system parameters $s=1$,
$\sigma_{L}=12.5\lambda$, and $D=20$. We have also shaded the regions
$\theta<\theta_{\rm min}$ and $\theta>\theta_{\rm max}$, which are likely
excluded from detection. Interestingly, integrating these intensities within
the angular spread $\Delta\theta=[0.2,0.35]$ yields detected coherent and
incoherent intensities of the same order already for $s\sim 1$. This is in
contrast to the prediction of the M-B model, which assumes that one is
detecting on-axis and where the coherent emission should dominate (see
$\theta\sim 0$ region in Fig. 9a, where $I_{\rm coh}(\theta)\gg I_{\rm
incoh}(\theta)$).

####  IV.3.2 Steady-state atomic emission $\gamma_{\rm at}^{\rm st}$

We use the previous assumptions to discuss the measurements of the steady-
state atomic emission $\gamma_{\rm at}^{\rm st}$ from Ref. Ferioli  _et al._
(2023). The data are shown in Fig. 9b as a function of $s$ and varying $D$.
Compared to what the M-B model predicts in Fig. 8a, the large coherent maxima
are mostly absent. Moreover, for drivings $s>1$, the atomic emission appears
to stabilize, suggesting saturation is achieved much faster than in Fig. 8a.

Interestingly, these features can be qualitatively reproduced with our
modified M-B model. Following the procedure used in Fig. 9a, we use Eq. (29)
to calculate the coherent and incoherent atomic intensities. Then, we
calculate the out-of-axis atomic emission by integrating the resulting
intensities within the interval of captured angles
$\Delta\theta\sim[0.2,0.35]$ rad. The result is shown in Fig. 9b with solid
lines, where this simple model accounting for the different transverse
profiles of coherent/incoherent emission seems sufficient to produce
consistency with the experimental results. According to this interpretation,
the incoherent component dominates the atomic emission already for $s\sim 1$.

We note that measuring an absolute value for $\gamma_{\rm at}$ is challenging,
as it requires precisely estimating the coupling efficiency to the fiber
detection mode, the filtering of the input field, or the mechanical
displacement of the cloud at large intensities, among other effects Ferioli
_et al._ (2023). For that reason, the experimental values for $\gamma_{\rm
at}$ in Figs. 9b have been effectively rescaled by an arbitrary global factor.

Due to the number of experimental uncertainties and the simplified manner in
which we try to incorporate three-dimensional effects into the 1D M-B model,
we emphasize that our model should not be taken as a full, quantitatively
accurate description of the three-dimensional propagation effects. However, we
believe that it should qualitatively capture the key physics and identifies an
important subtlety of measuring the forward atomic emission.

####  IV.3.3 Second-order photon correlations $g^{(2)}(0)$

Similarly, Fig. 9c shows the measurements of the $g^{(2)}(0)$ (circles) from
Ref. Ferioli  _et al._ (2023). In the plot, we also include the original M-B
model prediction (light blue curve). The correlations measured in the
magnetized regime ($s<1$) are slightly bunched, while $g^{(2)}(0)\approx 1$ in
the M-B model. Moreover, in the unmagnetized regime ($s>1$), the bunching in
the experimental data is much higher than what the M-B model predicts.

Again, these discrepancies can possibly be explained by assuming the detection
scheme. Starting from the mean-field Eq. (40), we substitute the coherent and
incoherent components by the integrals of Eqs. (29) over the same range of
angles $\Delta\theta=[0.2,0.35]$. Plotting this result in Figure 9c (dark blue
curve), we see that Eq. (40) predicts an increase in the bunching, which stems
from the additional purely incoherent light having chaotic correlations
$g^{(2)}(0)=2$. The details of the blue curve are highly sensitive to choices
of parameters in our model, but nonetheless this curve provides a plausible
picture of the origin of the correlations seen in the experiment.

##  V Conclusion

Here, we have investigated the properties of a driven, elongated atomic
ensemble using a model based on the Maxwell-Bloch (M-B) equations. At the
mean-field level, we identified non-analytic behavior in observables in the
thermodynamic limit ($D\rightarrow\infty$), such as the steady-state
transmission and magnetization, resembling that of the driven-dissipative
Dicke phase transition. This non-analyticity in free space, however, stems
from a lack of spatial homogeneity, which in turn results from propagation
effects absent in a cavity setup. Physically, the ensemble splits into
saturated and unsaturated regions, a phenomenon closer to what we term “phase
separation” rather than a conventional phase transition.

Under specific conditions, such as low optical depth or large driving
strength, the propagation effects become negligible and the mean-field
equations of the M-B model become equivalent to those of the Cavity model.
This mapping provides a microscopic justification for the similarities between
the measurements in Ref. Ferioli  _et al._ (2023) and the driven Dicke model.
Nevertheless, we have shown that the M-B model accurately reproduces most of
the experimental data.

Moving forward, our work holds direct relevance for state-of-the-art
experimental setups Ferioli  _et al._ (2023); Lechner  _et al._ (2023).
Inspired by the mapping between M-B and Cavity models, it would be interesting
to explore the conditions under which cavity phenomena, such as superradiant
lasing Jäger  _et al._ (2021); Shankar  _et al._ (2021); Ferioli  _et al._
(2021b), spin squeezing Puri and Lawande (1979); González-Tudela and Porras
(2013) or time crystalline behaviour Iemini  _et al._ (2018), manifest in
systems with propagation effects like free space atomic ensembles or chiral
waveguide QED setups. Specifically, it would be useful to quantify the
efficiency of such applications – e.g. scaling with optical depth – before the
phase separation and other effects limit their performance. Similarly, it
would also be interesting to extend the discussion of this mapping beyond
mean-field arguments Kusmierek  _et al._ (2023); Robicheaux and Suresh (2021).

During the completion of our work, we became aware of a related recent work
that reaches similar conclusions as us using a complementary theoretical
treatment Argawal  _et al._ (2024).

##  VI Acknowledgements

The authors acknowledge stimulating discussions with Fernando de Iemini, Jamir
Marino, Francis Robicheaux, and Ana Maria Rey. D.G. acknowledges support from
the Secretaria d’Universitats i Recerca de la Generalitat de Catalunya and the
European Social Fund (2020 FI B 00196). L.B. acknowledges support from the
European Union’s Horizon Europe research and innovation program under Grant
Agreement No. 101113690 (PASQuanS2.1). E.S. acknowledges support from the
Israel Science Foundation (ISF), the Center for New Scientists at the Weizmann
Institute of Science, and the Council for Higher Education (Israel). D.E.C.
acknowledges support from the European Union, under European Research Council
grant agreement No 101002107 (NEWSPIN), FET-Open grant agreement No 899275
(DAALI) and EIC Pathfinder Grant No 101115420 (PANDA); the Government of Spain
under Severo Ochoa Grant CEX2019-000910-S (MCIN/AEI/10.13039/501100011033);
Generalitat de Catalunya (CERCA program and AGAUR Project No. 2021 SGR 01442);
Fundació Cellex, and Fundació Mir-Puig. This work is also supported by the
Agence Nationale de la Recherche (ANR-22-PETQ-0004 France 2030, project
QuBitAF), and the European Research Council (Advanced grant No.
101018511-ATARAXIA) and the Horizon Europe programme HORIZON-
CL4-2022-QUANTUM-02-SGA via the project 101113690 (PASQuanS2.1).

## References

  * Hepp and Lieb (1973) K. Hepp and E. H. Lieb, Phys. Rev. A 8, 2517 (1973).
  * Agarwal  _et al._ (1977) G. S. Agarwal, A. C. Brown, L. M. Narducci, and G. Vetri, Phys. Rev. A 15, 1613 (1977).
  * Ritsch  _et al._ (2013) H. Ritsch, P. Domokos, F. Brennecke, and T. Esslinger, Rev. Mod. Phys. 85, 553 (2013).
  * Norcia  _et al._ (2016) M. A. Norcia, M. N. Winchester, J. R. K. Cline, and J. K. Thompson, Science Advances 2, e1601231 (2016).
  * Carmichael and Walls (1977) H. J. Carmichael and D. F. Walls, Journal of Physics B: Atomic and Molecular Physics 10, L685 (1977).
  * Carmichael (1980) H. J. Carmichael, Journal of Physics B: Atomic and Molecular Physics 13, 3551 (1980).
  * Schneider and Milburn (2002) S. Schneider and G. J. Milburn, Phys. Rev. A 65, 042107 (2002).
  * González-Tudela and Porras (2013) A. González-Tudela and D. Porras, Phys. Rev. Lett. 110, 080502 (2013).
  * Kirton and Keeling (2017) P. Kirton and J. Keeling, Phys. Rev. Lett. 118, 123602 (2017).
  * Hannukainen and Larson (2018) J. Hannukainen and J. Larson, Phys. Rev. A 98, 042113 (2018).
  * Kirton  _et al._ (2019) P. Kirton, M. M. Roses, J. Keeling, and E. G. Dalla Torre, Advanced Quantum Technologies 2, 1800043 (2019).
  * Barberena  _et al._ (2019) D. Barberena, R. J. Lewis-Swan, J. K. Thompson, and A. M. Rey, Phys. Rev. A 99, 053411 (2019).
  * Somech and Shahmoon (2022) O. Somech and E. Shahmoon, “Quantum entangled states of a classically radiating macroscopic spin,” (2022), arXiv:2204.05455 [quant-ph] .
  * Sánchez Muñoz  _et al._ (2019) C. Sánchez Muñoz, B. Buča, J. Tindall, A. González-Tudela, D. Jaksch, and D. Porras, Phys. Rev. A 100, 042113 (2019).
  * Drummond and Carmichael (1978) P. Drummond and H. Carmichael, Optics Communications 27, 160 (1978).
  * Iemini  _et al._ (2018) F. Iemini, A. Russomanno, J. Keeling, M. Schirò, M. Dalmonte, and R. Fazio, Phys. Rev. Lett. 121, 035301 (2018).
  * Zhu  _et al._ (2019) B. Zhu, J. Marino, N. Y. Yao, M. D. Lukin, and E. A. Demler, New Journal of Physics 21, 073028 (2019).
  * Puri and Lawande (1979) R. Puri and S. Lawande, Physics Letters A 72, 200 (1979).
  * Walls  _et al._ (1978) D. F. Walls, P. D. Drummond, S. S. Hassan, and H. J. Carmichael, Progress of Theoretical Physics Supplement 64, 307 (1978).
  * Gegg  _et al._ (2018) M. Gegg, A. Carmele, A. Knorr, and M. Richter, New Journal of Physics 20, 013006 (2018).
  * Ferioli  _et al._ (2023) G. Ferioli, A. Glicenstein, I. Ferrier-Barbut, and A. Browaeys, Nature Physics 19, 1345 (2023).
  * Phillips  _et al._ (2001) D. F. Phillips, A. Fleischhauer, A. Mair, R. L. Walsworth, and M. D. Lukin, Phys. Rev. Lett. 86, 783 (2001).
  * Liu  _et al._ (2001) C. Liu, Z. Dutton, C. H. Behroozi, and L. V. Hau, Nature 409, 490 (2001).
  * Fleischhauer  _et al._ (2005) M. Fleischhauer, A. Imamoglu, and J. P. Marangos, Rev. Mod. Phys. 77, 633 (2005).
  * Firstenberg  _et al._ (2013) O. Firstenberg, T. Peyronel, Q.-Y. Liang, A. V. Gorshkov, M. D. Lukin, and V. Vuletić, Nature 502, 71 (2013).
  * McCall and Hahn (1967) S. L. McCall and E. L. Hahn, Phys. Rev. Lett. 18, 908 (1967).
  * Arecchi and Bonifacio (1965) F. Arecchi and R. Bonifacio, IEEE Journal of Quantum Electronics 1, 169 (1965).
  * Haake  _et al._ (1979) F. Haake, H. King, G. Schröder, J. Haus, and R. Glauber, Phys. Rev. A 20, 2047 (1979).
  * Polder  _et al._ (1979) D. Polder, M. F. H. Schuurmans, and Q. H. F. Vrehen, Phys. Rev. A 19, 1192 (1979).
  * Gross and Haroche (1982) M. Gross and S. Haroche, Physics Reports 93, 301 (1982).
  * Castin and Mølmer (1995) Y. Castin and K. Mølmer, Phys. Rev. A 51, R3426 (1995).
  * Hammerer  _et al._ (2010) K. Hammerer, A. S. Sørensen, and E. S. Polzik, Rev. Mod. Phys. 82, 1041 (2010).
  * Gorshkov  _et al._ (2011) A. V. Gorshkov, J. Otterbach, M. Fleischhauer, T. Pohl, and M. D. Lukin, Phys. Rev. Lett. 107, 133602 (2011).
  * Argawal  _et al._ (2024) S. Argawal, E. Chaparro, D. Barberena, A. Piñeiro-Orioli, G. Ferioli, S. Pancaldi, I. Ferrier-Barbut, A. Browaeys, and A. M. Rey, “private communication (article in preparation),” (2024).
  * Dalla Torre  _et al._ (2013) E. G. Dalla Torre, J. Otterbach, E. Demler, V. Vuletić, and M. D. Lukin, Phys. Rev. Lett. 110, 120402 (2013).
  * Sarkar and Satchell (1987) S. Sarkar and J. S. Satchell, Journal of Physics A: Mathematical and General 20, 2147 (1987).
  * Shammah  _et al._ (2018) N. Shammah, S. Ahmed, N. Lambert, S. De Liberato, and F. Nori, Phys. Rev. A 98, 063815 (2018).
  * Xu  _et al._ (2013) M. Xu, D. A. Tieri, and M. J. Holland, Phys. Rev. A 87, 062101 (2013).
  * Chase and Geremia (2008) B. A. Chase and J. M. Geremia, Phys. Rev. A 78, 052101 (2008).
  * Calajó  _et al._ (2019) G. Calajó, Y.-L. L. Fang, H. U. Baranger, and F. Ciccarello, Phys. Rev. Lett. 122, 073601 (2019).
  * Shen and Fan (2009) J.-T. Shen and S. Fan, Phys. Rev. A 79, 023837 (2009).
  * Shen and Fan (2005) J.-T. Shen and S. Fan, Phys. Rev. Lett. 95, 213001 (2005).
  * Calajó and Chang (2022) G. Calajó and D. E. Chang, Phys. Rev. Res. 4, 023026 (2022).
  * Mahmoodian  _et al._ (2020) S. Mahmoodian, G. Calajó, D. E. Chang, K. Hammerer, and A. S. Sørensen, Phys. Rev. X 10, 031011 (2020).
  * Shi  _et al._ (2015) T. Shi, D. E. Chang, and J. I. Cirac, Phys. Rev. A 92, 053834 (2015).
  * Caneva  _et al._ (2015) T. Caneva, M. T. Manzoni, T. Shi, J. S. Douglas, J. I. Cirac, and D. E. Chang, New Journal of Physics 17, 113001 (2015).
  * Asenjo-Garcia  _et al._ (2017) A. Asenjo-Garcia, M. Moreno-Cardoner, A. Albrecht, H. J. Kimble, and D. E. Chang, Phys. Rev. X 7, 031024 (2017).
  * Meystre and Sargent (2007) P. Meystre and M. Sargent,  _Elements of Quantum Optics_ (Springer Berlin Heidelberg, 2007).
  * Lagendijk and van Tiggelen (1996) A. Lagendijk and B. A. van Tiggelen, Physics Reports 270, 143 (1996).
  * Schilder  _et al._ (2017) N. J. Schilder, C. Sauvan, Y. R. P. Sortais, A. Browaeys, and J.-J. Greffet, Phys. Rev. A 96, 013825 (2017).
  * Birkl  _et al._ (1995) G. Birkl, M. Gatzke, I. H. Deutsch, S. L. Rolston, and W. D. Phillips, Phys. Rev. Lett. 75, 2823 (1995).
  * Bettles  _et al._ (2015) R. J. Bettles, S. A. Gardiner, and C. S. Adams, Phys. Rev. A 92, 063822 (2015).
  * Shahmoon  _et al._ (2017) E. Shahmoon, D. S. Wild, M. D. Lukin, and S. F. Yelin, Phys. Rev. Lett. 118, 113601 (2017).
  * Rui  _et al._ (2020) J. Rui, D. Wei, A. Rubio-Abadal, S. Hollerith, J. Zeiher, D. M. Stamper-Kurn, C. Gross, and I. Bloch, Nature 583, 369 (2020).
  * Manzoni  _et al._ (2018) M. T. Manzoni, M. Moreno-Cardoner, A. Asenjo-Garcia, J. V. Porto, A. V. Gorshkov, and D. E. Chang, New Journal of Physics 20, 083048 (2018).
  * Bekenstein  _et al._ (2020) R. Bekenstein, I. Pikovski, H. Pichler, E. Shahmoon, S. F. Yelin, and M. D. Lukin, Nature Physics 16, 676 (2020).
  * Moreno-Cardoner  _et al._ (2021) M. Moreno-Cardoner, D. Goncalves, and D. E. Chang, Phys. Rev. Lett. 127, 263602 (2021).
  * Glicenstein  _et al._ (2020) A. Glicenstein, G. Ferioli, N. Šibalić, L. Brossard, I. Ferrier-Barbut, and A. Browaeys, Phys. Rev. Lett. 124, 253602 (2020).
  * Ferioli  _et al._ (2021a) G. Ferioli, A. Glicenstein, L. Henriet, I. Ferrier-Barbut, and A. Browaeys, Phys. Rev. X 11, 021031 (2021a).
  * Jennewein  _et al._ (2018) S. Jennewein, L. Brossard, Y. R. P. Sortais, A. Browaeys, P. Cheinet, J. Robert, and P. Pillet, Phys. Rev. A 97, 053816 (2018).
  * Jennewein  _et al._ (2016) S. Jennewein, M. Besbes, N. J. Schilder, S. D. Jenkins, C. Sauvan, J. Ruostekoski, J.-J. Greffet, Y. R. P. Sortais, and A. Browaeys, Phys. Rev. Lett. 116, 233601 (2016).
  * Ji  _et al._ (2023) L. Ji, Y. He, Q. Cai, Z. Fang, Y. Wang, L. Qiu, L. Zhou, S. Wu, S. Grava, and D. E. Chang, “Superradiant detection of microscopic optical dipolar interactions,” (2023), arXiv:2101.10779 [physics.atom-ph] .
  * Andreoli  _et al._ (2021) F. Andreoli, M. J. Gullans, A. A. High, A. Browaeys, and D. E. Chang, Phys. Rev. X 11, 011026 (2021).
  * Scully and Zubairy (1997) M. O. Scully and M. S. Zubairy,  _Quantum Optics_ (Cambridge University Press, 1997).
  * Kusmierek  _et al._ (2023) K. J. Kusmierek, S. Mahmoodian, M. Cordier, J. Hinney, A. Rauschenbeutel, M. Schemmer, P. Schneeweiss, J. Volz, and K. Hammerer, SciPost Phys. Core 6, 041 (2023).
  * Tey  _et al._ (2009) M. K. Tey, G. Maslennikov, T. C. H. Liew, S. A. Aljunid, F. Huber, B. Chng, Z. Chen, V. Scarani, and C. Kurtsiefer, New Journal of Physics 11, 043011 (2009).
  * Goncalves  _et al._ (2021) D. Goncalves, M. W. Mitchell, and D. E. Chang, Phys. Rev. A 104, 013724 (2021).
  * Corless and Jeffrey (2002) R. Corless and D. Jeffrey, Artificial Intelligence - AI , 76 (2002).
  * Veyron  _et al._ (2022) R. Veyron, V. Mancois, J.-B. Gerent, G. Baclet, P. Bouyer, and S. Bernon, Phys. Rev. Res. 4, 033033 (2022), arXiv:2110.12505 [quant-ph] .
  * Kim and Kwak (2020) G. Y. Kim and C. H. Kwak, Journal of Optics 22, 095502 (2020).
  * Allen and Eberly (1987) L. Allen and J. H. Eberly,  _Optical resonance and two-level atoms_ (Dover, 1987).
  * Tanji-Suzuki  _et al._ (2011) H. Tanji-Suzuki, I. D. Leroux, M. H. Schleier-Smith, M. Cetina, A. T. Grier, J. Simon, and V. Vuletić, in  _Interaction between Atomic Ensembles and Optical Resonators: Classical Description_, Advances In Atomic, Molecular, and Optical Physics, Vol. 60, edited by E. Arimondo, P. Berman, and C. Lin (Academic Press, 2011) pp. 201–237.
  * Sutherland and Robicheaux (2017) R. T. Sutherland and F. Robicheaux, Phys. Rev. A 95, 033839 (2017).
  * Dicke (1954) R. H. Dicke, Phys. Rev. 93, 99 (1954).
  * Mahmoodian  _et al._ (2018) S. Mahmoodian, M. Čepulkovskis, S. Das, P. Lodahl, K. Hammerer, and A. S. Sørensen, Phys. Rev. Lett. 121, 143601 (2018).
  * Lechner  _et al._ (2023) D. Lechner, R. Pennetta, M. Blaha, P. Schneeweiss, A. Rauschenbeutel, and J. Volz, Phys. Rev. Lett. 131, 103603 (2023).
  * Jäger  _et al._ (2021) S. B. Jäger, H. Liu, J. Cooper, T. L. Nicholson, and M. J. Holland, Phys. Rev. A 104, 033711 (2021).
  * Shankar  _et al._ (2021) A. Shankar, J. T. Reilly, S. B. Jäger, and M. J. Holland, Phys. Rev. Lett. 127, 073603 (2021).
  * Ferioli  _et al._ (2021b) G. Ferioli, A. Glicenstein, F. Robicheaux, R. T. Sutherland, A. Browaeys, and I. Ferrier-Barbut, Phys. Rev. Lett. 127, 243602 (2021b).
  * Robicheaux and Suresh (2021) F. Robicheaux and D. A. Suresh, Phys. Rev. A 104, 023702 (2021).

##  Appendix A Time dynamics of atomic emission

Here, we complement the discussion in Sec. IVB in the main text on the atomic
emission $\gamma_{\rm at}$ by discussing its dynamics at early times.
Specifically, we numerically solve the time evolution of the discrete Eqs.
(17) within mean-field approximation and starting with all atoms in their
ground state. Then, we calculate the atomic emission at different times using
Eq. (36). The results below are obtained with just the original M-B model,
i.e. without any additional correction for the measurement being off-axis
(like in Sec. IV-C).

In Fig. 10a, we plot $\gamma_{\rm at}(t)$ as a function of time for a
particular realization ($s=5$, $D=26$). At very early times ($t\Gamma<1$), the
atomic emission exhibits a prominent coherent peak, which is rapidly damped
into a steady-state plateau. Since the time duration of the peak is short,
dipole-dipole interactions cannot significantly influence the dynamics.
Therefore, this coherent peak emerges solely from a combination of Rabi
oscillations and constructive interference.

Figure 10a also shows the experimental measurements of $\gamma_{\rm at}(t)$
for the same realization, obtained with the setup in Fig. 7. Here, the
experimental values have been effectively rescaled by an arbitrary global
factor as in Fig. 9 in the main text for the same reasons. In other words,
Fig. 10a shows a good qualitative agreement between the experiment and the M-B
model, despite the possibility of their peak heights being different. The
discrepancies at later times are probably related to the atomic emission being
measured out of the optical axis, and they could be corrected by modifying the
M-B model following the discussion in Sec. IV-C.

We now explore the properties of the coherent maxima as a function of optical
depth $D$ and saturation parameter $s$. In Figs. 10b,c, we plot the peak
height $\gamma_{\rm max}=\max[\gamma_{\rm at}(t)]$ and damping rate $\tau_{\rm
at}$ as a function of the saturation parameter $s$. The damping rate is
obtained by fitting the function $\gamma_{\rm at}(t)=\gamma_{\rm
max}\exp(-t/\tau_{\rm at})$ to the interval right after the intensity maximum.
Overall, their scalings with the driving intensity are well captured by the
M-B model, despite not including any additional correction for the light being
measured off-axis.

##  Appendix B Dipole-dipole interaction matrix elements

Here, we present a complementary calculation for the mapping between the M-B
and Cavity models at large driving strengths. Specifically, we study
analytically how the coherent dipole-dipole interaction $\hat{H}_{\rm dd}$ in
Eq. (16b) couples different angular momentum states. To do so, we calculate
how $\hat{H}_{\rm dd}$ acts on the total angular momentum basis $|j,m\rangle$.
We note this is not trivial since $\hat{H}_{\rm dd}$ contains only local
operators that break the degeneracy of $|j,m\rangle$.

Figure 10: (a) Atomic emission $\gamma_{\rm at}(t)$ exhibits a large coherent
peak at early times that rapidly relaxes into a plateau. The solid line is a
numerical simulation of the M-B model, while the points are experimental
measurements using the setup from Ref. Ferioli  _et al._ (2023). We have
rescaled the experimental data as indicated in the main text. Here, $s=5$ and
$D=26$. (b,c) Peak height $\max[\gamma_{\rm at}(t)]$ and damping rate
$\tau_{\rm at}$ as a function of the saturation parameter $s$ for
$D=[26,30,38]$ (blue, gray and green, respectively). We again use solid lines
for the M-B model and points for the experimental data. The data in (b) has
been rescaled as indicated in the main text.

Since each term in $\hat{H}_{\rm dd}$ acts on a pair of atoms at a time, it is
convenient to split the atomic ensemble into two subsystems: an arbitrary two-
atom pair $(n,l)$ and the remaining $N-2$ atoms. We then expand a generic $N$
atom state $|j,m\rangle$ into the basis constructed from the tensor product
between the basis of each subsystem, i.e.
$|j_{nl},m_{nl}\rangle\otimes|j_{N-2},m_{N-2}\rangle$, such that

$|j,m\rangle=\sum_{\\{j_{N-2},m_{N-2}\\}\in\mathcal{C}}|0,0;j,m\rangle\ +\\\
+C^{1,0,j_{N-2},m}_{J,m}|1,0;j_{N-2},m_{N-2}\rangle+\\\
+C^{1,1,j_{N-2},m_{N-2}}_{j,m}|1,1;j_{N-2},m_{N-2}\rangle+\\\
+C^{1,-1,j_{N-2},m_{N-2}}_{j,m}|1,-1;j_{N-2},m_{N-2}\rangle,$ (42)

where the sum is restricted over the $j_{N-2}$ and $m_{N-2}$ that satisfy the
angular-momentum composition rules $\mathcal{C}$. In Eq. (42), we have defined
the Clebsch-Gordan coefficient
$C^{j_{nl},m_{nl},j_{N-2},m_{N-2}}_{j,m}\equiv\langle
j_{nl},m_{nl};j_{N-2},m_{N-2}|j,m\rangle$, and used the relations

$C^{0,0,j^{\prime},m^{\prime}}_{j,m}=\delta_{j,j^{\prime}}\delta_{m,m^{\prime}},$
(43a)
$C^{1,0,j^{\prime},m^{\prime}}_{j,m}=C^{1,0,j^{\prime},m^{\prime}}_{j,m}\delta_{m,m^{\prime}},$
(43b)

Here, we have simplified the notation by removing the label of the atom pair
$(n,l)$.

Having defined the Clebsch-Gordan coefficients between the basis, we now
calculate how the element ${\hat{H}_{{\rm
dd},nl}\equiv\sigma^{+}_{n}\sigma^{-}_{l}-\sigma^{+}_{l}\sigma^{-}_{n}}$,
appearing in the coherent dipole-dipole interaction ${\hat{H}_{\rm
dd}=-i\hbar\Gamma_{\rm 1D}/4\sum_{n>l}\hat{H}_{\rm dd,nl}}$, acts on the
angular momentum states. For simplicity, we focus only on the maximal angular
momentum states $j=N/2$, obtaining

$\hat{H}_{\rm dd}|N/2,m\rangle=i\hbar\Gamma_{\rm 1D}F(m)|\psi\rangle.$ (44)

Here, we have defined

$F(m)|\psi\rangle=\frac{C^{1,0,\left(N/2-1\right),m}_{N/2,m}}{4}\sum_{l>m}|0,0;(N/2-1),m\rangle_{n,l},$
(45)

where we exploited the permutational symmetry of the $j=N/2$ states and
imposed $\langle\psi|\psi\rangle=1$.

Using the Clebsh-Gordon coefficients defined above, one can show that
$\hat{H}_{\rm dd}$ couples the state $|N/2,m\rangle$ with the corresponding
unique state $|\psi\rangle$ in the $j=N/2-1$ manifold with strength

$F(m)=\sqrt{\frac{(N/2+m)\left(N/2-m\right)\left(N+1\right)}{24}}.$ (46)

Equation (46) provides the scaling of the coherent dipole-dipole matrix
elements with system size. In the thermodynamic limit $N\rightarrow\infty$,
the coupling scales linearly with atom number $F\propto N$ at the extremes of
the $j=N/2$ manifold ($|m|\approx N/2$), while $F\propto N^{3/2}$ at the
center ($m\approx 0$). Next, we identify two distinct regimes depending on how
the driving strength and collective decay compare to these scalings, as we
show in Fig. 11.

For weak drivings ($\Omega\lesssim\Gamma\sqrt{N}$), the effect of coherent
dipole-dipole interactions is not negligible, as $\hat{H}_{\rm dd}$ mixes
states from different angular momentum manifolds even for small number of
excitations. This suggests that propagation effects matter, and is consistent
with both our purity discussion based on the lateral decay rates (Sec. II.3)
and the discrepancy between the Cavity and M-B models in Sec. III.4.

On the other hand, for strong drivings ($\Omega\sim\Gamma N$), the effect of
$\hat{H}_{\rm dd}$ becomes negligible for experimentally realistic
cooperativities ($\Gamma\gg\Gamma_{\rm 1D}$). This result complements the
intuition developed in Sec. III.4, where the Cavity and M-B models converge as
propagation effects vanish.

Figure 11: Matrix elements of the coherent driving Hamiltonian $\hat{H}_{\rm
drive}$ (blue arrows), of the collective dissipation
$\mathcal{L}_{\Downarrow}$ (red arrows) and of the coherent dipole-dipole
interactions $\hat{H}_{\rm dd}$ (violet arrows) for Dicke states with $m=0$
and $m=-N/2+1$.