# Directional Superradiance in a Driven Ultracold Atomic Gas in Free-Space

Sanaa Agarwal, Edwin Chaparro, Diego Barberena, A. Piñeiro Orioli, G. Ferioli, S. Pancaldi, Igor Ferrier-Barbut, A. Browaeys, A.M. Rey

► To cite this version:

Sanaa Agarwal, Edwin Chaparro, Diego Barberena, A. Piñeiro Orioli, G. Ferioli, et al., Directional Superradiance in a Driven Ultracold Atomic Gas in Free-Space. PRX Quantum, 2024, 5 (4), pp.040335. 10.1103/PRXQuantum.5.040335. hal-05348473

HAL Id: hal-05348473

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

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
numerical capabilities [48] and many open questions remain. One particular exception is the limit when only permutationally symmetric states, also known as Dicke states [49] are populated. This situation arises naturally in optical cavities, where a single cavity mode talks to all atoms independent of their location in the array. In this limit, the theoretical treatment is significantly simpler and has been a focus of theoretical investigations for decades. One particular case is the so called Cooperative Resonance Fluorescence (CRF) [50–53] or collective atomic emission [54] that treats the behavior of a group of two-level atoms coupled identically to a single radiation mode. As first proposed by Dicke [49], it leads to a collective decay mechanism, superradiance, which, as was later theoretically shown, can be stabilized by driving the system below a critical drive strength. In the thermodynamic limit of a large number of atoms (or a large cooperation number), the competition between a coherent drive and decay gives rise to a non-equilibrium second-order phase transition. Below a threshold drive strength, the collective dipole reaches a highly pure steady-state characterized by a collective Bloch vector below the equator, pointing at a polar angle at which the superradiant decay and the external drive compensate each other. On the other hand, above a critical drive, collective dissipation is not enough to stabilize the strongly driven system and at the mean-field (MF) level, the system remains oscillating. Beyond-MF effects dampen the oscillations via phase diffusion and the system becomes a highly mixed state with a distribution centered around the equator.

A recent experiment [1] reported signatures of the above collective Dicke transition while interrogating a pencil-shaped cloud of N Rb atoms optically excited by a laser propagating along its main axis. The experiment observed clear manifestations of the two non-equilibrium phases depending on the ratio between the drive's single atom Rabi frequency, $Ω$, and the collective dissipation rate, $ΓN_{eff}$, characterized by the single particle decay rate, $Γ$, and an effective atom number, $N_{eff}$, which accounted for the finite extent of the cloud's diffraction mode, $N_{eff} ∝ N$. Above a critical drive $Ω > Ω_C$, a scaling of the photon emission rate consistent with $N_{eff}^2$ was observed. This scaling was modified below $Ω_C$. As the system crossed the critical point, the characteristics of the superradiant light changed as well.

Regardless of these clear signatures, as explained in that work [1], the applicability of a fully collective model scaled by an effective atom number to describe the light scattering of an elongated sample in free-space is highly unexpected. A justification of its validity starting from a microscopic model remains an open question. In this paper we solve this issue by performing a detailed study of the light scattering properties of pencil-shaped disordered arrays of two-level atoms by directly solving a master equation. Our analysis accounts for the spatial extent of the cloud and the dipole-dipole interactions across the array, including spatial fluctuations, elastic dipole-dipole interactions, and single particle decay. Starting from a

mean-field description of the master equation, expected to be valid in the weak and strong excitation limits, and complementing it with a Moving Average Cluster Expansion method (MACE-MF) and the cumulant method, we reproduce the experimental observations in all reported parameter regimes. Extended calculations over a larger atom number window reveal a modified scaling, which is not proportional to an effective atom number ($N_{eff}$) but to the square root of the total number of atoms ($\sqrt{N}$) in the array.

Our study demonstrates that the observed collective behavior is not valid over a broad parameter regime (see Fig. 1) since it misses important key features of free-space emission such as single atom spontaneous emission and the frequency shifts arising from inhomogeneous elastic dipolar interactions. By combining the latter mechanisms with collective decay plus drive, we are able to qualitatively reproduce all the features observed in the mean-field calculations that account for microscopic details. In the future, it will be interesting to push experiments to more dense regimes where the mean-field model becomes invalid and a genuine quantum treatment would be required in order to properly capture the role of photon-mediated dipolar interactions in driven ultra-cold atoms.

In the spirit of quantum simulation, our analysis takes advantage of state-of-the-art experimental capabilities in regimes challenging for theory, uses them to develop new theoretical insights that shed light on long standing problems, and makes predictions to further inspire experimental work.

In Sec. II, we introduce the dipolar model, which describes our system at the microscopic level. We further discuss the extremely dilute and dense limiting cases, which are analytically solvable and useful for understanding the role of interactions in different driving regimes. In Sec. III, we study the steady-state of the dipolar model using mean-field numerics in the weak and strong driving regimes, and we compare its emergent properties with the CRF and non-interacting models. Using these insights, we propose a phenomenological model in Sec. IV, which accounts for the inhomogeneity of dipolar interactions, is analytically tractable at the mean-field level, and qualitatively reproduces the dipolar model. In Sec. V, we compare our steady-state and dynamics results with the experimental data [1], and find fair agreement for most observables. Lastly, in Sec. VI, we discuss some concluding remarks and future directions.

## II. QUASI-1D GAS OF DIPOLAR-INTERACTING ATOMS

### A. Dipolar model

We consider a 3D elongated (pencil-shaped) gas of N point-like atoms fixed at their Gaussian-sampled positions ${\vec{r}_k}$, as shown in Fig. 2. The atoms have linewidth
<figure><img src="image_2.png" /><figcaption>FIG. 2. Schematic: A pencil-shaped ultracold gas of frozen two-level atoms interacting via photon-mediated interactions, with elastic (Re$G_{ij}$) and inelastic (Im$G_{ij}$) components. A continuous laser drive excites the atoms on-resonance with Rabi frequency $\Omega$, wavevector $\vec{k}_L = 2\pi/\lambda \hat{x}$, and polarization $\hat{e}_L = \hat{y}$ (perpendicular to the quantization axis, $\hat{z}$). Atoms spontaneously emit photons into free-space at rate $\Gamma$.</figcaption></figure>

(single atom decay rate) $\Gamma$ and are driven by a coherent laser drive with Rabi strength $\Omega = |\vec{E} \cdot \vec{d}|/\hbar$ ($\vec{E}$ is the electric field amplitude of the laser and $\vec{d}$ is the transition dipole moment of the atom), detuning $\Delta = \omega_L - \omega_0$ from the atomic transition (frequency $\omega_0$), wavevector $\vec{k}_L$, and polarisation $\hat{e}_L$. The quantization axis $\hat{e}_0$ is along the $\hat{z}$-axis and the axial direction of the gas is along the $\hat{x}$-axis. As in Ref. [1], we consider the case where $\hat{e}_L = \hat{e}_y = \hat{y}$, i.e., the driving laser is linearly polarized and excites the $\sigma_+$-transition between the hyperfine levels $|F=2, m_F=2\rangle \leftrightarrow |F=3, m_F=3\rangle$. The atoms interact via induced dipole-dipole interactions mediated by the vacuum electromagnetic modes [2].

We assume that the atomic cloud is cold enough that we can neglect motion of the atoms, treat them as frozen during the dynamics, and only focus on the internal state dynamics spanned by the two relevant levels in each atom ($|g\rangle \equiv |F=2, m_F=2\rangle$ ground and $|e\rangle \equiv |F=3, m_F=3\rangle$ excited), which define a spin-1/2 system. Under these conditions, which are similar to the system realized in Ref. [1], the state of an atom can be described by spin-1/2 Pauli operators: $\hat{\sigma}_k^z = |e_k\rangle\langle e_k| - |g_k\rangle\langle g_k|$, $\hat{\sigma}_k^+ = |e_k\rangle\langle g_k|$, $\hat{\sigma}_k^- = |g_k\rangle\langle e_k|$, for an atom $k$. The collective spin operators are denoted as $\hat{S}_\alpha = \sum_{k=1}^N \hat{\sigma}_k^\alpha / 2$ for $\alpha \in \{x, y, z\}$ and $\hat{S}_\pm = \sum_{k=1}^N \hat{\sigma}_k^\pm$.

The dynamics of the system is governed by the master equation obtained upon adiabatically eliminating the photonic degrees of freedom, $\hat{\rho} = -i[\hat{H}, \hat{\rho}] + \mathcal{L}(\hat{\rho})$, where $\hat{H}$ is the Hamiltonian, which leads to unitary evolution, and $\mathcal{L}(\hat{\rho})$ is the Lindbladian super-operator, which accounts for all the dissipative processes. The Hamiltonian for the system is given as $\hat{H} = \hat{H}_0 + \hat{H}_{\text{int}}$, where (setting $\hbar = 1$)

$$ \hat{H}_0 = -\frac{\Omega}{2} \sum_k (e^{i(\vec{k}_L \cdot \vec{r}_k - \pi/2)} \hat{\sigma}_k^+ + \text{h.c.}) - \Delta \hat{S}_z, \quad (1) $$

is the single-particle laser drive and

$$ \hat{H}_{\text{int}} = -\sum_{j \neq k} \mathcal{R}_{kj} \hat{\sigma}_k^+ \hat{\sigma}_j^-, \quad (2) $$

accounts for the dipole-dipole interaction. It sums over pairwise exchange processes among two different atoms in the array with the interaction strength given by the free-space electromagnetic Green's tensor $G(\mathbf{r}) = (3\Gamma/4) [(1 - \hat{r} \otimes \hat{r}) \frac{e^{ik_0 r}}{k_0 r} + (1 - 3\hat{r} \otimes \hat{r}) \left( \frac{ie^{ik_0 r}}{(k_0 r)^2} - \frac{e^{ik_0 r}}{(k_0 r)^3} \right)]$, where $\mathbf{r}$ is the vector connecting the two interacting atoms and $k_0 = 2\pi/\lambda$ ($\lambda$ is the atomic transition wavelength). The elastic part is determined by the real part of $G(\mathbf{r})$, $\mathcal{R}_{kj} = \hat{e}_+^* \cdot \text{Re}[G(\vec{r}_{kj})] \cdot \hat{e}_+$. Self-interactions are set to zero: $\mathcal{R}_{kk} = 0$. The Lindbladian for the system is expressed as

$$ \mathcal{L}(\hat{\rho}) = \sum_{j,k} \mathcal{I}_{kj} (2\hat{\sigma}_j^- \hat{\rho} \hat{\sigma}_k^+ - \{\hat{\sigma}_k^+ \hat{\sigma}_j^-, \hat{\rho}\}), \quad (3) $$

where $\mathcal{I}_{kj} = \hat{e}_+^* \cdot \text{Im}[G(\vec{r}_{kj})] \cdot \hat{e}_+$ is the inelastic dipolar interaction coefficient and $\mathcal{I}_{kk} = \Gamma/2$ is the spontaneous emission decay rate. We define the total dipolar interaction coefficient as $\mathcal{G}_{kj} = \mathcal{R}_{kj} + i\mathcal{I}_{kj}$.

For simplicity, we define a “spiral” basis, in which we absorb the drive phase in the coherences, as

$$ \hat{\sigma}_j^\pm = \hat{\sigma}_j^\pm e^{\pm i(\vec{k}_L \cdot \vec{r}_j - \pi/2)}. \quad (4) $$

Then, the spiral collective spin operators are $\hat{S}_\pm = \sum_{j=1}^N \hat{\sigma}_j^\pm e^{\pm i(\vec{k}_L \cdot \vec{r}_j - \pi/2)}$, $\hat{S}_x = (\hat{S}_+ + \hat{S}_-)/2$, $\hat{S}_y = -i(\hat{S}_+ - \hat{S}_-)/2$, and $\hat{S}_z = \hat{S}_z$. Accordingly, we also define new interaction coefficients $\tilde{\mathcal{G}}_{kj} = \mathcal{G}_{kj}e^{-i\vec{k}_L \cdot \vec{r}_{kj}}$, $\tilde{\mathcal{R}}_{kj} = \mathcal{R}_{kj}e^{-i\vec{k}_L \cdot \vec{r}_{kj}}$, and $\tilde{\mathcal{I}}_{kj} = \mathcal{I}_{kj}e^{-i\vec{k}_L \cdot \vec{r}_{kj}}$. In the spiral basis, the Hamiltonian and the Lindbladian are given as

$$ \hat{H}_0 = -\Omega \hat{S}_x - \Delta \hat{S}_z, \quad (5) $$

$$ \hat{H}_{\text{int}} = -\sum_{j \neq k} \tilde{\mathcal{R}}_{kj} \hat{\sigma}_k^+ \hat{\sigma}_j^-, \quad (6) $$

$$ \mathcal{L}(\hat{\rho}) = \sum_{j,k} \tilde{\mathcal{I}}_{kj} \left( 2 \hat{\sigma}_j^- \hat{\rho} \hat{\sigma}_k^+ - \{\hat{\sigma}_k^+ \hat{\sigma}_j^-, \hat{\rho}\} \right). \quad (7) $$

Thus, the drive acts collectively along the spiral-$\hat{x}$ direction in the Bloch sphere. Hereafter, we will work in the spiral basis and also set $\Delta = 0$, unless otherwise mentioned.

Similar to the experimental protocol in Ref. [1], we initialize all the atoms in the ground-state $|g\rangle^{\otimes N}$ and continuously drive them on resonance with the atomic transition, i.e., $\Delta = 0$. The excitation of atoms by the drive is counteracted by the free-space single particle and collective emission, which generates damping, allowing the system to eventually reach its steady-state. To characterize the system across a wide-range of drive strengths, $\Omega$, we look at the collective spin observables in the steady-state, namely, the atomic inversion $\langle \hat{S}_z \rangle$, the absolute
value of the spiral atomic coherence $|\langle\hat{S}_{+}\rangle|$, and its real part $|\langle\hat{S}_{x}\rangle|$.

The intensity operator describing the photon emission from the atomic sample along a direction $\vec{k}$ can be written in terms of the spin operators as [5, 55, 56],

$$ \hat{\mathcal{I}}(\vec{k}) \equiv I_0(\vec{k}) \sum_{i,j} \hat{\sigma}_i^+ \hat{\sigma}_j^- e^{i\vec{k} \cdot (\vec{r}_i - \vec{r}_j)}, \quad (8) $$

where $I_0(\vec{k})$ is a proportionality factor that accounts for the geometry of the dipolar emission pattern. In the spiral basis, the expectation value of the intensity can be expressed as $I(\vec{k}) \equiv \langle\hat{\mathcal{I}}(\vec{k})\rangle \equiv I_0(\vec{k}) \sum_{i,j} (\hat{\sigma}_i^+ \hat{\sigma}_j^-) e^{i(\vec{k} - \vec{k}_L) \cdot (\vec{r}_i - \vec{r}_j)}$. The intensity along the forward direction, $\vec{k} = \vec{k}_L$, features very interesting properties and will be the focus of this study. This direction is special because it is the direction along which the driving laser imposes coherence and therefore, the intensity can be enhanced due to constructive interference.

We also look into the equal-time two photon correlation function, defined as [56]

