# Cooperative effects in dense cold atomic gases including magnetic dipole
interactions

N. S. Baßler  Department of Physics, Friedrich-Alexander Universität Erlangen-
Nürnberg (FAU), Staudtstraße 7, D-91058 Erlangen, Germany  Max Planck
Institute for the Science of Light, Staudtstraße 2, D-91058 Erlangen, Germany
I. Varma  Institut für Physik, Johannes Gutenberg-Universität Mainz, 55122
Mainz, Germany  M. Proske  Institut für Physik, Johannes Gutenberg-Universität
Mainz, 55122 Mainz, Germany  P. Windpassinger  Institut für Physik, Johannes
Gutenberg-Universität Mainz, 55122 Mainz, Germany  K. P. Schmidt  Department
of Physics, Friedrich-Alexander Universität Erlangen-Nürnberg (FAU),
Staudtstraße 7, D-91058 Erlangen, Germany  C. Genes  Department of Physics,
Friedrich-Alexander Universität Erlangen-Nürnberg (FAU), Staudtstraße 7,
D-91058 Erlangen, Germany  Max Planck Institute for the Science of Light,
Staudtstraße 2, D-91058 Erlangen, Germany

###### Abstract

We theoretically investigate cooperative effects in cold atomic gases
exhibiting both electric and magnetic dipole-dipole interactions, such as
occurring for example in clouds of dysprosium atoms. We distinguish between
the quantum degenerate case, where we take a many body physics approach and
the quantum non-degenerate case, where we use the formalism of open system
dynamics. For quantum non-degenerate gases, we illustrate the emergence of
tailorable spin models in the high-excitation limit. In the low-excitation
limit, we provide analytical and numerical results detailing the effect of
magnetic interactions on the directionality of scattered light and
characterize sub- and superradiant effects. For quantum degenerate gases, we
study the interplay between sub- and superradiance effects and the fermionic
or bosonic quantum statistics nature of the ensemble.

##  I Introduction

Scattering of light off atomic gases necessarily involves aspects of quantum
cooperativity, arising from the common, hybrid reaction of closely positioned
and mutually coupled quantum emitters to the external stimulation. The optical
response can be very complex, as it strongly depends on the gas density and
temperature, the type of atoms comprising the gas, as well as on the strength
of the driving field. For low density and high temperature, an independent
scattering regime can be obtained where the gas response can be deduced from
the single atom response Labeyrie _et al._ (1999, 2003). For higher density,
weak excitation, and still high temperature, cooperative aspects such as
super- and subradiance (spontaneous emission rates larger or smaller than that
of an isolated particle) are emerging Dicke (1954); Kwong _et al._ (2015);
Bromley _et al._ (2016); Inouye _et al._ (1999, 1999); Weiss _et al._ (2019);
Jennewein _et al._ (2018); de Oliveira _et al._ (2014); Roof _et al._ (2016);
Araújo _et al._ (2016); Guerin _et al._ (2016); Gross and Haroche (1982);
Scully _et al._ (2006); A. _et al._ (2010); Keaveney _et al._ (2012);
Javanainen _et al._ (2014); Kaiser (2009); Zhu _et al._ (2016); Javanainen and
Ruostekoski (2016); Pellegrino _et al._ (2014); Jenkins _et al._ (2016). High
driving powers lead then to non-linear optical effects, as the atomic
transitions can saturate and atoms become fundamentally nonlinear elements. In
all three described cases, a simple quantum optics approach suffices, based on
the open quantum system formalism, where the atoms are treated as pseudo-spins
$1/2$ (with transitions between ground and excited electronic orbitals)
responding to an external stimulation.

At low temperature, the motional wavepackets of the atoms comprising the gas
can overlap, leading to a change in the theoretical framework, which must
necessarily include a many body formulation to the problem Lewenstein _et al._
(1994); Ruostekoski and Javanainen (1997); Javanainen _et al._ (1999); Lee _et
al._ (2016). The two distinct limits, of classical versus quantum degenerate
gases, are illustrated in Fig. 1.

Figure 1: (a) An atomic gas is illuminated by a laser with frequency
$k_{\ell}$ while a detector in the far-field at angles $\theta$ and $\phi$
measures the intensity of the scattered light. (b) In a first case, the de
Broglie wavelength of each atom is much smaller than the average inter-
particle separation, allowing one to treat the system as a classical thermal
gas far from quantum degeneracy. The optical response can then be obtained in
a simple fashion as the response of many coupled dipoles to external coherent
stimulation. (c) In the opposite limit, quantum degeneracy is achieved by
lowering the temperature and consequently producing particles with overlapping
deBroglie wavepackets. To compute the optical response, a quantum many body
approach to the problem is necessary.

We analyze here a few distinct regimes and tailor our theoretical approach to
each individual case, in order to deduce radiative emission properties. In a
first step, for quantum non-degenerate gases, we extend the formulation of
cooperative response to include magnetic interactions on top of the standard
electric dipole-dipole exchanges, as for example strongly present in
experiments with dysprosium atoms Childs (1970); Petersen _et al._ (2020). In
the weak driving limit, the results show slight modifications in the
superradiant response. In the strong drive limit, we show the emergence of a
spin Hamiltonian with tunable parameters. In the next step, we assume an
externally imposed potential and analyze the response of the quantum
degenerate gas to external light drive in the case of both bosonic and
fermionic statistics. Interesting aspects occur, as bosonic statistics implies
superradiant behavior, similar to the standard Dicke superradiance example,
even in the absence of any particle-particle interactions. In the fermionic
statistics case, subradiance instead characterizes the emission properties of
the gas.

The manuscript is organized as follows. In Sec. II we introduce the light-
matter interaction model at the Hamiltonian level for both degenerate and non-
degenerate gases including both, electric and magnetic, dipole-dipole
interactions. In Sec. III we reduce our treatment to the non-degenerate case,
where we first analyze the strong driving limit and show the emergence of a
many particle spin Hamiltonian with interactions of tunable strength. We then
exemplify the magnetic dipole interactions effect onto the light scattered in
the weak driving regime and find small deviations from a purely electrically
interacting gas. In Sec. IV we analyze the light scattered from a quantum
degenerate gas and find scaling laws for both bosonic and fermionic
statistics. In Sec. V we describe an experimental platform based on thermal
clouds of atomic dysprosium where the effects described above could possibly
be tested. We conclude in Sec. VI.

##  II Model

We consider $\mathcal{N}$ atoms (positioned at $\bm{R}_{i}$ with index $i$
running from $1$ to $\mathcal{N}$) trapped in an external potential. Optical
addressing of each atom (both by classical fields and the quantum
electromagnetic vacuum modes) is achieved by coupling to its single valence
electron, which we denote by its position $\bm{r}_{i}$ measured with respect
to the position of the nucleus. The atomic cloud is illuminated by a laser
with Rabi frequency $\Omega$, wavevector $k_{\ell}$ and frequency
$\omega_{\ell}=ck_{\ell}$ propagating in the $x$ direction as depicted in Fig.
1a). Two distinct regimes emerge: i) a classical limit where atomic motion is
described by a thermal distribution of velocities and ii) a quantum limit
where atomic motion is quantized and atoms become indistinguishable, thus the
quantum degenerate case. The distinction between the two cases is illustrated
in Fig. 1b) and Fig. 1c) where the overlap of individual atomic de Broglie
wavepackets indicates the criterion for the transition between the two limits.
In both cases we will be interested in spectroscopic quantities obtained from
a detector positioned in the far field regime with angles $\theta$ and $\phi$
with respect to the incoming laser.

For simplicity of presentation we proceed with a set of approximations which
are usually performed in treating light-matter interactions. We perform the
dipole approximation which assumes that the size of the electronic orbital is
negligible with respect to any relevant optical transition wavelength
$\lambda$. In addition, we assume non-overlapping electronic orbitals between
neighboring atoms. The electric and magnetic fields are quantized in a
fictitious box of volume $\mathcal{V}$ and expressed in terms of photon
creation and annihilation bosonic operators $\hat{a}_{\bm{k},\bm{\epsilon}}$
with frequencies $\omega_{k}=c|\bm{k}|$ and commutators
$\left[\hat{a}_{\bm{k},\bm{\epsilon}_{\bm{k}}},\hat{a}_{\bm{k}^{\prime},\bm{\epsilon}_{\bm{k}}^{\prime}}^{\dagger}\right]=\delta_{\bm{k},\bm{k}^{\prime}}\delta_{\bm{\epsilon}_{\bm{k}},\bm{\epsilon}_{\bm{k}}^{\prime}}$
(where $\epsilon_{\bm{k}}$ is the polarization for the mode $\bm{k}$). The
free space quantization yields the electric and magnetic field operators

$\displaystyle\hat{\bm{E}}$
$\displaystyle=\mathrm{i}\sum_{\bm{k},\bm{\epsilon}}g_{k}\omega_{k}\epsilon_{\bm{k}}\left(\hat{a}_{\bm{k},\bm{\epsilon}}\mathrm{e}^{\mathrm{i}\bm{k}\bm{R}}-\hat{a}_{\bm{k},\bm{\epsilon}}^{\dagger}\mathrm{e}^{-\mathrm{i}\bm{k}\bm{R}}\right)$
(1) $\displaystyle\hat{\bm{B}}$
$\displaystyle=\mathrm{i}\sum_{\bm{k},\bm{\epsilon}}\left(\bm{k}\times\epsilon_{\bm{k}}\right)g_{k}\left(\hat{a}_{\bm{k},\bm{\epsilon}}\mathrm{e}^{\mathrm{i}\bm{k}\bm{R}}-\hat{a}_{\bm{k},\bm{\epsilon}}^{\dagger}\mathrm{e}^{-\mathrm{i}\bm{k}\bm{R}}\right),$
(2)

with $g_{k}=1/\sqrt{2\omega_{k}\mathcal{V}\epsilon_{0}}$ being photon coupling
strength. From here on the treatment is distinct for the two assumed limits of
either a gas with a classical distribution of velocities or a quantum
degenerate gas of indistinguishable atoms.

###  II.1 Semiclassical gas approach

We follow the standard quantum optics approach where the internal electronic
dynamics of each atom is treated in terms of Pauli matrices. To this end, we
restrict the dynamics of the electron to a ground state $\lvert g\rangle_{i}$
and a single excited state $\lvert e\rangle_{i}$, separated by frequency
$\omega$. Notice that in practice, this assumption is a very good
approximation for optically pumped atoms. An example based on dysprosium atoms
is shown in Sec. V where the ground state is represented by a magnetic
sublevel with $J=8$ and $m_{J}=-8$ and the excited state with $J^{\prime}=9$
and $m_{J}^{\prime}=-9$.

The Pauli matrices represent transitions in the electronic degrees of freedom
and are defined as $\sigma_{i}=\lvert g\rangle_{i}\langle e\rvert_{i}$. The
coupling to electromagnetic waves occurs via either the transition electric
$\hat{\bm{d}}_{i}^{\phantom{\dagger}}=\bm{d}_{i}\sigma_{i}^{\phantom{\dagger}}+\bm{d}^{*}_{i}\sigma_{i}^{\dagger}$
and magnetic
$\hat{\bm{\mu}}_{i}^{\phantom{\dagger}}=\bm{\mu}_{i}\sigma_{i}^{\phantom{\dagger}}+\bm{\mu}^{*}_{i}\sigma_{i}^{\dagger}$
dipole operators or via the static magnetic dipole operator
$\hat{\bm{\mu}}_{s,i}^{\phantom{\dagger}}=\bm{\mu}_{i,e}\sigma_{i}^{\dagger}\sigma+\bm{\mu}_{i,g}\sigma_{i}\sigma^{\dagger}$.
The static components of the magnetic dipole are computed within the
respective electronic state $\bm{\mu}_{i,e}=-g\mu_{B}\langle
e_{i}\rvert\hat{\bm{L}}_{i}\lvert e_{i}\rangle$ and
$\bm{\mu}_{i,g}=-g\mu_{B}\langle g_{i}\rvert\hat{\bm{L}}_{i}\lvert
g_{i}\rangle$ where $\bm{L}_{i}$ is the angular momentum operator for atom
$i$, $\mu_{B}$ is the Bohr magneton and $g$ is the Landé factor. The
transition dipole matrix elements are computed between orbitals
$\bm{d}_{i}=-e\left\langle
e_{i}\middle|\hat{{\bm{r}}_{i}}\middle|g_{i}\right\rangle$ and
$\bm{\mu}_{i}=-g\mu_{B}\langle g_{i}\rvert\hat{\bm{L}}\lvert e_{i}\rangle$.
Due to the fixed parity of hydrogen-like orbitals the static electric dipole
moment must vanish which is not the case for the static magnetic dipole
moment. The resulting Hamiltonian

