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

D. Goncalves, L. Bombieri, G. Ferioli, S. Pancaldi, Igor Ferrier-Barbut, A.

Browaeys, E. Shahmoon, D.E. Chang

► To cite this version:

D. Goncalves, L. Bombieri, G. Ferioli, S. Pancaldi, Igor Ferrier-Barbut, et al. Driven-dissipative phase separation in free-space atomic ensembles. PRX Quantum, 2024, 6 (2), pp.020303. 10.1103/PRXQuantum.6.020303. [https://arxiv.org/abs/2003.05348]

HAL Id: hal-05348482

https://hal.science/hal-05348482v1

Submitted on 5 Nov 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Driven-dissipative phase separation in free-space atomic ensembles

D. Goncalves,¹ L. Bombieri,¹,²,³ G. Ferioli,⁴ S. Pancaldi,⁴

I. Ferrier-Barbut,⁴ A. Browaeys,⁴ E. Shahmoon,⁵ and D. E. Chang¹,⁶

¹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

⁴Université Paris-Saclay, Institut d'Optique Graduate School,
CNRS, Laboratoire Charles Fabry, 91127, Palaiseau, France.

⁵Department of Chemical & Biological Physics, Weizmann Institute of Science, Rehovot 7610001, Israel.

⁶ICREA-Institució Catalana de Recerca i Estudis Avançats, 08015 Barcelona, Spain.

(Dated: March 25, 2024)

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 [1–4], 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 [5–13], governed by the parameter $\beta = 2\Omega/N\Gamma_{1D}$, with $\Omega$ and $\Gamma_{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 [12, 14] and so-called time crystalline behavior [15–17]. The magnetized phase exhibits spin squeezing [8, 18], 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$ [5, 6, 15, 19, 20].

Surprisingly, a recent experiment [21] 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

<figure><img src="image_2.png" /><figcaption>FIG. 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_{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_{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.</figcaption></figure>
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 [22–25].

Here, we examine the behavior of a driven elongated atomic ensemble in free space to elucidate the possible ex- istence 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 [26–33]. 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 dissipa- tion, but features an extra coherent (dipole-dipole) in- teraction 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 ex- hibits what we term a “phase separation”, where different spatial regions of the ensemble acquire different proper- ties (saturated vs. magnetized) due to the spatial degree of freedom. Finally, we compare our M-B model with the experimental results of Ref. [21], obtaining a good agreement.

We highlight the complementary work from Ref. [34], which directly models a three-dimensional atomic ensem- ble, 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 fea- tures. Assuming that the cavity mode is lossy, and start- ing from the full atom-light master equation, one can trace out the cavity field within a Born-Markov approxi- mation to obtain a master equation for the atomic density

matrix $\hat{\rho}$ that reads [2, 5–7, 15]

$$ \dot{\hat{\rho}} = -\frac{i}{\hbar} [\hat{H}_{\text{drive}}, \hat{\rho}] + \mathcal{L}_{\downarrow}[\hat{\rho}]. \quad (1) $$

The Hamiltonian term

$$ \hat{H}_{\text{drive}} = \frac{\hbar \Omega}{2} (\hat{S}_{+} + \hat{S}_{-}), \quad (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 col- lective 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_{1D}}{2} \left( 2\hat{S}_{-}\hat{\rho}\hat{S}_{+} - \{\hat{\rho}, \hat{S}_{+}\hat{S}_{-}\} \right) \quad (3) $$

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

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

$$ \hat{S}^2 |j, m\rangle = j(j+1) |j, m\rangle, \quad (4a) $$

$$ \hat{S}_z |j, m\rangle = m |j, m\rangle, \quad (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.

### A. Driven-dissipative phase transition

As discussed in the introduction, the competition be- tween the coherent driving $\Omega$ and the collective dissi- pation $\Gamma_{1D}$ in Eq. (1) yields a driven-dissipative phase transition [5–13, 15]. As a function of the ratio $\beta = 2\Omega/N\Gamma_{1D}$, the steady-state solution of the density ma- trix 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 fo- cus 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 col- lective spin operators in the steady-state read

$$ \Omega \langle \hat{S}_z \rangle = \Gamma_{1D} \langle \hat{S}_z \rangle \langle \hat{S}_y \rangle, \quad (5a) $$

$$ \Omega \langle \hat{S}_y \rangle = \Gamma_{1D} (\langle \hat{S}^2 \rangle - \langle \hat{S}_z \rangle^2), \quad (5b) $$

where $\langle \hat{S}_x \rangle = 0$. To solve Eqs. (5), one usually ex- ploits the conservation of angular momentum to substi- tute $\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}, \quad \beta < 1. \quad (6) $$
<figure><img src="image_3.png" /><figcaption>FIG. 2. (a) Normalized steady-state magnetization $\langle \hat{S}_z \rangle/N$ as a function of $\beta = 2\Omega/N\Gamma_{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_{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_{j,m} \to \Gamma_{j,m}^\leftarrow$ and $\Gamma_{j,m} \to \Gamma_{j,m}^\rightarrow$ that modify $j$. (c) Ratio between the angular momentum-changing decay rates $\Gamma_{j,m}^\to/\Gamma_{j,m}^\leftarrow$ 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_{1D}$ (squares, points and triangles).</figcaption></figure>

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 [8, 10, 12], exact solution [15, 18], or an analysis based on approximate eigenstates [13]. Additional hallmarks of the phase transition are the steady-state density matrix going from a pure to completely mixed state [8, 18] and from squeezed to unsqueezed [8, 35] at the same critical point.

The unmagnetized phase ($\beta > 1$) also features super-radiant 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_{1D}\langle \hat{S}^+\hat{S}^- \rangle = \Gamma_{1D}\text{Tr}(\hat{S}^+\hat{S}^-)/(N+1) \propto N^2$. Notably, this super-radiance 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 [14], an effect which has been termed as a boundary time crystal [15–17].

## B. 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} [\hat{H}_{\text{drive}}, \hat{\rho}] + \mathcal{L}_{\downarrow}[\hat{\rho}] + \mathcal{L}_{\downarrow}[\hat{\rho}], \quad (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). \quad (7b) $$

Equations (7) will be referred to as the Cavity model, where atoms decay collectively into the cavity mode at a rate $\Gamma_{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 $O(N^3)$ elements instead of $O(2^N)$ [36–39].

## C. 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', m-1 | \dot{\hat{\rho}} | j', m-1 \rangle = \sum_{j,m} \Gamma_{j,m}^{\alpha} \langle j, m | \hat{\rho} | j, m \rangle + \dots \quad (8) $$

This reveals the rates $\Gamma_{j,m}^\alpha$ at which population is transferred from the state $|j, m\rangle$ into the states $|j', m-1\rangle$ via
spontaneous emission, where $\alpha = \{\leftarrow, \rightarrow, \downarrow\}$ indicates whether the angular momentum $j'$ 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_{j,m}^{\leftarrow}$ and $\Gamma_{j,m}^{\rightarrow}$, which change the value of $j$ by +1 and -1, respectively, and whose explicit form can be found in the literature to be [37]

$$ \Gamma_{j,m}^{\leftarrow} = \Gamma \frac{(j-m+1)(j-m+2)(\frac{N}{2}-j)}{2(j+1)(2j+1)}, \quad (9a) $$

$$ \Gamma_{j,m}^{\rightarrow} = \Gamma \frac{(j+m-1)(j+m)(j+1+\frac{N}{2})}{2j(2j+1)}. \quad (9b) $$

Figure 2c illustrates how the ratio between these decay rates depends on $j$ and $m$. In the red region, $\Gamma_{j,m}^{\leftarrow} \gg \Gamma_{j,m}^{\rightarrow}$ 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_{j,m}^{\leftarrow} = \Gamma_{j,m}^{\rightarrow}$, one finds a universal relation between angular momentum and magnetization in the steady state

$$ j = \frac{N}{2} \sqrt{1 - (1 + 2m/N)^2}, \quad (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_{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_{j,m}^{\leftarrow} < \Gamma_{j,m}^{\rightarrow}$. 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}}, \quad (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 \to \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$).

### D. 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_{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 \neq 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}, \quad \beta < \frac{1}{\sqrt{2}}, \quad (12) $$

are stable. Equation (12) indicates a first-order phase transition at a shifted critical point $\beta_c' = 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_c'$ prediction in Eq. (12). Again, these findings are supported by other methods such as a factorization ansatz [6, 15, 19] and numerics [20, 36].

### III. ENSEMBLE IN FREE SPACE

We now introduce our minimal model for the experiment reported in Ref. [21], 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 [26–33]. 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.

#### A. 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 Hamil-
tonian, including both atoms and photons, reads [40–44]

$$
\begin{equation}
\begin{split}
\hat{H} = {}& -i\hbar c \int dz [\hat{a}^{\dagger}(z) \partial_z \hat{a}(z)] + \\
& \qquad + \hbar \sqrt{\frac{c\Gamma_{1D}}{2}} \sum_{n=1}^{N} [\hat{a}^{\dagger}(z_n) \hat{\sigma}_{n}^{-} + \text{h.c.}].
\end{split}
\tag{13}
\end{equation}
$$

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')] = \delta(z - z')$. 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_{1D}$.

Next, we use Eq. (13) to calculate the Heisenberg equa-
tions 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_{1D}}{2c}} \sum_{n=1}^{N} \hat{\sigma}_{n}^{-} \delta(z - z_{n}), \quad (14a)
$$

$$
\frac{\partial}{\partial t} \hat{\sigma}_{n}^{-} = i \sqrt{\frac{c \Gamma_{1D}}{2}} \hat{a}(z_n) \hat{\sigma}_{n}^{z}, \quad (14\text{b})
$$

$$
\frac{\partial}{\partial t} \hat{\sigma}_n^z = i \sqrt{2c\Gamma_{1D}} [\hat{a}^\dagger(z_n) \hat{\sigma}_n^- - \text{h.c.}] . \quad (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 prop-
agation occurs in a time-scale much shorter than the
atomic dynamics, as Γ₁D ≫ 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 approx-
imation [30, 45, 46]. Integrating over space yields the
input-output relation

$$
\hat{a}(z) = a_{\text{in}} - i \sqrt{\frac{\Gamma_{1\text{D}}}{2c}} \sum_{n=1}^{N} \hat{\sigma}_{n}^{-} \theta(z - z_{n}), \quad (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_{\text{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 descrip-
tion is possible (see, e.g., Refs. [47, 48] for more details).

In addition to the propagating mode of interest, we as-
sume that the atomic interaction with the remaining free-
space modes can be modeled as independent emission Γ,
as in the Cavity model of Sec. II B. The resulting master
equation for the atomic density matrix $\hat{\rho}$ reads [43]

$$
\dot{\hat{\rho}} = -\frac{i}{\hbar} [(\hat{H}_{\text{drive}} + \hat{H}_{\text{dd}}), \hat{\rho}] + \mathcal{L}_{\downarrow}^{\prime}[\hat{\rho}] + \mathcal{L}_{\downarrow}[\hat{\rho}], \quad (16\text{a})
$$

$$
\hat{H}_{\text{dd}} = -\frac{i\hbar\Gamma_{1\text{D}}}{4} \sum_{n>m} [\hat{\sigma}_n^+ \hat{\sigma}_m^- - \hat{\sigma}_m^+ \hat{\sigma}_n^-], \quad (16\text{b})
$$

$$
\mathcal{L}'_{\downarrow}[\hat{\rho}] = \frac{\Gamma_{1\mathrm{D}}}{4} \left[ 2\hat{S}^{-} \hat{\rho} \hat{S}^{+} - \{\hat{\rho}, \hat{S}^{+} \hat{S}^{-}\} \right], \quad (16\mathrm{c})
$$

where $\hat{H}_{\text{drive}} = \hbar \sqrt{c\Gamma_{1D}/2} \sum_n (a_{\text{in}}^*(z_n) \hat{\sigma}_n^- + \text{h.c.})$ describes the driving by the input field and $\hat{H}_{\text{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}_{\text{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 [3]. However, the $\hat{H}_{\text{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 prop-
agation in atomic ensembles. As the emission into the
modes not of explicit interest is considered to be indepen-
dent, it neglects effects associated with wave interference
and multiple scattering of light in these modes. For ex-
ample, this model cannot capture the “speckled” pattern
of the scattered far-field intensity when the atomic posi-
tions are disordered [49, 50], nor can it capture Bragg
diffraction [51] or the perfect reflection of light from or-
dered arrays [52–54], which has recently been realized to
be a valuable resource for applications [55–57]. 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 en-
ergies and decay rates [47, 58–61], fundamental dephas-
ing mechanisms [62], or limitations to the ensemble re-
fRACTIVE index [63].

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 in-dependent emission. We also again point to the complementary work of Ref. [34], where the three-dimensional nature of the problem is treated in more detail, and which arrives at similar conclusions.

## B. 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 [64]

$$ \partial_z \hat{E}(z) = -i \sqrt{\frac{\Gamma_{1D}}{2}} \frac{N}{L} \hat{\sigma}^-(z), \quad (17a) $$

$$ \dot{\hat{\sigma}}^-(z) = -\frac{\Gamma}{2} \hat{\sigma}^-(z) + i \sqrt{\frac{\Gamma_{1D}}{2}} \hat{E}(z) \hat{\sigma}^z(z), \quad (17b) $$

$$ \dot{\hat{\sigma}}^z(z) = -\Gamma(\hat{\sigma}^z(z) + 1) + i\sqrt{2\Gamma_{1D}}(\hat{E}^\dagger(z)\hat{\sigma}^-(z) - \text{h.c.}), \quad (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 [30, 65]. Taking $\Gamma \gg \Gamma_{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)}, \quad (18) $$

where $\eta = \Gamma_{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 [66, 67],

$$ \eta = \frac{3}{8\pi} \frac{\lambda^2}{A}, \quad (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 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_{in}$, the amplitude of the field within a 1D ensemble spanning $z \in [0, L]$ obeys

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

Here, $W(x)$ is the Wright Omega function [65, 68], 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 [69, 70].

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_{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)}. \quad (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_{in}$ remains roughly constant throughout the ensemble, yielding $T \approx 1$.

## C. 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 \to \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 \ge 1, \end{cases} \qquad (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_{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) dz$ in the ensemble by substituting Eq. (20) into Eq. (18). In the thermodynamic limit, one has

$$\bar{s}_{\infty}^{z}(s) \approx \begin{cases} \frac{s-1}{2} & \text{for } s < 1, \\ 0 & \text{for } s \ge 1. \end{cases} \qquad (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.

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) dz + \frac{1}{2L} \int_{z_c}^{L} \sigma^{z}(z) dz = \frac{s-1}{2}, \quad (24)$$

which is exactly Eq. (23). This is in stark contrast with the magnetized phase existing in the Cavity and Dicke

<figure><img src="image_4.png" /><figcaption>FIG. 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.</figcaption></figure>

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 \to 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_{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.
## D. 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_{in}$ and substitute the result in Eq. (17b). Then, we use the relation $\Omega = \sqrt{2\Gamma_{1D}}E_{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}_i^\alpha\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 \hat{S}_{-}^{\text{Cav}} \rangle \approx -\frac{\Gamma N}{2} \langle \hat{\sigma}^{-} \rangle + \frac{i\Omega N}{2} \langle \hat{\sigma}^{z} \rangle + \frac{\Gamma_{1D} N^2}{2} \langle \hat{\sigma}^{z} \rangle \langle \hat{\sigma}^{-} \rangle, \quad (25a) $$

$$ \langle \hat{S}_{-}^{\text{MB}} \rangle \approx -\frac{\Gamma N}{2} \langle \hat{\sigma}^{-} \rangle + \frac{i\Omega N}{2} \langle \hat{\sigma}^{z} \rangle + \frac{\Gamma_{1D} N^2}{4} \langle \hat{\sigma}^{z} \rangle \langle \hat{\sigma}^{-} \rangle. \quad (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}_{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}_{dd} \rangle \approx \frac{\Gamma_{1D} N^2}{8} (\langle \hat{\sigma}^+ \rangle \langle \hat{\sigma}^- \rangle - \langle \hat{\sigma}^+ \rangle \langle \hat{\sigma}^- \rangle) = 0. \quad (26) $$

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

<figure><img src="image_5.png" /><figcaption>FIG. 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.</figcaption></figure>

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}_n^\pm\rangle$ vanish and the chiral dipole-dipole interactions become negligible $\langle\hat{H}_{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}_{dd}$ acts on the angular momentum states in the strong driving limit.

## IV. MODELING OF EXPERIMENTAL DATA

We now reproduce some key experimental results from Ref. [21] 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$ [21]. 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 [30]. 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_{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, [71]

$$P = \int_{\Delta\Theta} d\Omega_{\mathbf{k}} I_N(\mathbf{k}), \qquad (27)$$

where $I_N(\mathbf{k})$ is the intensity emitted in the direction of the wavevector $\mathbf{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} d\Omega_{\mathbf{k}} I_1(\mathbf{k}) \equiv \Gamma \langle \hat{\sigma}^{+} \hat{\sigma}^{-} \rangle, \qquad (28)$$

where $I_1(\mathbf{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_{in}z_n}$, where $k_{in}$ is the input field wavevector. For $N \gg 1$, the intensity $I_N(\mathbf{k})$ emitted in the direction $\mathbf{k}$ can be separated into coherent and incoherent components as [21, 72, 73]

$$I_{\text{coh}}(\mathbf{k}) = I_1(\mathbf{k}) N^2 |\langle \hat{\sigma}^- \rangle|^2 F(\mathbf{k}, \mathbf{k}_{\text{in}}), \qquad (29a)$$

$$I_{\text{incoh}}(\mathbf{k}) = I_1(\mathbf{k}) N \langle \hat{\sigma}^+ \hat{\sigma}^- \rangle, \qquad (29b)$$

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

$$F(\mathbf{k}, \mathbf{k}_{\text{in}}) = \frac{1}{N^2} \sum_{n \neq m} e^{ik_{\text{in}}(1-\cos\theta)(z_n-z_m)}. \qquad (30)$$

Physically, $F(\mathbf{k}, \mathbf{k}_{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 $\mathbf{k}_{in}$. Indeed, assuming $z_n$ follow a

<figure><img src="image_6.png" /><figcaption>FIG. 5. Form factor $F(\mathbf{k}, \mathbf{k}_{in})$ as a function of the angle $\theta$ between $\mathbf{k}$ and the input field $\mathbf{k}_{in}$, and calculated for typical ensemble sizes of the experiment of Ref. [21]. 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.</figcaption></figure>

Gaussian distribution $\rho(z; \sigma_L)$ with variance $\sigma_L^2$ and zero mean, Eq. (30) in the continuous limit becomes [21]

$$F(\theta) \approx \left| \int_{-\infty}^{\infty} dz \, \rho(z; \sigma_L) e^{ik_{\text{in}}(1-\cos\theta)z} \right|^2 = \\ = \exp[-\sigma_L^2 k_{\text{in}}^2 (1 - \cos\theta)^2]. \qquad (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_{ax} \sim 20-25\lambda$ in Ref. [21], consistently with the discussion surrounding the Appendix C of Ref. [21].

We then estimate the effective single-atom cooperativity $\eta_{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_{\text{eff}} = \frac{\int_{k_z>0} d\Omega_{\mathbf{k}} F(\mathbf{k}, \mathbf{k}_{\text{in}}) I_1(\mathbf{k})}{\int_{k_z>0} d\Omega_{\mathbf{k}} I_1(\mathbf{k})}. \qquad (32)$$

Intuitively, $\eta_{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$) [74]. On the other hand, $\eta_{eff} = 0$ when atoms can have arbitrary relative phases and do not constructively interfere. In the case of the experiment [21], where atoms were circularly polarized $I_1(\mathbf{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_{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. [21], 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. [21], using
the relation $D = 2\eta_{\text{eff}}N$. For the saturation parameter $s$, we use the relation

$$ \beta = \frac{2\Omega}{N\eta_{\text{eff}}\Gamma} \equiv \sqrt{\frac{8s}{D}}, \qquad (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.

### A. 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. [21]. 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. [21]. 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. [21], and thus we cannot directly compare theory and experiment.

Nevertheless, one can infer an approximated value for $T$ from the measurements of $\Omega_{\text{eff}}$ in Ref. [21]. There, $\Omega_{\text{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_{\text{eff}} \approx \Omega - i\Gamma_{\text{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_{\text{eff}}/\Omega)^2$.

Comparing the transmittance from the M-B model and the one inferred from the measurements of $\Omega_{\text{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><img src="image_7.png" /><figcaption>FIG. 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. [21], 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. [21]. 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_{\text{eff}}/\Omega)^2$ from Ref. [21], 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$.</figcaption></figure>

### B. Atomic emission

We now discuss the intensity of the light emitted by the atomic cloud in the forward (transmitted) direction, $\gamma_{\text{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_{\text{at}} = \frac{\Gamma_{\text{1D}}}{2} \langle \hat{S}^+ \hat{S}^- \rangle. \qquad (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 [21]. 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_{\text{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_{\text{at}}$ within the ideal M-B model.
<figure><img src="image_8.png" /><figcaption>FIG. 7. Sketch of the experimental setup in Ref. [21] for the detection of the atomic field intensity $\gamma_{at}$. A spatial filter (SF) blocks the input field component (dark red) by eliminating modes with $k \approx k_{in}$. The measured atomic emission comprises modes slightly out of the optical axis $k \not\approx k_{in}$ that still fall within the numerical aperture of the lens, defining the solid angle $\Delta\Theta$. The measured $\gamma_{at}$ thus contains both coherent (light red) and incoherent (lighter red) atomic emission.</figcaption></figure>

## 1. Product state ansatz

To calculate quantities like atomic emission $\gamma_{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 [65]

$$ \hat{\rho} = \bigotimes_{n=1}^{N} \hat{\rho}_n, \qquad \hat{\rho}_n = \frac{1}{2} \begin{pmatrix} 1 - \langle \hat{\sigma}_n^z \rangle & 2 \langle \hat{\sigma}_n^- \rangle \\ 2 \langle \hat{\sigma}_n^- \rangle^* & 1 + \langle \hat{\sigma}_n^z \rangle \end{pmatrix}. \quad (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_{at} = \frac{\Gamma_{1D}}{2} \sum_{n=1}^{N} \frac{\langle \hat{\sigma}_n^z \rangle + 1}{2} + \frac{\Gamma_{1D}}{2} \sum_{n \neq m} \langle \hat{\sigma}_n^+ \rangle \langle \hat{\sigma}_m^- \rangle. \quad (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).

## 2. Steady-state atomic emission

In the following, we describe the properties of the steady-state atomic emission $\gamma_{at}^{st}$ in Eq. (36), while we

discuss its early-time dynamics in Appendix A. Figure 8a shows $\gamma_{at}^{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 iluminated atoms can radiate, such that $\gamma_{at}^{st} \propto \Gamma N z_c / L \propto \Gamma Ds \propto |E_{in}|^2$ independently of system size. Conversely, deep in the unmagnetized regime ($s \gg 1$), atomic coherence vanishes and $\gamma_{at}^{st} = \Gamma_{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_{at}^{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_{m,j}^{out}$ from $|j, m\rangle$ to any other state $|j', m-1\rangle$ [37], 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_{j,m}^{\text{out}}(n) \approx \Gamma_{1D} N^2 (n - 2n^2) + \Gamma N n, \quad n \le \frac{1}{2}, \quad (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 D, we map Eq. (37) to the M-B model by replacing $\Gamma$ by $\Gamma_{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_{at}^{st}(\bar{n}) \approx \frac{\Gamma_{1D}}{2} \left( \frac{N^2}{2} (\bar{n} - 2\bar{n}^2) + N\bar{n} \right), \quad (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
<figure><img src="image_9.png" /><figcaption>FIG. 8. (a) Steady-state atomic emission $\gamma_{at}^{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 from the mapping to the Cavity model [Eq. (38)] for the same optical depths. Inset: Scaling of the steady-state maxima $\max[\gamma_{at}^{st}]$ 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$.</figcaption></figure>

prediction from Eq. (38) breaks down at larger optical depths, where the mapping between the Cavity and M-B models is no longer valid.

### 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}. \qquad (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 [13]. 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 [75].

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}_n^z \rangle$ and $\langle \hat{\sigma}_m^+ \rangle \langle \hat{\sigma}_n^+ \rangle$ with $m \neq n$, which we connect to the coherent $P_{\text{coh}}$ and incoherent $P_{\text{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_{\text{incoh}}^2 + 4P_{\text{incoh}}P_{\text{coh}} + P_{\text{coh}}^2}{(P_{\text{coh}} + P_{\text{incoh}})^2}, \qquad (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_{\text{incoh}} \approx \Gamma_{1D}N\bar{n}/2$ and $P_{\text{coh}} \approx \Gamma_{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{[N\bar{n} + f(\bar{n})]^2 - 2N^2\bar{n}^2}{[f(\bar{n})]^2}, \qquad (41)$$

where the function $f(\bar{n}) = 2\gamma_{at}^{st}(\bar{n})/\Gamma_{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
<figure><img src="image_10.png" /></figure>

FIG. 9. (a) Coherent $I_{\text{coh}}(\theta)$ and incoherent $I_{\text{incoh}}(\theta)$ intensity components of the atomic emission as a function of the angle $\theta$ between the detected mode $\mathbf{k}$ and the optical axis $\mathbf{k}_{\text{in}}$ defined by the input field. The grey region $\theta \le \theta_{\text{min}}$ represents the range of angles blocked by the spatial filter SF, while the one for $\theta \ge \theta_{\text{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. [21], 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. [21], 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.

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.

### C. Deviations from the M-B model

The measurements [21] of the steady-state atomic emission $\gamma_{\text{at}}^{\text{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_{\text{at}}^{\text{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. [21].

#### 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_{\text{at}}$ within a solid angle $\Delta\Theta = 2\pi\Delta\theta$. The range of $\theta_{\text{min}} < \theta < \theta_{\text{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 [21], the maximal captured angle is approximately $\theta_{\text{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_{\text{min}} \approx \lambda/\pi w_0 = 0.05$ rad [21]. We note that $\theta_{\text{min}}$ was most likely larger than this value, to ensure that no input light reached the detector. While the specific $\theta_{\text{min}}$ is unknown, we estimate that it could have been as much as $\theta_{\text{min}} \lesssim 0.2$ 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_{\text{min}}$ and $\theta > \theta_{\text{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_{\text{coh}}(\theta) \gg I_{\text{incoh}}(\theta)$).
## 2. Steady-state atomic emission $\gamma_{at}^{st}$

We use the previous assumptions to discuss the measurements of the steady-state atomic emission $\gamma_{at}^{st}$ from Ref. [21]. 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_{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 [21]. For that reason, the experimental values for $\gamma_{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.

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

Similarly, Fig. 9c shows the measurements of the $g^{(2)}(0)$ (circles) from Ref. [21]. 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 \to \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. [21] 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 [21, 76]. 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 [77–79], spin squeezing [8, 18] or time crystalline behaviour [16], 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 [65, 80].

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 [34].

# 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).

[1] K. Hepp and E. H. Lieb, Phys. Rev. A **8**, 2517 (1973).

[2] G. S. Agarwal, A. C. Brown, L. M. Narducci, and G. Vetri, Phys. Rev. A **15**, 1613 (1977).

[3] H. Ritsch, P. Domokos, F. Brennecke, and T. Esslinger, Rev. Mod. Phys. **85**, 553 (2013).

[4] M. A. Norcia, M. N. Winchester, J. R. K. Cline, and J. K. Thompson, Science Advances **2**, e1601231 (2016).

[5] H. J. Carmichael and D. F. Walls, Journal of Physics B: Atomic and Molecular Physics **10**, L685 (1977).

[6] H. J. Carmichael, Journal of Physics B: Atomic and Molecular Physics **13**, 3551 (1980).

[7] S. Schneider and G. J. Milburn, Phys. Rev. A **65**, 042107 (2002).

[8] A. González-Tudela and D. Porras, Phys. Rev. Lett. **110**, 080502 (2013).

[9] P. Kirton and J. Keeling, Phys. Rev. Lett. **118**, 123602 (2017).

[10] J. Hannukainen and J. Larson, Phys. Rev. A **98**, 042113 (2018).

[11] P. Kirton, M. M. Roses, J. Keeling, and E. G. Dalla Torre, Advanced Quantum Technologies **2**, 1800043 (2019).

[12] D. Barberena, R. J. Lewis-Swan, J. K. Thompson, and A. M. Rey, Phys. Rev. A **99**, 053411 (2019).

[13] O. Somech and E. Shahmoon, “Quantum entangled states of a classically radiating macroscopic spin,” (2022), arXiv:2204.05455 [quant-ph].

[14] 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).

[15] P. Drummond and H. Carmichael, Optics Communications **27**, 160 (1978).

[16] F. Iemini, A. Russomanno, J. Keeling, M. Schirò, M. Dalmonte, and R. Fazio, Phys. Rev. Lett. **121**, 035301 (2018).

[17] B. Zhu, J. Marino, N. Y. Yao, M. D. Lukin, and E. A. Demler, New Journal of Physics **21**, 073028 (2019).

[18] R. Puri and S. Lawande, Physics Letters A **72**, 200 (1979).

[19] D. F. Walls, P. D. Drummond, S. S. Hassan, and H. J. Carmichael, Progress of Theoretical Physics Supplement **64**, 307 (1978).

[20] M. Gegg, A. Carmele, A. Knorr, and M. Richter, New Journal of Physics **20**, 013006 (2018).

[21] G. Ferioli, A. Glicenstein, I. Ferrier-Barbut, and A. Browaeys, Nature Physics **19**, 1345 (2023).

[22] D. F. Phillips, A. Fleischhauer, A. Mair, R. L. Walsworth, and M. D. Lukin, Phys. Rev. Lett. **86**, 783 (2001).

[23] C. Liu, Z. Dutton, C. H. Behroozi, and L. V. Hau, Nature **409**, 490 (2001).

[24] M. Fleischhauer, A. Imamoglu, and J. P. Marangos, Rev. Mod. Phys. **77**, 633 (2005).

[25] O. Firstenberg, T. Peyronel, Q.-Y. Liang, A. V. Gorshkov, M. D. Lukin, and V. Vuletić, Nature **502**, 71 (2013).

[26] S. L. McCall and E. L. Hahn, Phys. Rev. Lett. **18**, 908 (1967).

[27] F. Arecchi and R. Bonifacio, IEEE Journal of Quantum Electronics **1**, 169 (1965).

[28] F. Haake, H. King, G. Schröder, J. Haus, and R. Glauber, Phys. Rev. A **20**, 2047 (1979).

[29] D. Polder, M. F. H. Schuurmans, and Q. H. F. Vrehen, Phys. Rev. A **19**, 1192 (1979).

[30] M. Gross and S. Haroche, Physics Reports **93**, 301 (1982).

[31] Y. Castin and K. Mølmer, Phys. Rev. A **51**, R3426 (1995).

[32] K. Hammerer, A. S. Sørensen, and E. S. Polzik, Rev. Mod. Phys. **82**, 1041 (2010).

[33] A. V. Gorshkov, J. Otterbach, M. Fleischhauer, T. Pohl, and M. D. Lukin, Phys. Rev. Lett. **107**, 133602 (2011).

[34] 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).

[35] E. G. Dalla Torre, J. Otterbach, E. Demler, V. Vuletić, and M. D. Lukin, Phys. Rev. Lett. **110**, 120402 (2013).

[36] S. Sarkar and J. S. Satchell, Journal of Physics A: Mathematical and General **20**, 2147 (1987).

[37] N. Shammah, S. Ahmed, N. Lambert, S. De Liberato, and F. Nori, Phys. Rev. A **98**, 063815 (2018).

[38] M. Xu, D. A. Tieri, and M. J. Holland, Phys. Rev. A **87**, 062101 (2013).

[39] B. A. Chase and J. M. Geremia, Phys. Rev. A **78**, 052101 (2008).

[40] G. Calajó, Y.-L. L. Fang, H. U. Baranger, and F. Ciccarello, Phys. Rev. Lett. **122**, 073601 (2019).

[41] J.-T. Shen and S. Fan, Phys. Rev. A **79**, 023837 (2009).

[42] J.-T. Shen and S. Fan, Phys. Rev. Lett. **95**, 213001 (2005).

[43] G. Calajó and D. E. Chang, Phys. Rev. Res. **4**, 023026 (2022).

[44] S. Mahmoodian, G. Calajó, D. E. Chang, K. Hammerer, and A. S. Sørensen, Phys. Rev. X **10**, 031011 (2020).
[45] T. Shi, D. E. Chang, and J. I. Cirac, Phys. Rev. A **92**, 053834 (2015).

[46] T. Caneva, M. T. Manzoni, T. Shi, J. S. Douglas, J. I. Cirac, and D. E. Chang, New Journal of Physics **17**, 113001 (2015).

[47] A. Asenjo-Garcia, M. Moreno-Cardoner, A. Albrecht, H. J. Kimble, and D. E. Chang, Phys. Rev. X **7**, 031024 (2017).

[48] P. Meystre and M. Sargent, Elements of Quantum Optics (Springer Berlin Heidelberg, 2007).

[49] A. Lagendijk and B. A. van Tiggelen, Physics Reports **270**, 143 (1996).

[50] N. J. Schilder, C. Sauvan, Y. R. P. Sortais, A. Browaeys, and J.-J. Greffet, Phys. Rev. A **96**, 013825 (2017).

[51] G. Birkl, M. Gatzke, I. H. Deutsch, S. L. Rolston, and W. D. Phillips, Phys. Rev. Lett. **75**, 2823 (1995).

[52] R. J. Bettles, S. A. Gardiner, and C. S. Adams, Phys. Rev. A **92**, 063822 (2015).

[53] E. Shahmoon, D. S. Wild, M. D. Lukin, and S. F. Yelin, Phys. Rev. Lett. **118**, 113601 (2017).

[54] J. Rui, D. Wei, A. Rubio-Abadal, S. Hollerith, J. Zeiher, D. M. Stamper-Kurn, C. Gross, and I. Bloch, Nature **583**, 369 (2020).

[55] 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).

[56] R. Bekenstein, I. Pikovski, H. Pichler, E. Shahmoon, S. F. Yelin, and M. D. Lukin, Nature Physics **16**, 676 (2020).

[57] M. Moreno-Cardoner, D. Goncalves, and D. E. Chang, Phys. Rev. Lett. **127**, 263602 (2021).

[58] A. Glicenstein, G. Ferioli, N. Šibalić, L. Brossard, I. Ferrier-Barbut, and A. Browaeys, Phys. Rev. Lett. **124**, 253602 (2020).

[59] G. Ferioli, A. Glicenstein, L. Henriet, I. Ferrier-Barbut, and A. Browaeys, Phys. Rev. X **11**, 021031 (2021).

[60] S. Jennewein, L. Brossard, Y. R. P. Sortais, A. Browaeys, P. Cheinet, J. Robert, and P. Pillet, Phys. Rev. A **97**, 053816 (2018).

[61] 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).

[62] 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].

[63] F. Andreoli, M. J. Gullans, A. A. High, A. Browaeys, and D. E. Chang, Phys. Rev. X **11**, 011026 (2021).

[64] M. O. Scully and M. S. Zubairy, *Quantum Optics* (Cambridge University Press, 1997).

[65] K. J. Kusmierek, S. Mahmoodian, M. Cordier, J. Hney, A. Rauschenbeutel, M. Schemmer, P. Schneeweiss, J. Volz, and K. Hammerer, SciPost Phys. Core **6**, 041 (2023).

[66] M. K. Tey, G. Maslennikov, T. C. H. Liew, S. A. Aljunid, F. Huber, B. Chng, Z. Chen, V. Scarani, and C. Kurssiefer, New Journal of Physics **11**, 043011 (2009).

[67] D. Goncalves, M. W. Mitchell, and D. E. Chang, Phys. Rev. A **104**, 013724 (2021).

[68] R. Corless and D. Jeffrey, Artificial Intelligence - AI, 76 (2002).

[69] 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].

[70] G. Y. Kim and C. H. Kwak, Journal of Optics **22**, 095502 (2020).

[71] L. Allen and J. H. Eberly, *Optical resonance and two-level atoms* (Dover, 1987).

[72] 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.

[73] R. T. Sutherland and F. Robicheaux, Phys. Rev. A **95**, 033839 (2017).

[74] R. H. Dicke, Phys. Rev. **93**, 99 (1954).

[75] S. Mahmoodian, M. Čepulkovskis, S. Das, P. Lodahl, K. Hammerer, and A. S. Sørensen, Phys. Rev. Lett. **121**, 143601 (2018).

[76] D. Lechner, R. Pennetta, M. Blaha, P. Schneeweiss, A. Rauschenbeutel, and J. Volz, Phys. Rev. Lett. **131**, 103603 (2023).

[77] S. B. Jäger, H. Liu, J. Cooper, T. L. Nicholson, and M. J. Holland, Phys. Rev. A **104**, 033711 (2021).

[78] A. Shankar, J. T. Reilly, S. B. Jäger, and M. J. Holland, Phys. Rev. Lett. **127**, 073603 (2021).

[79] G. Ferioli, A. Glicenstein, F. Robicheaux, R. T. Sutherland, A. Browaeys, and I. Ferrier-Barbut, Phys. Rev. Lett. **127**, 243602 (2021).

[80] 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_{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_{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_{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
<figure><img src="image_11.png" /><figcaption>FIG. 10. (a) Atomic emission $\gamma_{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. [21]. We have rescaled the experimental data as indicated in the main text. Here, $s = 5$ and $D = 26$. (b,c) Peak height $\max[\gamma_{at}(t)]$ and damping rate $\tau_{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.</figcaption></figure>

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_{max} = \max[\gamma_{at}(t)]$ and damping rate $\tau_{at}$ as a function of the saturation parameter $s$. The damping rate is obtained by fitting the function $\gamma_{at}(t) = \gamma_{max} \exp(-t/\tau_{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}_{dd}$ in Eq. (16b) couples different angular momentum states. To do so, we calculate how $\hat{H}_{dd}$ acts on the total angular momentum basis $|j, m\rangle$. We note this is not trivial since $\hat{H}_{dd}$ contains only local operators that break the degeneracy of $|j, m\rangle$.

Since each term in $\hat{H}_{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 C} |0, 0; j, m\rangle + C_{j,m}^{1,0,j_{N-2},m} |1, 0; j_{N-2}, m_{N-2}\rangle + C_{j,m}^{1,1,j_{N-2},m_{N-2}} |1, 1; j_{N-2}, m_{N-2}\rangle + C_{j,m}^{1,-1,j_{N-2},m_{N-2}} |1, -1; j_{N-2}, m_{N-2}\rangle, \quad (B1)$$

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

$$C_{j,m}^{0,0,j',m'} = \delta_{j,j'} \delta_{m,m'}, \quad (B2a)$$

$$C_{j,m}^{1,0,j',m'} = C_{j,m}^{1,0,j',m'} \delta_{m,m'}, \quad (B2b)$$

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}_{dd,nl} \equiv \sigma_n^+ \sigma_l^- - \sigma_l^+ \sigma_n^-$, appearing in the coherent dipole-dipole interaction $\hat{H}_{dd} = -i\hbar\Gamma_{1D}/4 \sum_{n>l} \hat{H}_{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}_{dd}|N/2, m\rangle = i\hbar\Gamma_{1D}F(m)|\psi\rangle. \quad (B3)$$

Here, we have defined

$$F(m)|\psi\rangle = \frac{C_{N/2,m}^{1,0,(N/2-1),m}}{4} \sum_{l>m} |0, 0; (N/2 - 1), m\rangle_{n,l}, \quad (B4)$$

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

Using the Clebsch-Gordon coefficients defined above, one can show that $\hat{H}_{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)(N/2 - m)(N + 1)}{24}}. \quad (B5)$$

Equation (B5) provides the scaling of the coherent dipole-dipole matrix elements with system size. In the thermodynamic limit $N \to \infty$, the coupling scales linearly with
<figure><img src="image_12.png" /><figcaption>FIG. 11. Matrix elements of the coherent driving Hamiltonian $\hat{H}_{\text{drive}}$ (blue arrows), of the collective dissipation $\mathcal{L}_{\downarrow}$ (red arrows) and of the coherent dipole-dipole interactions $\hat{H}_{\text{dd}}$ (violet arrows) for Dicke states with $m=0$ and $m=-N/2+1$.</figcaption></figure>

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}_{\text{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 C) and the discrepancy between the Cavity and M-B models in Sec. IIID.

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