$$ \tilde{g}_2(\vec{k}, \vec{k}') = \frac{\langle:\hat{\mathcal{I}}(\vec{k})\hat{\mathcal{I}}(\vec{k}'):\rangle}{I(\vec{k})I(\vec{k}')}, \quad (9) $$

which gives the likelihood of simultaneously emitting a photon along $\vec{k}'$ and another photon along the $\vec{k}$-direction. The :: in the above expression implies it needs to be evaluated using normal-ordered operators. In the experiment [1] in question, the detector was placed in the forward direction to measure the above mentioned observables. Even though the experiment did not measure the atomic coherence, we will include it in our study as it is useful for gaining an intuitive understanding of the system. In addition to the steady-state, we also compute the dynamics of the forward intensity, $I(\vec{k}_L, t)$, and the excitation fraction,

$$ n_e(t) = \frac{1}{2N} \sum_{j=1}^{N} (\langle\hat{\sigma}_j^z(t)\rangle + 1). \quad (10) $$

In the following, we use these observables to investigate distinct key features of the system as a function of the strength of the drive, from the weak to the strong driving regimes. To understand the role of interactions, we consider the cases of the extremely dense and dilute limits, which are analytically tractable and well-studied. Then, we study the moderately dense ensemble of the experiment (Fig. 2) and compare it with the extremely dilute and dense limits as a means to investigate the behaviors emergent from dipolar interactions. Similar to the experiment, we keep the size of the atomic cloud fixed at an rms axial length $l_{ax} = 20\lambda$ and radial length $l_{rad} = \lambda/2$, as shown in Fig. 2. By varying $N$ we can tune the density of the cloud and thereby, the strength of the dipole-dipole interactions between the atoms.

## B. Cooperative Resonance Fluorescence (CRF)

In the limit in which all the atoms are confined within a single wavelength, $|r_{0,N}| < \lambda$, we have $\mathcal{R}_{kj} \sim 1/(|r_{kj}|/\lambda)^3$, $\mathcal{I}_{kj} \to \Gamma/2$. If we choose to neglect the elastic part of the interactions ($\mathcal{R}_{kj} = 0$), and freeze the motion of the atoms, we can emulate the situation found in an optical cavity, where all the atoms interact via a single electromagnetic mode. When the cavity mode is set on resonance with the atomic transition, the elastic interactions are suppressed and the dominant interactions are the ones responsible for collective decay or superradiance, via which atoms emit collectively at an enhanced emission rate $\Gamma N$. In the presence of an additional drive, $\Omega$, the system reduces to the well-studied Cooperative-Resonance-Fluorescence (CRF) model, which describes a system of $N$ atoms driven with a resonant laser drive at a Rabi frequency $\Omega$, and subject to collective decay (superradiance) described by the jump operator $\sqrt{\Gamma} \hat{S}_{-}$ with $\hat{S}_{\pm} = \sum_i \hat{\sigma}_i^{\pm}$ [57]. Here we use the notation CRF to distinguish from the so-called Dicke model [58–60] which contains both the rotating and counter-rotating terms. In the CRF model, the counter-rotating terms are irrelevant and therefore neglected. The master equation governing the dynamics of the CRF model is given by:

$$ \partial_t \hat{\rho} = -i\Omega[\hat{S}_x, \hat{\rho}] + \frac{\Gamma}{2}(2\hat{S}_-\hat{\rho}\hat{S}_+ - \{\hat{S}_+\hat{S}_-, \hat{\rho}\}) \quad (11) $$

where the drive, $\Omega$, is applied along the x-direction of the Bloch-sphere.

The CRF exhibits a steady-state phase transition at a critical frequency $\Omega_{CRF}^{C} = N\Gamma/2$ [61], which delineates two distinct steady-state behaviors depending on the dimensionless parameter $\beta = \Omega/\Omega_{CRF}^{C}$. For $\beta < 1$, the system is in the “superradiant phase”, where the drive is balanced by the collective emission, canceling the total electric field experienced by the collective dipole. We use quotation marks to remark that this regime was initially referred to instead as “superfluorescent” by the authors [50, 52, 61]. In this phase, the collective dipole lies in the $y-z$ plane of the Bloch sphere at an angle $\theta$ from the south pole, given by $\sin\theta = \beta$. For $\beta > 1$, the system transitions into a highly mixed steady-state, known as the “normal” phase. The statistical mixture arises from the collective emission not being strong enough to compensate the excitation from the drive. As the atomic self-radiated field is not canceled by the drive, the dipoles undergo collective Rabi flopping. Although Rabi oscillations persist at long times at the mean-field level, quantum fluctuations lead to phase diffusion, which dampens the oscillations and destroys the coherences, although the system always remains in the collective manifold by construction. These two distinct steady-state behaviors are separated by a second-order phase transition at $\beta = 1$, as shown in Fig. 3a.

Here it is important to note that Ref. [1] used the opposite convention and denoted the $\beta > 1$ regime as the superradiant phase. This is inspired by Refs. [62, 63] and
the fact that Eq. (11) has a discrete symmetry characterized by its invariance under a mirror reflection $\hat{S}_x \to -\hat{S}_x$ followed by complex conjugation. This mirror symmetry of the Lindblad generator is spontaneously broken in the $\beta > 1$ phase and therefore under this view, this is the phase featuring spontaneous symmetry breaking. For the purpose of this work, this is just a different convention and does not affect our conclusions. We will stick to the more standard CRF convention that uses $\hat{S}_z$ as the physically motivated order parameter and denotes the $\beta < 1$ regime as the superradiant phase.

By doing a Holstein-Primakoff expansion around the mean-field steady-state in the superradiant regime as well as introducing a semi-classical method to average in the normal phase [64], we can express the steady-state analytically across all driving regimes by varying $\beta$ (with $\eta = \sqrt{\beta^2 - 1}$) as (see Appendix C for derivation):

$$ \langle \hat{S}_z \rangle = \begin{cases} -\frac{N}{2} \sqrt{1 - \beta^2}, & \beta < 1 \\ 0, & \beta \ge 1 \end{cases} \qquad (12) $$

$$ \langle \hat{S}_x \rangle = 0, \qquad (13) $$

$$ \langle \hat{S}_y \rangle = \begin{cases} \frac{N}{2} \beta, & \beta < 1 \\ \frac{N}{2} \left( \beta - \frac{\eta}{\beta \arctan\left(\frac{1}{\eta}\right)} \right), & \beta \ge 1 \end{cases} \qquad (14) $$

$$ \langle \hat{S}_+ \hat{S}_- \rangle = \begin{cases} \frac{N^2}{4} \beta^2, & \beta < 1 \\ \frac{N^2}{4} \left( \beta^2 - \frac{\eta}{\arctan\left(\frac{1}{\eta}\right)} \right), & \beta \ge 1 \end{cases} \qquad (15) $$

$$ g_2(0) = \begin{cases} 1, & \beta < 1 \\ \eta \left( \frac{\beta^4}{\eta} \arctan\left(\frac{1}{\eta}\right) - \beta^2 - \frac{2}{3} \right) \\ \times \arctan\left(\frac{1}{\eta}\right) \left( \beta^2 \arctan\left(\frac{1}{\eta}\right) - \eta \right)^2, & \beta \ge 1 \end{cases} \qquad (16) $$

For this phase transition, the order parameter is the population inversion, $\langle \hat{S}_z \rangle$, which is non-zero in the superradiant phase and zero in the normal phase (Eq. (12)). $\langle \hat{S}_z \rangle$ is continuous but features an abrupt change in its derivative at the critical point in the thermodynamic limit ($N \to \infty$), as shown in Fig. 3a. The finite-N steady-state in Fig. 3 is obtained using exact diagonalization.

The observable $\langle \hat{S}_x \rangle$ is conserved and zero in all regimes (Eq. (13)). This is a distinguishing feature of this model, which will become important when comparing with the other models in this paper. Thus, the atomic coherence is purely imaginary, $\langle \hat{S}_+ \rangle = i \langle \hat{S}_y \rangle$. In the superradiant phase, $\langle \hat{S}_y \rangle$ grows linearly with $\beta$. On the other hand, as $\beta$ is increased in the normal phase, the steady-state becomes more and more mixed, and $\langle \hat{S}_y \rangle$ goes to zero, as shown in Fig. 3b.

The photon emission rate (intensity), $I = \langle \hat{S}_+ \hat{S}_- \rangle$, at a fixed $\beta$, increases as $N^2\beta^2$ in the superradiant phase. In the normal phase, the intensity still scales as $N^2$, reflecting the collective nature of the system, as $\langle \hat{S} \cdot \hat{S} \rangle = N/2(N/2 + 1)$, where $\hat{S} = \{\hat{S}_x, \hat{S}_y, \hat{S}_z\}$. However, the normal phase intensity is independent of $\beta$ due

<figure><img src="image_3.png" /><figcaption>FIG. 3. Steady-state of the CRF model using exact diagonalization (extremely dense limit, $r_{ij} \to 0$, and artificially setting $R_{ij} = 0$): (a) atomic inversion $\langle \hat{S}_z \rangle/(N/2)$ shows the superradiant phase transition in the thermodynamic limit, (b) total coherence $|\langle \hat{S}_+ \rangle|/N = |\langle \hat{S}_y \rangle|/N$ driven by the laser drive, which is along the $\hat{S}_x$-direction, such that $\langle \hat{S}_x \rangle = 0$, (c) intensity scales as $N^2$ and remains constant as $\Omega$ is varied in the strong-drive regime (due to preservation of $\langle \hat{S}^2 \rangle = \langle \hat{S}_+ \hat{S}_- \rangle + \langle \hat{S}_z^2 \rangle = N(N+1)/4$), (d) two-photon correlation function $g_2(0)$ does not reach the maximally mixed state ($g_2(0) = 2$) due to preserved quantum correlations. $N \to \infty$ curves correspond to Eqs. (12)-(16).</figcaption></figure>

to the mixed nature of the state in the normal phase. In fact, for $\beta \gg 1$, the intensity reaches the asymptotic value of $I = N^2/6$, as shown in Fig. 3c. The coherent nature of the steady-state in the superradiant phase is also evident in the value of the two-photon correlation function (Fig. 3d), as $g_2(0) = 1$ for $\beta < 1$, in the large-N limit. In the normal phase, on the contrary, the system enters a regime where $g_2(0) > 1$, suggesting the mixed nature of the steady-state and the build-up of classical correlations.

### C. Non-interacting model

In the extremely dilute limit, the mean inter-atomic distance is large, $|r_{ij}| \gg \lambda$, and thus to a good approximation we can neglect the interactions in Eq. (6) and Eq. (7) by setting $R_{kj} = 0$ and $I_{k \neq j} = 0$. This gives a simplified master equation for the non-interacting model, $\hat{\rho} = -i[\hat{H}_0, \hat{\rho}] + \mathcal{L}_0(\hat{\rho})$, that describes an array of independent atoms coherently interacting with a classical laser drive while spontaneously emitting photons at a rate $\Gamma$. The Lindbladian, $\mathcal{L}_0(\hat{\rho}) = (\Gamma/2) \sum_j (2\hat{\sigma}_j^- \hat{\rho} \hat{\sigma}_j^+ - \{\hat{\sigma}_j^z, \hat{\rho}\}/2 - \hat{\rho})$, captures the single-particle spontaneous emission. The laser drive acts as an effective global magnetic field, $\vec{B}_L = (\Omega/2, 0, \Delta/2)$. The corresponding dynamics of the Bloch vector, $\langle \hat{\vec{\sigma}} \rangle = (\langle \hat{\sigma}^x \rangle, \langle \hat{\sigma}^y \rangle, \langle \hat{\sigma}^z \rangle)$, can be described by the well-known Bloch equations as

$$ \frac{d\langle \hat{\vec{\sigma}} \rangle}{dt} = 2\langle \hat{\vec{\sigma}} \rangle \times \vec{B}_L - \vec{f}(\langle \hat{\vec{\sigma}} \rangle), \qquad (17) $$
where $\vec{f}(\langle\hat{\vec{\sigma}}\rangle) = -(\Gamma/2)[\langle\hat{\vec{\sigma}}\rangle + (\langle\hat{\sigma}^z\rangle + 2)\hat{z}]$ accounts for the damping.

For the on-resonant case ($\Delta = 0$), the steady-state Bloch vector is obtained as

$$ \langle\hat{\sigma}^+\rangle = i\frac{\Omega}{\Gamma}\langle\hat{\sigma}^z\rangle, \quad \langle\hat{\sigma}^z\rangle = -R \equiv -\frac{1}{1+2(\Omega/\Gamma)^2}, \quad (18) $$

and the steady-state intensity is

$$ I(\vec{k}) = I_0(\vec{k}) \frac{(\Omega/\Gamma)^2}{1 + 2(\Omega/\Gamma)^2} \left[ N + \sum_{j, i \neq j} \frac{e^{i(\vec{k} - \vec{k}_L) \cdot (\vec{r}_i - \vec{r}_j)}}{1 + 2(\Omega/\Gamma)^2} \right], \quad (19) $$

which is dominated by the coherences in the weak-driving regime ($\Omega \ll \Gamma$) and by the incoherent single-particle-like term in the strong-driving regime ($\Omega \gg \Gamma$). In the large-$N$ limit, the steady-state intensity can be obtained analytically by converting the sum above to an integral (with $\vec{k} = (\cos\theta, \sin\theta \cos\phi, \sin\theta \sin\phi)$), as

$$ I(\hat{k}) \propto \frac{2 \left(\frac{\Omega}{\Gamma}\right)^4 N}{\left[1 + 2 \left(\frac{\Omega}{\Gamma}\right)^2\right]^2} \left[1 + \frac{N e^{-(2\pi)^2 (\sigma_{ax}^2 (1-\cos\theta)^2 + \sigma_{rad}^2 \sin^2\theta)}}{2(\Omega/\Gamma)^2}\right] \quad (20) $$

where $\sigma_{ax}$ and $\sigma_{rad}$ are the axial and radial extents of the cloud, respectively. As discussed earlier, when the laser-wavevector aligns with the direction of observation ($\vec{k} = \vec{k}_L \Rightarrow \theta = 0$), the imprinted phases interfere constructively. In the weak-drive limit, the forward emission is dominated by its coherent part, $I(\vec{k}_L) \approx I_0(\vec{k}_L)|\langle\hat{\sigma}^+\rangle|^2 N(N-1)$. This $N^2$-enhancement of the emission rate is only seen in the forward direction and vanishes exponentially fast with increasing $\theta$ when $\theta \ll 1$ (Eq. (A19)). Along other $\vec{k}$, the phases in Eq. (19) do not cancel and the intensity, dominated by its incoherent part, scales as $N$ far away from $\vec{k}_L$. It is worth noting that this $N^2$-enhancement is not due to quantum correlations and arises purely from the coherences of the atoms.

For comparison with the other models in this work, we define the notion of a “critical” drive strength, $\Omega_{\mathrm{C}}^{\mathrm{Non-int.}}$, that maximizes the “forward” intensity. There is, however, no phase transition in the non-interacting model. This definition of the critical drive is inspired by the CRF model, in which the total coherence and the total intensity peak at the critical drive strength and a phase transition is observed (see Sec. II B for details). For the non-interacting model, we set $d\langle\hat{\sigma}^+\rangle/d\Omega = 0$ to obtain the critical drive strength as $\Omega_{\mathrm{C}}^{\mathrm{Non-int.}} = \Gamma/\sqrt{2}$.

As shown in Fig. 4a, in the weak-driving limit ($\Omega \ll \Omega_{\mathrm{C}}^{\mathrm{Non-int.}}$), the forward intensity scales as $I(\vec{k}_L) \propto (N\Omega/\Gamma)^2$, as shown in Fig. 4c. Moreover, we find that $\tilde{g}_2(0) = 1$ (Eq. (A11)) in this limit, as shown in Fig. 4d, which reflects that the system is in a coherent state with the majority of the atoms in the ground state.

In the strong-driving limit ($\Omega \gg \Omega_{\mathrm{C}}^{\mathrm{Non-int.}}$), radiative decay disrupts the coherences and leads to a completely

<figure><img src="image_4.png" /><figcaption>FIG. 4. Steady-state in the non-interacting (dilute) limit across driving regimes: (a) atomic inversion $\langle\hat{S}_z\rangle/(N/2)$, (b) total coherence $|\langle\hat{S}_+\rangle|/N = |\langle\hat{S}_y\rangle|/N$ driven by the laser in the $\hat{S}_x$-direction such that $\langle\hat{S}_x\rangle = 0$, (c) “forward” intensity (along the axial direction) scales as $N^2$ in the weak-drive regime and as (e) $N$ in the strong-drive regime $(\langle\hat{S}_z\rangle/(N/2) \to 0)$, (d) two-photon correlation function $\tilde{g}_2(0)$ in the “forward” direction. (e), (f): Due to the constructive interference of coherences in the “forward” direction, an $N$-dependent stronger drive ($\Omega \gtrsim \Gamma\sqrt{N}$) is needed to reach the thermal state $(\langle\hat{S}_+\hat{S}_-\rangle \to N/2, \tilde{g}_2(0) \to 2)$.</figcaption></figure>

mixed steady-state ($R \to 0$ in Eq. (18)). In this mixed state, the Bloch vector of each spin-1/2 (atomic dipole) is reduced to a point at the center of the Bloch sphere such that $\langle\hat{\sigma}^z\rangle \to 0$ (Fig. 4a) and $|\langle\hat{\sigma}^+\rangle| \to 0$ (Fig. 4b). In the absence of a finite coherence, there is no constructive interference in the forward direction, i.e., $I(\vec{k}_L) \propto N/2$ (Fig. 4e). We obtain $\tilde{g}_2(0) = 2$ in this limit, as shown in Fig. 4d, since the system is in a thermal state. As shown in Fig. 4f, a very strong drive, $\Omega \gtrsim \Gamma\sqrt{N}$, is needed to fully reach this condition. This $N$-dependent drive scaling emerges due to the $N^2$ enhancement from the constructive interference of coherences, and is further explained in Appendix A.

Lastly, the on-resonant drive Hamiltonian commutes with $\hat{\sigma}^x$, so starting from $\langle\hat{\sigma}^x\rangle = 0$ in the initial state, it remains zero at all times, irrespective of $\Omega$, as shown in Fig. 4b.

### III. MEAN-FIELD DIPOLAR MODEL

#### A. General description

For a general distribution of atomic positions, the interacting system is free to explore an exponentially large Hilbert space ($2^N$) and its exact dynamics is not tractable. For simplicity, we take the mean-field (MF)
approximation here to characterize the steady-state of the system. In Appendix G, we compare MF with beyond-MF approximation methods to include the effects of correlations, namely the MACE-MF (see Appendix F 2) and cumulant approximation methods. We show that at the atomic densities considered in this paper, the MF approximation is valid when describing the observables under consideration.

We obtain the MF equations of motion (see Appendix B) and express them as the Bloch equation, same as Eq. (17), but here the components of the effective magnetic field, $\vec{B}_i^{\text{MF}}$ at atom i, are site-dependent due to the inclusion of the dipolar interactions as

$$B_{i,x}^{\text{MF}} = \frac{\Omega}{2} + \frac{1}{4} \sum_{j \neq i} \left( \text{Re} \tilde{\mathcal{G}}_{ij} \langle \hat{\sigma}_j^x \rangle + \text{Im} \tilde{\mathcal{G}}_{ij} \langle \hat{\sigma}_j^y \rangle \right), \quad (21)$$

$$B_{i,y}^{\text{MF}} = \frac{1}{4} \sum_{j \neq i} \left( \text{Re} \tilde{\mathcal{G}}_{ij} \langle \hat{\sigma}_j^y \rangle - \text{Im} \tilde{\mathcal{G}}_{ij} \langle \hat{\sigma}_j^x \rangle \right), \quad (22)$$

$$B_{i,z}^{\text{MF}} = \frac{\Delta}{2}. \quad (23)$$

**Short-time physics:** Although the MF treatment above significantly reduces the complexity of the problem, the non-linearity of these equations makes them difficult to solve analytically and apart from simple limiting cases, we need to solve them numerically. The early time dynamics ($\Gamma t \ll 1$) is one such regime where it is possible to gain analytical insights.

To understand the role of interactions at short-times, we include the inelastic interaction term in the Hamiltonian $\hat{H}$ (Eqs. (5), (6)) and obtain the non-Hermitian Hamiltonian in the spiral basis as $\hat{H}_{\text{NH}} = \hat{H} - i \sum_{k,j \neq k} \tilde{\mathcal{G}}_{kj} \hat{\sigma}_k^+ \hat{\sigma}_j^-$, where $\hat{H}_0$ describes the laser drive. Given that the system is initialized in a state where all the Bloch vectors are identical in the tilde basis, i.e., $\langle \hat{\tilde{\sigma}}_k \rangle \equiv \langle \hat{\tilde{\sigma}} \rangle$, at short times, each Bloch vector commutes with terms of the form $(\hat{\tilde{\sigma}}_k \cdot \hat{\tilde{\sigma}}_j)$, since $\sum_{k,j \neq k} \langle [\hat{\tilde{\sigma}}_k \cdot \hat{\tilde{\sigma}}_j, \hat{\tilde{\sigma}}_i] \rangle \sim \sum_{j \neq i} \langle \hat{\tilde{\sigma}}_j \rangle \times \langle \hat{\tilde{\sigma}}_i \rangle \sim 0$.

The emergence of a density shift can be elucidated by adding such terms, $(\hat{\tilde{\sigma}}_k \cdot \hat{\tilde{\sigma}}_j)$, to $\hat{H}_{\text{NH}}$ without altering the physics of the original system at short times as

$$\hat{H} = \hat{H}_{\text{NH}} + \frac{1}{4} \sum_{k,j \neq k} \text{Re} \tilde{\mathcal{G}}_{kj} (\hat{\tilde{\sigma}}_k \cdot \hat{\tilde{\sigma}}_j), \quad (24)$$

which allows us to rewrite the real part of the new Hamiltonian as $\text{Re}\hat{\mathcal{H}} = \hat{H}_0 + \frac{1}{4} \sum_{k,j \neq k} \text{Re}\tilde{\mathcal{G}}_{kj} \hat{\sigma}_k^z \hat{\sigma}_j^z$, and at the mean field as, $\text{Re}\hat{\mathcal{H}}^{\text{MF}} = - \sum_i \vec{B}_i^{\text{MF}} \cdot \hat{\tilde{\sigma}}_i$, where $\vec{B}_i^{\text{MF}} = \vec{B}_L - (\delta_i/2)\hat{z} = (\Omega/2, 0, \Delta/2 - \delta_i/2)$, and $\delta_i = \sum_{j \neq i} \text{Re}\tilde{\mathcal{G}}_{ij} \langle \sigma_j^z \rangle / 2$ is the interacting part of $\vec{B}_i^{\text{MF}}$. This acts as a self-adjusting magnetic field along $\hat{z}$ that depends on the atomic inversion, $\langle \sigma_j^z \rangle$. This effective field generated by other atoms in the array, leads to a precession of each atomic dipole i about the $\hat{z}$-axis of the Bloch sphere, generating what is known as a density shift or

a collective Lamb-shift [65]. As the number of atoms in the cloud is increased, with the spatial extent fixed, the cloud gets denser and the induced frequency shift at each atom gets larger. The frequency shift has two components – a homogeneous component, $\Delta_{\text{MF}}$, which is the average shift across atoms, and an inhomogeneous component, associated with the random distribution of the atoms in the array. The homogeneous component is non-zero when the cloud is not spherically symmetric, and it pushes the atoms out of resonance, suppressing the growth of coherences and the excitation fraction. The inhomogeneous component leads to dephasing. These two key features will be an important consideration when we derive a simplified model of the dipolar Hamiltonian in Sec. IV.

The average frequency shift can be measured in Rabi spectroscopy by scanning the detuning for the new resonance condition, i.e., $\Delta_{\text{MF}} = \sum_i \delta_i / N = \sum_{i,j \neq i} \text{Re}\tilde{\mathcal{G}}_{ij} \langle \sigma_j^z \rangle / (2N)$, where the RHS has been averaged over all atoms i [65]. The density shift can also be measured via Ramsey spectroscopy [66]. In this case, even at $\Delta = 0$, there is a residual precession which leads to a non-zero $\langle \hat{S}_x \rangle$ during the dynamics of the system even when the system is initialized with $\langle \hat{S}_x \rangle = 0$. This feature distinguishes the dipolar model from the CRF and non-interacting models, where $\langle \hat{S}_x \rangle = 0$ at all times.

**Weak-drive limit:** Another simplified case is when the laser drive is very weak, $\Omega \ll \Gamma$, such that individual dipoles remain close to the south-pole of the Bloch sphere, i.e., $\langle \hat{\sigma}_j^z \rangle \approx -1$, even in the steady-state. In this regime, correlations are suppressed by factors of $\Omega/\Gamma$ and the system can be described almost exactly using mean-field theory by setting $\langle \hat{\sigma}_j^z \rangle = -1$ at all times and only considering the dynamics of the coherences $\langle \hat{\sigma}_j^\dagger \rangle$. This regime has been intensively studied for dilute samples of dipolar-interacting atomic gases [13, 18, 19, 21, 22, 25, 29–32, 37].

In the weak-drive limit ($\langle \hat{\sigma}_j^z \rangle = -1$) of our pencil-shaped Gaussian cloud, the average frequency shift is non-zero and as discussed above, can be obtained at short times as $\Delta_{\text{MF}} = - \sum_{j \neq i}^{N} \text{Re}\tilde{\mathcal{G}}_{ij}/(2N)$. In the far-field limit ($|r_{ij}| \gg \lambda \Rightarrow \tilde{\mathcal{G}}_{ij} \propto 1/|r_{ij}|$), we find that $\sum_{j \neq i}^{N} \text{Re}\tilde{\mathcal{G}}_{ij}/(2N) \propto N$, with a small ($\ll \Gamma$) proportionality constant that depends on the cloud extent and can be estimated numerically (see Appendix D for details).

In the weak-drive limit, the mean-field steady-state coherence, can also be obtained for a dilute gas by treating the interactions perturbatively as [18, 31]

$$\langle \hat{\sigma}_j^- \rangle = \frac{i\Omega/2}{i\Delta + \Gamma/2} \left[ 1 - \frac{i}{i\Delta + \Gamma/2} \sum_{k \neq j} \tilde{\mathcal{G}}_{jk} \right], \quad (25)$$
and the forward intensity is obtained as

$$ \frac{I(\vec{k}_{\text{L}})}{I_0(\vec{k}_{\text{L}})} = \frac{N(N-1)\Omega^2/4}{\Delta^2 + \Gamma^2/4} \left[ 1 - \frac{\sum_{j,k \neq j} (\Gamma \text{Im}\tilde{G}_{jk} - 2\Delta \text{Re}\tilde{G}_{jk})}{(\Delta^2 + \Gamma^2/4)N} \right] \quad (26) $$

The steady-state density shift and the linewidth broadening can be obtained from the lineshape as $\Delta = -\sum_{j,k \neq j} \text{Re}\tilde{G}_{jk}/N$ and $\bar{\Gamma} = 2\sum_{j,k \neq j} \text{Im}\tilde{G}_{jk}/N$, respectively [18, 31]. On resonance ($\Delta = 0$), the x-component of the coherence is proportional to the density shift, $\langle \hat{S}_x \rangle/N = 2\Omega\bar{\Delta}/\Gamma^2$ and the forward intensity is reduced by the inhomogeneous broadening as

$$ I(\vec{k}_{\text{L}}) \propto N(N-1)(\Omega/\Gamma)^2(1-2\bar{\Gamma}/\Gamma). \quad (27) $$

However, in a highly dense gas, multi-photon scattering processes become significant, making the steady-state more complicated and beyond the perturbative limit [31].

## B. Dynamical Regimes in the Steady State

Outside the special limits discussed above, we resort to numerical solutions for obtaining the steady-state of the system. We find that our system exhibits different dynamical behaviors characterized by distinct light emission properties, as we vary the laser driving strength, $\Omega$. Here, we draw parallels and distinctions among the steady-state properties of the dipolar, the non-interacting, and the CRF models, across different driving regimes.

### 1. Inversion, forward intensity, and two-photon correlation function

In Fig. 5a, we show the atomic inversion of the dipolar model across a range of $\Omega$. Unlike the CRF model, the dipolar model does not exhibit a sharp phase transition and the system crosses over smoothly from the weak to the strong driving regime. Nevertheless, we find it useful to define a “critical” driving strength, $\Omega_{\text{C}}^{\text{Dipolar}}$, for the dipolar model, as the drive strength at which the forward intensity, $\langle \hat{S}_{+} \hat{S}_{-} \rangle$, peaks for a given $N$. The intensity peak separates the dynamical behaviors in the different driving regimes based on their emission properties, similar to the CRF model. We obtain $\Omega_{\text{C}}^{\text{Dipolar}} = c\Gamma\sqrt{N}$ ($c \approx 0.08$) from the numerical steady-state of the system. By plotting the steady-state atomic inversion as a function of a normalized drive strength, $\Omega/\Omega_{\text{C}}^{\text{Dipolar}}$, we observe that the curves for different $N$ appear to collapse, as shown in Fig. 5a.

As $N$ decreases, the cloud becomes dilute and the dipolar interaction coefficients become less important. At very small $N$, the steady-state of the non-interacting

<figure><img src="image_5.png" /><figcaption>FIG. 5. Steady-state of the MF Dipolar model: (a) atomic inversion $\langle \hat{S}_z \rangle/(N/2)$, (b) contrast $|\langle \hat{S}_+ \rangle|/N$, (c) “forward” intensity (along the axial direction), (d) two-photon correlation function $\tilde{g}_2(0)$ in the “forward” direction. (a)-(d) The dipolar model (solid) converges to the non-interacting model (black dashed) in the dilute limit (small $N$). The x-axis of the non-interacting curve is scaled as $\Omega/\Omega_{\text{C}}^{\text{Non-int.}} \sim \mathcal{O}(1)$ in (a)-(c) and as $\Omega/\Omega_{\text{C}}^{\text{Dipolar}} \propto 1/\sqrt{N}$ in (d). In (a), (b), and (d), the non-interacting curve is independent of $N$; in (c), the non-interacting curve is for $N = 100$.</figcaption></figure>

model is smoothly recovered for all the observables, as shown in Fig. 5.

In Fig. 5c of the dipolar model, we show that the variation of the forward intensity $\propto \langle \hat{S}_+ \hat{S}_- \rangle$ curve with the normalized drive strength is very similar to that of the non-interacting model (Fig. 4c). The main difference between these two models is evident in the $N$-scaling of $\langle \hat{S}_+ \hat{S}_- \rangle$. In the weak-driving regime of the CRF and non-interacting models, where phase matching leads to a perfect $N^2$ enhancement, we show that the forward intensity curves for different $N$ collapse when divided by $N^2$, in Fig. 3c and Fig. 4c, respectively. But for the dipolar model, we find that the curves do not collapse when divided by $N^2$. At weak drive intensities, the site-dependent shifts due to the dipolar interactions imprint random phases on the coherences and thereby suppress the $N^2$ enhancement of $\langle \hat{S}_+ \hat{S}_- \rangle$ [31, 37]. Hence, the $N$-scaling of $\langle \hat{S}_+ \hat{S}_- \rangle$ decreases as $N$ increases, as later shown in Fig. 8b.

To further understand the emission properties of the steady-state, we look at the forward two-photon correlation function $\tilde{g}_2(0) = \frac{\langle \hat{S}_+ \hat{S}_+ \hat{S}_- \hat{S}_- \rangle}{\langle \hat{S}_+ \hat{S}_- \rangle^2}$. As shown in Fig. 5d, in the weak driving regime ($\Omega \ll \Omega_{\text{C}}^{\text{Dipolar}}$), we have $\tilde{g}_2(0) \approx 1$, which corresponds to a coherence state. This is expected because the steady-state of the system resembles a coherent state on the Bloch sphere with a finite coherence $|\langle \hat{S}_+ \rangle|/N$, as seen in Fig. 5b and later validated analytically using a simpler model (Sec. IV). Given that the CRF, non-interacting, and dipolar models all feature
almost pure Gaussian-like steady-states in this regime, $\tilde{g}_2(0)$ is not able to distinguish their subtle differences and all the three models have $\tilde{g}_2(0) \sim 1$ in the weak excitation limit (see Figs. 3d, 4d, and 5d).

When the drive strength is above the “critical” drive strength, i.e., $\Omega \gg \Omega_C^{\text{Dipolar}}$, the drive quickly dephases the array. As a consequence, dipolar interactions become subdominant compared to the rapid rotation induced by the drive. Consistently, the steady-state starts to become fully mixed, with suppressed coherences, signaled by the fact that $\tilde{g}_2(0) \sim 2$ approaches its thermal value (see Fig. 5d).

This is in strike contrast to the CRF case, where the collective nature of the master equation enforces the preservation of the collective nature of the state even in the large driving limit ($\Omega \gg \Omega_C^{\text{CRF}}$). Therefore, the intensity keeps scaling as $N^2$ (Eq. (15)), as shown in Fig. 3c and $\tilde{g}_2(0)$ remains always below the thermal value and saturates as $\tilde{g}_2(0) \to 1.2$ (see Fig. 3d). These key distinctions imply that a free-space atomic cloud does not behave as a CRF model in the strong drive regime.

## 2. Steady-state frequency shift and contrast

In Fig. 6a, we see that the steady-state of the dipolar model has $\langle \hat{S}_x \rangle \neq 0$, unlike the non-interacting (Fig. 4b) and CRF (Fig. 3b) models where the *x*-component of the collective coherence is zero. The resonant laser drive commutes with $\langle \hat{S}_x \rangle$ and the finite $\langle \hat{S}_x \rangle$ arises from the dispersive (elastic) part of the dipolar interactions, which leads to frequency shifts, as discussed in Sec. III A. The latter is most prominent in the weak-drive limit, where the inversion remains close to its initial minimum value. In this limit, the fractional *x*-component of the coherence, $|\langle \hat{S}_x \rangle|/|\langle \hat{S}_+ \rangle|$, increases with *N*, which suggests that denser clouds have larger collective shifts.

Beyond the critical drive, as shown in Fig. 5b, the contrast, $|\langle \hat{S}_+ \rangle| = \sqrt{\langle \hat{S}_x \rangle^2 + \langle \hat{S}_y \rangle^2}$, starts decaying as the state begins to get mixed and the Bloch vector length gets smaller. This behavior continues with increasing $\Omega/\Gamma$ until eventually the system’s Bloch vector reduces to the one expected for the non-interacting system. The approach to the non-interacting regime as the system reaches the strong drive regime is also signaled by the dominant $\langle \hat{S}_{y,z} \rangle$ spin projections, and the reduced $\langle \hat{S}_x \rangle$ one. This in turn leads to a reduced value of $|\langle \hat{S}_x \rangle|/|\langle \hat{S}_+ \rangle|$, as shown in Fig. 6a.

In Fig. 6b, we look at the lineshape of the steady-state intensity in the weak-driving regime ($\Omega/\Gamma = 0.1$) by varying the laser detuning, $\Delta$, of the MF dipolar model. We find that increasing *N* leads to a bimodal distribution with a peak at $\Delta = 0$ and a second resonance at $\Delta_{\text{peak}}$. Note the second peak arises from dipolar interactions, and scales linearly with *N*, as shown in Fig. 6c.

The linear *N*-scaling is consistent with previous work

<figure><img src="image_6.png" /><figcaption>FIG. 6. (a) The *x*-component of the coherence, $|\langle \hat{S}_x \rangle|/|\langle \hat{S}_+ \rangle| \neq 0$, arising from frequency shifts due to the elastic dipolar interactions. For the non-interacting model, $\langle \hat{S}_x \rangle = 0$ (black dashed). (b) Steady-state “forward” intensity for a range of laser detunings, $\Delta$: MF dipolar model numerics (dots), interpolated data (solid lines). (c) Steady-state frequency shift computed as the detuning at which the maximum intensity is reached (circles) which follows a linear dependence with *N* as can see from the linear best fit (solid line). Data from numerics of the MF dipolar model in the weak-drive regime ($\Omega/\Gamma = 0.1$), averaged over 10 realizations.</figcaption></figure>

done in the weak drive limit [18, 31, 37, 67, 68] that observed a frequency shift proportional to the atomic density. In our quasi-1D configuration along the $\hat{x}$-axis, the atomic density scales as, $\rho \sim N/(a_{ho}^2 L)$, where $L = 2l_{ax}$ is the axial length of the pencil-shaped cloud and $a_{ho} \ll x_{ij}$ sets the radial confinement. The existence of a prominent second peak can be naively understood from the fact that the inter-atomic distances are determined by the axial spacing, i.e., $k_0|\vec{r}_{ij}| = k_0\sqrt{x_{ij}^2 + y_{ij}^2 + z_{ij}^2} \approx k_0|x_{ij}|$. As such, the pairwise interactions between atoms along the $\hat{x}$-axis contribute constructively to the shift. Mathematically, this is because the phases imprinted by the laser cancel the ones of the Green’s function $\tilde{G}_{i>j} \propto \exp(ik_0|\vec{r}_{ij}| - i\vec{k}_L \cdot \vec{r}_{ij}) \sim 1$ (see Appendix B for more details). Hence, by aligning the laser wavevector along the elongated geometry of the cloud, one can induce a non-zero global frequency shift even in a disordered configuration [24].

As discussed earlier in Sec. III A, the steady-state density shift of the lineshape in the weak-drive limit for a dilute gas is $\Delta = -\sum_{k \neq j} \text{Re}\tilde{G}_{jk}/N$. For simplicity, we consider the far-field regime ($\propto 1/r$) of the interactions, which is the dominant term in the dilute case, and obtain $\tilde{\Delta}$ in the large-*N* limit by integrating over the pencil-shaped distribution of the cloud, as $\tilde{\Delta} = -c_R N\Gamma$, where $c_R \approx 0.003$ is a constant depending on the spatial extent of the cloud. Similarly, we obtain the linewidth broadening as $\tilde{\Gamma} = c_I N\Gamma$ with $c_I \approx 0.004$. In the regime where $2\tilde{\Gamma} > \Gamma$, i.e., $N > 250$, the dilute gas description is no longer valid on resonance as, under this assumption, the predicted forward intensity becomes unphysical (Eq. (27)), and the perturbative treatment fails. Thus, even though the prediction of $\tilde{\Delta}$ is very close to $\Delta_{\text{peak}}$, the weak-drive physics in our system is beyond the first-
order expansion in interactions.

Of course, our analysis neglects motional effects in
thermal samples, laser forces, and dipolar forces. As
shown in Ref. [31], motional decoherence can induce a
reduction of the density shift and wash out the double
peak structure compared to what is expected for frozen
atoms.

C. N-Scaling for Different Driving Strengths

<figure><img src="image_7.png" /><figcaption>FIG. 7. (a) Steady-state atomic inversion $\langle\hat{S}_z\rangle/(N/2)$ for the dipolar (solid), non-interacting (dotted), and CRF (dashed dotted) models. The $x$-axis is scaled with the critical drive $\Omega_C^m$ for each model, $m$. (b) $N$-scaling of $\Omega_C^m$ for the dipolar ($\sim\sqrt{N}$), non-interacting ($\sim O(1)$), and CRF ($\sim N$) models.</figcaption></figure>

To emphasize the collapse of different *N* curves and the
different critical properties across models, we compare
the steady-state atomic inversion across a wide range of
driving strengths for the dipolar model (solid), the CRF
model (dashed-dotted), and the non-interacting model
(dotted), in Fig. 7a. The *x*-axis of Fig. 7a has been
scaled differently for each model, corresponding to the
scaling of the “critical” driving strength $\Omega_C^m$, where $m \in$
{CRF model, Dipolar model, Non-int. (Non-interacting)
model}. The *N*-scaling of $\Omega_C^m$ for these models has been
shown in Fig. 7b. The *N*-scaling of the critical drive for
the dipolar model ($\propto\sqrt{N}$) clearly differs from that of
the non-interacting ($\sim O(1)$) and CRF ($\propto N$) models.

In Fig. 8, we show the *N*-scaling of the steady-state in-
tensity, i.e., we define $\alpha$ such that $I(\vec{k}_L) \propto \langle \hat{S}_+ \hat{S}_- \rangle \propto N^\alpha$
at fixed $\Omega/\Gamma$, across a wide-range of driving strengths.
In Fig. 8a, we choose the *x*-axis to be the bare driving
strength $\Omega/\Gamma$, to compare our predictions with the mea-
surements of a recent experiment [1]. In the weak-driving
regime, we find $\alpha < 1$ for the dipolar and the CRF mod-
els. Usually, $\alpha < 1$ is understood as a signature of sub-
radiance but in our case this is simply an artifact due
to the *N*-scaling being calculated at a fixed bare driv-
ing strength, $\Omega/\Gamma$. At a fixed value of $\Omega/\Gamma$, the system
crosses over from the strong ($\Omega \gg \Omega_C^{\text{Dipolar}}$, $\Omega \gg \Omega_C^{\text{CRF}}$)
to the weak ($\Omega \ll \Omega_C^{\text{Dipolar}}$, $\Omega \ll \Omega_C^{\text{CRF}}$) driving regime
as *N* is increased, leading to an inaccurate *N*-scaling.
Similarly, we see a diverging value of $\alpha$, i.e., $\alpha > 2$,
for the mean-field dipolar model at intermediate driving

<figure><img src="image_8.png" /><figcaption>FIG. 8. *N*-scaling exponent of steady-state “forward” intensity, $\langle \hat{S}_+ \hat{S}_- \rangle \propto N^\alpha$: (a) $\alpha$ computed at a fixed Rabi frequency, $\Omega/\Gamma$, similar to Ref. [1]. $\alpha$ is similar for the CRF ($N=10$) and dipolar model ($N=2000$) in the intermediate driving regime for a small range of $\Omega/\Gamma$, which was probed in Ref. [1]. At large $\Omega/\Gamma$, $\alpha=2$ for the CRF due to the preservation of $\langle \hat{S}^2 \rangle$, but the dipolar model becomes single-particle-like ($\alpha=1$) due to spontaneous emission. Due to finite-size effects, $\alpha < 2$ for $N=10$ (CRF). (b) $\alpha$ computed at a fixed $\Omega/\Omega_C^m$, the ratio of the Rabi frequency and the critical drive, for each model. MF dipolar (solid), non-interacting (black dashed), and CRF (green dashed-dotted).</figcaption></figure>

strengths in Fig. 8a, which is again an artifact of fixing
$\Omega/\Gamma$. These artifacts are not seen for the non-interacting
model because $\Omega_C^{\text{Non-int.}} \sim O(1)$ (in the large-$N$ limit),
so fixing $\Omega/\Gamma$ also fixes the regime for all $N$.

In Fig. 8a, we see that $\alpha$ of the dipolar model coincides
with that of the non-interacting model in the regime of
strong driving for all $N \le 2000$, i.e., $\Omega \gg \Omega_C^{\text{Dipolar}} =
c\Gamma\sqrt{N} \approx 3.6\Gamma$. This confirms our expectation that in the
strong driving regime, the dipolar model behaves like the
non-interacting model. In Fig. 8a, we also see the clear
difference between the CRF and the dipolar models in
the strong driving regime. We show the CRF model for
$N=10$ to compare with the experimental prediction [1],
where it is reported that a cloud of $N=2000$ atoms
can equivalently be described by the CRF model of a
reduced effective atom number, $N_{\text{eff}}=10$. Close to the
peak, the $\alpha$ of the dipolar model agrees with that of the
CRF model in a small range of $\Omega/\Gamma$. Even though, the
two models have very different $\alpha$ in general, this was
the regime probed by the experiment, partly justifying
the conclusions drawn in Ref. [1]. Reaching the genuine
strong drive limit requires large enough Rabi frequencies
which are not so accessible in current experiments.

In Fig. 8b, we calculate the *N*-scaling by properly fix-
ing the scaled driving strength, $\Omega/\Omega_C^{\text{Dipolar}}$, and as ex-
pected, we find no signatures of subradiance. Moreover,
the divergence of $\alpha$ seen in Fig. 8a vanishes. Instead at
weak driving, $\alpha$ decreases from 2 (perfect phase match-
ing) to close to 1 (randomized phases) for the dipolar
model as *N* is increased and the collective shifts get
stronger, suppressing the *N*²-enhancement. In the strong
driving regime, $\alpha \to 1$ for all *N*, as the system becomes
single-particle-like. For the CRF model, in the large-*N*
limit, $\alpha$ scales as *N*² in all driving regimes.
IV. MODIFIED-CRF MODEL

In this section, we propose a simplified theoretical model at the mean-field level to describe the emergent properties of our system. In the microscopic picture of the dipolar model, the dispersive part of the dipole-dipole interactions between atoms leads to a frequency shift in the transition frequency of each atom, as discussed in Sec. III A. This shift can equivalently be captured by an additional term $\sum_{i=1}^{N} \delta_i(t)\sigma_i^z$ in the CRF Hamiltonian as

$$ \hat{H}_{\text{Mod-CRF}} = -\Omega \hat{S}_x + \sum_{i=1}^{N} \delta_i(t) \sigma_i^z, \quad (28) $$

where the frequency shift of an atom $i$ is described as $\delta_i(t) = \sum_{j \neq i} \text{Re}\tilde{G}_{ij}\langle\hat{\sigma}_j^z(t)\rangle/4$, inspired from the short-time dynamics of the MF dipolar model. To distinguish the collective slow varying part of the interaction-induced shift from the inhomogeneous, fast time varying effects, we express the frequency shift as $\delta_i(t) = \bar{\delta}(t) + h_i(t)$ where, $\bar{\delta}(t) = \sum_{i,j \neq i} \text{Re}\tilde{G}_{ij}\langle\hat{\sigma}_j^z(t)\rangle/(4N)$ is the average shift across all the atoms in the cloud at any time $t$ and $h_i(t)$ is the fast time varying inhomogeneous component, which describes beyond short-time dynamics. The homogeneous component, $\bar{\delta}(t)$, is dominated by the shifts that add up along the axial direction due to the constructive interference of the laser-induced phases with the interaction-induced phases. To capture the constructive shifts in a simplified way, we assume that $\bar{\delta}(t)$ for our cloud can be described by a collective system, such as a translationally invariant atomic array, giving $\bar{\delta}(t) \approx \sum_i \text{Re}\tilde{G}_{i0}\langle\hat{S}_z(t)\rangle/(2N)$. The interaction-induced time varying inhomogeneous component, $h_i(t)$, accounts for the random and dynamic evolution of the atomic dipolar phases. $h_i(t)$ has a zero mean, is time-dependent, and is therefore not removable by a simple echo pulse. For simplicity, these conditions are roughly incorporated by considering the $h_i$ functions as stochastic white noise variables with spectral function: $\bar{h}_i(t)h_j(t') = (\gamma_d/2)\delta_{ij}\delta(t-t')$ [69].

The net dephasing arising from the $h_i(t)\hat{\sigma}_i^z$ term can be accounted for by local jump operators $\sqrt{\gamma_d}\hat{\sigma}_i^z$ in the master equation formulation, where $\gamma_d$ is a common dephasing rate for all the atoms, which we set to be proportional to the variance of the frequency shifts, which scales linearly with the density of the atomic ensemble [6, 18, 31]. Thus, we define the dephasing rate as $\gamma_d = c_d N \gamma_s$, where $c_d$ is a phenomenological constant that depends on the overall volume and geometry of the atomic cloud and $\gamma_s$ is the spontaneous emission decay rate.

Accounting for the overall elastic interactions, the collective superradiant emission from nearby atoms, and the single-particle spontaneous emission (since the purely collective behavior is only possible in the high density limit, not achieved in the experiment), we obtain the modified CRF model described by a master equation of the form $\hat{\rho} = -i[\hat{H}_{\text{Mod-CRF}}, \hat{\rho}] + \mathcal{L}_{\text{Mod-CRF}}(\hat{\rho})$, where

the Hamiltonian is

$$ \hat{H}_{\text{Mod-CRF}} = -\Omega \hat{S}_x + \chi \langle \hat{S}_z \rangle \hat{S}_z, \quad (29) $$

where $\chi = 2\bar{\delta}(t)/\langle\hat{S}_z(t)\rangle = \sum_{i \neq 0} \text{Re}\tilde{G}_{i0}/N$ is a constant. $\chi\langle\hat{S}_z(t)\rangle$ acts as a time-dependent global magnetic field, which describes the shearing of the collective Bloch vector via one-axis twisting (OAT) at the mean-field level. The Lindbladian is expressed as

$$ \mathcal{L}_{\text{Mod-CRF}}(\hat{\rho}) = \frac{\Gamma_D}{2} \left( 2\hat{S}_-\hat{\rho}\hat{S}_+ - \{\hat{S}_+\hat{S}_-, \hat{\rho}\} \right) \\ + \frac{\gamma_s}{2} \sum_j \left( 2\hat{\sigma}_j^- \hat{\rho}\hat{\sigma}_j^+ - \{\hat{\sigma}_j^z, \hat{\rho}\}/2 - \hat{\rho} \right) \\ + \gamma_d \sum_j \left( \hat{\sigma}_j^z \hat{\rho}\hat{\sigma}_j^z - \hat{\rho} \right), \quad (30) $$

where the first line is the collective dissipation in the CRF model with rate $\Gamma_D = \sum_{j \neq 0} \text{Im}\tilde{G}_{0j}/N$, the second line is the single-particle spontaneous emission with decay rate $\gamma_s$, and the third line is the single-particle dephasing with rate $\gamma_d$. Unlike the CRF model, this modified model does not preserve $\langle\hat{S}^2\rangle$. The case without dephasing ($\gamma_d = 0$) and OAT ($\chi = 0$) has been predicted to show a first-order phase transition and bistability [51, 61, 70] at a critical driving strength given by $\Omega_C \approx \Gamma_D N/(2\sqrt{2})$. Hence, spontaneous emission only changes the critical drive by a factor of $1/\sqrt{2}$ compared to the original CRF model ($\Omega_C^{\text{CRF}} = \Gamma_N N/2$). When we include the effect of dephasing ($\gamma_d \neq 0$) and dispersive interactions ($\chi \neq 0$), we find that the modified model still undergoes a first-order phase transition with bistability, but the critical point depends on the $N$-scaling of $\gamma_d$. We obtain the mean-field steady-state in the superradiant phase for $N \gg 1$ from the following self-consistent equations (see Appendix D for derivation) –

$$ \frac{\langle \hat{S}_z \rangle}{N/2} = -\frac{1}{2} \pm \frac{1}{2} \sqrt{1 - 8 \left(1 + \frac{4\gamma_d}{\gamma_s}\right) \left(\frac{\langle \hat{S}_x \rangle^2 + \langle \hat{S}_y \rangle^2}{N^2}\right)}, \\ \frac{\langle \hat{S}_x \rangle^2 + \langle \hat{S}_y \rangle^2}{N^2} = \frac{\Omega^2 - \Gamma_D \gamma_s N (1 + \langle \hat{S}_z \rangle / (N/2))}{N^2 (\Gamma_D^2 + \chi^2)} $$

The critical point is defined as the driving strength at which the solution above ceases to be valid and is obtained as (with $\gamma_d = c_d N \gamma_s / 4$)

$$ \Omega_C^{\text{Mod-CRF}} = \Gamma_D \sqrt{N} f \quad (31) $$

where $f \approx \frac{1}{2\sqrt{2}} \sqrt{\frac{1+(\chi/\Gamma_D)^2}{c_d} + \frac{4\gamma_s}{\Gamma_D}}$ is a constant.

In Fig. 9, we show the steady-state values of the Mod-CRF model obtained from mean-field numerics, where we have set $c_d = 0.002$. This specific value of $c_d$ is not special and we find qualitatively similar results for other values as long as $c_d \ll 1$. We have $\chi = \sum_{j \neq 0} \text{Re}\tilde{G}_{0j}/N \approx 0.003$ and $\Gamma_D = \sum_{j \neq 0} \text{Im}\tilde{G}_{0j}/N \approx 0.002$ for our pencil-shaped cloud (see Appendix D for details). The steady-state of
<figure><img src="image_9.png" /></figure>

FIG. 9. Steady-state of the modified CRF model from MF numerics: (a) atomic inversion $\langle\hat{S}_z\rangle/(N/2)$, (b) total coherence $|\langle\hat{S}_+\rangle|/N$, (c) intensity $\langle\hat{S}_+\hat{S}_-\rangle$, which largely scales as $\sim N$ due to dephasing, (d) two-photon correlation function $g_2(0)$. (e) The x-component of the coherence, $|\langle\hat{S}_x\rangle|/|\langle\hat{S}_+\rangle| \neq 0$, arising from the mean-field frequency shift. In the large $N$ limit, $\Omega_C^{\text{Mod-CRF}} \sim \Gamma\sqrt{N}$ is the critical drive for this model (Eq. (31)).

the Mod-CRF model shows qualitatively similar behaviors to the MF dipolar model (Fig. 5). The critical drive strength for the Mod-CRF model, $\Omega_C^{\text{Mod-CRF}}$, scales as $\sim \sqrt{N}$, similar to $\Omega_C^{\text{Dipolar}}$. This $\sqrt{N}$-scaling is also recovered in the numerics as the scaled numerical curves (x-axis scaling: $\Omega/\Omega_C^{\text{Mod-CRF}}$) of the atomic inversion become indistinguishable at large $N$ (Fig. 9a). Similarly, the scaled numerical curves of the intensity (Fig. 9c), $g_2(0)$ (Fig. 9d), and the atomic coherences (Figs. 9b and 9e) also look indistinguishable. In the weak-drive regime, the real-part of the atomic coherence, $\langle\hat{S}_x\rangle$, is non-zero and increases with $N$ due to collective shifts (Fig. 9e). In the strong-driving regime ($\Omega \gg \Omega_C^{\text{Mod-CRF}}$), the intensity of the Mod-CRF becomes single-particle-like, i.e., $\langle\hat{S}_+\hat{S}_-\rangle = N/2$, as the coherences go to zero (Fig. 9b). In this regime, the atomic inversion goes to zero as well (Fig. 9a), and $g_2(0)$ reaches its thermal value of 2 (Fig. 9d). These characteristics describe a departure from the CRF model and qualitatively resemble the properties of the MF dipolar model (Fig. 5, Fig. 6a).

V. COMPARISONS WITH EXPERIMENT

In this section we compare our mean-field numerics for the dipolar model with experimental results [1]. The experimental system and measurement protocols are all described in Ref. [1]. The participating states correspond to the $\sigma_+$-polarised atomic transition between the two levels $|5S_{1/2}, F = 2, m_F = 2\rangle \rightarrow |5P_{3/2}, F = 3, m_F = 3\rangle$ in $^{87}$Rb with linewidth $\Gamma = 2\pi \times 6$MHz. The Clebsch-Gordan coefficient for this transition is 1. This system is

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

FIG. 10. Excitation fraction $n_e$ vs time for fixed $\Omega/\Gamma = 4.5$: MF dipolar numerics for $N_{\text{MF}}$ atoms (red thick) on top of experimental data (blue traces) for $N_{\text{exp}}$ taken from Ref. [1]. The shaded red regions show the numerics data in the specified $N_{\text{MF}}$ ranges. The dashed black line is the solution of the non-interacting model.

the same as the one described in Fig. 2.

In Fig. 10, we show the dynamics of the total excitation fraction $n_e(t) = \sum_i (\langle\sigma_z^i(t)\rangle + 1)/(2N)$. To do this, we evolve the system under a continuous drive with Rabi frequency $\Omega = 4.5\Gamma$ from the ground state to the steady-state. We find that the mean-field numerics (red lines) agree more or less with the experiment (blue traces), for different values of $N$. The shaded red region shows $n_e$ within a range of $N$ values to account for experimental uncertainties. Both the numerics and the experimental data feature a dephasing of the Rabi oscillations with increasing particle number, which is not seen in the non-interacting model (black dashed).

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

FIG. 11. $\langle\hat{S}_z\rangle/(N/2)$: Steady-state atomic inversion. (a) MF Dipolar numerics (solid lines) and experimental data (symbols) from Ref. [1]. (b) MF data with x-axis scaling $\sim \Omega/\Omega_C^{\text{Dipolar}}$, showing collapse of curves for different $N$.

In Fig. 11a, we look at the steady-state atomic inversion and find very good agreement of our numerics (continuous lines) with the experimental data (symbols) for $N = 300$. An experimental data point in Fig. 11a for $N = 1500$ at $\Omega/\Gamma \approx 4.5$ ($\langle\hat{S}_z\rangle/(N/2) \approx -0.45$) is not consistent with the steady-state in Fig. 10c ($n_e(t \rightarrow \infty) \approx 0.5 \Rightarrow \langle\hat{S}_z\rangle/(N/2) \approx 0$). This could be due to calibration errors or other experimental systematics. Other than these quantitative differences, it is clear that there is a qualitative agreement with the theory in the trends observed in the experiment. In Fig. 11b, we plot the
<figure><img src="image_12.png" /></figure>

FIG. 12. (a) Steady-state intensity of the MF dipolar model in the forward direction ($\theta = 0$, dashed), averaged over a solid angle ($\Delta\theta = \theta_{out} - \theta_{in}$) around the forward direction with $\theta_{in} = \pi/22$ and $\theta_{out} = \pi/7 \pm \pi/50$ (solid lines, shaded region), and experimental data (dots) from Ref. [1]. $I_N$ is a fitting factor for the MF dipolar model chosen to match the theory peak with the experimental values. (b) MF data ($\theta_{in} = \pi/22$ and $\theta_{out} = \pi/7$) with x-axis scaling $\sim \Omega/\Omega_C^{\text{Dipolar}}$ shows collapse of curves, as expected.

numerical data with our scaling of the Rabi frequency, $\Omega/\Omega_C^{\text{Dipolar}} \sim \Omega/(0.08\Gamma\sqrt{N})$, and the data collapses to a single curve, as previously discussed.

In Fig. 12a, we compare the steady-state intensity from MF numerics with the experimental data [1] and find fair agreement. Here, we look at the steady-state intensity, $I_0(\vec{k}) \sum_{ij} \langle\hat{\sigma}_i^+ \hat{\sigma}_j^-\rangle e^{i\vec{k} \cdot \vec{r}_{ij}}$, strictly in the forward direction ($\vec{k} = \vec{k}_L$, dashed) and averaged over a solid angle around the forward direction (shaded region). Here, $\hat{k} = (\cos\theta, \sin\theta \cos\phi, \sin\theta \sin\phi)$, is the direction of observation, with $\theta$ the angle from the $\hat{x}$-axis (forward direction) and $\phi$ the azimuthal angle in the $y-z$ plane from $\hat{y}$. For the circularly-polarized transition, $\hat{e}_q = \hat{e}_+ = -(\hat{x} + i\hat{y})/\sqrt{2}$, and we get $I_0(\vec{k}) \propto (1 + \sin^2\theta \sin^2\phi)/2$ (see Appendix A 3). As some of the light in the forward direction is filtered out in the experiment to remove the laser light, we average the intensity over an annular region [31, 71], such that the averaged intensity is expressed as $\int_{\theta_{in}}^{\theta_{out}} d\theta \int_0^{2\pi} d\phi \sin\theta I_0(\vec{k}) \sum_{ij} \langle\hat{\sigma}_+ \hat{\sigma}_-\rangle e^{i\vec{k} \cdot \vec{r}_{ij}}/A$, where $\theta_{out} = \pi/7 \pm \pi/50$, $\theta_{in} = \pi/22$, and $A = \int_{\theta_{in}}^{\theta_{out}} d\theta \int_0^{2\pi} d\phi \sin\theta$ is the normalization factor. The values of $\theta_{out}$ and $\theta_{in}$ are extremely difficult to determine experimentally and hence, we consider them as fitting parameters, while making sure their values remain within a reasonable range.

Furthermore, we have included direct contributions from the probe light to the intensity, considering that the filter may not be perfect (see Appendix E for details). We use a $\sim 4\%$ leakage of probe light intensity, which is also a fitting parameter and its value is consistent across Figs. 12, 13, and 14. The experimental data, which is the photon rate in arbitrary units, is scaled with a fitting factor, which was used in Ref. [1] to compare their data with the CRF model. To compare our results with the experimental data, we scale the intensity with a

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

FIG. 13. Dynamics of intensity for the MF dipolar model in the forward direction (black dashed), averaged over a solid angle ($\Delta\theta = \theta_{out} - \theta_{in}$) about the forward direction with $\theta_{in} = \pi/22$ and $\theta_{out} = \pi/7 \pm \pi/50$ (blue solid, blue shaded region), experimental data (red traces) [72]. The intensity $I(t)$ has been scaled by its steady-state value, $I_{ss}$, for each dataset. (a)-(c) weak and intermediate driving regimes ($\Omega \le \Omega_C^{\text{Dipolar}}$), (d)-(f) strong-driving regime ($\Omega > \Omega_C^{\text{Dipolar}}$).

fitting factor $I_N$, which is obtained by matching the peak of the experimental curve with the peak of the numerics for each $N$.

In Fig. 13, we compare the dynamics of the forward intensity, $I(t)$, measured in the experiment (red traces) [72] with the intensity from the MF numerics in the strictly forward direction (black dashed) and the MF numerics data averaged over a solid angle about the forward direction (blue solid) with $\theta_{in} = \pi/22$ and a range of $\theta_{out} = \pi/7 \pm \pi/50$ (blue shaded region). This is the same range of solid angle values discussed earlier for the steady-state intensity. To properly compare the bare photon rate from the experiment with the numerical values of the intensity, we scale the dynamical intensity $I(t)$ by its steady-state value $I_{ss}$, thus making the scaled intensity, $I(t)/I_{ss}$, independent of the distance between the atomic cloud and the detector, which is difficult to determine exactly for the experimental setup. In the weak ($\Omega < \Omega_C^{\text{Dipolar}}$) and intermediate ($\Omega \sim \Omega_C^{\text{Dipolar}}$) driving regimes, the experiment agrees well with the purely forward intensity of the dipolar model (black dashed) for different $N$, as shown in Figs. 13(a)-(c). In Figs. 13(d)-(f), we see that in the strong-driving regime ($\Omega > \Omega_C^{\text{Dipolar}}$), the experiment agrees well with the dipolar model intensity when averaged over the solid angle (blue shaded region) for different $N$. Thus, the solid angle $(\theta_{in}, \theta_{out})$, for which the numerics and the experiment agree, depends on the drive strength, $\Omega/\Gamma$. One possible explanation is the fact that the radiative force of the laser pushes the atomic cloud closer to the filter and the detector, changing the solid angle of the detected light. Of course, the radiative force may also lead to other motional effects such as dephasing, which we have not considered here.

In Fig. 14, we look at the steady-state equal-time two-point photon correlation function $\hat{g}_2(0)$, as defined in Eq. (9), in the strictly forward direction ($\theta = 0$, black
<figure><img src="image_14.png" /><figcaption>FIG. 14. $\tilde{g}_2(0)$ in the steady-state for $N = 1400$: in the forward direction ($\theta = 0$, black dashed), averaged over a solid angle ($\Delta\theta = \theta_{\text{out}} - \theta_{\text{in}}$) about the forward direction with $\theta_{\text{in}} = \pi/22$ and $\theta_{\text{out}} = \pi/16 \pm \pi/100$ (green shaded region), experimental data (red dots) from Ref. [1].</figcaption></figure>

dashed) and averaged over an annular region around the forward direction, same as we did for the intensity earlier, with the same $\theta_{\text{in}} = \pi/22$ but with a smaller $\theta_{\text{out}} = \pi/16 \pm \pi/100$ (green shaded region). Again, we find fair agreement between the experimental (red dots) [1] and numerical (green shaded region) values, when averaged over the solid angle. We see some discrepancies at intermediate $\Omega/\Gamma$ values, which could arise from the dephasing of coherences due to atomic motion in the experiment. The best-fit $\theta_{\text{out}}$ values are not identical to the ones used for the intensity plots but such changes could occur between different experimental measurements. Here, we plot the numerical values for the non-interacting model instead of the dipolar model, as we find it is quicker to integrate over the solid angle analytically over the non-interacting model's steady-state, than to do a discrete sum in this case. We have already shown that the $\tilde{g}_2(0)$ values for the non-interacting model coincide with those of the MF dipolar model in Fig. 5d, so we expect these values to be valid for our system. For further details about the analytical averaging of $\tilde{g}_2(0)$ for the non-interacting model, see Appendix A.

## VI. CONCLUSIONS

We have studied the age-old, albeit complicated, problem of characterizing the emission properties of a spatially-extended ensemble of driven interacting atoms in free-space. What makes this problem complicated is the lack of symmetries such as translational invariance in periodic arrays and permutational invariance in optical cavities. In the absence of symmetries, an exact microscopic treatment of the system is impossible due to its exponentially large Hilbert space ($\sim 2^N$). Nevertheless, we are able to greatly reduce the complexity of the problem by identifying two distinct driving regimes (weak and strong), where a mean-field (MF) treatment is almost exact and the problem becomes at least numerically tractable with just 3N degrees of freedom. Furthermore, we find that our MF numerics are able to qualita-

tively capture the physics across all driving regimes and beyond-MF methods (MACE-MF and cumulant) lead to negligible corrections.

We find a crossover between non-equilibrium phases in the steady-state of our system as a function of the driving strength $\Omega/\Gamma$, akin to the superradiant phase transition in the CRF model. We find that the inhomogeneity of free-space dipolar interactions plays a key role in making our system strikingly different from the CRF model. This is most prominent in the N-scaling of the critical order parameter, which is $\propto \sqrt{N}$ in our system and differs from that in the permutationally invariant CRF model ($\propto N$). Moreover, the strongly-driven phase in our model is completely mixed and single-particle-like, whereas the “normal” phase of the CRF has collective quantum correlations.

At the cloud densities considered, our system is governed by the far-field dipolar interaction ($\sim 1/r$), which is long-range in quasi-1D. The collective effects of the long-range interactions are tempered by the dephasing arising from inhomogeneous frequency shifts and free-space emission. Using these physical insights, we propose a heuristic modified CRF model capable of describing our system, which includes single-particle dephasing, shearing, and spontaneous emission. This modified CRF is able to qualitatively reproduce the emergent features of our complicated microscopic dipolar model and is analytically solvable at the mean-field level.

More importantly, our results are able to reproduce most of the experimental findings from a recent work [1]. Our work bridges the theoretical gap between the widely studied CRF model and the spatially-extended inhomogeneous atomic ensembles accessible in current experimental setups. While previous work has shed some light on the properties of dilute atomic ensembles, our work extends this knowledge to moderately dense ensembles where collective effects are relevant beyond just frequency shifts and linewidth broadening.

For future work, it would be useful to measure (using Ramsey spectroscopy) the collective frequency shift arising from dipolar interactions, which leads to a non-zero $\langle \hat{S}_x \rangle$ in the dipolar model, as opposed to the CRF and non-interacting models, where $\langle \hat{S}_x \rangle = 0$. This is a smoking gun signature of many-body effects in this system. The next step would be studying highly dense atomic ensembles, where the near-field elastic dipolar interactions ($\sim 1/r^3$) dominate the physics and the MF treatment breaks down. The effects of atomic motion would become pertinent in the presence of strong light-matter interactions in dense ensembles. Another avenue worth exploring is the case of multilevel atoms, which includes the internal level structure of alkaline-earth(-like) atoms due to hyperfine splitting. Even in the weak-driving regime, the multilevel system can have quantum correlations in the ground-state manifold [73] and is distinct from the semiclassical two-level picture. While there has been some recent work on the emission properties of multilevel arrays
in free-space, it is largely confined to the early-time dynamics [74] and inverted arrays [75]. Characterizing the steady-state properties of multilevel ensembles is crucial for current experiments and remains an unsolved problem.

During the completion of our work, we became aware of a related recent work that reaches similar conclusions as we using a complementary theoretical treatment [71].

## ACKNOWLEDGMENTS

We thank Jarrod T. Reilly, Eric Yilun Song, and James K. Thompson for insightful discussions and feedback on our manuscript. We wish to acknowledge support by the VBFF, AFOSR grant FA9550-18-1-0319, by the NSF JILA-PFC PHY-2317149, QLCI-OMA-2016244, by the U.S. Department of Energy, Office of Science, National Quantum Information Science Research Centers Quantum Systems Accelerator, and by NIST.

## Appendix A: Non-interacting model

Here we consider the non-interacting system and its steady-state solution. The master equation describing this system is $\hat{\rho} = -i[\hat{H}_0, \hat{\rho}] + \mathcal{L}_0(\hat{\rho})$, where $\hat{H}_0 = -(\Omega/2) \sum_k (e^{i\vec{k}_L \cdot \vec{r}_k} \sigma_k^+ + \text{h.c.})$ is the Hamiltonian and $\mathcal{L}_0(\hat{\rho}) = (\Gamma/2) \sum_j (2\sigma_j^- \hat{\rho} \sigma_j^+ - \{\sigma_j^z, \hat{\rho}\}/2 - \hat{\rho})$ is the Lindbladian. The atomic equations of motion reduce to the standard optical Bloch equations:

$$ \langle \dot{\sigma}_k^z \rangle = -\Gamma (\langle \sigma_k^z \rangle + 1) + i\Omega (e^{i\vec{k}_L \cdot \vec{r}_k} \langle \sigma_k^+ \rangle - e^{-i\vec{k}_L \cdot \vec{r}_k} \langle \sigma_k^- \rangle), \quad (A1) $$

$$ \langle \dot{\sigma}_k^+ \rangle = -\frac{\Gamma}{2} \langle \sigma_k^+ \rangle + i \frac{\Omega}{2} e^{-i \vec{k}_L \cdot \vec{r}_k} \langle \sigma_k^z \rangle. \quad (A2) $$

and the steady-state can be obtained by setting $\langle \dot{\sigma}_k^z \rangle = 0$ and $\langle \dot{\sigma}_k^+ \rangle = 0$ for all $k$, as

$$ \langle \sigma_k^z \rangle = - \frac{\Gamma^2}{\Gamma^2 + 2\Omega^2} \equiv -R \quad (A3) $$

$$ \langle \sigma_k^+ \rangle = -i \frac{\Omega}{\Gamma} e^{-i \vec{k}_L \cdot \vec{r}_k} R. \quad (A4) $$

The steady-state intensity in an observed direction $\vec{k}$ is

$$ \frac{I(\vec{k})}{I_0(\vec{k})} = \sum_j \left( \frac{\langle \sigma_j^z \rangle + 1}{2} \right) + \sum_{i \neq j} \langle \sigma_i^+ \rangle \langle \sigma_j^- \rangle e^{i\vec{k} \cdot (\vec{r}_i - \vec{r}_j)} \\ = \frac{(1-R)}{2} N + \frac{\Omega^2}{\Gamma^2} R^2 \sum_{i \neq j} e^{i(\vec{k}-\vec{k}_L) \cdot (\vec{r}_i - \vec{r}_j)}. \quad (A5) $$

The two photon correlation function in an observed direction $\vec{k}$ is

$$ \tilde{g}_2(0)(\vec{k}) = \frac{\langle \tilde{S}^+(\vec{k}) \tilde{S}^+(\vec{k}) \tilde{S}^-(\vec{k}) \tilde{S}^-(\vec{k}) \rangle}{\langle \tilde{S}^+(\vec{k}) \tilde{S}^-(\vec{k}) \rangle^2}, \quad (A6) $$

where $\tilde{S}_{\pm}(\vec{k}) = \sum_i \sigma_i^{\pm} e^{\pm i\vec{k} \cdot \vec{r}_i}$. For non-interacting particles, the correlations factor for different atoms as $\langle \sigma_{\alpha}^i \sigma_{\beta}^j \rangle = \langle \sigma_{\alpha}^i \rangle \langle \sigma_{\beta}^j \rangle$. For the same atom, we use the commutation relations of the Pauli operators as $\langle \sigma_{\alpha}^i \sigma_{\beta}^j \rangle = \delta_{\alpha,\beta} \mathbb{1} + 2i\epsilon^{\alpha\beta\gamma} \langle \sigma_{\gamma}^i \rangle$, where $\alpha, \beta, \gamma \in \{x, y, z\}$. Using this scheme, the two photon correlation function along the measurement direction can be expressed as

$$ \tilde{g}_2(0)(\vec{k}_L) = \left( \frac{1}{2} \sum_{\langle ij \rangle} (\langle \sigma_i^z \rangle + 1) (\langle \sigma_j^z \rangle + 1) \right. \\ \left. + 2 \sum_{\langle ijk \rangle} (\langle \sigma_i^z \rangle + 1) \langle \tilde{\sigma}_j^+ \rangle \langle \tilde{\sigma}_k^- \rangle + \sum_{\langle ijkl \rangle} (\langle \tilde{\sigma}_i^+ \rangle \langle \tilde{\sigma}_j^+ \rangle \langle \tilde{\sigma}_k^- \rangle \langle \tilde{\sigma}_l^- \rangle) \right) \\ \times \left( \frac{1}{2} \sum_i (\langle \sigma_i^z \rangle + 1) + \sum_{\langle ij \rangle} (\langle \tilde{\sigma}_i^+ \rangle \langle \tilde{\sigma}_j^- \rangle) \right)^{-2}, \quad (A7) $$

where we have suppressed the notation to denote $\tilde{\sigma}_j^+ \equiv \sigma_j^+ e^{i\vec{k}_L \cdot \vec{r}_j}$ and $\langle ij \dots \rangle$ denotes sum over unlike indices, i.e., $i \neq j \neq \dots$.

## 1. Weak-driving regime

In the weak-driving regime ($\Omega \ll \Gamma$), the steady-state solution can be expressed as

$$ \langle \sigma_k^z \rangle = - \left( 1 - \frac{2\Omega^2}{\Gamma^2} \right) \quad (A8) $$

$$ \langle \sigma_k^+ \rangle = i \frac{\Omega}{\Gamma} e^{-i \vec{k}_L \cdot \vec{r}_k}, \quad (A9) $$

and the intensity along the laser wavevector is

$$ \frac{I(\vec{k}_L)}{I_0(\vec{k}_L)} = \frac{\Omega^2}{\Gamma^2} N \left( N - \frac{1}{2} \right) \approx \frac{\Omega^2}{\Gamma^2} N^2, \quad (A10) $$

which scales as $N^2$ in the large-$N$ limit due to the constructive interference of coherences along the measurement direction. In this regime, the two-photon correlation function can be expressed as

$$ \tilde{g}_2(0)(\vec{k}_L) = 1 - \frac{2}{N} + \frac{1}{N^2} \Rightarrow \lim_{N \to \infty} \tilde{g}_2(0)(\vec{k}_L) = 1, \quad (A11) $$

which is consistent with the value for a coherent state. In the extreme weak-driving limit, the atoms are very weakly-excited such that $\langle \sigma_j^z \rangle \to -1$, which along with Eq. (A7) gives $\tilde{g}_2(0)(\vec{k}_L) \to 1$.

## 2. Strong-driving regime

In the strong-driving regime ($\Omega \gg \Gamma$), the steady-state solution can be expressed as

$$ \langle \sigma_k^z \rangle = - \frac{\Gamma^2}{2\Omega^2} \quad (A12) $$
$$ \langle \sigma_k^+ \rangle = i \frac{\Gamma}{2\Omega} e^{-i\vec{k}_L \cdot \vec{r}_k}, \qquad (A13) $$

and the intensity along the laser wavevector given by

$$ \frac{I(\vec{k}_L)}{I_0(\vec{k}_L)} = \frac{N}{2} + \frac{\Gamma^2}{4\Omega^2} N(N-2). \qquad (A14) $$

When the drive is in the regime $\frac{\Gamma^2}{4\Omega^2} < 1/N \Rightarrow \Omega > \Gamma\sqrt{N}/2$ for all N, the intensity scales as N. As the drive gets extremely large, $\Omega \gg N\Gamma/2$, the state gets fully mixed and $I(\vec{k}_L) \to \frac{N}{2}$. In this limit, we get

$$ \tilde{g}_2(0) = 2 - \frac{2}{N} \Rightarrow \lim_{N \to \infty} \tilde{g}_2(0) = 2, \qquad (A15) $$

consistent with the fact that the system is described by a thermal state. In the extreme strong-driving limit, the system is maximally mixed such that $\langle \sigma_j^z \rangle \to 0$ and $\langle \sigma_j^{+-} \rangle \to 0$. Plugging this into Eq. (A7) also gives $\tilde{g}_2(0)(\vec{k}_L) \to 2$.

### 3. Averaging over a solid angle

For our pencil-shaped cloud (Fig. 2), the atomic positions are distributed in a Gaussian distribution,

$$ \rho(x, y, z) = \frac{N \exp \left( -\frac{x^2}{2\sigma_{\text{ax}}^2} - \frac{y^2+z^2}{2\sigma_{\text{rad}}^2} \right)}{(2\pi)^{3/2} \sigma_{\text{ax}} \sigma_{\text{rad}}^2}, \qquad (A16) $$

where $\sigma_{\text{ax}} = 20\lambda$ and $\sigma_{\text{rad}} = \lambda/2$ are the axial and radial standard deviations of the cloud, respectively. Then, in the large-$N$ limit, we can replace the sum over atoms $i$ and $j$ in the intensity in Eq. (A5) by an integral over the positions, weighted by their distribution $\rho(x, y, z)$, as

$$ \frac{I(\vec{k})}{I_0(\vec{k})} \approx \frac{(1-R)}{2} N - \frac{\Omega^2}{\Gamma^2} R^2 N + \frac{\Omega^2}{\Gamma^2} R^2 \int d^3 r \int d^3 r' \rho(\vec{r}) \rho(\vec{r}') e^{i(\vec{k} - \vec{k}_L) \cdot (\vec{r} - \vec{r}')} $$

where in the first line we have subtracted the contribution from the $i=j$ term in the integral. We define $\vec{q} = \vec{k} - \vec{k}_L = (q_x, q_y, q_z)$. Then, it is easy to do obtain the integral in Cartesian coordinates as

$$ A = \left| \int_{-\infty}^{\infty} dx \int_{-\infty}^{\infty} dy \int_{-\infty}^{\infty} dz \rho(x, y, z) e^{i(q_x x + q_y y + q_z z)} \right|^2 \\ = N^2 \exp\left(-\left(\sigma_{\text{ax}}^2 q_x^2 + \sigma_{\text{rad}}^2 q_y^2 + \sigma_{\text{rad}}^2 q_z^2\right)\right). \qquad (A17) $$

Now, we substitute $\vec{k}_L = (2\pi/\lambda)\hat{x}$ and $\vec{k} = (2\pi/\lambda)(\cos\theta, \sin\theta\cos\phi, \sin\theta\sin\phi)$, with $\theta$ the angle from the $\hat{x}$-axis (forward direction) and $\phi$ the azimuthal angle in the $y-z$ plane from $\hat{y}$. Then, we obtain the intensity as

$$ \frac{I(\theta, \phi)}{I_0(\vec{k})} = \frac{2(\Omega/\Gamma)^4}{(1 + 2(\Omega/\Gamma)^2)^2} N \left[ 1 \right] $$

$$ + \frac{N \exp\left(-\left(2\pi\right)^2\left(\sigma_{\text{ax}}^2\left(1 - \cos\theta\right)^2 + \sigma_{\text{rad}}^2\sin^2\theta\right)\right)}{2(\Omega/\Gamma)^2} \right]. \qquad (A18) $$

Close to the forward direction $\theta \equiv \delta\theta \approx 0$, the intensity can be expressed as

$$ \frac{I(\delta\theta, \phi)}{I_0(\vec{k})} \approx \frac{2(\Omega/\Gamma)^4}{(1 + 2(\Omega/\Gamma)^2)^2} N \left[ 1 + \frac{N e^{-(2\pi)^2(\sigma_{\text{rad}}^2\delta\theta^2 + \sigma_{\text{ax}}^2\delta\theta^4/4)}}{2(\Omega/\Gamma)^2} \right], \qquad (A19) $$

such that the $N^2$-enhancement falls off exponentially fast as the observation direction deviates from the forward direction.

We need to multiply the intensity above (Eq. (A18)) by the geometric factor $I_0(\hat{k})$. Due to the large distance of the detector from the cloud, $I_0(\hat{k})$ is the factor associated with the far-field part of the Electromagnetic Green's tensor and appears in the expression of the dipolar intensity pattern as $I_0(\hat{k}) = \hat{e}_q \cdot (1 - \hat{k} \otimes \hat{k}) \cdot (1 - \hat{k} \otimes \hat{k}) \cdot \hat{e}_q^*$, where $\hat{e}_q$ is the polarization of light associated with the atomic transition corresponding to fluorescence. For circularly polarized light, we get $I_0(\hat{k}) = (1 + \sin^2\theta \sin^2\phi)/2$.

Then, the intensity can be averaged over an annular region, which has $\theta_{\text{in}} = \pi/22$ and $\theta_{\text{out}} = \pi/7$ as the inner and outer boundaries, to obtain $\langle I \rangle \equiv \int_0^{2\pi} d\phi \int_{\theta_{\text{in}}}^{\theta_{\text{out}}} d\theta I(\theta, \phi) \sin\theta(1 + \sin^2\theta \sin^2\phi)/(2A) \approx \frac{2(\Omega/\Gamma)^4}{(1+2(\Omega/\Gamma)^2)^2} N [c_A^I + \frac{c_A^{II}N}{2(\Omega/\Gamma)^2}]$, where $A = \int_0^{2\pi} d\phi \int_{\theta_{\text{in}}}^{\theta_{\text{out}}} d\theta \sin\theta$. It can be seen that the extra factor of $c_A^I/c_A^I \approx 0.022$, due to the averaging over a finite solid angle around the forward direction, reduces the relative strength of the coherent emission and shifts the peak of the intensity to lower $\Omega/\Gamma$. We have substituted $\sigma_{\text{ax}} = \lambda/2$ and $\sigma_{\text{rad}} = 20\lambda$ above, corresponding to our system (Fig. 2), to get a numerical value.

Similarly, we can also integrate over the atomic positions, weighted by the Gaussian distribution, to obtain the analytical expression for the steady-state two-photon correlation function, $g_2(\vec{k}, \vec{k}')$, defined in Eq. (9), in terms of $\hat{k} \equiv (\theta, \phi)$ and $\hat{k}' \equiv (\theta', \phi')$.

## Appendix B: Mean-field Dipolar Model

The mean-field equations of motion for the dipolar model can be obtained from the master equation by factoring the multi-atom correlations as $\langle \hat{A}_k \hat{B}_j \rangle = \langle \hat{A}_k \rangle \langle \hat{B}_j \rangle$, where $j \neq k$. For an atom $k$, this treatment leads to:

$$ \langle \hat{\sigma}_k^z \rangle = -\Gamma (\langle \hat{\sigma}_k^z \rangle + 1) + 2i \left( \frac{\Omega}{2} e^{i\vec{k}_L \cdot \vec{r}_k} + \sum_{\substack{j=1 \\ j \neq k}}^{N} \mathcal{G}_{kj} \langle \sigma_j^- \rangle \right) \langle \hat{\sigma}_k^+ \rangle $$
$$ \begin{aligned} & -2i \left( \frac{\Omega}{2} e^{-i\vec{k}_L \cdot \vec{r}_k} + \sum_{\substack{j=1 \\ j \neq k}}^{N} \mathcal{G}_{kj}^* \langle \sigma_j^+ \rangle \right) \langle \sigma_k^- \rangle, \quad (B1) \\ & \langle \dot{\sigma}_k^+ \rangle = -\frac{\Gamma}{2} \langle \sigma_k^+ \rangle + i \left( \frac{\Omega}{2} e^{-i\vec{k}_L \cdot \vec{r}_k} + \sum_{\substack{j=1 \\ j \neq k}}^{N} \mathcal{G}_{kj}^* \langle \sigma_j^+ \rangle \right) \langle \sigma_k^z \rangle, \quad (B2) \end{aligned} $$

where $\mathcal{G}_{kj} = \mathcal{R}_{kj} + i\mathcal{I}_{kj} = \mathcal{G}_{jk}$ is the dipolar interaction coefficient and is symmetric with respect to the indices. We can see from the equations above that the dipolar interaction term acts like an effective time-dependent complex drive with a strength that depends on the coherence of other atoms in the array. The elastic and inelastic dipolar coefficients are the real and imaginary parts, respectively, of the free-space electromagnetic Green's function,

$$ \begin{aligned} \mathcal{R}_{kj} &= (3\Gamma/4) \left[ (1 - \cos^2 \theta) \frac{\cos(k_0 r)}{k_0 r} \right. \\ & \qquad \left. + (1 - 3 \cos^2 \theta) \left( -\frac{\sin(k_0 r)}{(k_0 r)^2} - \frac{\cos(k_0 r)}{(k_0 r)^3} \right) \right] \\ \mathcal{I}_{kj} &= (3\Gamma/4) \left[ (1 - \cos^2 \theta) \frac{\sin(k_0 r)}{k_0 r} \right. \\ & \qquad \left. + (1 - 3 \cos^2 \theta) \left( \frac{\cos(k_0 r)}{(k_0 r)^2} - \frac{\sin(k_0 r)}{(k_0 r)^3} \right) \right] \end{aligned} $$

where $r = |\vec{r}_{kj}|$, $\cos \theta = \hat{r}_{kj} \cdot \hat{e}_q$, $q = 0, \pm 1$, $\theta$ is the angle between the polarization of the atomic transition ($\hat{e}_q$, orientation of the transition “dipole”) and the interatomic distance. For the MF dipolar model at the densities considered in this paper, $k_0 \vec{r} \ge 1$ and the physics is dominated by the $1/r$ (far-field) term. The far-field interaction coefficients are somewhat simpler –

$$ \begin{aligned} \mathcal{R}_{kj} &\approx (3\Gamma/4)(1 - \cos^2 \theta) \frac{\cos(k_0 r)}{k_0 r} \\ \mathcal{I}_{kj} &\approx (3\Gamma/4)(1 - \cos^2 \theta) \frac{\sin(k_0 r)}{k_0 r} \end{aligned} $$

and we have $\mathcal{G}_{kj} = \mathcal{R}_{kj} + i\mathcal{I}_{kj} \approx (3\Gamma/4)(1 - \cos^2 \theta)e^{ik_0 r}/(k_0 r)$. By gauging away the phase of the laser, we can rewrite the MF equations as –

$$ \begin{aligned} \langle \dot{\sigma}_k^z \rangle = & -\Gamma (\langle \sigma_k^z \rangle + 1) + 2i \left( \frac{\Omega}{2} + \sum_{\substack{j=1 \\ j \neq k}}^{N} \tilde{\mathcal{G}}_{kj} \langle \tilde{\sigma}_j^- \rangle \right) \langle \tilde{\sigma}_k^+ \rangle \\ & - 2i \left( \frac{\Omega}{2} + \sum_{\substack{j=1 \\ j \neq k}}^{N} \tilde{\mathcal{G}}_{kj}^* \langle \tilde{\sigma}_j^+ \rangle \right) \langle \tilde{\sigma}_k^- \rangle, \quad (B3) \end{aligned} $$

$$ \langle \dot{\tilde{\sigma}}_k^+ \rangle = -\frac{\Gamma}{2} \langle \tilde{\sigma}_k^+ \rangle + i \left( \frac{\Omega}{2} + \sum_{\substack{j=1 \\ j \neq k}}^{N} \tilde{\mathcal{G}}_{kj}^* \langle \tilde{\sigma}_j^+ \rangle \right) \langle \tilde{\sigma}_k^z \rangle, \quad (B4) $$

where we have included the laser phase in the interaction coefficient as $\tilde{\mathcal{G}}_{kj} = \mathcal{G}_{kj} e^{-i\vec{k}_L \cdot \vec{r}_{kj}} \approx (3\Gamma/4)(1 -

$\cos^2 \theta)e^{i(k_0 r - \vec{k}_L \cdot \vec{r}_{kj})}/(k_0 r)$, in the far-field limit. This additional phase leads to constructive interference of interaction terms in the direction of the laser wavevector, i.e., $k_0 r - \vec{k}_L \cdot \vec{r}_{kj} = 0$. For our setup, $\vec{k}_L = 2\pi/\lambda \hat{x} \Rightarrow \vec{k}_L \cdot \vec{r}_{kj} = k_0(x_k - x_j)$ and in the quasi-1D gas, we have $|\vec{r}_{kj}| \approx |x_k - x_j|$ and $\cos^2 \theta \approx |\hat{x} \cdot \hat{e}_q|^2$ ($q = 0, \pm 1$) for the majority of atomic pairs. For a linearly polarized atomic transition ($\hat{e}_0 = \hat{z}$), $|\hat{x} \cdot \hat{e}_q|^2 = 0$ and for a circularly polarized transition ($\hat{e}_+ = -(\hat{x} + i\hat{y})/\sqrt{2}$), $|\hat{x} \cdot \hat{e}_q|^2 = 1/2$. Then, $\tilde{\mathcal{G}}_{kj,k>j} \approx (3\Gamma/8)/(k_0 r)$ and $\tilde{\mathcal{G}}_{kj,k<j} \approx (3\Gamma/8)e^{2ik_0 r}/(k_0 r)$, where we have assumed that the atomic indices are sorted in increasing order of position along the x-axis. Unlike in periodic arrays with special lattice spacing such that $e^{2ik_0 r} = 1$, in a disordered array the $e^{2ik_0 r}$ phases in $\tilde{\mathcal{G}}_{kj,k<j}$ terms would get washed out in comparison to the $\tilde{\mathcal{G}}_{kj,k>j}$ terms. In the MF equations of motion, we have constructive interference of terms $\sim \tilde{\mathcal{G}}_{kj,k>j} \langle \tilde{\sigma}_j^- \rangle \langle \tilde{\sigma}_k^+ \rangle$. This implies that due to the phase matching of dipolar interactions with the laser drive, effects will constructively add up from an atom k further along in the path of the laser absorbing a photon emitted from an atom j earlier in the path of the laser.

## Appendix C: CRF model analysis

Consider the CRF master equation given by:

$$ \frac{\partial \hat{\rho}}{\partial t} = -i\Omega [\hat{S}_x, \hat{\rho}] + \Gamma \left( \hat{S}_- \hat{\rho} \hat{S}_+ - \frac{1}{2} \{ \hat{S}_+ \hat{S}_-, \hat{\rho} \} \right) \quad (C1) $$

Eq. (C1) can be rewritten in the following form:

$$ \frac{\partial \hat{\rho}}{\partial t} = \Gamma \left( \hat{O} \hat{\rho} \hat{O}^{\dagger} - \frac{1}{2} \{ \hat{O}^{\dagger} \hat{O}, \hat{\rho} \} \right), \quad (C2) $$

where

$$ \hat{O} = \hat{S}_- + i\beta \frac{N}{2} = \hat{S}_x - i\hat{S}_y + i\beta \frac{N}{2}, \quad (C3) $$

where $\beta$ is the order parameter and is given by $\beta = \frac{\Omega}{\Omega_c}$ with $\Omega_c = \frac{N\Gamma}{2}$. From the factorization of Eq. (C2) we can infer the properties of the steady-state [76].

For $\beta < 1$ the system is in the superradiant phase, to analyze the stable state within the polarized phase, we employ a Holstein-Primakoff expansion centered around the polarization direction. In mean field theory, the Bloch vector stabilizes at an angle $\theta$, where $\sin \theta = \beta$. This leads us to consider a rotated coordinate system.

$$ \hat{S}_x = \hat{S}'_x, \quad (C4) $$

$$ \hat{S}_y = \hat{S}'_y \cos \theta - \hat{S}'_z \sin \theta, \quad (C5) $$

$$ \hat{S}_z = \hat{S}'_z \cos \theta + \hat{S}'_y \sin \theta, \quad (C6) $$

such that the Bloch vector is aligned along $-z'$. We then do a Holstein-Primakoff (HP) expansion about this
direction, with the (lowest order) replacements $\hat{S}'_z \approx -\frac{N}{2}$, $\hat{S}'_x \approx \hat{x}\sqrt{\frac{N}{2}}$, $\hat{S}'_y \approx -\hat{p}\sqrt{\frac{N}{2}}$ with $\hat{x} = \frac{\hat{a}+\hat{a}^\dagger}{\sqrt{2}}$ and $\hat{p} = \frac{\hat{a}-\hat{a}^\dagger}{\sqrt{2i}}$ which satisfies $[\hat{x}, \hat{p}] = i$. The operator $\hat{O}$ in Eq. (C3) in the HP expansion will have the form:

$$ \hat{O} = \hat{S}_x - i\hat{S}_y + i\beta \frac{N}{2} \approx \hat{x} + i\hat{p} \cos \theta. \quad (C7) $$

For $\beta < 1$ the state is in the superradiant phase and will remain coherent, this implies that the density matrix of the steady-state satisfies $\hat{\rho} = |0_D\rangle\langle 0_D|$ and $\hat{O}|0_D\rangle = 0$, that is, $\hat{O}$ corresponds to an annihilation operator with unique dark state $|0_D\rangle$. To normalize properly $\hat{O}$ we recall that the annihilation operator $\hat{a}_D = k\hat{O}$ satisfies the commutation relation $[\hat{a}_D, \hat{a}_D^\dagger] = 1$, replacing we find that $k^2[\hat{x} + i\hat{p}\cos\theta, \hat{x} - i\hat{p}\cos\theta] = k^2\cos\theta(-2i^2) = 2k^2\cos\theta = 1$, so $k = 1/\sqrt{2\cos\theta}$ and the annihilation operator reads as $\hat{a}_D = (\hat{x} + i\hat{p}\cos\theta)/\sqrt{2\cos\theta}$. In these terms we can define $\hat{x}_D = \frac{\hat{a}_D+\hat{a}_D^\dagger}{\sqrt{2}}$ and $\hat{p}_D = \frac{\hat{a}_D-\hat{a}_D^\dagger}{\sqrt{2i}}$, so we can find that $\hat{x} = \sqrt{\cos\theta}\hat{x}_D$ and $\hat{p} = \hat{p}_D/\sqrt{\cos\theta}$. These operators will allow to calculate expected values in the superradiant regime.

Our aim is to find expressions that represent the behavior of expected value of operators of interest at the limit of a large number of particles. Replacing the operators of Eq. (C6) expanded in HP we find that:

$$ \hat{S}_x \approx \hat{x}_D \sqrt{\frac{N}{2} \cos \theta}, \quad (C8) $$

$$ \hat{S}_y \approx \frac{N}{2} \sin \theta - \hat{p}_D \sqrt{\frac{N}{2} \cos \theta}, \quad (C9) $$

$$ \hat{S}_z \approx -\frac{N}{2} \cos \theta - \hat{p}_D \sqrt{\frac{N}{2 \cos \theta}} \sin \theta. \quad (C10) $$

From the mean field analysis we can obtain that $\langle \hat{S}_x \rangle = 0$, $\langle \hat{S}_y \rangle = \frac{N}{2} \sin \theta$ and $\langle \hat{S}_z \rangle = -\frac{N}{2} \cos \theta$ [76]. For the expected value $\langle \hat{S}_+ \hat{S}_- \rangle$ we can calculate using the HP expansion:

$$ \begin{aligned} \hat{S}_+ = \hat{S}_x + i\hat{S}_y &\approx \hat{x}_D \sqrt{\frac{N}{2} \cos \theta} + i \left( \frac{N}{2} \sin \theta - \hat{p}_D \sqrt{\frac{N}{2} \cos \theta} \right) \\ &= i \frac{N}{2} \sin \theta + \sqrt{\frac{N}{2} \cos \theta} (\hat{x}_D - i\hat{p}_D) \quad (C11) \\ &= i \frac{N}{2} \sin \theta + \sqrt{N \cos \theta} \hat{a}_D^\dagger. \end{aligned} $$

So, $\hat{S}_+\hat{S}_-$ expanded as a function of the HP operators and its expected value is:

$$ \begin{aligned} \hat{S}_+\hat{S}_- &\approx \left(i\frac{N}{2}\sin\theta + \sqrt{N\cos\theta}\hat{a}_D^\dagger\right)\left(-i\frac{N}{2}\sin\theta + \sqrt{N\cos\theta}\hat{a}_D\right), \\ \langle \hat{S}_+\hat{S}_- \rangle &= \langle 0_D | \hat{S}_+\hat{S}_- | 0_D \rangle \approx i\frac{N}{2}\sin\theta\left(-i\frac{N}{2}\sin\theta\right) \\ &= \left(\frac{N}{2}\sin\theta\right)^2. \end{aligned} \quad (C12) $$

We used the fact that $\langle \hat{a}_D^\dagger \rangle = \langle \hat{a}_D \rangle = \langle \hat{a}_D^\dagger \hat{a}_D \rangle = 0$. This formula is valid for $\theta < \pi/2$ where the HP approximation holds.

Now we estimate the expected value $\langle \hat{S}_+\hat{S}_+\hat{S}_-\hat{S}_- \rangle$ as:

$$ \begin{aligned} \langle \hat{S}_+ \hat{S}_+ \hat{S}_- \hat{S}_- \rangle &= \langle 0_D | \hat{S}_+ \hat{S}_+ \hat{S}_- \hat{S}_- | 0_D \rangle, \\ \hat{S}_- | 0_D \rangle &\approx -i \frac{N}{2} \sin \theta | 0_D \rangle, \\ \langle \hat{S}_+ \hat{S}_+ \hat{S}_- \hat{S}_- \rangle &\approx \left(i \frac{N}{2} \sin \theta\right) \left(-i \frac{N}{2} \sin \theta\right) \langle 0_D | \hat{S}_+ \hat{S}_- | 0_D \rangle \\ &= \left(\frac{N}{2} \sin \theta\right)^4. \end{aligned} \quad (C13) $$

In this case we find that $g_2(0) = \frac{\langle \hat{S}_+ \hat{S}_+ \hat{S}_- \hat{S}_- \rangle}{\langle \hat{S}_+ \hat{S}_- \rangle^2}$ at the HP limit is $g_2(0) \approx \frac{\left(\frac{N}{2} \sin \theta\right)^4}{\left(\left(\frac{N}{2} \sin \theta\right)^2\right)^2} = 1$ for the superradiant phase. Regarding the HP approximation, this tends to improve at the limit of large number of particles such that it converges to the exact result [76], so we can expect that for $\beta < 1$ the previous result is valid at $N \to \infty$ and in these conditions $g_2(0) = 1$.

For $\beta > 1$ we are in the normal phase, in that case the steady-state is highly mixed. In general, the steady-state can be written formally as:

$$ \hat{\rho}_{ss} = \mathcal{N} \frac{1}{\frac{\hat{S}_-}{N/2} + i\beta} \cdot \frac{1}{\frac{\hat{S}_+}{N/2} - i\beta} \quad (C14) $$

where $N$ is a normalization constant and $\beta = 2\Omega/(\Gamma N)$. In the normal phase ($\beta > 1$), the properties of the steady-state are derived through semiclassical analysis. We parameterize the phase space using angles $\theta$ and $\phi$ as follows $(\hat{S}_x, \hat{S}_y, \hat{S}_z) \to \frac{N}{2}(\sin\theta \cos\phi, \sin\theta \sin\phi, \cos\theta)$ so $\hat{S}_+ \to \frac{N}{2}\sin\theta e^{i\phi}$, and operator traces are replaced by integrals over the sphere, with the measure $\frac{N}{4\pi}\sin\theta d\theta d\phi$. The normalization constant $N$ is determined through this process, and a justification for the semiclassical approximation is provided in [64]. To obtain $\mathcal{N}$ we recall the normalization condition Tr($\hat{\rho}_{ss}$) = 1, so we require to satisfy:

$$ 1 = \mathcal{N} \int_0^\pi \int_0^{2\pi} \frac{1}{(\sin\theta e^{-i\phi} + i\beta)(\sin\theta e^{i\phi} - i\beta)} \frac{N}{4\pi} \sin\theta d\theta d\phi. \quad (C15) $$

This integral can be done using calculus of residues around the unit circle $z = e^{i\phi}$ in the counterclockwise direction and we must be aware that $\sin\theta/\beta < 1$ in the normal phase, which will be relevant for determining the suitable residues inside the integration curve. In particular, using this substitution in (C15) leads to:

$$ 1 = -\mathcal{N} \frac{N}{4\pi\beta} \int_0^\pi \left( \oint_{|z|=1} \frac{dz}{\left(z - i \frac{\sin\theta}{\beta}\right) \left(z - i \frac{\beta}{\sin\theta}\right)} \right) d\theta. \quad (C16) $$
In (C16) we can identify that the residue $z_{\text{res},1} = i \frac{\sin \theta}{\beta}$ belongs to the interior of the integration curve whereas $z_{\text{res},2} = i \frac{\beta}{\sin \theta}$ does not belong because for all $\theta$ we have $\sin \theta < \beta$ in the normal phase. By the Cauchy's integral formula we deduce that:

$$1 = \mathcal{N} \frac{N}{2} \int_{0}^{\pi} \frac{\sin \theta}{\beta^2 - \sin^2 \theta} d\theta. \qquad (C17)$$

From here we find that $\mathcal{N} = \frac{\sqrt{\beta^2 - 1}}{N \arctan((\sqrt{\beta^2 - 1})^{-1})}$.

The estimation of the expected values can be performed using a similar integration process such that $\hat{O} = \hat{O}(\theta, \phi)$ and $\eta = \sqrt{\beta^2 - 1}$ as:

$$\langle \hat{O} \rangle \approx \frac{\eta}{4\pi \arctan(1/\eta)} \int_0^\pi \int_0^{2\pi} \frac{\hat{O}(\theta, \phi) \sin \theta d\phi d\theta}{(\sin \theta e^{-i\phi} + i\beta)(\sin \theta e^{i\phi} - i\beta)}. \qquad (C18)$$

For example in the case of $\langle \hat{S}_x^2 \rangle$ we can observe that $\hat{S}_x^2 \rightarrow (\frac{N}{2})^2 \sin^2 \theta \cos^2 \phi$ so in the normal phase we have that:

$$\langle \hat{S}_x^2 \rangle \approx \frac{\eta}{4\pi \arctan(1/\eta)} \int_0^\pi \int_0^{2\pi} \frac{(\frac{N}{2})^2 \sin^3 \theta \cos^2 \phi d\phi d\theta}{(\sin \theta e^{-i\phi} + i\beta)(\sin \theta e^{i\phi} - i\beta)} \\ = \left(\frac{N}{2}\right)^2 \frac{\eta}{3\beta^2 \arctan(1/\eta)}. \qquad (C19)$$

In particular for calculating $g_2(0)$ we require to estimate $\langle \hat{S}_+ \hat{S}_+ \hat{S}_- \hat{S}_- \rangle$. Given that we are estimating this quantity at the limit of large $N$ then we expand the operator in the semiclassical approximation at the leading order of $N$. This imply that $\hat{S}_+ \hat{S}_+ \hat{S}_- \hat{S}_- = \hat{S}_+^2 \hat{S}_-^2 \rightarrow (\frac{N}{2} \sin \theta e^{i\phi})^2 (\frac{N}{2} \sin \theta e^{-i\phi})^2 + O(N^3) \approx (\frac{N}{2})^4 \sin^4 \theta$ for large $N$. In this case the mathematical expression for finding $g_2(0)$ at the leading order is given by:

$$m_n = \frac{\eta}{4\pi \arctan(1/\eta)} \int_0^\pi \int_0^{2\pi} \frac{(\frac{N}{2})^n \sin^{n+1} \theta d\phi d\theta}{(\sin \theta e^{-i\phi} + i\beta)(\sin \theta e^{i\phi} - i\beta)}, \qquad (C20)$$

$$g_2(0) \approx \frac{m_4}{m_2^2} \\ \approx \left( \frac{\beta^4}{\eta} \arctan\left(\frac{1}{\eta}\right) - \beta^2 - \frac{2}{3} \right) \left( \beta^2 \arctan\left(\frac{1}{\eta}\right) - \eta \right)^2 \\ \times \eta \arctan\left(\frac{1}{\eta}\right). \qquad (C21)$$

The remanent integrals used to estimate the quantities of interest in Eqs. 13-16 are the following:

$$\langle \hat{S}_y \rangle \approx \frac{\eta}{4\pi \arctan\left(\frac{1}{\eta}\right)} \int_0^\pi \int_0^{2\pi} \frac{\frac{N}{2} \sin^2 \theta \sin \phi d\phi d\theta}{(\sin \theta e^{-i\phi} + i\beta)(\sin \theta e^{i\phi} - i\beta)} \\ = \frac{N}{2} \left( \beta - \frac{\eta}{\beta \arctan(1/\eta)} \right), \qquad (C22)$$

$$\langle \hat{S}_+ \hat{S}_- \rangle \approx \langle \hat{S}_x^2 + \hat{S}_y^2 \rangle \\ = \frac{\eta}{4\pi \arctan\left(\frac{1}{\eta}\right)} \int_0^\pi \int_0^{2\pi} \frac{(\frac{N}{2})^2 \sin^3 \theta d\phi d\theta}{(\sin \theta e^{-i\phi} + i\beta)(\sin \theta e^{i\phi} - i\beta)} \\ = \frac{N^2}{4} \left( \beta^2 - \frac{\eta}{\arctan\left(\frac{1}{\eta}\right)} \right). \qquad (C23)$$

## Appendix D: Modified CRF model

The main steady-state features of the dipolar model can be recovered from a much simpler model, which can be solved analytically, the Cooperative Resonance Fluorescence model with one-axis twisting, dephasing, and spontaneous emission. Even though this model has a bistable solution in the superradiant phase. For our chosen initial condition $|g\rangle^{\otimes N}$, we find that the system reaches a steady state very similar to the dipolar model. To see this, first we obtain the mean field steady state, which can be derived from the equations of motion of this model, given by:

$$\langle \hat{S}_z \rangle = -\Omega \langle S_y \rangle - \gamma_s \left( \frac{N}{2} + \langle S_z \rangle \right) - \Gamma_D (\langle S_x \rangle^2 + \langle S_y \rangle^2), \\ \langle \hat{S}_y \rangle = \Omega \langle S_z \rangle + \Gamma_D \langle S_z \rangle \langle S_y \rangle - \frac{(\gamma_s + 4\gamma_d)}{2} \langle S_y \rangle \\ + 2\chi \langle S_x \rangle \langle S_z \rangle, \\ \langle \hat{S}_x \rangle = \Gamma_D \langle S_z \rangle \langle S_x \rangle - \frac{(\gamma_s + 4\gamma_d)}{2} \langle S_x \rangle - 2\chi \langle S_y \rangle \langle S_z \rangle. \qquad (D1)$$

where, as usual, we have used $\langle \sigma_j^\alpha \sigma_k^\beta \rangle = \langle \sigma_j^\alpha \rangle \langle \sigma_k^\beta \rangle$ for $j \neq k$ and $\langle \sigma_k^\alpha \sigma_k^\beta \rangle = 2i\epsilon_{\alpha\beta\gamma}(\langle \sigma_k^\gamma \rangle)$. For the cooperative emission and one-axis twisting (OAT) terms, we have $\langle \sigma_j^\alpha \sigma_k^\beta \rangle = \langle \sigma_j^\alpha \rangle \langle \sigma_k^\beta \rangle$ for all $j, k$ [61]. We can see from the equations above that for $\chi = 0$, $\langle \hat{S}_x \rangle = \Gamma_D \langle S_z \rangle \langle S_x \rangle - \frac{(\gamma_s + 4\gamma_d)}{2} \langle S_x \rangle \Rightarrow \langle \hat{S}_x(t) \rangle = \langle S_x(0) \rangle \int dt \exp \left[ \Gamma_D \langle S_z(t) \rangle - \frac{(\gamma_s + 4\gamma_d)}{2} \right]$. Hence, for an initial state with $\langle \hat{S}_x(0) \rangle = 0$ and $\chi = 0$, the value of $\langle \hat{S}_x(t) \rangle$ remains zero for the full dynamics of the system. This is because the laser drive is along $\hat{S}_x$, so it commutes with and does not alter $\hat{S}_x$. The other terms are dissipative and they destroy the coherence, $\hat{S}_x$. However, the OAT ($\chi \neq 0$) term causes shearing of the collective Bloch vector about the z-axis, which leads to $\hat{S}_x \neq 0$.

For simplicity, we will here on use the notation $z \equiv \langle S_z \rangle/(N/2)$ and we define new variables $(r, \phi)$ to express the collective coherence as $\langle S_+ \rangle/N = re^{i\phi}$, such that $r = \sqrt{\langle S_x \rangle^2 + \langle S_y \rangle^2}/N$ is the contrast and $\phi$ is the phase in the XY-plane of the collective Bloch sphere [77]. Then, the MF equations can be expressed in terms of $(z, r, \phi)$ as

$$\dot{z} = -2\Omega r \sin \phi - \gamma_s(1+z) - 2\Gamma_D N r^2, \qquad (D2)$$
$$ \dot{r} = \frac{z}{2}(\Omega \sin \phi + \Gamma_D Nr) - \frac{(\gamma_s + 4\gamma_d)}{2}r, \quad (D3) $$

$$ \dot{\phi} = \frac{z}{2r} \Omega \cos \phi + \frac{\chi N z}{2}. \quad (D4) $$

The steady state can be obtained by setting $\dot{z} = \dot{r} = \dot{\phi} = 0$. Then, we multiply $4r/z$ to Eq. (D3) to get

$$ -2\Omega r \sin \phi - 2\Gamma_D N r^2 = -2(\gamma_s + 4\gamma_d) \frac{r^2}{z} \quad (D5) $$

We plug Eq. (D5) into Eq. (D2) to get the steady state as

$$ \begin{aligned} & -2(\gamma_s + 4\gamma_d)\frac{r^2}{z} - \gamma_s(1+z) = 0 \\ \Rightarrow \quad & z = -\frac{1}{2} \pm \frac{1}{2}\sqrt{1 - 8\left(1 + \frac{4\gamma_d}{\gamma_s}\right)r^2}. \end{aligned} \quad (D6) $$

There is a phase-transition when the term inside the square-root reaches zero and the critical value $r_c$ can be obtained as

$$ r_c = \frac{1}{2\sqrt{2(1 + 4\gamma_d/\gamma_s)}}, \quad z_c = -1/2. \quad (D7) $$

Now, we want to obtain the critical driving strength. We can square and rewrite Eq. (D2) and Eq. (D4) as

$$ 4\Omega^2 r^2 \sin^2 \phi = [2\Gamma_D N r^2 + \gamma_s (1+z)^2]^2, \quad (D8) $$

$$ 4\Omega^2 r^2 \cos^2 \phi = 4\chi^2 N^2 r^4. \quad (D9) $$

We add Eq. (D8) and Eq. (D9), and keep upto leading order terms in $N$ to obtain

$$ \Omega = \pm \sqrt{(\chi^2 + \Gamma_D^2) N^2 r^2 + \Gamma_D \gamma_s N (1+z)}. \quad (D10) $$

The sign of $\Omega$ determines the direction of the drive along $\pm \hat{S}_x$. For $\chi = 0$, we get back the same MF equations when we flip the signs of $\Omega$ and $\phi$ simultaneously. Thus, the two drive directions $\pm \hat{S}_x$ are physically equivalent. However, a non-zero $\chi$ breaks this symmetry. Nevertheless, we get qualitatively similar steady states for the negative and positive values of $\Omega$, so we will only consider the positive value for illustration here. Now, we plug $z_c$ and $r_c$ from Eq. (D7) into Eq. (D10) to get the critical drive strength as

$$ \Omega_C^{\text{Mod-CRF}} = \frac{\Gamma_D}{2\sqrt{2}} \sqrt{\frac{1 + (\chi/\gamma_d)^2}{1 + (4\gamma_d/\gamma_s)} N^2 + \frac{4\gamma_s}{\Gamma_D} N}. \quad (D11) $$

If we consider the case without dephasing and OAT, i.e., $\gamma_d = 0$ and $\chi = 0$, we recover $\Omega_C \approx \Gamma_D N/(2\sqrt{2})$ [51, 61, 70]. Thus, the inclusion of spontaneous emission in the CRF model only modifies the critical driving strength by a constant factor without altering its $N$-scaling. As discussed in the main text, accounting for the frequency shifts from dipolar interactions is important. For that we set the dephasing rate to be $4\gamma_d/\gamma_s = c_d N$

(c_d=constant). Then, we take the large-$N$ limit, to obtain the critical point as

$$ \Omega_C \approx \frac{\Gamma_D}{2\sqrt{2}} \sqrt{N} \sqrt{\frac{1 + (\chi/\Gamma_D)^2}{c_d} + \frac{4\gamma_s}{\Gamma_D}}. \quad (D12) $$

<figure><img src="image_15.png" /><figcaption>FIG. 15. (a) Dynamics of the excitation fraction $n_e$ for fixed $\Omega/\Gamma \equiv \Omega/\gamma_s = 5.1$ for the CRF model with dephasing ($c_d = 0.002$), OAT, and spontaneous emission.</figcaption></figure>

In Fig. 15, we show the dynamics of the excitation fraction $n_e(t)$, using the mean-field equations for the modified CRF model (Eq. (D1)), with $c_d = 0.002$. This specific value of $c_d$ is not special and we find qualitatively similar results for other values as long as $c_d \ll 1$. We have taken the continuum limit to obtain $\chi = \frac{1}{N} \int d^3r \text{Re}\tilde{G}(\vec{r})\rho(\vec{r}) \approx 0.003$ and $\Gamma_D = \frac{1}{N} \int d^3r \text{Im}\tilde{G}(\vec{r})\rho(\vec{r}) \approx 0.002$, where $\rho(\vec{r})$ (Eq. (A16)) is the pencil-shape distribution of the cloud. Similar to the MF dipolar model, we find that for fixed $\Omega/\gamma_s$, increasing $N$, leads to larger dephasing and to a suppression of $n_e$.

## Appendix E: Contribution of laser light to the intensity

We consider a coordinate system where the centre of the atomic cloud is the reference point for the origin, $C = (0, 0, 0)$. All the other vectors/positions are defined with respect to this point. Then, the position of an atom $j$ (with respect to $C$) is given as $\vec{r}_j$. The position of the laser source is $\vec{d}_L = (-d, 0, 0)$. The position of a point on the detector is $\vec{R} = R(\cos \theta, \sin \theta \cos \phi, \sin \theta \sin \phi)$, where $\theta$ is the angle from the $x$-axis and $\phi$ is the azimuthal angle in the $y-z$ plane.

The total electric field and the intensity at a point $\vec{R}$ on the detector are given as

$$ \vec{E}_{\text{total}}^{+}(\vec{R}) = \vec{E}_{\text{drive}}^{+}(\vec{R}) + \langle \vec{E}_{\text{atoms}}^{+}(\vec{R}) \rangle, \quad (E1) $$

$$ I_{\text{total}}(\vec{R}) = \vec{E}_{\text{total}}^{-}(\vec{R}) \cdot \vec{E}_{\text{total}}^{+}(\vec{R}), \quad (E2) $$

where the electric field of the driving laser is

$$ \vec{E}_{\text{drive}}^{+}(\vec{R}) = E(\vec{R}) \hat{e}_{L} = (\hbar\Omega/|\vec{d}|) \hat{e}_{L} e^{i\vec{k}_{L} \cdot \vec{R}'} \quad (E3) $$

where $\vec{R}' = \vec{R} - \vec{d}_L = (R \cos \theta + d, R \sin \theta \cos \phi, R \sin \theta \sin \phi)$ is the distance vector
from the laser source to the detector, the laser polariza-
tion is $\hat{e}_L = \hat{e}_y$, and its wavevector is $\vec{k}_L = (2\pi/\lambda)\hat{x}$ for
our system. The electric field of the atomic dipoles is

$$
\vec{E}_{\text{atoms}}^{+}(\vec{R}) = \frac{\mu_0 \omega_0^3}{3\pi \Gamma c} \sum_j \mathbf{G}(\vec{R} - \vec{r}_j) \cdot \hat{e}_+ |\vec{d}| \sigma_j^{-} \quad (\text{E4})
$$

where $\vec{R} - \vec{r}_j$ is the distance vector from atom $j$ to the
detector.

Then, considering that the detector is very far from
the atomic cloud, it is safe to assume that $R \gg |\vec{r}_j|$ and
the electric field of the atoms is dominated by the far-
field limit terms ($R \gg \lambda \Rightarrow \mathbf{G} \propto 1/(k_0 R)$). Under these
assumptions, the total intensity (laser + atomic emission)
is obtained as

$$
I_{\text{total}}(\vec{R}) = \left( \frac{\hbar \Gamma}{k_0 |\vec{R}| |\vec{d}|} \right)^2 \left[ \frac{\Omega^2}{\Gamma^2} \right. \\
\left. \qquad + \frac{3\Omega}{2\sqrt{2\Gamma}} \sum_j \left( \text{Im}\langle\tilde{\sigma}_j^-\rangle (1 - \sin^2 \theta \cos^2 \phi) \right. \right. \\
\left. \left. \qquad + \text{Re}\langle\tilde{\sigma}_j^-\rangle \frac{\sin(2\theta)}{2} \cos \phi \right) \right. \\
\left. \qquad + \frac{9}{16} \sum_{j,k} \frac{(1 + \sin^2 \theta \cos^2 \phi)}{2} \langle\sigma_j^+ \sigma_k^-\rangle e^{i\vec{k} \cdot \vec{r}_{jk}} \right], \tag{E5}
$$

where $\vec{k} = (2\pi/\lambda)\hat{R}$, and we have defined

$$
\begin{align*}
\tilde{\sigma}_j^- &= \sigma_j^- \exp(i\vec{k} \cdot (\vec{R} - \vec{r}_j) - i\vec{k}_L \cdot \vec{R}') \\
&= \sigma_j^- \exp(i\vec{k} \cdot (\vec{R} - \vec{r}_j) - i\vec{k}_L \cdot (\vec{R} - \vec{d}_L)) \\
&= \sigma_j^- \exp(-i(\vec{k} \cdot \vec{r}_j - \vec{k}_L \cdot \vec{d}_L) - i(\vec{k}_L - \vec{k}) \cdot \vec{R}) \quad (\text{E6})
\end{align*}
$$

The spiral basis was defined to absorb the laser-induced
phase in the Hamiltonian. Keeping the same definition,
in the current coordinate system, the spiral basis is de-
fined as $\tilde{\sigma}_j^- = i\sigma_j^- e^{-i\vec{k}_L \cdot (\vec{r}_j - \vec{d}_L)}$, where the extra phase
factor of $i$ comes from the overlap of the circular polar-
ization of the transition with the laser polarization. So
we can rewrite Eq. (E6) in terms of the spiral basis op-
erator as

$$
\tilde{\sigma}_j^- = -i\tilde{\sigma}_j^- \exp(-i(\vec{k} - \vec{k}_L) \cdot (\vec{r}_j - \vec{R})) \quad (\text{E7})
$$

In the forward direction $\vec{k} = \vec{k}_L$, the phases cancel, and
we have $\tilde{\sigma}_j^- = -i\tilde{\sigma}_j^-$. For simplicity, we assume $(\vec{k} -
\vec{k}_L) \cdot \vec{R} \approx 2\pi$ as it is difficult to determine its exact value
experimentally. Since we are averaging over $\vec{k}$ to obtain
the final numerical values, we expect that changing the
value of $|\vec{R}|$ will not change the result significantly.

For simplicity, when including the contribution of the
probe light in estimating $\tilde{g}_2(0)$, we have neglected the
term accounting for the interference between the probe
and the atomic electric fields.

**Appendix F: Numerical Approximations**

1. Mean-field methods

We use mean-field theory to obtain the dynamics and
steady-state of the dipolar model. Under mean-field the-
ory, the connected part of the atom-atom correlations is
assumed to be negligible. In this case all correlations can
be expressed as products of classical expectation values,
as shown in Fig. 16. To benchmark the results obtained
from mean-field theory, we use the cumulant method,
which goes beyond the mean-field approximation by in-
cluding two-point connected correlations (see Fig. 16 for
an intuitive picture). The MF equations of motion for
the dipolar model are given in Appendix B.

<table>
  <thead>
    <tr>
      <th>Mean-Field (MF)</th>
      <th>Cumulant</th>
      <th>MACE-MF</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <td>⟨σ<sup>2</sup><sub>j</sub>σ<sup>2</sup><sub>i</sub>⟩ ≈ ⟨σ<sup>2</sup><sub>j</sub>σ<sup>2</sup><sub>i</sub>⟩⟨σ<sup>2</sup><sub>i</sub>⟩ + ⟨σ<sup>2</sup><sub>i</sub>σ<sup>2</sup><sub>j</sub>⟩⟨σ<sup>2</sup><sub>j</sub>⟩ + ⟨σ<sup>2</sup><sub>i</sub>σ<sup>2</sup><sub>i</sub>⟩⟨σ<sup>2</sup><sub>j</sub>⟩<br>-2⟨σ<sup>2</sup><sub>j</sub>⟩⟨σ<sup>2</sup><sub>i</sub>⟩⟨σ<sup>2</sup><sub>i</sub>⟩</td>
      <td>≈</td>
      <td>H<sub>MACE-MF</sub> = H<sub>MACE</sub> + &Sigma;<sub>i&isin;C<sub>i</sub></sub> &lt;[σ&#773;<sup>2</sup>]<sub>i&isin;C<sub>i</sub></sub>&gt;, σ&#773;<sub>i</sub></td>
    </tr>
    <tr>
      <td>σ&#773;<sub>j</sub>σ&#773;<sub>i</sub> ≈ σ<sup>2</sup><sub>j</sub>σ<sup>2</sup><sub>i</sub></td>
      <td>≈</td>
      <td>C<sub>i</sub>σ&#773;<sub>i</sub> ≡ C<sub>i</sub>σ&#773;<sub>i</sub></td>
    </tr>
  </tbody>
</table>

FIG. 16. Numerical approximations. Mean-Field (Left): Uses non-overlapping clusters, factoring many-body correlations into single-particle values. Cumulant Approximation (Center): At second order, factors 3-body correlations into two-body and single-body correlations. MACE-MF (Right): Treats a cluster and its surrounding particles as a single cluster influenced by an external time-dependent magnetic field.

2. MACE-MF Method

We introduce an improved cluster technique which we
call MACE-MF (Moving Average Cluster Expansion +
Mean-Field). The method combines dynamics consider-
ing both local clusters and mean-field interactions with
external particles. This method is based on a cluster ap-
proximation, called MACE [78], developed to deal with
dipolar interactions in dilute arrays in the past. It ex-
actly solves for the dynamics inside clusters which are
not rigid but instead adjusted depending of the observ-
able in consideration. For single point observables σ̃k,
it chooses a cluster that contains the strongest coupled
particles to the k-particle in consideration. It disregards
the influence of particles outside the cluster. In fact, at
a cluster size of 1, MACE essentially reduces to single-
particle dynamics.

To account for interaction effects with the particles
outside a given cluster, we include them via a mean-
field approach. In this way the dipole-dipole interaction
with the outside particles is mapped as an effective time-
dependent magnetic field acting upon the atoms in the
cluster as shown in Fig. 16. So by combining the exact
dynamics within the cluster with the mean-field contri-
butions we aim to improve the description of the non-equilibrium dynamics of the system.

The clusters are chosen and adapted to optimize the single or two particle observables in consideration by including in the cluster the particles that feature the strongest coupling constants to the particle or particles the observable is acting on.

We generalize it to also include the effect of the external particles to each cluster via the mean-field approximation. This means that the interaction of a given cluster with the surrounding particles is reduced to an effective magnetic field acting on the particles in the cluster. We denote this approach as MACE-Mean-Field (MACEMF). Over all, the method aims to include the best of both the MF and the MACE solutions and improve upon them. In the limit of one particle per cluster the MACEMF reduces to the MF approximation, and if we neglect the MF couplings between clusters then the MACEMF reduces to the MACE approximation.

For the j particle in the array, the effective Hamiltonian and Lindblandian associated to the cluster of the particle j read as:

$$H_{\text{CMF}}^{j} = -\Omega \sum_{i \in C_j} \left( e^{i \vec{k}_L \cdot \vec{r}_i} (e_L \cdot e_0^*) \sigma_i^+ + \text{h.c.} \right) \quad \text{(F1a)}$$

$$- \sum_{i,k \in C_j} \mathcal{R}_{ik} \sigma_i^+ \sigma_k^- - \sum_{\substack{i \in C_j, \\ k \notin C_j}} \mathcal{R}_{ik} (\sigma_i^+ \langle \sigma_k^- \rangle + \text{h.c.}),$$

$$L[\rho] = \sum_{i,k} \mathcal{I}_{ik} (2\sigma_k^- \rho \sigma_i^+ - \{\sigma_i^+ \sigma_k^- , \rho\}) \quad \text{(F1b)}$$

$$+ \sum_{\substack{i \in C_j, \\ k \notin C_j}} \mathcal{I}_{ik} (\langle \sigma_k^- \rangle [\rho_j, \sigma_i^+] - \langle \sigma_i^+ \rangle [\rho_j, \sigma_i^-]).$$

where $C_j$ denotes the set of particles that are contained in the j-particle cluster and $\rho_j$ is the density matrix that describes this cluster.

The cluster $C_j$ is selected by choosing the particles that have the strongest interactions with the j-particle (largest $|G_{jk}|$).

Identifying the second part of the lindblandian F1b as an effective element of the hamiltonian we can rewrite:

$$H_{\text{CMF}}^{j} = -\sum_{i \in C_j} \left( \left( \Omega e^{i \vec{k}_L \cdot \vec{r}_i} (e_L \cdot e_0^*) + \sum_{k \notin C_j} \mathcal{R}_{ik} \langle \sigma_k^- \rangle \right) \sigma_i^+ + \text{h.c.} \right) \quad \text{(F2a)}$$

$$- \sum_{i,k \in C_j} \mathcal{R}_{ik} \sigma_i^+ \sigma_k^- - i \sum_{\substack{i \in C_j, \\ k \notin C_j}} \mathcal{I}_{ik} (\langle \sigma_k^- \rangle \sigma_i^+ - \langle \sigma_i^+ \rangle \sigma_k^-),$$

$$L[\rho] = \sum_{i,k \in C_j} \mathcal{I}_{ik} (2\sigma_k^- \rho_j \sigma_i^+ - \{\sigma_i^+ \sigma_k^- , \rho_j\}). \quad \text{(F2b)}$$

Rearranging further the hamiltonian F2a we obtain that:

$$H_{\text{CMF}}^{j} = -\sum_{i \in C_j} \left( \left( \Omega e^{i \vec{k}_L \cdot \vec{r}_i} (e_L \cdot e_0^*) + \sum_{k \notin C_j} \mathcal{G}_{ik} \langle \sigma_k^- \rangle \right) \sigma_i^+ + \text{h.c.} \right) \quad \text{(F3)}$$

$$- \sum_{i,k \in C_j} \mathcal{R}_{ik} \sigma_i^+ \sigma_k^-$$

For estimating second-order correlations of the form $\langle \sigma_j^\alpha \sigma_k^\beta \rangle$, where $\alpha, \beta = x, y, z$ and $j \neq k$, we consider the following cases (recall that by default $i \in C_i$ for all $i$):

$$\langle \sigma_j^\alpha \sigma_k^\beta \rangle \approx \begin{cases} \frac{\langle \sigma_j^\alpha \sigma_k^\beta \rangle_{C_j} + \langle \sigma_j^\alpha \sigma_k^\beta \rangle_{C_k}}{2}, & \text{if } k \in C_j \text{ and } j \in C_k, \\ \frac{\langle \sigma_j^\alpha \sigma_k^\beta \rangle_{C_j} + \langle \sigma_j^\alpha \rangle \langle \sigma_k^\beta \rangle}{\langle \sigma_j^\alpha \rangle \langle \sigma_k^\beta \rangle}, & \text{if } k \in C_j \text{ and } j \notin C_k, \\ \frac{\langle \sigma_j^\alpha \rangle \langle \sigma_k^\beta \rangle + \langle \sigma_j^\alpha \sigma_k^\beta \rangle_{C_k}}{2}, & \text{if } k \notin C_j \text{ and } j \in C_k, \\ \frac{\langle \sigma_j^\alpha \rangle \langle \sigma_k^\beta \rangle}{\langle \sigma_j^\alpha \rangle \langle \sigma_k^\beta \rangle}, & \text{if } k \notin C_j \text{ and } j \notin C_k \end{cases} \quad \text{(F4)}$$

Here $\langle \mathcal{O} \rangle_{C_j} = \text{Tr}(\rho_j \mathcal{O})$, which represents the expected value of the operator $\mathcal{O}$ respect the cluster $C_j$. Similar decomposition can be used to estimate higher-order correlations.

## Appendix G: Validity of Mean-field theory

In this section, we benchmark the steady-state properties of the MF dipolar model by comparing with the results of beyond mean-field methods such as MACE-MF, exact diagonalization (ED), and the cumulant approximation. We average the data over ~ 10 realizations for the MF, MACE-MF, and ED results.

<figure><img src="image_16.png" /><figcaption>FIG. 17. Comparing MACE-MF C6, i.e., a cluster-size of 6 atoms, with MF in the steady-state for $N = 2000$</figcaption></figure>

The MACE-MF includes short-range beyond-MF correlations by treating them exactly (using ED) within a given cluster size. In Fig. 17, we show that the MF steady-state atomic inversion and forward intensity for $N = 2000$ agree extremely well with those obtained from MACE-MF with cluster size 6, across a range of $\Omega/\Gamma$ spanning the different dynamical phases. We find that the MACE-MF results do not change appreciably as we increase cluster size beyond 6. The agreement between MF and MACE-MF suggests that the effect of short-range interactions is minimal in the steady-state and gets
washed out at these densities. At higher densities, the inter-atomic distances would be shorter and 1/r³ inter- actions would dominate, causing MF to break down.

<figure><img src="image_17.png" /><figcaption>FIG. 18. Comparing Cumulant (single realization) with MF (averaged over 10 realizations) in the steady-state for $N = 800$</figcaption></figure>

The cumulant method includes long-range beyond-MF two-point correlations. In Fig. 18, we show that the MF steady-state atomic inversion and forward intensity for $N = 800$ agree extremely well with those obtained from cumulant, across a range of $\Omega/\Gamma$ in the intermediate driving regime, where we expect the effect of correlations to be most apparent, compared to the weak and strong driving regimes. From the excellent agreement between MF and cumulant, we can see that the effect of long-range correlations is negligible in the steady-state of our system. Cumulant numerics with higher densities and larger-$N$ take extremely long runtimes, restricting us up to $N = 800$ here. However, as we already see good agreement with MACE-MF at $N = 2000$, we do not expect the cumulant results to change at least qualitatively up to $N = 2000$.

<figure><img src="image_18.png" /><figcaption>FIG. 19. Comparing ED with MF for the two-photon correlation function in the steady-state $\tilde{g}_2(0)$</figcaption></figure>

In Fig. 19, we compare the two-photon correlation function $\tilde{g}_2(0)$ in the steady-state of the MF dipolar model with that obtained from ED. We choose small $N$ for this comparison as the Hilbert space grows exponentially with $N$ and ED becomes numerically intractable for larger $N$. We keep the OD fixed while reducing $N$ so that the interaction strength is still similar to the $N = 2000$ system. We find that at small $\Omega/\Gamma$, there is a lot of variance in $\tilde{g}_2(0)$ across realizations and within these error bars, there seems to be a fair agreement between ED and MF. At large $\Omega/\Gamma$, the system becomes single-particle-like and error bars get smaller as interactions do not play a role and the change in atomic positions across realizations stops affecting the physics. Hence, there is a much better agreement between ED and MF in this regime. Due to finite $N$ corrections, $\tilde{g}_2(0)$ does not go from exactly 1 to 2, as the drive strength is increased.

[1] G. Ferioli, A. Glicenstein, I. Ferrier-Barbut, and A. Browaeys, Nat. Phys. pp. 1–5 (2023), ISSN 1745-2481, publisher: Nature Publishing Group, URL https://www.nature.com/articles/s41567-023-02064-w.

[2] M. Gross and S. Haroche, Physics Reports 93, 301 (1982), ISSN 0370-1573, URL https://www.sciencedirect.com/science/article/pii/S0370157382901028.

[3] K. Hammerer, A. S. Sørensen, and E. S. Polzik, Rev. Mod. Phys. 82, 1041 (2010), URL https://link.aps.org/doi/10.1103/RevModPhys.82.1041.

[4] M. J. Stephen, The Journal of Chemical Physics 40, 669 (1964), ISSN 0021-9606, https://pubs.aip.org/aip/jcp/article-pdf/40/3/669/11344002/669_1_online.pdf, URL https://doi.org/10.1063/1.1725188.

[5] R. H. Lehmberg, Phys. Rev. A 2, 883 (1970), URL https://link.aps.org/doi/10.1103/PhysRevA.2.883.

[6] R. Friedberg, S. Hartmann, and J. Manassah, Physics Reports 7, 101 (1973), ISSN 0370-1573, URL https://www.sciencedirect.com/science/article/pii/S037015737390001X.

[7] D. F. V. James, Phys. Rev. A 47, 1336 (1993), URL https://link.aps.org/doi/10.1103/PhysRevA.47.1336.

[8] S. D. Jenkins and J. Ruostekoski, Phys. Rev. A 86, 031602 (2012), URL https://link.aps.org/doi/10.1103/PhysRevA.86.031602.

[9] J. Javanainen and J. Ruostekoski, Opt. Express 24, 993 (2016), ISSN 1094-4087, URL https://opg.optica.org/abstract.cfm?URI=oe-24-2-993.

[10] D. Porras and J. I. Cirac, Phys. Rev. A 78, 053816 (2008), URL https://link.aps.org/doi/10.1103/PhysRevA.78.053816.

[11] A. A. Svidzinsky, J.-T. Chang, and M. O. Scully, Phys. Rev. Lett. 100, 160504 (2008), URL https://link.aps.org/doi/10.1103/PhysRevLett.100.160504.

[12] T. Ido, T. H. Loftus, M. M. Boyd, A. D. Ludlow, K. W. Holman, and J. Ye, Phys. Rev. Lett. 94, 153001 (2005), URL https://link.aps.org/doi/10.1103/PhysRevLett.94.153001.

[13] M. O. Scully, Phys. Rev. Lett. 102, 143601 (2009), URL https://link.aps.org/doi/10.1103/PhysRevLett.102.143601.

[14] R. Röhlsberger, K. Schlage, B. Sahoo, S. Couet, and R. Rüffer, Science 328, 1248 (2010), https://www.science.org/doi/pdf/10.1126/science.1187770, URL https://www.science.org/doi/abs/10.1126/science.1187770.

[15] J. Keaveney, A. Sargsyan, U. Krohn, I. G. Hughes,
D. Sarkisyan, and C. S. Adams, Phys. Rev. Lett. **108**, 173601 (2012), URL https://link.aps.org/doi/10.1103/PhysRevLett.108.173601.

[16] Z. Meir, O. Schwartz, E. Shahmoon, D. Oron, and R. Ozeri, Phys. Rev. Lett. **113**, 193002 (2014), URL https://link.aps.org/doi/10.1103/PhysRevLett.113.193002.

[17] J. Javanainen, J. Ruostekoski, Y. Li, and S.-M. Yoo, Phys. Rev. Lett. **112**, 113603 (2014), URL https://link.aps.org/doi/10.1103/PhysRevLett.112.113603.

[18] S. L. Bromley, B. Zhu, M. Bishof, X. Zhang, T. Bothwell, J. Schachenmayer, T. L. Nicholson, R. Kaiser, S. F. Yelin, M. D. Lukin, et al., Nat Commun **7**, 11039 (2016), ISSN 2041-1723, number: 1 Publisher: Nature Publishing Group, URL https://www.nature.com/articles/ncomms11039.

[19] S. J. Roof, K. J. Kemp, M. D. Havey, and I. M. Sokolov, Phys. Rev. Lett. **117**, 073003 (2016), URL https://link.aps.org/doi/10.1103/PhysRevLett.117.073003.

[20] J. Rui, D. Wei, A. Rubio-Abadal, S. Hollerith, J. Zeiher, D. M. Stamper-Kurn, C. Gross, and I. Bloch, Nature **583**, 369 (2020), ISSN 1476-4687, number: 7816 Publisher: Nature Publishing Group, URL https://www.nature.com/articles/s41586-020-2463-x.

[21] H. Ma and S. F. Yelin, *Collective lamb shift and modified linewidth of an interacting atomic gas* (2024), 2305.01865.

[22] N. O. Gjonbalaj, S. Ostermann, and S. F. Yelin, Phys. Rev. A **109**, 013720 (2024), URL https://link.aps.org/doi/10.1103/PhysRevA.109.013720.

[23] S. Jennewein, L. Brossard, Y. R. P. Sortais, A. Browaeys, P. Cheinet, J. Robert, and P. Pillet, Phys. Rev. A **97**, 053816 (2018), URL https://link.aps.org/doi/10.1103/PhysRevA.97.053816.

[24] A. Glicenstein, G. Ferioli, N. Šibalić, L. Brossard, I. Ferrier-Barbut, and A. Browaeys, Phys. Rev. Lett. **124**, 253602 (2020), URL https://link.aps.org/doi/10.1103/PhysRevLett.124.253602.

[25] M. O. Scully, E. S. Fry, C. H. R. Ooi, and K. Wódkiewicz, Phys. Rev. Lett. **96**, 010501 (2006), URL https://link.aps.org/doi/10.1103/PhysRevLett.96.010501.

[26] M. O. Scully, Laser Phys. **17**, 635 (2007), ISSN 1555-6611, URL https://doi.org/10.1134/S1054660X07050064.

[27] P. R. Berman, *American Journal of Physics* **78**, 1323 (2010), ISSN 0002-9505, URL https://pubs.aip.org/aapt/ajp/article-pdf/78/12/1323/12837154/1323.1.online.pdf, URL https://doi.org/10.1119/1.3475686.

[28] R. A. de Oliveira, M. S. Mendes, W. S. Martins, P. L. Saldanha, J. W. R. Tabosa, and D. Felinto, Phys. Rev. A **90**, 023848 (2014), URL https://link.aps.org/doi/10.1103/PhysRevA.90.023848.

[29] M. O. Araújo, I. Krešić, R. Kaiser, and W. Guerin, Phys. Rev. Lett. **117**, 073002 (2016), URL https://link.aps.org/doi/10.1103/PhysRevLett.117.073002.

[30] R. T. Sutherland and F. Robicheaux, Phys. Rev. A **93**, 023407 (2016), URL https://link.aps.org/doi/10.1103/PhysRevA.93.023407.

[31] B. Zhu, J. Cooper, J. Ye, and A. M. Rey, Phys. Rev. A **94**, 023612 (2016), URL https://link.aps.org/doi/10.1103/PhysRevA.94.023612.

[32] W. Guerin, in *Advances in Atomic, Molecular, and Optical Physics*, edited by L. F. DiMauro, H. Perrin, and S. F. Yelin (Academic Press, 2023), vol. 72 of Advances

*In Atomic, Molecular, and Optical Physics*, pp. 253-296, URL https://www.sciencedirect.com/science/article/pii/S1049250X23000022.

[33] J. Pellegrino, R. Bourgain, S. Jennewein, Y. R. P. Sortais, A. Browaeys, S. D. Jenkins, and J. Ruostekoski, Phys. Rev. Lett. **113**, 133602 (2014), URL https://link.aps.org/doi/10.1103/PhysRevLett.113.133602.

[34] 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), URL https://link.aps.org/doi/10.1103/PhysRevLett.116.233601.

[35] D. Pavolini, A. Crubellier, P. Pillet, L. Cabaret, and S. Liberman, Phys. Rev. Lett. **54**, 1917 (1985), URL https://link.aps.org/doi/10.1103/PhysRevLett.54.1917.

[36] T. Bienaimé, N. Piovella, and R. Kaiser, Phys. Rev. Lett. **108**, 123602 (2012), URL https://link.aps.org/doi/10.1103/PhysRevLett.108.123602.

[37] T. Bienaimé, R. Bachelard, N. Piovella, and R. Kaiser, Fortschritte der Physik **61**, 377 (2013), https://onlinelibrary.wiley.com/doi/pdf/10.1002/prop.201200089, URL https://onlinelibrary.wiley.com/doi/abs/10.1002/prop.201200089.

[38] M. O. Scully, Phys. Rev. Lett. **115**, 243602 (2015), URL https://link.aps.org/doi/10.1103/PhysRevLett.115.243602.

[39] W. Guerin, M. O. Araújo, and R. Kaiser, Phys. Rev. Lett. **116**, 083601 (2016), URL https://link.aps.org/doi/10.1103/PhysRevLett.116.083601.

[40] O. Rubies-Bigorda, S. Ostermann, and S. F. Yelin, Phys. Rev. A **107**, L051701 (2023), URL https://link.aps.org/doi/10.1103/PhysRevA.107.L051701.

[41] G. Ferioli, A. Glicenstein, L. Henriet, I. Ferrier-Barbut, and A. Browaeys, Phys. Rev. X **11**, 021031 (2021), URL https://link.aps.org/doi/10.1103/PhysRevX.11.021031.

[42] M.-T. Rouabah, M. Samoylova, R. Bachelard, P. W. Courteille, R. Kaiser, and N. Piovella, J. Opt. Soc. Am. A **31**, 1031 (2014), URL https://opg.optica.org/josaa/abstract.cfm?URI=josaa-31-5-1031.

[43] D. Plankensteiner, L. Ostermann, H. Ritsch, and C. Genes, Sci Rep **5**, 16231 (2015), ISSN 2045-2322, number: 1 Publisher: Nature Publishing Group, URL https://www.nature.com/articles/srep16231.

[44] R. J. Bettles, S. A. Gardiner, and C. S. Adams, Phys. Rev. A **94**, 043844 (2016), URL https://link.aps.org/doi/10.1103/PhysRevA.94.043844.

[45] G. Facchinetti, S. D. Jenkins, and J. Ruostekoski, Phys. Rev. Lett. **117**, 243601 (2016), URL https://link.aps.org/doi/10.1103/PhysRevLett.117.243601.

[46] R. T. Sutherland and F. Robicheaux, Phys. Rev. A **94**, 013847 (2016), URL https://link.aps.org/doi/10.1103/PhysRevA.94.013847.

[47] H. H. Jen, M.-S. Chang, and Y.-C. Chen, Phys. Rev. A **94**, 013803 (2016), URL https://link.aps.org/doi/10.1103/PhysRevA.94.013803.

[48] C. D. Mink and M. Fleischhauer, SciPost Phys. **15**, 233 (2023), URL https://scipost.org/10.21468/SciPostPhys.15.6.233.

[49] R. H. Dicke, Phys. Rev. **93**, 99 (1954), URL https://link.aps.org/doi/10.1103/PhysRev.93.99.

[50] H. J. Carmichael and D. F. Walls, Journal of Physics B: Atomic and Molecular Physics **10**, L685 (1977), URL
https://dx.doi.org/10.1088/0022-3700/10/18/002.

[51] D. F. Walls, P. D. Drummond, S. S. Hassan, and H. J. Carmichael, Progress of Theoretical Physics Supplement **64**, 307 (1978), ISSN 0375-9687, https://academic.oup.com/ptps/article-pdf/doi/10.1143/PTPS.64.307/5292058/64-307.pdf, URL https://doi.org/10.1143/PTPS.64.307.

[52] P. Drummond and H. Carmichael, Optics Communications **27**, 160 (1978), ISSN 00304018, URL https://linkinghub.elsevier.com/retrieve/pii/0030401878901980.

[53] D. F. Walls, Journal of Physics B: Atomic and Molecular Physics **13**, 2001 (1980), URL https://dx.doi.org/10.1088/0022-3700/13/10/008.

[54] L. M. Narducci, D. H. Feng, R. Gilmore, and G. S. Agarwal, Phys. Rev. A **18**, 1571 (1978), URL https://link.aps.org/doi/10.1103/PhysRevA.18.1571.

[55] B. R. Mollow, Phys. Rev. **188**, 1969 (1969), URL https://link.aps.org/doi/10.1103/PhysRev.188.1969.

[56] R. Loudon, *The quantum theory of light* (OUP Oxford, 2000).

[57] F. H. Bonifacio, P. Schwendimann, Phys. Rev. A **4**, 302 (1971), ISSN 2469-9934, URL https://doi.org/10.1103/PhysRevA.4.302.

[58] K. Hepp and E. H. Lieb, Annals of Physics **76**, 360 (1973), ISSN 0003-4916, URL https://www.sciencedirect.com/science/article/pii/0003491673900390.

[59] Y. K. Wang and F. T. Hioe, Phys. Rev. A **7**, 831 (1973), URL https://link.aps.org/doi/10.1103/PhysRevA.7.831.

[60] M. M. Roses and E. G. Dalla Torre, PLOS ONE **15**, 1 (2020), URL https://doi.org/10.1371/journal.pone.0235197.

[61] H. J. Carmichael, J. Phys. B: At. Mol. Phys. **13**, 3551 (1980), ISSN 0022-3700, URL https://iopscience.iop.org/article/10.1088/0022-3700/13/18/009.

[62] J. Hannukainen and J. Larson, Phys. Rev. A **98**, 042113 (2018), URL https://link.aps.org/doi/10.1103/PhysRevA.98.042113.

[63] V. Link, K. Luoma, and W. T. Strunz, Phys. Rev. A **99**, 062120 (2019), URL https://link.aps.org/doi/10.1103/PhysRevA.99.062120.

[64] D. Barberena and A. M. Rey, *Relevant averages in normal phase (in preparation)* (2023).

[65] D. E. Chang, J. Ye, and M. D. Lukin, Phys. Rev. A **69**, 023810 (2004), URL https://link.aps.org/doi/10.1103/PhysRevA.69.023810.

[66] R. B. Hutson, W. R. Milner, L. Yan, J. Ye, and C. Sanner, Science **383**, 384 (2024), URL https://doi.org/10.1126/science.adh4477.

[67] Y. C. O Morice and J. Dalibard, Phys. Rev. A **51**, 3896 (1995), URL https://link.aps.org/doi/10.1103/PhysRevLett.128.093001.

[68] L. Chomaz, L. Corman, T. Yefsah, R. Desbuquois, and J. Dalibard, New Journal of Physics **14**, 055001 (2012), URL https://dx.doi.org/10.1088/1367-2630/14/5/055001.

[69] A. M. Rey, L. Jiang, and M. D. Lukin, Phys. Rev. A **76**, 053617 (2007), URL https://link.aps.org/doi/10.1103/PhysRevA.76.053617.

[70] D. Roberts and A. A. Clerk, Phys. Rev. Lett. **131**, 190403 (2023), URL https://link.aps.org/doi/10.1103/PhysRevLett.131.190403.

[71] D. Goncalves, L. Bombieri, G. Ferioli, S. Pancaldi, I. Ferrier-Barbut, A. Browaeys, E. Shamoon, and D. E. Chang, *private communication (article in preparation)* (2024).

[72] G. Ferioli, S. Pancaldi, I. Ferrier-Barbut, and A. Browaeys, *private communication* (2024).

[73] A. Asenjo-Garcia, H. J. Kimble, and D. E. Chang, Proceedings of the National Academy of Sciences **116**, 25503 (2019), https://www.pnas.org/doi/pdf/10.1073/pnas.1911467116, URL https://www.pnas.org/doi/abs/10.1073/pnas.1911467116.

[74] T. Bilitewski, A. P. n. Orioli, C. Sanner, L. Sonderhouse, R. B. Hutson, L. Yan, W. R. Milner, J. Ye, and A. M. Rey, Phys. Rev. Lett. **128**, 093001 (2022), URL https://link.aps.org/doi/10.1103/PhysRevLett.128.093001.

[75] S. J. Masson, J. P. Covey, S. Will, and A. Asenjo-Garcia, *Dicke superradiance in ordered arrays of multilevel atoms* (2023), 2304.00093.

[76] D. Barberena, R. J. Lewis-Swan, J. K. Thompson, and A. M. Rey, Phys. Rev. A **99**, 053411 (2019), ISSN 2469-9926, 2469-9934, URL https://link.aps.org/doi/10.1103/PhysRevA.99.053411.

[77] K. Tucker, D. Barberena, R. J. Lewis-Swan, J. K. Thompson, J. G. Restrepo, and A. M. Rey, Phys. Rev. A **102**, 051701 (2020), URL https://link.aps.org/doi/10.1103/PhysRevA.102.051701.

[78] K. R. Hazzard, B. Gadway, M. Foss-Feig, B. Yan, S. A. Moses, J. P. Covey, N. Y. Yao, M. D. Lukin, J. Ye, D. S. Jin, et al., Phys. Rev. Lett. **113**, 195302 (2014), publisher: American Physical Society, URL https://link.aps.org/doi/10.1103/PhysRevLett.113.195302.