$\mathcal{H}=\mathcal{H}_{\text{at}}+\mathcal{H}_{\text{em}}+\mathcal{H}_{\text{el+mag}}+\mathcal{H}_{\text{drive}},$
(3)

is a sum over the free Hamiltonians of the atoms $\mathcal{H}_{\text{at}}$,
the electromagnetic vacuum Hamiltonian $\mathcal{H}_{\text{em}}$, the electric
and magnetic dipole coupling to the radiation field
$\mathcal{H}_{\text{el+mag}}$, and the semiclassical drive
$\mathcal{H}_{\text{drive}}$.
The first two terms are explicitly written as

$\mathcal{H}_{\text{at}}+\mathcal{H}_{\text{em}}=\omega\sum_{i}\sigma_{i}^{\dagger}\sigma_{i}^{\phantom{\dagger}}+\sum_{\bm{k},\bm{\epsilon}}\omega_{k}\hat{a}_{\bm{k},\bm{\epsilon}_{\bm{k}}}^{\dagger}\hat{a}_{\bm{k},\bm{\epsilon}_{\bm{k}}}^{\phantom{\dagger}},$
(4)

where we have set the zero of the energy at the ground electronic state level
and ignored the zero point energy of the vacuum modes. The coupling of the
atomic system to the electric and magnetic quantum fields is

$\mathcal{H}_{\text{el+mag}}=\sum_{i}\hat{\bm{d}}_{i}\cdot\hat{\bm{E}}(\bm{R}_{i})+\left(\hat{\bm{\mu}}_{t,i}+\hat{\bm{\mu}}_{s,i}\right)\cdot\hat{\bm{B}}(\bm{R}_{i}).$
(5)

The semiclassical drive is characterized by the Rabi frequency $\Omega$ and is
written as

$\mathcal{H}_{\text{drive}}=\Omega\sum_{i}\left(\mathrm{e}^{-\mathrm{i}{\bm{k}}_{\ell}\bm{R}_{i}}\sigma_{i}+\mathrm{e}^{\mathrm{i}{\bm{k}}_{\ell}\bm{R}_{i}}\sigma_{i}^{\dagger}\right),$
(6)

where ${\bm{k}}_{\ell}=k_{\ell}\hat{e}_{x}$ and $\hat{e}_{x}$ is the unit
vector in the x-direction. The Hamiltonian above is the starting point for the
models analyzed in Sec. III.

###  II.2 Quantum degenerate gas approach

In the opposite limit of a quantum degenerate gas, it is more convenient to
introduce a two-species model where field operators
$\Psi^{\dagger}_{e,g}(\bm{R})$ create atoms at some position $\bm{R}$ with
their electron either in the ground or in the excited state. The dipole
moments can now be written in terms of field operators as the action of the
Pauli matrices is now expressed by combinations of the field operators. For
example, the operator
$\hat{\Sigma}(\bm{R})=\Psi^{\dagger}_{g}(\bm{R})\Psi^{\phantom{\dagger}}_{e}(\bm{R})$
is the field-theoretical equivalent of the matrix operator $\sigma_{i}$ where
the particle location, which was previously denoted by an index, is now
denoted by a continuous variable due to the replacement of localized
scatterers by fields. The meaning is that an atom in the excited state is
destroyed while another atom in the ground state is created in exactly the
same position. The commutation relations are the standard ones
$\left[\Psi^{\phantom{\dagger}}_{\alpha}(\bm{R}),\Psi^{\dagger}_{\beta}(\bm{R})\right]_{\zeta}=\delta(\bm{R}-\bm{R}^{\prime})\delta_{\alpha,\beta}$,
where $\alpha,\beta\in\\{e,g\\}$ are species indices while $\zeta$ specifies
commutation or anticommutation relations depending on the bosonic or fermionic
nature of the gas. The single particle motional Hamiltonian is written as

$\mathcal{H}_{0}=-\bm{\nabla}^{2}/2M+V_{\text{ext}}(\bm{R}),$ (7)

where $V_{\text{ext}}(\bm{R})$ is an externally applied potential (assumed
quadratic in the following) that allows for the use of a simpler notation
using the trap basis states. Notice that, with a state-independent choice for
$V_{\text{ext}}(\bm{R})$, the motional wavefunctions are the same for the
excited type and ground state type atoms.

In second quantization, the total system Hamiltonian is obtained as an
integration of the Hamiltonian density. As a next step, the field operators
can be expanded in a conveniently chosen basis. In particular, we will
consider the trap basis defined by the eigenvectors
$\mathcal{H}_{0}\phi_{\bm{n}}=\omega_{\bm{n}}\phi_{\bm{n}}$. Particle creation
and annihilation operators can then be defined as

$\begin{split}\hat{g}_{\bm{n}}&=\int\mathop{}\\!\mathrm{d}{\bm{R}}\left[\phi_{\bm{n}}(\bm{R})\right]^{*}\Psi_{g}(\bm{R})\\\
\hat{e}_{\bm{n}}&=\int\mathop{}\\!\mathrm{d}{\bm{R}}\left[\phi_{\bm{n}}(\bm{R})\right]^{*}\Psi_{e}(\bm{R}).\end{split}$
(8)

The operators can be interpreted in the following way: when
$\hat{g}_{\bm{n}}^{\phantom{\dagger}}$ is applied to the vacuum, it creates an
atom in the trap state $\bm{n}$ and in the electronic ground state. Similarly,
$\hat{e}_{\bm{n}}^{\phantom{\dagger}}$ creates an atom in the trap state
$\bm{n}$ and in the electronic excited state. The fermionic or bosonic nature
of the atoms is consistently taken into account by the commutation relations
of these operators. The transition from free particles to trapped ones is
performed by the tuning of the single particle trapping potential which in
turn affects the shape of the trap basis eigenvectors. With this, the free
Hamiltonian of the atoms can be written as

$\mathcal{H}_{\text{at}}=\sum_{\bm{n}}\omega_{\bm{n}}\hat{g}^{\dagger}_{\bm{n}}\hat{g}^{\phantom{\dagger}}_{\bm{n}}+\sum_{\bm{n}}(\omega_{0}+\omega_{\bm{n}})\hat{e}^{\dagger}_{\bm{n}}\hat{e}^{\phantom{\dagger}}_{\bm{n}}.$
(9)

The electromagnetic modes Hamiltonian $\mathcal{H}_{\text{em}}$ is the same as
before. The coupling of the atomic system to the electric and magnetic quantum
fields is described by the following Hamiltonian

$\begin{split}\mathcal{H}_{\text{el+mag}}&=\sum_{\bm{n},\bm{m}}\left(\hat{\Sigma}_{\bm{n},\bm{m}}^{\phantom{\dagger}}+\hat{\Sigma}_{\bm{n},\bm{m}}^{\dagger}\right)\left\langle\bm{n}\middle|\bm{d}\cdot\bm{E}(\bm{R})\middle|\bm{m}\right\rangle\\\
&+\sum_{\bm{n},\bm{m}}\left(\hat{\Sigma}_{\bm{n},\bm{m}}^{\phantom{\dagger}}+\hat{\Sigma}_{\bm{n},\bm{m}}^{\dagger}\right)\left\langle\bm{n}\middle|\bm{\mu}\cdot\bm{B}(\bm{R})\middle|\bm{m}\right\rangle\\\
&+\sum_{\bm{n},\bm{m}}\left(\hat{g}^{\dagger}_{\bm{n}}\hat{g}^{\phantom{\dagger}}_{\bm{m}}\bm{\mu}_{g}+\hat{e}^{\dagger}_{\bm{n}}\hat{e}^{\phantom{\dagger}}_{\bm{m}}\bm{\mu}_{e}\right)\cdot\left\langle\bm{n}\middle|\bm{B}(\bm{R})\middle|\bm{m}\right\rangle,\end{split}$
(10)

where
$\hat{\Sigma}_{\bm{n},\bm{m}}^{\phantom{\dagger}}=\hat{g}^{\dagger}_{\bm{n}}\hat{e}^{\phantom{\dagger}}_{\bm{m}}$
are ladder operators destroying an excited atom in trap state $\bm{m}$ and
creating a ground state atom in trap state $\bm{n}$. The semiclassical drive
is now written as

$\mathcal{H}_{\text{drive}}=\Omega\sum_{\bm{n},\bm{m}}{\hat{\Sigma}}_{\bm{n},\bm{m}}^{\phantom{\dagger}}\eta_{\bm{n},\bm{m}}+\hat{\Sigma}_{\bm{n},\bm{m}}^{\dagger}\eta_{\bm{n},\bm{m}},$
(11)

where we have defined the Franck-Condon factors
$\eta_{\bm{n},\bm{m}}=\left\langle\bm{n}\middle|\mathrm{e}^{-\mathrm{i}\bm{k}_{\ell}\bm{R}}\middle|\bm{m}\right\rangle=\int\mathop{}\\!\mathrm{d}{\bm{R}}\left[\phi_{\bm{n}}(\bm{R})\right]^{*}{\mathrm{e}^{-\mathrm{i}\bm{k}_{\ell}\bm{R}}}\phi_{\bm{m}}(\bm{R})$.
For tight trapping conditions, where the localization of the atoms is on a
level much smaller than the wavelength $2\pi/k_{\ell}$, the exponential can be
approximated with unity and the Franck-Condon factors are equal to
$\delta_{\bm{n}\bm{m}}$. In addition, to take into account the first-order
correction, a Lamb-Dicke limit approximation can be made and only the matrix
elements of the linear term $\mathrm{i}\bm{k}_{\ell}\bm{R}$ would need to be
considered.

##  III The classical motion limit

Let us now focus on the semiclassical case, which assumes classical dynamics
for the atomic motion, while the electronic dynamics is described in a quantum
fashion in terms of Pauli matrices. As a function of the drive intensity, the
emergent physics can be quite distinct. In a first case, where a high-
intensity drive is assumed, high-excitation levels can be reached, with many
atoms being excited at the same time. This is the standard regime for Dicke
superradiance, i.e. the quick burst of spontaneous emission from an initial
fully excited ensemble of closely spaced atoms. The other limit we consider is
the weak-excitation limit, where there are hardly any excitations present in
the system, rendering it possible to perform a transformation to a fully
classical coupled dipole model. Despite its simplicity, the weak-excitation
limit gives insights in the emergence of cooperative effects and their role in
modifying directional scattering of light. In particular, we emphasize the
role of magnetic dipole-dipole interactions and their competition with the
widely studied electric counterpart.

Figure 2: Illustration of interactions present in the effective Hamiltonian
Eq. (12). On the left is the excitation hopping due to the electronic dipole
interaction between sites i and j. Initial occupations are indicated in filled
open circles and the final occupations in dotted open circles. As indicated by
the filled circle, the initially excited state on site $i$ is changed to the
ground state by $\sigma_{i}$. Similarly, the ground state on site j, indicated
by a filled circle is excited due to $\sigma_{j}^{\dagger}$. On the right side
is the density-density interaction due to the static magnetic dipole
interaction with the same color coding. Here, either the ground or excited
state can be occupied on $i$ and $j$, leading to different interactions
indicated by arrows between them and annotated with the appropriate
interaction terms.

###  III.1 Effective spin Hamiltonian

From the Hamiltonian listed in Eq. (3), one can derive a master equation for
the evolution of the $\mathcal{N}$ electronic systems under the approximation
of a frozen gas. Such approximation can hold for low enough temperatures,
where motion evolution is slow compared to the time taken by any radiative
processes. The derivation is based on the elimination of the photonic degrees
of freedom Lewenstein _et al._ (1994); James (1993); Lehmberg (1970);
Ruostekoski and Javanainen (1997); Javanainen _et al._ (1999), as outlined in
detail in Apps. A, B, C, D. The effective Hamiltonian describing dynamics in
the reduced subspace of dimension $2^{\mathcal{N}}$ of the electronic degrees
of freedom reads

$\begin{split}\mathcal{H}_{\text{eff}}&=\Delta\sum_{i}\sigma_{i}^{\dagger}\sigma_{i}+\sum_{i\neq
j}\left(g^{(d)}_{ji}+g^{(\mu)}_{ji}\right)\sigma_{i}^{\dagger}\sigma_{j}^{\phantom{\dagger}}\\\
&+\Omega\sum_{i}\left(\mathrm{e}^{-\mathrm{i}\bm{R}_{i}{\bm{k}}_{\ell}}\sigma_{i}+\mathrm{e}^{\mathrm{i}\bm{R}_{i}{\bm{k}}_{\ell}}\sigma_{i}^{\dagger}\right)\\\
&+\frac{1}{8}\sum_{i\neq
j}\Omega^{e,e}_{i,j}(1+\sigma_{i}^{z})(1+\sigma_{j}^{z})\\\
&+\frac{1}{8}\sum_{i\neq
j}\Omega^{g,g}_{i,j}(1-\sigma_{i}^{z})(1-\sigma_{j}^{z})\\\
&+\frac{1}{4}\sum_{i\neq
j}\Omega^{g,e}_{i,j}(1-\sigma_{i}^{z})(1+\sigma_{j}^{z}).\end{split}$ (12)

This describes the unitary part of the interaction corresponding to a full-
fledged spin-1/2 XXZ model. Here the detuning is defined as
$\Delta=\omega-\omega_{\ell}$. The second term gives the usual electric and
magnetic dipole-dipole interactions allowing for the hopping of excitations
within the whole ensemble, describing the XY part. The last three terms
represent the contribution from the static magnetic interactions and lead to
frequency shifts conditioned on the occupancy of the pair of atoms involved in
the interaction, which is typically referred to as an effective Ising
interaction. The coherent photon exchange via electric dipole-dipole
interactions $g^{(d)}_{ji}$ has been widely studied Bromley _et al._ (2016).
The magnetic transition dipole-dipole couplings $g^{(\mu)}_{ji}$ generally can
be ignored as they are of very small magnitude compared to the electric ones.
However, other important couplings occur

$\Omega^{\alpha,\beta}_{i,j}=\frac{\mu_{0}}{4\pi
R_{ij}^{3}}\left[3\frac{(\bm{\mu}_{i,\alpha}\cdot\bm{R}_{ij})(\bm{\mu}_{j,\beta}\cdot\bm{R}_{ij})}{R_{ij}^{2}}-\bm{\mu}_{i,\alpha}\cdot\bm{\mu}_{j,\beta}\right]`$
(13)

with $\alpha,\beta\in\\{e,g\\}$, which describe density-density interactions
between atoms owing to the static magnetic dipoles, illustrated in Fig. 2 on
the right side.

With the elimination of the electromagnetic vacuum, the system is
characterized by open system dynamics where the collective dissipation is
included in the master equation

$\frac{\mathop{}\\!\mathrm{d}\rho}{\mathop{}\\!\mathrm{d}t}=-\mathrm{i}\left[\mathcal{H}_{\text{eff}},\rho\right]+\mathcal{L}_{\mu}[\rho]+\mathcal{L}_{d}[\rho].$
(14)

We assume standard Lindblad form for the loss terms which we write as follows

$\mathcal{L}_{\alpha}[\rho]=\sum_{i,j=1}^{\mathcal{N}}f^{(\alpha)}_{ij}\left(2\sigma_{i}\rho\sigma_{j}^{\dagger}-\left\\{\sigma_{j}^{\dagger}\sigma_{i},\rho\right\\}\right).$
(15)

The independent radiative loss rates are
$\Gamma=k_{0}^{3}d^{2}/(3\pi\hbar\epsilon_{0})$ and
$\Gamma_{\mu}=k_{0}^{3}\mu^{2}/3\pi\hbar\epsilon_{0}$ stemming from electric
and magnetic contributions. Moreover, collective dissipation at rates
$f^{(d)}_{ij}$ (electric) and $f^{(\mu)}_{ij}$ (magnetic) are also present,
with exact expressions listed in App. D.

###  III.2 The weak-excitation limit: coupled dipoles model

The Hamiltonian in Eq. (12) is of general validity. However, by considering
the weak-excitation limit, the evolution of the system is restricted to a very
small subspace where the analytical description of the dynamics can be greatly
simplified. In order to do this, we linearize the time evolution of the system
by replacing all population operators with $-1$ and factorizing all two
operator correlations. We gather all expectation values of atomic coherences
in a single vector
$\bm{v}=(\left\langle\sigma_{1}\right\rangle,...,\left\langle\sigma_{\mathcal{N}}\right\rangle)^{\top}$
and write an effective first-order differential equation

$\dot{\bm{v}}(t)=-i\mathcal{M}\bm{v}(t)+\Omega\bm{v}_{\text{drive}},$ (16)

where the drive vector incorporates all the drive phases
$\bm{v}_{\text{drive}}=(\mathrm{e}^{\mathrm{i}{\bm{k}}_{\ell}\bm{R}_{1}},...,\mathrm{e}^{\mathrm{i}{\bm{k}}_{\ell}\bm{R}_{\mathcal{N}}})^{\top}$.
The weak driving approximation $\Omega\ll\lvert\bm{f}^{(d)}\rvert$ assumes
that the Rabi frequency is much weaker than the dissipative part of the
dipole-dipole interaction. Notice also that the equation above can be solved
both in steady state to derive spectroscopic features of the ensemble as well
as in the time domain (by imposing a time dependence on $\Omega$). The matrix
$\mathcal{M}$ incorporates both coherent and dissipative cooperative behavior

$\mathcal{M}_{jj^{\prime}}=\left[\Delta+\delta\omega_{j}\right]\delta_{jj^{\prime}}-\left[g_{jj^{\prime}}^{(d)}+\mathrm{i}f_{jj^{\prime}}^{(d)}\right].$
(17)

We used the low-excitation condition in order to approximate
$\left\langle\sigma^{z}_{j}\sigma_{j^{\prime}}^{\phantom{\dagger}}\right\rangle\approx-\left\langle\sigma_{j^{\prime}}^{\phantom{\dagger}}\right\rangle$.
This is always true if $j=j^{\prime}$, but only approximately valid under the
condition that very few excitations are present in the system such that, on
average, each site has much lower than unit population in the excited state.
This leads to the definition of a local frequency shift

$\begin{split}\delta\omega_{j}&=\sum_{j^{\prime}\neq
j}\left(\Omega_{j^{\prime}j}^{e,g}-\Omega_{j^{\prime}j}^{g,g}\right)\\\
&=\frac{(\mu_{e}-\mu_{g})\mu_{g}\mu_{0}}{4\pi}\sum_{j^{\prime}\neq
j}\frac{3z_{ij}^{2}/R_{ij}^{2}-1}{R_{ij}^{3}},\end{split}$ (18)

where $\delta\omega_{j}$ denotes the total shift acquired owing to magnetic
interactions. We have assumed that all magnetic dipole moments in the ground
state are $\mu_{g}\hat{e}_{z}$ and all magnetic moments in the excited state
are $\mu_{e}\hat{e}_{z}$. This leads to the definition of the magnetic
interaction rate

$\Omega_{\mu}=\frac{\mu_{0}(\mu_{e}-\mu_{g})\mu_{g}}{4\pi}.$ (19)

Additionally, we have neglected contributions from the magnetic transition
dipole moments couplings as they are much smaller than the electrically
mediated couplings. We then consider the coherent part of the far-field
intensity radiated by a system defined by these dipoles

$I(\bm{r}_{s})\propto\left|\sum_{i}(1-z_{s}^{2})\mathrm{e}^{-\mathrm{i}\bm{k}_{s}\bm{R}_{i}}v_{i}\right|^{2}$
(20)

where $\bm{r}_{s}$ indicates the detection angle,
$\bm{k}_{s}\parallel\bm{r}_{s}$ and $\lvert\bm{k}_{s}\rvert=k_{0}$ is the wave
vector for the propagation in detection direction.

Figure 3: Numerical results based on solution of Eq. 16 in steady state. We
consider multiple configurations of a cloud where the particles are
distributed according to a Gaussian probability distribution with optical
depth $\tau=3\mathcal{N}/(2k_{0}^{2}R^{2})$, simulating averaging over atomic
motion. The number of configurations is $\mathcal{N}_{c}=1000$ and
$\mathcal{N}=1000$. The optical depth is $\tau=15$. $\Omega_{\mu}$ quantifies
the strength of the magnetic interaction. In (a) we show how the forward
enhancement is modified by the magnetic interaction by plotting the intensity
variation versus the in-plane angle $\theta$ on the single-particle resonance
$\Delta=0$. There is a reduction in the forward enhancement for
$\Omega_{\mu}=0.03\Gamma\lambda^{-3}$ ($\Omega_{\mu}k_{0}^{-3}=8\Gamma$)
compared to no magnetic interaction. In (b) we investigate the lineshape of
the cloud by changing the laser detuning $\Delta$ and compare
$\Omega_{\mu}=0.03\Gamma\lambda^{-3}$ ($\Omega_{\mu}k_{0}^{-3}=8\Gamma$) with
$\Omega_{\mu}=0$ at an in-plane detection angle of $\theta=\pi/10$. The double
peak structure in the non-magnetic case becomes asymmetric when the magnetic
interaction is turned on. Changing the sign flips the effect on the two peaks
(not shown). In (c) and (d) we investigate the effect of the magnetic
interaction on the cooperative properties of the sample. To this end we
determine the collective lineshift $\bar{\Delta}$ and the collective linewidth
$\bar{\Gamma}$ in the forward direction $\theta=0$ by calculating the
lineshapes and then extracting these parameters from a Lorentzian fit. The
lineshift and linewidth enhancement are shown for different optical densities
$\tau$ with the magnetic interaction strength on the x-axis.

In steady state, assuming time-independent driving, the linear system obtained
from setting the left-hand side of Eq.(16) to zero has to be solved. In
practice, this equation needs to be solved for many different configurations.
Since we assume the cloud to be thermal, we assume this particle distribution
with respect to the center of the cloud to be Gaussian. In Fig. 3(a), we can
see the straightforward effect of a reduction of the collective motion and
consequently, the forward scattering enhancement due to frequency disorder.
This is simply due to the fact that motional disorder increases the coupling
between the driven mode and modes which do not emit in the laser direction.
The lineshape of the light scattered from a cloud, for two values of the
magnetic interaction strength $\Omega_{\mu}$, is shown in Fig. 3(b). For the
dense ensemble we can already observe a double peak structure. The magnetic
field shifts the spectral weight of the peaks to yield an asymmetric double
peak structure which can be interpreted as a spectral shift of the eigenvalues
due to the frequency disorder.

In Fig. 3(c) and (d) we perform studies of the lineshift and linewidth
modification due to the presence of the magnetic interaction for different
optical densities $\tau$. From the lineshift we can estimate the total
symmetric contribution to the frequency shift of the ensemble, which causes a
lineshift that depends on the sign of the magnetic dipole-dipole interaction.
From the linewidth enhancement one can see the effect of the asymmetric
contribution to the frequency disorder. From $\tau=20$ it is fairly clear that
there are two different regimes. For large magnetic interaction $\Omega_{\mu}$
the large frequency disorder simply leads to an inhomogeneous linewidth
enhancement. For smaller $\Omega_{\mu}\approx 0.005\lambda^{-3}\Gamma$
($\Omega_{\mu}k_{0}^{-3}\approx\Gamma_{0}$) the linewidth is actually
decreased. One explanation for this is that the frequency disorder facilitates
scattering into states which are naturally more subradiant than the originally
driven state without disorder.

##  IV Quantum degenerate case

We consider now the case of a quantum degenerate gas and investigate the
effect of the quantum statistical properties of the particles. As a first
step, we follow a procedure which is standard for the derivation of the Dicke
superradiance phenomenon with inverted two-level systems. This means, we will
consider a non-interacting bosonic/fermionic trapped gas and will ask how does
the emitted intensity scale as a function of the quantum statistics and as a
function of trap parameters.

We start with a quantum degenerate gas in contact with a thermal reservoir at
controllable temperature $T$. Tuning $T$ allows one to go from the quantum
degenerate case to a completely classical, thermal gas regime (for high
temperature). We will compute the intensity $I(\bm{k}_{s})$ of angle-resolved
spontaneous emission from the system with $\bm{k}_{s}$ being the wave vector
in direction of the detector and $|\bm{k}_{s}|=k_{0}$. This intensity is
defined as the correlation of the Fourier transforms of the polarization
operators

$\begin{split}I(\bm{k}_{s})\propto\left\langle\hat{\Sigma}^{\dagger}(\bm{k}_{s})\hat{\Sigma}(\bm{k}_{s})\right\rangle,\end{split}$
(21)

where the Fourier transformation is defined as follows
$\hat{\Sigma}(\bm{k}_{s})=\int\mathop{}\\!\mathrm{d}{\bm{R}}\mathrm{e}^{-\mathrm{i}\bm{k}_{s}\bm{R}}\hat{\Sigma}(\bm{R})$.
Here $\left\langle\right\rangle$ stands for the ensemble average. We take the
state in the following to be thermal in the excited and ground state so that
the average becomes a thermal one where excited and ground state factorize.

For non-interacting scatterers in thermal equilibrium, the expectation value
of a generic observable $\hat{A}$ can be written as

$\left\langle\hat{A}\right\rangle=\frac{1}{Z}\sum_{\bm{N}}\mathrm{e}^{-\beta\omega(\bm{N})}\left\langle\bm{N}\middle|\hat{A}\middle|\bm{N}\right\rangle,$
(22)

where
$\bm{N}=(\mathcal{N}_{g_{0}},\mathcal{N}_{e_{0}},\mathcal{N}_{g_{1}},\mathcal{N}_{e_{1}},\ldots)$
defines a microstate where the subindices $0,1,2,\ldots$ represent higher
dimensional tuples of indices. This notation simply indicates that there are
$\mathcal{N}_{\alpha_{i}}$ particles in the state with quantum numbers
$\alpha_{i}$. Denoted by $\lvert\bm{N}\rangle$ is the quantum state
corresponding to these quantum numbers and $\omega(\bm{N})$ is the total
energy of the state. In the case of harmonically (frequency $\omega_{t}$)
trapped atoms, the total energy decomposes into summands

$\omega(\bm{N})=\omega_{0}\sum_{n}\mathcal{N}_{e_{n}}+\omega_{t}\sum_{n}n\left(\mathcal{N}_{g_{n}}+\mathcal{N}_{e_{n}}\right).$
(23)

Let us first remark that, at moderate temperatures, owing to the fact that
$\omega_{0}\ll\omega_{t}$, one can safely neglect the thermally induced
occupancy of the excited electronic state, meaning that the first term
containing $\omega_{0}$ can be ignored in the exponent
$\mathrm{e}^{-\beta\omega(\bm{N})}$.

In order to estimate the angle resolved scattered intensity we will make use
of the trap basis. We utilize the decomposition of the field operators in this
basis $\Psi_{g}(\bm{R})=\sum_{\bm{n}}\phi_{\bm{n}}(\bm{R})\hat{g}_{\bm{n}}$
and $\Psi_{e}(\bm{R})=\sum_{\bm{n}}\phi_{\bm{n}}(\bm{R})\hat{e}_{\bm{n}}$ and
expand the two position correlation function as follows

$\displaystyle\left\langle\hat{\Sigma}^{\dagger}(\bm{R})\hat{\Sigma}(\bm{R}^{\prime})\right\rangle$
$\displaystyle=\sum_{\begin{subarray}{c}\bm{n}_{1},\bm{n}_{2}\\\
\bm{m}_{1},\bm{m}_{2}\end{subarray}}\phi^{*}_{\bm{m}_{1}}(\bm{R})\phi_{\bm{n}_{1}}(\bm{R})\phi_{\bm{m}_{2}}(\bm{R}^{\prime})\phi^{*}_{\bm{n}_{2}}(\bm{R}^{\prime})$
$\displaystyle\phantom{\sum_{\bm{n}_{1},\bm{n}_{2}}}\times\left\langle\hat{e}_{\bm{m}_{1}}^{\dagger}\hat{g}_{\bm{n}_{1}}^{\phantom{\dagger}}\hat{g}_{\bm{n}_{2}}^{\dagger}\hat{e}_{\bm{m}_{2}}^{\phantom{\dagger}}\right\rangle.$
(24)

In order to compute the expression in Eq. (21) we apply the Fourier transform
definition and identify the Franck-Condon factors stemming from the
integration over space

$\eta_{\bm{n},\bm{m}}=\int\mathop{}\\!\mathrm{d}{\bm{R}}\mathrm{e}^{-\mathrm{i}\bm{R}\bm{k}_{s}}\phi_{\bm{n}}^{*}(\bm{R})\phi_{\bm{m}}^{\phantom{\dagger}}(\bm{R}).$
(25)

Next, we determine the thermal expectation values. In doing so, we ignore the
contribution in the exponential stemming from $\omega_{0}$, corresponding to
the energy of the electronic excited state, as its thermal occupancy is
negligible even at room temperature. The result for the four-operator
correlation becomes

$\left\langle\hat{e}_{\bm{m}_{1}}^{\dagger}\hat{g}_{\bm{n}_{1}}^{\phantom{\dagger}}\hat{g}_{\bm{n}_{2}}^{\dagger}\hat{e}_{\bm{m}_{2}}\right\rangle^{\phantom{\dagger}}=\delta_{\bm{n}_{1}\bm{n}_{2}}\delta_{\bm{m}_{1}\bm{m}_{2}}(\zeta
p_{g_{\bm{n}}}+1)p_{e_{\bm{m}}},$ (26)

where $\sum_{\bm{n}}p_{g_{\bm{n}}}=\mathcal{N}_{g}$ and
$\sum_{\bm{n}}p_{e_{\bm{n}}}=\mathcal{N}_{e}$ are the number of excitations in
the ground and excited electronic bands, respectively, and $p_{g_{\bm{n}}}$
and $p_{e_{\bm{n}}}$ are the corresponding number distributions indexed by
trap state $\bm{n}$. The parameter $\zeta$ is defined with two values:
$\zeta=1$ for bosons and $\zeta=-1$ for fermions. These results can be
combined to find a general expression for the scaling of the radiated
intensity

$I(\bm{k}_{s})\propto\mathcal{N}_{e}+\zeta\sum_{\bm{n},\bm{m}}\lvert\eta_{\bm{n},\bm{m}}(\bm{k}_{s})\rvert^{2}p_{g_{\bm{n}}}p_{e_{\bm{m}}}.$
(27)

The first term in the expression of $I(\bm{k}_{s})$ comes from independent
emission, while the effect of quantum statistics and correlations comes from
the second term
$\sum_{n,m}\lvert\eta_{n,m}(\bm{k}_{s})\rvert^{2}p_{g_{\bm{n}}}p_{e_{\bm{m}}}$.
Some insight can be obtained from limiting cases. The infinite-temperature
case is trivial as the expected emission for either bosons or fermions is
fully independent and is proportional to $\mathcal{N}_{e}$, as expected from a
classical thermal gas, i.e.
$I_{\text{bos/ferm}}(\bm{k}_{s},T=\infty)\propto\mathcal{N}_{e}$ where all
particles emit independently. The zero-temperature case is much more
intriguing as we obtain the following expressions:

$\displaystyle I_{\text{bos}}(\bm{k}_{s},T=0)$
$\displaystyle\propto\mathrm{e}^{-\left(\bm{k}_{s}\cdot\bm{r}_{\text{zpm}}\right)^{2}}\mathcal{N}_{e}\left(\mathcal{N}_{g}+1\right),$
(28a) $\displaystyle I_{\text{ferm}}(\bm{k}_{s},T=0)$
$\displaystyle\propto\mathcal{N}_{e}-\sum_{\bm{n},\bm{m}}^{\text{Filling}}\lvert\eta_{\bm{n},\bm{m}}\rvert^{2},$
(28b)

for bosons and fermions respectively.

Let us first discuss the emission from a gas of bosons condensed onto the
ground state of the external trap (which is the case at zero temperature). The
Franck-Condon factor $\eta_{0,0}(\bm{k}_{s})$ in this case is easily estimated
and it includes the zero point motion $\bm{r}_{\text{zpm}}$ which indicates
the extent of the wavefunction in the ground state of the harmonic trap. For
deep trapping conditions, the Franck Condon factor is close to unity and the
emission shows the same features as observed in standard Dicke superradiance
in quantum optics. For sake of comparison, let us consider a fixed number of
particles $\mathcal{N}=\mathcal{N}_{e}+\mathcal{N}_{g}$. For an ensemble of
particles at the same spot, as assumed in the Dicke model, it is usually
convenient to use the set of collective Dicke states $\lvert J,m\rangle$, with
$J=\mathcal{N}/2$ and $m=(\mathcal{N}_{g}-\mathcal{N}_{e})/2$, situated on the
surface of the Bloch sphere. Superradiance refers then to a situation where
the emission of the collective state situated close to the equator shows a
proportionality to $J^{2}$. By considering deeply trapped particles instead,
the same kind of quick emission emerges as a fundamental property of the zero-
temperature bosonic wavefunction. This is to be expected, as both models
assume indistinguishable particles: in Dicke superradiance, this emerges from
the condition that all emitters are placed in the same spot while for the
bosonic cloud this is achieved by the symmetry of the ground state. The
scaling of the emitted intensity with the number of excited atoms is seen in
Fig. 4.

In the fermionic case, the summation limit “Filling” in the equation above
indicates that indices run from $(0,0)$ to the indices of the state at the
Fermi energy level. The summation in Eq. (28b) can be simplified in the case
of a deep trap where $\eta_{\bm{n},\bm{m}}\approx\delta_{\bm{n},\bm{m}}$ so
that
$I_{\text{ferm}}(\bm{k}_{s},T=0)\propto\mathcal{N}_{e}-\min(\mathcal{N}_{e},\mathcal{N}_{g})$.
This can be understood as a dynamic suppression of spontaneous emission into
occupied states due to the Pauli exclusion principle. For example, notice that
below the equator of the Bloch sphere all states are fully subradiant. This is
illustrated in Fig. 4 as the red, dashed curve and compared in Fig.4 against
the case of distinguishable particles (also reached for high temperatures)
where a linear scaling with the number of $\mathcal{N}_{e}$ is obtained.

Figure 4: Comparison between scattered intensity from distinguishable,
fermionic and bosonic systems in the deep trap limit where the equations Eq.
(28) apply, for $\mathcal{N}=4$. The intensity is shown as a function of the
excited particle fraction $\mathcal{N}_{e}/\mathcal{N}$. The blue dotted line
shows bosonic particles, the red dashed line fermionic particles and the solid
black line distinguishable particles (which is the case for a high temperature
ensemble).

##  V Experimental realization

In order to experimentally observe the effects discussed in the previous
sections, we have developed and built an apparatus to cool and trap neutral
dysprosium Mühlbauer _et al._ (2018). Dysprosium belongs to the group of
lanthanide elements whose characteristic open f-shell electron configuration
[Xe]4f106s2 with spin S = 2, orbital angular momentum L = 6, and total angular
momentum J = 8 gives rise to its high magnetic moment of 10 Bohr magnetons
($\mu\sim 10\mu_{B}$). Compared to alkali atoms, whose magnetic moments are on
the order of only $\mu\sim 1\mu_{B}$, the magnetic dipole-dipole interactions
in ultra-cold gases of dysprosium is about 100 times stronger, making it an
ideal candidate for investigating its contribution to cooperative effects. In
addition, there is an almost equal abundance of stable bosonic and fermionic
isotopes allowing for the creation of both types of quantum degenerate gases
Lu _et al._ (2011, 2012).

The general experimental scheme for producing laser cooled samples of
lanthanide atoms is implemented as follows Frisch _et al._ (2012); Maier _et
al._ (2014): A strong optical transition in the blue spectrum range is used to
precool the atoms in a Zeeman slower (ZS) before capturing them in a narrow
line magneto optical trap (MOT). Operating a MOT on transitions with natural
linewidths on the order of 100 kHz is required to reach Doppler temperatures
below 5$\mu$K, which allows one to directly transfer the atoms from the MOT
into an optical dipole trap (ODT). In our setup, we employ a strong J = 8
$\rightarrow$ J′ = 9 transition at 421 nm, with a natural linewidth of
$\Gamma_{421}=2\pi\cdot 32$ MHz for precooling. A thermal beam of atomic
dysprosium, transversally cooled on this broad transition, reaches the ZS with
initial velocities of several hundred meters per second. The atoms are then
longitudinally decelerated to a velocity of about 24 m/s in the spin-flip
configuration ZS before entering the main chamber, where they are captured in
a 6 beam 3D MOT setup. The MOT transition, on the other hand, is a closed,
narrow linewidth J = 8 $\rightarrow$ J′ = 9 transition at 626 nm with a
natural linewidth of $\Gamma_{626}=2\pi\cdot 136$ kHz which corresponds to a
Doppler temperature of 3.2 $\mu$K.

Our experimental setup consists of a high numerical aperture (NA) science cell
connected to an ultra-high vacuum (UHV) chamber (MOT) and a dedicated laser
system capable of generating tunable and frequency stabilized laser radiation.
Consequently, the laser frequency can be readily adapted to capture either
bosons or fermions without any optical or mechanical adjustments. Fig. 5 shows
a CAD render of the vacuum system. To reduce the influence of undesired
external magnetic fields, the vacuum system is assembled exclusively from non-
magnetic stainless steel, titanium and ceramic glass. Additionally, a
commercial, three axis magnetic field compensation system (Stefan Mayer
Instruments GmbH: MR-3) is used to drive a 3D coil system and compensate low
frequency magnetic field disturbances.

Figure 5: CAD render of the vacuum system. Dy atoms reach the ZS with a
velocity of few hundred meters per second. They are longitudinally slowed in
the ZS with resonant laser light at the 421 nm before being captured in a
narrow line 3D MOT at 626 nm. An achromatic lens in conjunction with an air-
bearing translation stage is used to focus down the optical transport beam at
the position of the MOT and create a deep ODT. The focal spot of this beam is
then moved from the MOT chamber to the center of the science cell by moving
the translation stage, thereby transferring the atoms alongside. In order to
retain a substantial number of atoms in the trap after this sequence, the
process of transport and re-trapping of atoms needs to be done on a timescale
that is shorter than the lifetime of the ODT.

The next step for creating dense samples of atomic dysprosium is to transfer
atoms to the science cell attached to the MOT chamber. The cell offers a high
optical access (with nine optical view ports), is designed in-house and made
out of the machinable glass-ceramic MACOR (Corning Inc.). A deep ODT at the
MOT position will be created, using a diode pumped solid-state (DPSS) laser
(Coherent: Mephisto MOPA 55 W) in order to transport the atoms. This will be
achieved by using an achromatic lens to focus down the transport beam to a
beam waist of $\omega_{0}\sim 66\mu m$ and two mirrors mounted on an air-
bearing translation stage (Aerotech: ABL1500-300) which enables reproducible
and precise movement of the focal spot of the beam, thereby, transporting the
atoms alongside Gustavson _et al._ (2001). The atoms are re-trapped inside the
science cell using a custom made multi-lens micro objective, designed to offer
diffraction-limited performance at the trapping wavelength of 1064 nmAlt
(2002); Li _et al._ (2018). A five-lens configuration consisting solely of
commercial singlets was chosen to achieve a high NA of 0.53 and a working
distance of 22 mm. Simulation and optimization of objective parameters like
lens curvature, relative spacing, and thickness was done using a ray tracing
software (Opticstudio) in order to minimize optical aberrations and achieve
diffraction-limited performance. The beam waist at the focus of the micro-
objective was measured using a piezo-controlled knife edge and found to be
$\omega_{0h}=5.94\pm 1.18\mu$m and $\omega_{0h}=6.99\pm 0.67\mu$m, in the
horizontal and vertical direction, respectively. We estimate that an ODT
created with this objective should provide a sufficiently tight confinement
for a few hundred to a few thousand dysprosium atoms to reach densities on the
order of $\sim 10^{13}$ atoms/cm3. This condition would allow then us to
explore dynamics in the high-density regime where the inter-particle/inter-
emitter spacing is lower than the wavelength of the exciting transition,
giving rise to effects like Dicke sub- and superradiance.

##  VI Conclusions

We have considered the emission and scattering properties of a dense atomic
cloud first by assuming distinguishable, classically moving particles and then
moving onto the quantum degenerate case, where quantum statistics play an
important role. In the case where atomic motion is treated classically, we
made use of an open quantum system approach where we have shown that tunable
spin Hamiltonians can be designed. In addition, under the assumption of weak
external drive, and in particular for atoms where magnetic interactions are
strong, we have characterized their effects onto the properties of scattered
light. For quantum degenerate gases, in the non-interacting regime, we have
quantified collective emission properties for both bosons and fermions. We
expect that our findings are of direct relevance for dense atomic dysprosium
clouds for which we describe an explicit experimental set up. In the future,
we will explicitly include the interaction of atoms via the electromagnetic
vacuum in order to derive an open system dynamics approach tailored to lossy
quantum degenerate gases.

## Acknowledgments

We acknowledge financial support from the Max Planck Society and from the
Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) – Project-ID
429529648 – TRR 306 QuCoLiMa (”Quantum Cooperativity of Light and Matter”).

## References

  * Labeyrie _et al._ (1999) G. Labeyrie, F. de Tomasi, J.-C. Bernard, C. A. Müller, C. Miniatura, and R. Kaiser, “Coherent backscattering of light by cold atoms,” Phys. Rev. Lett. 83, 5266–5269 (1999).
  * Labeyrie _et al._ (2003) G. Labeyrie, D. Delande, C. A. Müller, C. Miniatura, and R. Kaiser, “Coherent backscattering of light by an inhomogeneous cloud of cold atoms,” Phys. Rev. A 67, 033814 (2003).
  * Dicke (1954) R. H. Dicke, “Coherence in spontaneous radiation processes,” Phys. Rev. 93, 99–110 (1954).
  * Kwong _et al._ (2015) C. C. Kwong, T. Yang, D. Delande, R. Pierrat, and D. Wilkowski, “Cooperative emission of a pulse train in an optically thick scattering medium,” Phys. Rev. Lett. 115, 223601 (2015).
  * Bromley _et al._ (2016) S. L. Bromley, B. Zhu, M. Bishof, X. Zhang, T. Bothwell, J. Schachenmayer, T. L. Nicholson, R. Kaiser, S. F. Yelin, M. D. Lukin, A. M. Rey, and J. Ye, “Collective atomic scattering and motional effects in a dense coherent medium,” Nature Communications 7, 11039 (2016).
  * Inouye _et al._ (1999) S. Inouye, A. P. Chikkatur, D. M. Stamper-Kurn, J. Stenger, D. E. Pritchard, and W. Ketterle, “Superradiant rayleigh scattering from a bose-einstein condensate,” Science 285, 571–574 (1999).
  * Weiss _et al._ (2019) P. Weiss, A. Cipris, M. O. Araújo, R. Kaiser, and W. Guerin, “Robustness of dicke subradiance against thermal decoherence,” Phys. Rev. A 100, 033833 (2019).
  * Jennewein _et al._ (2018) S. Jennewein, L. Brossard, Y. R. P. Sortais, A. Browaeys, P. Cheinet, J. Robert, and P. Pillet, “Coherent scattering of near-resonant light by a dense, microscopic cloud of cold two-level atoms: Experiment versus theory,” Phys. Rev. A 97, 053816 (2018).
  * de Oliveira _et al._ (2014) R. A. de Oliveira, M. S. Mendes, W. S. Martins, P. L. Saldanha, J. W. R. Tabosa, and D. Felinto, “Single-photon superradiance in cold atoms,” Phys. Rev. A 90, 023848 (2014).
  * Roof _et al._ (2016) S. J. Roof, K. J. Kemp, M. D. Havey, and I. M. Sokolov, “Observation of single-photon superradiance and the cooperative lamb shift in an extended sample of cold atoms,” Phys. Rev. Lett. 117, 073003 (2016).
  * Araújo _et al._ (2016) M. O. Araújo, I. Krešić, R. Kaiser, and W. Guerin, “Superradiance in a large and dilute cloud of cold atoms in the linear-optics regime,” Phys. Rev. Lett. 117, 073002 (2016).
  * Guerin _et al._ (2016) W. Guerin, M. O. Araújo, and R. Kaiser, “Subradiance in a large cloud of cold atoms,” Phys. Rev. Lett. 116, 083601 (2016).
  * Gross and Haroche (1982) M. Gross and S. Haroche, “Superradiance: An essay on the theory of collective spontaneous emission,” Physics Reports 93, 301–396 (1982).
  * Scully _et al._ (2006) M. O. Scully, E. S. Fry, C. H. R. Ooi, and K. Wódkiewicz, “Directed spontaneous emission from an extended ensemble of $n$ atoms: Timing is everything,” Phys. Rev. Lett. 96, 010501 (2006).
  * A. _et al._ (2010) S. A. A., J.-T. Chang, and M. O. Scully, “Cooperative spontaneous emission of $n$ atoms: Many-body eigenstates, the effect of virtual lamb shift processes, and analogy with radiation of $n$ classical oscillators,” Phys. Rev. A 81, 053821 (2010).
  * Keaveney _et al._ (2012) J. Keaveney, A. Sargsyan, U. Krohn, I. G. Hughes, D. Sarkisyan, and C. S. Adams, “Cooperative lamb shift in an atomic vapor layer of nanometer thickness,” Phys. Rev. Lett. 108, 173601 (2012).
  * Javanainen _et al._ (2014) J. Javanainen, J. Ruostekoski, Y. Li, and S.-M. Yoo, “Shifts of a resonance line in a dense atomic sample,” Phys. Rev. Lett. 112, 113603 (2014).
  * Kaiser (2009) R. Kaiser, “Quantum multiple scattering,” Journal of Modern Optics 56, 2082–2088 (2009), https://doi.org/10.1080/09500340903082663 .
  * Zhu _et al._ (2016) B. Zhu, J. Cooper, J. Ye, and A. M. Rey, “Light scattering from dense cold atomic media,” Phys. Rev. A 94, 023612 (2016).
  * Javanainen and Ruostekoski (2016) J. Javanainen and J. Ruostekoski, “Light propagation beyond the mean-field theory of standard optics,” Opt. Express 24, 993–1001 (2016).
  * Pellegrino _et al._ (2014) J. Pellegrino, R. Bourgain, S. Jennewein, Y. R. P. Sortais, A. Browaeys, S. D. Jenkins, and J. Ruostekoski, “Observation of suppression of light scattering induced by dipole-dipole interactions in a cold-atom ensemble,” Phys. Rev. Lett. 113, 133602 (2014).
  * Jenkins _et al._ (2016) S. D. Jenkins, J. Ruostekoski, J. Javanainen, R. Bourgain, S. Jennewein, Y. R. P. Sortais, and A. Browaeys, “Optical resonance shifts in the fluorescence of thermal and cold atomic gases,” Phys. Rev. Lett. 116, 183601 (2016).
  * Lewenstein _et al._ (1994) M. Lewenstein, L. You, J. Cooper, and K. Burnett, “Quantum field theory of atoms interacting with photons: Foundations,” Phys. Rev. A 50, 2207–2231 (1994).
  * Ruostekoski and Javanainen (1997) J. Ruostekoski and J. Javanainen, “Quantum field theory of cooperative atom response: Low light intensity,” Phys. Rev. A 55, 513–526 (1997).
  * Javanainen _et al._ (1999) J. Javanainen, J. Ruostekoski, B. Vestergaard, and M. R. Francis, “One-dimensional modeling of light propagation in dense and degenerate samples,” Phys. Rev. A 59, 649–666 (1999).
  * Lee _et al._ (2016) M. D. Lee, S. D. Jenkins, and J. Ruostekoski, “Stochastic methods for light propagation and recurrent scattering in saturated and nonsaturated atomic ensembles,” Phys. Rev. A 93, 063803 (2016).
  * Childs (1970) W. J. Childs, “Hyperfine structure of ${}^{5}i_{8,7}$ atomic states of ${\mathrm{dy}}^{161,163}$ and the ground-state nuclear moments,” Phys. Rev. A 2, 1692–1701 (1970).
  * Petersen _et al._ (2020) N. Petersen, M. Trümper, and P. Windpassinger, “Spectroscopy of the 1001-nm transition in atomic dysprosium,” Phys. Rev. A 101, 042502 (2020).
  * James (1993) D. F. V. James, “Frequency shifts in spontaneous emission from two interacting atoms,” Phys. Rev. A 47, 1336–1346 (1993).
  * Lehmberg (1970) R. H. Lehmberg, “Radiation from an $n$-atom system. i. general formalism,” Phys. Rev. A 2, 883–888 (1970).
  * Mühlbauer _et al._ (2018) F. Mühlbauer, N. Petersen, C. Baumgärtner, L. Maske, and P. Windpassinger, “Systematic optimization of laser cooling of dysprosium,” Applied Physics B 124, 120 (2018).
  * Lu _et al._ (2011) M. Lu, N. Q. Burdick, S. H. Youn, and B. L. Lev, “Strongly dipolar bose-einstein condensate of dysprosium,” Phys. Rev. Lett. 107, 190401 (2011).
  * Lu _et al._ (2012) M. Lu, N. Q. Burdick, and B. L. Lev, “Quantum degenerate dipolar fermi gas,” Phys. Rev. Lett. 108, 215301 (2012).
  * Frisch _et al._ (2012) A. Frisch, K. Aikawa, M. Mark, A. Rietzler, J. Schindler, E. Zupanič, R. Grimm, and F. Ferlaino, “Narrow-line magneto-optical trap for erbium,” Phys. Rev. A 85, 051401 (2012).
  * Maier _et al._ (2014) T. Maier, H. Kadau, M. Schmitt, A. Griesmaier, and T. Pfau, “Narrow-line magneto-optical trap for dysprosium atoms,” Opt. Lett. 39, 3138–3141 (2014).
  * Gustavson _et al._ (2001) T. L. Gustavson, A. P. Chikkatur, A. E. Leanhardt, A. Görlitz, S. Gupta, D. E. Pritchard, and W. Ketterle, “Transport of bose-einstein condensates with optical tweezers,” Phys. Rev. Lett. 88, 020401 (2001).
  * Alt (2002) W. Alt, “An objective lens for efficient fluorescence detection of single atoms,” Optik 113, 142–144 (2002).
  * Li _et al._ (2018) X. Li, F. Zhou, M. Ke, P. Xu, X.-D. He, J. Wang, and M.-S. Zhan, “High-resolution ex vacuo objective for cold atom experiments,” Appl. Opt. 57, 7584–7590 (2018).

##  Appendix A Magnetic interaction terms in light matter Hamiltonian from
first principles

We will be using the following definitions for the gauge, electric and
magnetic field operators

$\displaystyle\hat{\bm{A}}$
$\displaystyle=\sum_{\bm{k},\bm{\epsilon}}g_{k}\epsilon_{\bm{k}}\left(\hat{a}_{\bm{k},\bm{\epsilon}}\mathrm{e}^{\mathrm{i}\bm{k}\bm{R}}+\hat{a}_{\bm{k},\bm{\epsilon}}^{\dagger}\mathrm{e}^{-\mathrm{i}\bm{k}\bm{R}}\right)$
(29) $\displaystyle\hat{\bm{E}}$
$\displaystyle=\mathrm{i}\sum_{\bm{k},\bm{\epsilon}}g_{k}\omega_{k}\epsilon_{\bm{k}}\left(\hat{a}_{\bm{k},\bm{\epsilon}}\mathrm{e}^{\mathrm{i}\bm{k}\bm{R}}-\hat{a}_{\bm{k},\bm{\epsilon}}^{\dagger}\mathrm{e}^{-\mathrm{i}\bm{k}\bm{R}}\right)$
(30) $\displaystyle\hat{\bm{B}}$
$\displaystyle=\mathrm{i}\sum_{\bm{k},\bm{\epsilon}}\left(\bm{k}\times\epsilon_{\bm{k}}\right)g_{k}\left(\hat{a}_{\bm{k},\bm{\epsilon}}\mathrm{e}^{\mathrm{i}\bm{k}\bm{R}}-\hat{a}_{\bm{k},\bm{\epsilon}}^{\dagger}\mathrm{e}^{-\mathrm{i}\bm{k}\bm{R}}\right).$
(31)

with $g_{k}=1/\sqrt{2\omega_{k}\mathcal{V}\epsilon_{0}}$ photon coupling
strength with photonic operators as defined in the main text.

If we omit terms where the square of the gauge field appears and the term
where the atom core momentum is coupled to the gauge field (amounting to a
Born-Oppenheimer and subsequently a rotating wave approximation), the
simplified Hamiltonian is composed of three parts, written in the Coulomb
gauge

$\displaystyle\mathcal{H}$
$\displaystyle\approx\mathcal{H}_{0}+\mathcal{H}_{\text{af}}+\mathcal{H}_{\text{dipole}}$
(32) $\displaystyle\mathcal{H}_{0}$
$\displaystyle=\sum_{i}^{N}\frac{\hat{\bm{P}}_{i}^{2}}{2M}+\frac{\hat{\bm{p}}_{i}^{2}}{2m}+V_{ec}(\bm{R}_{i},\bm{r}_{i})$
(33) $\displaystyle\mathcal{H}_{af}$
$\displaystyle=\frac{1}{m}\sum_{i}^{N}\hat{\bm{A}}(\bm{r}_{i})\cdot\hat{\bm{p}}_{i}$
(34) $\displaystyle\mathcal{H}_{\text{dipole}}$
$\displaystyle=\frac{1}{2}\sum_{i\neq
j}\sum\frac{\bm{d}_{i}\cdot\bm{d}_{j}-3(\bm{R}_{ij}\bm{d}_{i})(\bm{R}_{ij}\cdot\bm{d}_{j})}{R^{3}_{ij}}.$
(35)

where we have already performed a dipole approximation for the Coulomb
interaction between different atoms leading to the explicit form of
$\mathcal{H}_{\text{dipole}}$. The Hamiltonian $\mathcal{H}_{af}$ is the
canonical coupling between the electronic degree of freedom and the transverse
modes of the light-field. We have not completed the dipolar approximation in
the gauge field which still contains the position of the electron $\bm{r}_{i}$
which is a quantum operator in this description. The Hamiltonian
$\mathcal{H}_{0}$ contains the electron kinetic energy, the center of mass
kinetic energy and the interaction between the electron and the core. The
electron kinetic energy and the core-electron potential give rise to the level
structure of the atom. We consider a Born-Oppenheimer picture where
$\bm{R}_{i}$ represents the center of mass momentum.

Instead of making the usual dipolar approximation for the gauge field, we
consider one order higher in the Taylor approximation of the electron position
around the center of mass position
$\mathrm{e}^{\mathrm{i}\bm{k}\bm{r}_{i}}\approx\mathrm{e}^{\mathrm{i}\bm{k}\bm{R}_{i}}\left[1+\mathrm{i}\bm{k}\bm{x}_{i}\right]$
with the distance vector $\bm{x}_{i}=\bm{r}_{i}-\bm{R}_{i}$. The atom field
Hamiltonian then splits into two contributions

$\begin{split}\mathcal{H}_{af}&=\frac{1}{m}\sum_{i}^{N}\hat{\bm{p}}_{i}\cdot\hat{\bm{A}}(\bm{R}_{i})+\mathrm{i}\sum_{\bm{k},\bm{\epsilon}}g_{k}\left(\hat{\bm{p}}_{i}\cdot\epsilon_{\bm{k}}\right)(\bm{k}\bm{x}_{i})\left(\hat{a}_{\bm{k},\bm{\epsilon}}\mathrm{e}^{\mathrm{i}\bm{k}\bm{R}_{i}}-\hat{a}_{\bm{k},\bm{\epsilon}}^{\dagger}\mathrm{e}^{-\mathrm{i}\bm{k}\bm{R}_{i}}\right).\end{split}$
(36)

The first term is the simple dipole approximation which leads to the coupled
dipole model after some additional approximations and a gauge transformation.
The second term is of higher order and contains the magnetic and quadrupole
electric moment as we will see. In order to see that the double dot product
indeed yields the magnetic and quadrupole term, it needs to be rearranged

$\begin{split}(\hat{\bm{p}}_{i}\cdot\bm{\epsilon}_{\bm{k}})(\bm{k}\cdot\bm{x}_{i})&=(\hat{\bm{p}}_{i}\cdot\bm{k})(\bm{x}_{i}\cdot\bm{\epsilon}_{\bm{k}})+(\bm{x}_{i}\times\hat{\bm{p}}_{i})(\bm{k}\times\epsilon_{\bm{k}})\\\
&=\frac{1}{2}\left((\hat{\bm{p}}_{i}\cdot\bm{\epsilon}_{\bm{k}})(\bm{k}\cdot\bm{x}_{i})+(\hat{\bm{p}}_{i}\cdot\bm{k})(\bm{x}_{i}\cdot\bm{\epsilon}_{\bm{k}})\right)+\frac{1}{2}\hat{\bm{L}}(\bm{k}\times\epsilon_{\bm{k}})\end{split}$
(37)

The first term is the electronic quadrupole term whereas the second term
contains the angular momentum implying that it is indeed the magnetic
coupling. Summarizing the calculations above, the terms that must be added to
the Hamiltonian if the dipole approximation is performed to higher order are

$\begin{split}\mathcal{H}_{m}&=\sum_{i}^{N}\frac{\bm{\hat{L}}_{i}}{2m}\mathrm{i}\sum_{\bm{k},\bm{\epsilon}}g_{k}(\bm{k}\times\bm{\epsilon}_{\bm{k}})\left(\hat{a}_{\bm{k},\bm{\epsilon}}\mathrm{e}^{\mathrm{i}\bm{k}\bm{R}_{i}}-\hat{a}_{\bm{k},\bm{\epsilon}}^{\dagger}\mathrm{e}^{-\mathrm{i}\bm{k}\bm{R}_{i}}\right)=\sum_{i}^{N}\hat{\bm{\mu}}_{i}\hat{\bm{B}}(\bm{R}_{i})\\\
\mathcal{H}_{q}&=\sum_{i,\bm{k},\bm{\epsilon}}\frac{\mathrm{i}\lvert\bm{k}\rvert}{2m}(p_{\bm{k}}r_{\bm{\epsilon}}+p_{\bm{\epsilon}}r_{\bm{k}})g_{k}\left(\hat{a}_{\bm{k},\bm{\epsilon}}\mathrm{e}^{\mathrm{i}\bm{k}\bm{R}_{i}}-\hat{a}_{\bm{k},\bm{\epsilon}}^{\dagger}\mathrm{e}^{-\mathrm{i}\bm{k}\bm{R}_{i}}\right)\end{split}$
(38)

Here we avoid complications with spin degrees of freedom by making the
association $\hat{\bm{\mu}}_{i}=-\frac{\bm{\hat{L}}_{i}}{2m}$ and noting that
spin degrees of freedom need to be explicitly inserted into the Hamiltonian
which ultimately change the definition of $\hat{\bm{\mu}}_{i}$. Since we go
into a two level picture this is however already taken care of and we need not
consider these complications in detail.

##  Appendix B Representation of the Hamiltonian in the two-level
approximation

We perform a two-level approximation in which we assume there exist two
relevant levels per atom, the ground state $\lvert g_{i}\rangle$ and the
excited state $\lvert e_{i}\rangle$ separated by a frequency $\omega_{0}$.
Taking the projection of the momentum and magnetic moment operator onto this
subspace, we can explicilty write the Hamiltonian using the matrix elements
within this subspace. One relevant point is here that the momentum operator
has no diagonal components whereas the magnetic moment operator has no such
symmetry restriction.

Omitting the quadrupole term, the atom-field Hamiltonian then becomes

$\begin{split}\mathcal{H}_{\text{af}}&=\mathrm{i}\sum_{i}^{N}\left(\bm{d}_{i}\sigma_{i}-\bm{d}^{*}_{i}\sigma_{i}^{\dagger}\right)\hat{\bm{A}}(\bm{R}_{i})+\left(\bm{\mu}_{i}\sigma_{i}+\bm{\mu}^{*}_{i}\sigma_{i}^{\dagger}+\bm{\mu}_{i,e}\sigma^{\dagger}\sigma+\bm{\mu}_{i,g}\sigma\sigma^{\dagger}\right)\hat{\bm{B}}(\bm{R}_{i}).\end{split}$
(39)

with $\bm{d}_{i}=\left\langle
g_{i}\middle|\hat{x}_{i}\middle|e_{i}\right\rangle$ and
$\bm{\mu}_{i,e}=\left\langle
g_{i}\middle|\hat{\bm{\mu}}_{i}\middle|e_{i}\right\rangle$ for the transition
dipoles and $\bm{\mu}_{i,e}=\left\langle
e_{i}\middle|\hat{\bm{\mu}}_{i}\middle|e_{i}\right\rangle$,
$\bm{\mu}_{i,g}=\left\langle
g_{i}\middle|\hat{\bm{\mu}}_{i}\middle|g_{i}\right\rangle$ for the static
magnetic dipole moment. The sigma matrices are transition operators of the
type $\sigma_{i}=\lvert g_{i}\rangle\langle e_{i}\rvert$.

##  Appendix C Unitary transformation into dipolar gauge

The next step in deriving the light-matter Hamiltonian is the gauge
transformation

$U=\exp\left(-\mathrm{i}\sum_{i}^{N}\int_{\bm{r}_{i}}^{\bm{R}_{i}}\mathop{}\\!\mathrm{d}{\bm{r}}\hat{\bm{A}}(\bm{r})\right)$
(40)

We note here that the gauge field operator commutes with the magnetic field
operator $\left[\hat{\bm{A}}(\bm{r}),\hat{\bm{B}}(\bm{r}^{\prime})\right]=0$
so that $H_{m}$ is not affected by this transformation. The effect on the
dipole electronic transition and the static dipole
$\mathcal{H}_{\text{dipole}}$ is as usual so that, after also performing a
rotating wave approximation, the Hamiltonian is

$\begin{split}\mathcal{H}&=\omega\sum_{i}\sigma_{i}^{\dagger}\sigma_{i}+\sum_{\bm{k},\bm{\epsilon}}\omega_{k}\hat{a}_{\bm{k},\bm{\epsilon}}^{\dagger}\hat{a}_{\bm{k},\bm{\epsilon}}+\sum_{i}\hat{\bm{d}}_{i}\cdot\hat{\bm{E}}(\bm{R}_{i})+\left(\hat{\bm{\mu}}_{t,i}+\hat{\bm{\mu}}_{s,i}\right)\cdot\hat{\bm{B}}(\bm{R}_{i}).\end{split}$
(41)

with definitions
$\hat{\bm{d}}_{j}=\bm{d}\sigma_{j}+\bm{d}^{*}\sigma_{j}^{\dagger}$,
$\hat{\bm{\mu}}_{t,i}=\bm{\mu}_{i}\sigma_{i}+\bm{\mu}^{*}_{i}\sigma_{i}^{\dagger}$
and
$\hat{\bm{\mu}}_{s,i}=\bm{\mu}_{i,e}\sigma^{\dagger}\sigma+\bm{\mu}_{i,g}\sigma\sigma^{\dagger}$.

##  Appendix D Elimination of photonic degrees of freedom from Hamiltonian

The last step in determining the effective system dynamics is to perform an
adiabatic elimination of the light modes for this Hamiltonian Eq. (41). The
first step is to determine the equation of motion for the photonic operators

$\begin{split}\frac{\mathop{}\\!\mathrm{d}\hat{a}_{\bm{k},\bm{\epsilon}}}{\mathop{}\\!\mathrm{d}t}&=-\mathrm{i}\omega_{k}\hat{a}_{\bm{k},\bm{\epsilon}}+\sum_{i}g_{k}\omega_{k}(\hat{\bm{d}}_{i}\cdot\bm{\epsilon}_{\bm{k}})\mathrm{e}^{-\mathrm{i}\bm{k}\hat{\bm{R}}_{i}}+\sum_{i}g_{k}\left(\left[\hat{\bm{\mu}}_{t,i}+\hat{\bm{\mu}}_{s,i}\right]\cdot\left[\bm{k}\times\bm{\epsilon}_{\bm{k}}\right]\right)\mathrm{e}^{-\mathrm{i}\bm{k}\hat{\bm{R}}_{i}}\end{split}$
(42)

which can then be integrated to yield

$\begin{split}\hat{a}_{\bm{k},\bm{\epsilon}}(t)&=\hat{a}_{\bm{k},\bm{\epsilon}}(0)\mathrm{e}^{-\mathrm{i}\omega_{k}t}+\int_{0}^{t}\mathop{}\\!\mathrm{d}{s}\mathrm{e}^{-\mathrm{i}\omega_{k}(t-s)}\sum_{i}g_{k}\omega_{k}(\hat{\bm{d}}_{i}(s)\cdot\bm{\epsilon}_{\bm{k}})\mathrm{e}^{-\mathrm{i}\bm{k}\hat{\bm{R}}_{i}}\\\
&+\int_{0}^{t}\mathop{}\\!\mathrm{d}{s}\mathrm{e}^{-\mathrm{i}\omega_{k}(t-s)}\sum_{i}g_{k}\left(\left[\hat{\bm{\mu}}_{t,i}(s)+\hat{\bm{\mu}}_{s,i}(s)\right]\cdot\left[\bm{k}\times\bm{\epsilon}_{\bm{k}}\right]\right)\mathrm{e}^{-\mathrm{i}\bm{k}\hat{\bm{R}}_{i}}.\end{split}$
(43)

A simplification for the calculation that follows is to consider negative and
positive frequency components of the magnetic and electric field in particular

$\begin{split}\hat{\bm{B}}^{\pm}&=\pm\mathrm{i}\sum_{\bm{k},\bm{\epsilon}}\left(\bm{k}\times\epsilon_{\bm{k}}\right)g_{k}\hat{a}_{\bm{k},\bm{\epsilon}}\mathrm{e}^{\pm\mathrm{i}\bm{k}\bm{R}}\\\
\hat{\bm{E}}^{\pm}&=\pm\mathrm{i}\sum_{\bm{k},\bm{\epsilon}}g_{k}\omega_{k}\epsilon_{\bm{k}}\hat{a}_{\bm{k},\bm{\epsilon}}\mathrm{e}^{\pm\mathrm{i}\bm{k}\bm{R}}\end{split}$
(44)

so that we only have to consider the equation of motion for one of the
photonic operators. To shorten the notation we also introduce the total
magnetic dipole operator
$\hat{\bm{\mu}}_{i}(t)=\hat{\bm{\mu}}_{t,i}(t)+\hat{\bm{\mu}}_{s,i}(t)$
leading to expressions for the positive field components of the electric and
magnetic field

$\begin{split}\hat{\bm{E}}^{+}(\bm{R}_{i})&=\hat{\bm{E}}_{0}^{+}(\bm{R}_{i})+\mathrm{i}\sum_{\bm{k},\bm{\epsilon},j}(g_{k}\omega_{k})^{2}\mathrm{e}^{\mathrm{i}\bm{k}\bm{R}_{ij}}\epsilon_{\bm{k}}\int_{0}^{t}\mathop{}\\!\mathrm{d}{s}\mathrm{e}^{-\mathrm{i}\omega_{k}(t-s)}(\hat{\bm{d}}_{j}(s)\cdot\bm{\epsilon}_{\bm{k}})\\\
&+\mathrm{i}\sum_{\bm{k},\bm{\epsilon},j}g^{2}_{k}\omega_{k}\mathrm{e}^{\mathrm{i}\bm{k}\bm{R}_{ij}}\epsilon_{\bm{k}}\int_{0}^{t}\mathop{}\\!\mathrm{d}{s}\mathrm{e}^{-\mathrm{i}\omega_{k}(t-s)}\left(\hat{\bm{\mu}}_{i}(s)\cdot\left[\bm{k}\times\bm{\epsilon}_{\bm{k}}\right]\right)\\\
\hat{\bm{B}}^{+}(\bm{R}_{i})&=\hat{\bm{B}}_{0}^{+}(\bm{R}_{i})+\mathrm{i}\sum_{\bm{k},\bm{\epsilon},j}\left(\bm{k}\times\epsilon_{\bm{k}}\right)g^{2}_{k}\omega_{k}\mathrm{e}^{\mathrm{i}\bm{k}\bm{R}_{ij}}\int_{0}^{t}\mathop{}\\!\mathrm{d}{s}\mathrm{e}^{-\mathrm{i}\omega_{k}(t-s)}(\hat{\bm{d}}_{j}(s)\cdot\bm{\epsilon}_{\bm{k}})\\\
&+\mathrm{i}\sum_{\bm{k},\bm{\epsilon},j}\left(\bm{k}\times\epsilon_{\bm{k}}\right)g^{2}_{k}\mathrm{e}^{\mathrm{i}\bm{k}\bm{R}_{ij}}\int_{0}^{t}\mathop{}\\!\mathrm{d}{s}\mathrm{e}^{-\mathrm{i}\omega_{k}(t-s)}\left(\hat{\bm{\mu}}_{i}(s)\cdot\left[\bm{k}\times\bm{\epsilon}_{\bm{k}}\right]\right)\end{split}$
(45)

where the zero subscript indicates that the photonic operators in the field
expression were taken at time zero. Thus, both fields have similar terms.
First, the vacuum time evolution, then the interaction with the proper dipole
moment and then a cross term between the magnetic and electric components of
the interaction.

Using the usual results for sums over polarization vectors since they form a
orthonormal basis with the normalized wavevector, the cross terms can be
combined into

$\begin{split}\sum_{\bm{\epsilon}}\left((\hat{\bm{d}}_{i}(t)\cdot\bm{\epsilon}_{\bm{k}})\hat{\bm{\mu}}_{i}(s)\cdot\left[\bm{k}\times\bm{\epsilon}_{\bm{k}}\right]+\hat{\bm{\mu}}_{i}(t)\cdot\left[\bm{k}\times\bm{\epsilon}_{\bm{k}}\right](\hat{\bm{d}}_{j}(s)\cdot\bm{\epsilon}_{\bm{k}})\right)=\left(\hat{\bm{\mu}}_{i}(t)\times\hat{\bm{d}}_{i}(s)-\hat{\bm{d}}_{i}(t)\times\hat{\bm{\mu}}_{i}(s)\right)\cdot\bm{k}\end{split}$
(46)

which is linear in $\bm{k}$, such that the negative components of the sum of
$\hat{\bm{E}}^{+}(\bm{R}_{i})$ cancel the positive wave vectors in
$\hat{\bm{E}}^{-}(\bm{R}_{i})$ since the dipole operator is hermitian and the
polarization has been chosen real. The same goes for the magnetic field,
canceling the cross terms.

After the cross-terms are cancelled, the total magnetic and electric field
operators can be obtained from the positive frequency components by adding the
hermitian conjugate so that after performing the usual summation over the
polarization vectors

$\displaystyle\hat{\bm{E}}(\bm{R}_{i})$
$\displaystyle=\hat{\bm{E}}_{0}(\bm{R}_{i})+\mathrm{i}\sum_{\bm{k},j}(g_{k}\omega_{k})^{2}\mathrm{e}^{\mathrm{i}\bm{k}\bm{R}_{ij}}\int_{0}^{t}\mathop{}\\!\mathrm{d}{s}\mathrm{e}^{-\mathrm{i}\omega_{k}(t-s)}\cdot\left[\hat{\bm{d}}_{j}(s)-\frac{\bm{k}(\hat{\bm{d}}_{j}(s)\cdot\bm{k})}{k^{2}}\right]+\mathrm{h.c}$
(47) $\displaystyle\hat{\bm{B}}(\bm{R}_{i})$
$\displaystyle=\hat{\bm{B}}_{0}(\bm{R}_{i})+\mathrm{i}\sum_{\bm{k},j}(g_{k}\omega_{k})^{2}\mathrm{e}^{\mathrm{i}\bm{k}\bm{R}_{ij}}\int_{0}^{t}\mathop{}\\!\mathrm{d}{s}\mathrm{e}^{-\mathrm{i}\omega_{k}(t-s)}\cdot\left[\hat{\bm{\mu}}_{i}(s)-\frac{\bm{k}(\hat{\bm{\mu}}_{i}(s)\cdot\bm{k})}{k^{2}}\right]+\mathrm{h.c}$
(48)

so that the magnetic and electric fields have the same structure. Now, the
continuum limit is taken for the electromagnetic field modes
$\frac{1}{V}\sum_{\bm{k}}\rightarrow\frac{1}{(2\pi)^{3}}\int\mathop{}\\!\mathrm{d}{\bm{k}}$
and the angular integral is performed in spherical coordinates so that,
defining the dipolar Green’s function

$\begin{split}\mathbf{G}_{k}(\mathbf{R})&=\left(\mathds{1}+\frac{1}{k^{2}}\nabla\otimes\nabla\right)\frac{\mathrm{e}^{\mathrm{i}kR}}{4\pi
R}\\\ &=\frac{\mathrm{e}^{\mathrm{i}kR}}{4\pi
k^{2}}\left[\left(\frac{k^{2}}{R}+\frac{\mathrm{i}k}{R^{2}}-\frac{1}{R^{3}}\right)\mathds{1}+\left(-\frac{k^{2}}{R}-\frac{3\mathrm{i}k}{R^{2}}+\frac{3}{R^{3}}\right)\frac{\mathbf{R}\otimes\mathbf{R}}{R^{2}}\right]\end{split}$
(49)

with the usual shorthand $\bm{G}(\bm{R})=\bm{G}_{k_{0}}(\bm{R})$ and the
decomposition
$\bm{G}_{k}(\bm{R})=\bm{\Omega}_{k}(\bm{R})-\mathrm{i}\bm{\Gamma}_{k}(\bm{R})$.
The expressions for the electric and magnetic field become

$\displaystyle\hat{\bm{E}}(\bm{R}_{i})$
$\displaystyle=\hat{\bm{E}}_{0}(\bm{R}_{i})+\frac{\mathrm{i}}{\pi\epsilon_{0}}\sum_{j}\int_{0}^{\infty}\mathop{}\\!\mathrm{d}{k}\int_{0}^{t}\mathop{}\\!\mathrm{d}{s}k^{2}\bm{\Gamma}_{\bm{k}}(\bm{R}_{ij})\cdot\hat{\bm{d}}_{j}(s)\mathrm{e}^{-\mathrm{i}\omega_{k}(t-s)}+\mathrm{h.c}$
(50) $\displaystyle\hat{\bm{B}}(\bm{R}_{i})$
$\displaystyle=\hat{\bm{B}}_{0}(\bm{R}_{i})+\frac{\mathrm{i}}{\pi\epsilon_{0}}\sum_{j}\int_{0}^{\infty}\mathop{}\\!\mathrm{d}{k}\int_{0}^{t}\mathop{}\\!\mathrm{d}{s}k^{2}\bm{\Gamma}_{\bm{k}}(\bm{R}_{ij})\cdot\hat{\bm{\mu}}_{j}(s)\mathrm{e}^{-\mathrm{i}\omega_{k}(t-s)}+\mathrm{h.c}.$
(51)

The key is now that this elimination cannot be performed exactly and one must
make a perturbative Ansatz by plugging the time evolution of the dipole
operators without light-matter interactions into these integrals. The free-
space time evolution of the dipole operators is explicitly

$\begin{split}\hat{\bm{d}}_{j}(s)&=\bm{d}\sigma_{j}(s)+\bm{d}^{*}\sigma_{j}^{\dagger}=\bm{d}\sigma_{j}(t)\mathrm{e}^{-\mathrm{i}\omega(s-t)}+\bm{d}^{*}\sigma_{j}(t)^{\dagger}\mathrm{e}^{\mathrm{i}\omega(s-t)}\\\
\hat{\bm{\mu}}_{t,j}(s)&=\bm{\mu}_{j}\sigma_{j}(s)+\bm{\mu}^{*}_{j}\sigma_{j}^{\dagger}(s)=\bm{\mu}_{j}\sigma_{j}(t)\mathrm{e}^{-\mathrm{i}\omega(s-t)}+\bm{\mu}^{*}_{j}\sigma_{j}^{\dagger}(t)\mathrm{e}^{\mathrm{i}\omega(s-t)}\\\
\hat{\bm{\mu}}_{s,j}(s)&=\bm{\mu}_{j,e}\sigma_{j}^{\dagger}(s)\sigma_{j}(s)+\bm{\mu}_{j,g}\sigma_{j}(s)\sigma_{j}^{\dagger}(s)=\hat{\bm{\mu}}_{s,j}(t).\end{split}$
(52)

This implies that the fundamental difference between the static and transition
dipole moments is given by the frequency they rotate at. All these time
evolutions will lead to integrations of the type

$\zeta(\omega)=\mathrm{i}\int_{0}^{\infty}\mathrm{e}^{\mathrm{i}\omega
t}\mathop{}\\!\mathrm{d}{t}=\mathcal{P}\frac{1}{\omega}+\mathrm{i}\pi\delta(\omega)$
(53)

so that, splitting the magnetic field into a transition part
$\hat{\bm{B}}=\hat{\bm{B}}_{t}+\hat{\bm{B}}_{s}$ while keeping the vacuum
contribution in the transition part

$\begin{split}\hat{\bm{E}}(\bm{R}_{i})&=\hat{\bm{E}}_{0}(\bm{R}_{i})+\frac{1}{\pi\epsilon_{0}}\sum_{j}\int_{0}^{\infty}\mathop{}\\!\mathrm{d}{k}k^{2}\bm{\Gamma}_{k}(\bm{R}_{ij})\cdot\bm{d}\left(\sigma_{j}\zeta(\omega_{k}+\omega)+\sigma_{j}^{\dagger}\zeta(\omega_{k}-\omega)\right)+\mathrm{h.c}\\\
\hat{\bm{B}}_{t}(\bm{R}_{i})&=\hat{\bm{B}}_{0}(\bm{R}_{i})+\frac{1}{\pi\epsilon_{0}}\sum_{j}\int_{0}^{\infty}\mathop{}\\!\mathrm{d}{k}k^{2}\bm{\Gamma}_{k}(\bm{R}_{ij})\cdot\bm{\mu}_{t,j}\left(\sigma_{j}\zeta(\omega_{k}+\omega)+\sigma_{j}^{\dagger}\zeta(\omega_{k}-\omega)\right)+\mathrm{h.c}\\\
\hat{\bm{B}}_{s}(\bm{R}_{i})&=\frac{1}{\pi\epsilon_{0}}\sum_{j}\int_{0}^{\infty}\mathop{}\\!\mathrm{d}{k}k^{2}\bm{\Gamma}_{k}(\bm{R}_{ij})\cdot\bm{\zeta}(\omega_{k})\left(\bm{\mu}_{j,e}\sigma_{j}^{\dagger}\sigma_{j}+\bm{\mu}_{j,g}\sigma_{j}\sigma_{j}^{\dagger}\right)+\mathrm{h.c}.\end{split}$
(54)

The previous statement about the rotational frequencies now becomes clear as
the transition fields contain zeta functions of the form
$\zeta(\omega_{k}\pm\omega)$ where the Green’s function around
$\bm{\Gamma}_{k_{0}}(\bm{R})$ matters. The static contributions however
contain a zeta function of the form $\zeta(\omega_{k})$ where
$\bm{\Gamma}_{k=0}(\bm{R})$ matters. Since the zeta function essentially only
evaluates the Kramers Kronig relations for the complex function
$\bm{G}(\bm{R})$, i.e. $(k^{\prime
2})\Omega_{k^{\prime}}(\bm{R})=\int\mathop{}\\!\mathrm{d}{k}\frac{k^{2}\Gamma_{k}(\bm{R})}{\omega^{\prime}-\omega}$.
Ultimately, this leads to the explicit representation of the fields, using the
fact that $\lim_{k\to 0}k^{2}G_{k}(\bm{R})=\frac{1}{4\pi
R_{ij}^{3}}\left(3\frac{\bm{R}_{ij}\bm{R}_{ij}}{R_{ij}R_{ij}}-1\right)$ so
that only the real part remains

$\begin{split}\hat{\bm{E}}(\bm{R}_{i})&=\hat{\bm{E}}_{0}(\bm{R}_{i})+\frac{k_{0}^{2}}{\epsilon_{0}}\sum_{j}\left[\left(\bm{\Gamma}(\bm{R}_{ij})-\mathrm{i}\bm{\Omega}(\bm{R}_{ij})\right)\sigma_{j}+\left(\bm{\Gamma}(\bm{R}_{ij})+\mathrm{i}\bm{\Omega}(\bm{R}_{ij})\right)\sigma_{j}^{\dagger}\right]\bm{d}\\\
\hat{\bm{B}}(\bm{R}_{i})&=\hat{\bm{B}}_{0}(\bm{R}_{i})+\frac{k_{0}^{2}}{\epsilon_{0}}\sum_{j}\left[\left(\bm{\Gamma}(\bm{R}_{ij})-\mathrm{i}\bm{\Omega}(\bm{R}_{ij})\right)\sigma_{j}+\left(\bm{\Gamma}(\bm{R}_{ij})+\mathrm{i}\bm{\Omega}(\bm{R}_{ij})\right)\sigma_{j}^{\dagger}\right]\bm{\mu}_{t}\\\
\hat{\bm{B}}_{s}(\bm{R}_{i})&=\frac{\mu_{0}}{4\pi}\sum_{j}\frac{1}{R_{ij}^{3}}\left(3\frac{\bm{R}_{ij}\bm{R}_{ij}}{R_{ij}R_{ij}}-1\right)\hat{\bm{\mu}_{s,j}}.\end{split}$
(55)

The next step is to plug these expressions into the equations of motion, to
use the rotating wave approximation and the definitions
$\Gamma=k_{0}^{3}d^{2}/(3\pi\hbar\epsilon_{0})$ and
$\Gamma_{\mu}=k_{0}^{3}\mu^{2}/3\pi\hbar\epsilon_{0}$ in order to define the
Green’s function for the matter problem. Indeed, the definitions
$f^{(\mu)}_{ij}=\frac{k_{0}^{2}}{\epsilon_{0}}\bm{\mu}_{i,t}\bm{\Omega}(\bm{R}_{ij})\bm{\mu}_{j,t}$,
$f^{(g)}_{ij}=\frac{k_{0}^{2}}{\epsilon_{0}}\bm{d}_{i}\bm{\Omega}(\bm{R}_{ij})\bm{d}_{j}$,
equivalent definitions for the imaginary part and Eq. (13) lead precisely to
Eq. (12), Eq. (14) and Eq. (15) after using the rotating wave approximation.

After adding a classical drive just as in the main text, the equations of
motion for the transition operators become

$\begin{split}\frac{\mathop{}\\!\mathrm{d}\sigma_{k}}{\mathop{}\\!\mathrm{d}t}&=-\mathrm{i}\left[\Delta+\delta\hat{\omega}_{k}\right]\sigma_{k}+\mathrm{i}\sigma_{k}^{z}\sum_{i\neq
k}\left(g_{ki}^{(d)}+g_{ki}^{(\mu)}\right)\sigma_{i}+\sigma_{k}^{z}\sum_{i}\left(f_{ki}^{(d)}+f_{ki}^{(\mu)}\right)\sigma_{i}-\mathrm{i}\Omega\mathrm{e}^{\mathrm{i}\bm{k}_{0}\bm{R}_{i}}\sigma_{k}^{z}+\sigma_{k,\text{in}}\end{split}$
(56)

with a frequency operator

$\delta\hat{\omega}_{k}=\frac{1}{2}\sum_{i\neq
k}\left[\left(\Omega_{ik}^{e,e}-\Omega_{ik}^{g,g}\right)+\left(\Omega_{ik}^{e,e}+\Omega_{ik}^{g,g}-2\Omega_{ik}^{e,g}\right)\sigma_{i}^{z}\right].$
(57)

This is a configuration and density-dependent Zeeman shift of the $k$-th
atom’s frequency.

##  Appendix E Weak-excitation limit

The weak-excitation limit, i.e.
$\left\langle\sigma^{z}_{j}\sigma_{j^{\prime}}\right\rangle\approx-\left\langle\sigma_{j^{\prime}}\right\rangle$,
the equations of motion for the expectation values of the transition operators
become

$\begin{split}\frac{\mathop{}\\!\mathrm{d}\left\langle\sigma_{k}\right\rangle}{\mathop{}\\!\mathrm{d}t}&=-\mathrm{i}\left[\Delta+\delta\omega_{k}\right]\left\langle\sigma_{k}\right\rangle-\mathrm{i}\sum_{i\neq
k}\left(g_{ki}^{(d)}+g_{ki}^{(\mu)}\right)\left\langle\sigma_{i}\right\rangle-\sum_{i}\left(f_{ki}^{(d)}+f_{ki}^{(\mu)}\right)\left\langle\sigma_{i}\right\rangle+\mathrm{i}\Omega\mathrm{e}^{\mathrm{i}\bm{k}_{0}\bm{R}_{i}}\end{split}$
(58)

which is written in matrix notation in Eq. (16).