#  Superradiant transfer of quantized orbital angular momentum between light
and atoms in a ring trap

N. Piovella  Dipartimento di Fisica ”Aldo Pontremoli”, Università degli Studi
di Milano, Via Celoria 16, I-20133 Milano, Italy  G.R.M. Robb  SUPA and
Department of Physics, University of Strathclyde,Glasgow G4 0NG, Scotland, UK
R. Bachelard  Université Côte d’Azur, CNRS, Institut de Physique de Nice,
06560 Valbonne, France  Instituto de Física de São Carlos, Universidade de São
Paulo, 13560-970 São Carlos, SP, Brazil

###### Abstract

The orbital angular momentum (OAM) from a laser beam can be coherently
transferred to a Bose-Einstein condensate in a ring trap, in quantized units
of $\hbar$. The light-matter coupling allows for the superradiant transfer of
the atoms between the discrete OAM states. Tuning the ring parameters and
winding number of the pump light, specific angular momentum states can be
populated. This in turn allows control of the emission to generate light with
OAM different from that of the pump, as the atomic ring imprints its
contribution on the scattered light.

Introduction – Photons possess both translational and spin-angular momentum,
where the latter is associated with different polarization states of the light
wave [1]. Light may also be prepared in states with well-defined orbital
angular momentum (OAM), such as in Laguerre-Gaussian (LG) beams [2]. Coherent
control of the interaction between light and matter facilitates the
macroscopic transfer of those quantities from one medium to another, opening
up exciting and novel possibilities for quantum information engineering, such
as storing and transferring light in ensembles of atoms [3].

The spatial modes of the matter field possessing orbital angular momentum are
associated with a discrete infinite-dimensional (because of the quantization
of the OAM) Hilbert space. Transferring orbital angular momentum from light to
a Bose-Einstein condensate provides a route to entanglement that involves many
orthogonal quantum states, rather than just two [4, 5]. Multi-dimensional
entangled states could be of considerable importance in the field of quantum
information, enabling, for example, a more efficient use of communication
channels in quantum cryptography [6].

Recently, it has been shown that cold atoms trapped in a ring trap and driven
by a laser beam carrying OAM are set in rotation and acquire bunching in phase
[7]. Depending on the pump winding number $\ell$ and on the ring radius,
different harmonics of the phase become bunched, similarly to the longitudinal
configuration of light-matter instabilities [8]. However, a classical
treatment of the azimuthal motion of the atoms means that their OAM is not
quantized.

We here investigate the collective coupling between the ring of atoms and
light within a quantized approach, so the OAM of the trapped condensate
evolves through states with integer numbers of $\hbar$. The light-mediated
dipole-dipole interactions between the atoms results in a “superradiant”
coupling with the light, in analogy with the coupling achieved in cigar-shaped
matter waves and linear-momenta states [9], yet here promoting the transfer of
the atoms toward different OAM states. The nonlinear nature of this light-atom
dynamics leads to the population of higher OAM states for the atoms and,
consequently, for the scattered photons. This population of higher OAM states
stems from the presence of higher harmonics in the potential experienced by
the atoms. We show that this in turn enables the creation of cat states and
the emission of light with OAM different from the pump.

Modelling the light-atom interaction – Let us consider that the atoms are
driven resonantly by an OAM laser beam in the mode $\mathrm{LG}_{0\ell}$, with
electric field

$E_{0}(\rho,\phi,z)=A_{0}R(\rho)e^{i(k_{0}z+\ell\phi-\omega_{0}t)}+\mathrm{c.c.},$
(1)

with $R(\rho)=(\sqrt{2}\rho/w)^{\ell}\exp(-\rho^{2}/w^{2})$ using cylindrical
coordinates, $w$ the beam waist and $\ell$ the pump winding number. The atoms
rotate in a ring trap in the plane perpendicular to the pump propagation axis
$z$ (see Fig. 1a), and they are described by the matter wave function
$\Psi(\phi,\tau)$, where $\phi$ refers to the azimuthal angle. Their dynamics
is given by [7, 10]

$\displaystyle i\frac{\partial\Psi(\phi,\tau)}{\partial\tau}$ $\displaystyle=$
$\displaystyle-\frac{\partial^{2}\Psi(\phi,\tau)}{\partial\phi^{2}}$ (2)
$\displaystyle+$
$\displaystyle\frac{\gamma}{2}\int_{0}^{2\pi}d\phi^{\prime}V(\phi-\phi^{\prime})|\Psi(\phi^{\prime},\tau)|^{2}\Psi(\phi,\tau),$

with $\tau=\omega_{\phi}t$ the dimensionless time,
$\omega_{\phi}=\hbar/2M\rho^{2}$ the angular recoil frequency, $M$ the atom
mass and $\rho$ the ring radius. The light-atom coupling is characterized by
the parameter
$\gamma=N(\Gamma/\omega_{\phi})(\Omega_{0}/2\Delta_{0})^{2}R^{2}(\rho)$, where
$N$ is the atom number, $\Omega_{0}=dA_{0}/\hbar$ the pump Rabi frequency,
$\Gamma$ the atomic decay rate for the two-level transition at frequency
$\omega_{a}$, $d$ the associated dipole moment, and
$\Delta_{0}=\omega_{a}-\omega_{0}$ the pump-atom detuning. The potential $V$
reads

$V(\varphi)=-\frac{\cos[2k_{0}\rho q(\varphi)-\ell\varphi]}{k_{0}\rho
q(\varphi)},$ (3)

where $q(\varphi)=\sqrt{\sin^{2}(\varphi/2)+\epsilon^{2}}$ and $\epsilon$ is a
cut-off parameter. The potential $V$ has a finite range due to the multimode
emission in vacuum [11], and it decreases as $1/r_{jm}$ with the distance
$r_{jm}$ between the pair of atoms $(j,m)$. Since the atoms are set on a ring
of radius $\rho$, we have $r_{jm}=2\rho|\sin[(\phi_{j}-\phi_{m})/2]|$. The
phenomenological cut-off $\epsilon$ in $q(\phi_{j}-\phi_{m})$ has been added
to avoid the singularity in the denominator and can be chosen arbitrarily
small.

The wave function $\Psi$ describes the atoms in a BEC state, with the
normalization condition $\int_{0}^{2\pi}|\Psi(\phi,\tau)|^{2}d\phi=1$. This
wave-function can be expanded on the basis of the angular momentum states
$|m\rangle$ as
$\Psi(\phi,\tau)=(1/\sqrt{2\pi})\sum_{m}c_{m}(\tau)\exp(im\phi)$, while the
potential is expanded as a Fourier series:
$V(\varphi)=\sum_{k}V_{k}\exp(ik\varphi)$, with
$V_{k}=(1/2\pi)\int_{0}^{2\pi}V(\varphi)\exp(-ik\varphi)d\varphi$. By
projecting Eq. (2) on the state $|m\rangle$, we obtain the dynamical equations
for the probability amplitudes $c_{m}(\tau)$:

$\displaystyle\frac{dc_{m}}{d\tau}$ $\displaystyle=$ $\displaystyle-
im^{2}c_{m}-i\frac{\gamma}{2}\sum_{k}V_{k}c_{m-k}\sum_{n}c^{*}_{n-k}c_{n},$
(4)

where the normalization condition now reads $\sum_{m}|c_{m}(\tau)|^{2}=1$. The
radiated intensity emitted by the atoms is
$I(\theta,\phi)=I_{1}R^{2}(\rho)N^{2}|M(\theta,\phi)|^{2}$, where [10]

$\displaystyle M(\theta,\phi)$ $\displaystyle=$
$\displaystyle\int_{0}^{2\pi}e^{-ik_{0}\rho\sin\theta\cos(\phi-\phi^{\prime})+i\ell\phi^{\prime}}|\Psi(\phi^{\prime})|^{2}\frac{d\phi^{\prime}}{2\pi}$
(5)

is the dimensionless electric field, $\theta$ is the polar angle,
$I_{1}=(\hbar\omega_{0}\Gamma/8\pi r^{2})(\Omega_{0}/2\Delta_{0})^{2}$ and $r$
is the distance from the detector. Introducing the $m$-th azimuthal bunching
coefficient

$\Phi_{m}(\tau)=\sum_{n=-\infty}^{+\infty}c_{n-m}^{*}(\tau)c_{n}(\tau),$ (6)

the electric field rewrites as the sum of different OAM components, composed
by photons with angular momentum $\hbar\ell^{\prime}=\hbar(\ell+m)$:

$M(\theta,\phi)=\sum_{m}(-i)^{\ell+m}J_{\ell+m}(k_{0}\rho\sin\theta)\Phi_{m}(\tau)e^{i(\ell+m)\phi},$
(7)

where $J_{n}(x)$ is the Bessel function of the first kind and of order $n$.
The rescaled intensity averaged over the azimuthal angle is
$\bar{I}(\theta)=\sum_{m}J_{\ell+m}^{2}(k_{0}\rho\sin\theta)|\Phi_{m}|^{2}$.

Linear stability analysis – Let us now discuss the atom dynamics. We observe
that the cubic nonlinearity in Eq. (4) couples state $m$ with the states $m\pm
k$ provided that a $k$-th Fourier component $V_{k}$ of the potential and a
$k$-th azimuthal bunching $\Phi_{k}$ are present. More quantitatively, the
dynamical growth of the azimuthal bunching at harmonic $m$ is obtained by a
linear stability analysis of Eq. (4) around the equilibrium point
$\Phi_{m}^{(0)}=\delta_{m,0}$, which describes uniformly distributed phases.
This leads to the complex eigenvalues [10]

$\lambda_{m}=\pm im\sqrt{m^{2}+\gamma V_{m}},$ (8)

where the $m$-th azimuthal bunching $\Phi_{m}$ grows exponentially as
$\exp[\mathrm{Re}(\lambda_{m})\tau]$. Hence, the exponential growth of the
azimuthal bunching $\Phi_{m}$ requires a finite imaginary part of the $m$-th
Fourier component of the potential, $V_{m}$. For $\ell=0$ the potential
$V(\varphi)$ is an even function of $\varphi$ and $V_{m}$ is real, hence the
system is unstable only in presence of OAM light, with $\ell\neq 0$.
Intuitively, the OAM pump exerts on the atoms a torque forcing them to rotate.
Since the atomic angular momentum is quantized, the atoms rotate with a
discrete angular momentum, whose value is determined by the anharmonicity of
the potential $V$.

Figure 1: (a) Scheme of the atomic ring pumped with OAM light, so the atoms
acquire phase bunching. (b) Growth rate $|\mathrm{Re}(\lambda_{m})|$ for the
modes $m$, as a function of the ring radius $k_{0}\rho$. The mode with the
strongest rate corresponds to $m\approx k\rho$. Simulations realized with
$\gamma=0.2$ and $\epsilon=0.1$.

When uniformly distributed in phase, all the atoms rotate with the same
angular velocity and no macroscopic rotation is visible. However, in presence
of an OAM beam, the combined effect of the pump and scattered beam creates a
self-consistent potential which modulates the azimuthal atomic phases. As a
result, atoms are coupled by the scattered photons, recoiling in phase as
angular momentum is transferred from the light. Because the beam-matter
coupling grows with the bunching in phase, an instability develops and both
the bunching in phase and the scattered intensity grow exponentially until the
moment when the atoms are transferred to a new OAM state $|m\rangle$. In the
case of a subwavelength trap, the ring of atoms essentially behaves as a
point-like dipole and the matter is transferred to the state $|\ell\rangle$,
where $\hbar\ell$ is the angular momentum of the pump photons. However,
increasing the ring radius $\rho$, the potential presents higher harmonics and
is able to transfer the atoms to OAM states with $m$ larger than $\ell$.

Figure 2: (a) Populations $N_{m}$ and (b) average angular velocity
$\langle\omega\rangle$, in units of $\omega_{\phi}$, as a function of the
dimensionless time $\tau$. (c) Dimensionless scattered intensity
$|M(\theta,\phi)|^{2}$ and (d) $\bar{I}(\theta)$ vs $\theta$ for
$k_{0}\rho=1$, $\ell=1$ and $\gamma=0.05$; The blue line is the component with
$\ell^{\prime}=0$, the red line is the component with $\ell^{\prime}=1$ and
the dashed black line is the total intensity. for $\gamma=0.05$, $k_{0}\rho=1$
and $\ell=1$. The inset in (b) presents the rate coefficients $g_{k}$ (in
a.u.) as a function of $k$, showing that the dominant transition is for $k=1$.

According to Eq. (8), the system admits two different regimes: a classical
regime when $\gamma|V_{m}|\gg m^{2}$ and thus $\lambda_{m}\approx\pm
im\sqrt{\gamma V_{m}}$, which is characterized by an exponential growth rate
proportional to $\sqrt{N}$, and a quantum regime, when $\gamma|V_{m}|\ll
m^{2}$ so that $\lambda_{m}\approx\pm i(m^{2}+\gamma V_{m}/2)$, which presents
a rate proportional to $N$ [12]. In the classical regime, both the states
$|m\rangle$ and $|-m\rangle$ become populated almost simultaneously.
Differently, in the quantum regime, the dynamical term $\exp(\pm im^{2}\tau)$
corresponds to a phase shift proportional to the rotational energy
$L_{z}^{2}/2M\rho^{2}=\hbar\omega_{\phi}m^{2}$, with orbital angular momentum
$L_{z}=\hbar m$. This energy shift, which stems from energy and angular
momentum conservation in the photon-atom scattering process, hinders the
occupation of one of the two states $|\pm m\rangle$ (which one is populated
depends on the sign of $\mathrm{Im}(V_{m})$). This leads to the scaling of the
growth rate as $N$ for the quantum regime [12].

The scattering of the pump light into different OAM states can be deduced from
the eigenspectrum, whose real part $|\mathrm{Re}(\lambda_{m})|$ is presented
in Fig. 1b for the quantum regime ($\gamma=0.2$), as a function of the radius
$k_{0}\rho$ and of $m$. The fact that for increasing $k_{0}\rho$ the maximum
occurs at values of $m$ which scales linearly with $k_{0}\rho$ implies that
larger ring are able to populate higher OAM states of the condensate. This is
a notable difference with the case of the one-dimensional quantum collective
atomic recoil lasing (CARL) [8, 13, 14], where only transitions between
adjacent linear momentum states are possible. This is due to the fact that
CARL is characterized by a sinusoidal potential at the light wavelength, which
results from the interference between the pump and the backward scattered
modes. In the ring trap, the presence of higher harmonics in potential (3)
allows the atoms to bunch at higher OAM states but also, as we shall see
later, to scatter photons of different OAM in different directions.

Rate equations – In the quantum regime, when $\gamma|V_{m}|\ll m^{2}$, the
atoms are transferred to OAM states $|m\rangle$ through a superradiant cascade
which involves, essentially, only two states at the time. Then, Eqs. (4) are
approximated by the following rate equations for the state populations
$N_{m}=|c_{m}|^{2}$:

$\frac{dN_{m}}{d\tau}=\left[\sum_{k=1}^{m}g_{k}N_{m-k}-\sum_{k=1}^{+\infty}g_{k}N_{m+k}\right]N_{m},$
(9)

where $g_{k}=\gamma|\mathrm{Im}(V_{k})|$ are the rate coefficients, and with
the normalization condition $\sum_{m}N_{m}=1$. These coupled equations
generalize the superradiant transition between two states $N_{0}$ and $N_{k}$
when only a single coefficient $g_{k}$ is present, and it solves as
$N_{0,k}(\tau)=(1/2)\\{1\mp\tanh[g_{k}(\tau-\tau_{0})/2]\\}$, with
$\tau_{0}=(1/g_{k})\ln[2/\sqrt{N_{k}(0)}]$ the delay time. In the superradiant
cascade the coherence of the matter wave function is preserved and the phases
of the complex amplitudes $c_{m}=\sqrt{N_{m}}e^{i\phi_{m}}$ are determined by
the equations

$\frac{d\phi_{m}}{d\tau}=-(m^{2}+\gamma
V_{0})-\left[\sum_{k=1}^{m}\alpha_{k}N_{m-k}+\sum_{k=1}^{+\infty}\alpha_{k}N_{m+k}\right],$
(10)

with $\alpha_{k}=(\gamma/2)\mathrm{Re}(V_{k})$.

For a small ring radius, $k_{0}\rho\ll 1$, the atoms are transferred
sequentially through the OAM states by steps of $\Delta m=\ell$. This is
illustrated in Fig. 2, which depicts the temporal evolution of the populations
$N_{m}$ for the first OAM states (see panel (a)) and the average angular
velocity $\langle\omega\rangle=\sum_{m}mN_{m}$ (see panel (b)). The angular
distribution of the scattered field intensity $|M(\theta,\phi)|^{2}$ is
plotted in Fig. 2(c) at the time when the bunching is maximum
($|\Phi_{1}|\approx 1/2$), while panel (d) shows the total average intensity
$\bar{I}(\theta)$ as a function of the polar angle $\theta$ (dashed line),
together with the main OAM components. In this case the field (5) is
approximated by $M(\theta,\phi)\approx-
iJ_{1}(k_{0}\rho\sin\theta)\exp(i\phi)+\Phi_{1}^{*}J_{0}(k_{0}\rho\sin\theta)$,
with $|\Phi_{1}|\approx 1/2$ and where we have neglected the term proportional
to $J_{2}(\sin\theta)$. The average intensity is $\bar{I}(\theta)\approx
J_{1}^{2}(k_{0}\rho\sin\theta)+|\Phi_{1}|^{2}J_{0}^{2}(k_{0}\rho\sin\theta)$.
In Fig. 2 we can observe the component $\ell^{\prime}=1$ (red line), the
component $\ell^{\prime}=0$ (blue line) and the total intensity (dashed black
line). The light with the same orbital angular momentum of the pump $\ell=1$
is emitted in the perpendicular direction, whereas the light with $\ell=0$ is
emitted forward and backward within a large diffraction angle. The total
intensity is almost isotropic, as expected for a small ring.

Populating higher OAM states – Atomic rings larger than the wavelength are
able to populate directly atomic states with large OAM, taking advantage of
the higher harmonics of the potential $V$. In Fig. 3 we illustrate how by
choosing parameters such that the potential favors the coupling to several
higher OAM states, the system can be cast in a macroscopic superposition. The
pump with $\ell=2$ transfers most of the atomic population from $m=0$ to
$m=6$, and to $m=5$ to a lesser extent, due to the set of coefficients $g_{k}$
among which $g_{6}$ and $g_{5}$ are highest (the parameters have been chosen
so $g_{1}=0$, so transitions with $\Delta m=1$ are suppressed). The first
transition sees a transfer of population from $m=0$ to $m=6$ ($N_{6}=91\%$)
and to $m=5$ ($N_{5}=9\%$). During the next transitions, the atoms are
transferred to $m=12$ and $m=11$, then to $m=18$ and $m=17$, and so on, since
the specific set of parameters favors transitions with $\Delta m=6$. Thus,
after two transition channels are open (that is, with $\Delta m=5$ and $6$), a
macroscopic superposition of two rotational states, is created. Throughout the
successive transitions toward higher phase bunching states, superpositions of
two states persist. The possibility that these superposition states are
Schrodinger cat states will be investigated elsewhere.

Figure 3: Populations $N_{m}$ as a function of the dimensionless time $\tau$.
Simulations realized for $k_{0}\rho=5.605$, $\ell=2$ and $\gamma=1$. The inset
presents the rate coefficients $g_{k}$ (in a.u.) as a function of $k$, showing
that the dominant transitions are for $k=5$ and $k=6$.

Tailoring the light emission – The dynamical evolution to higher OAM atomic
states in turn results in the emission of light whose angular momentum differs
from the pump one. This can be a tool either for reading the atomic state, or
for producing light with a specific OAM. In particular, one may choose the
ring radius such that no photons with the same OAM of the pump are emitted in
a particular direction, so light with a different OAM can be collected more
efficiently. As an example, we show in Fig. 4(a) how an atomic ring with
radius $k_{0}\rho=5$, along with a pump of winding number $\ell=2$, is
transferred to the OAM state $m=5$.

Figure 4: (a) Population of the state $m=0$ and $m=5$ as a function of the
dimensionless time $\tau$. The inset shows the coefficients $g_{k}$. (b)
Angular distribution of the radiation intensity $|M(\theta,\phi)|^{2}$, at
time $\tau=800$. (c) Dimensionless intensity $\bar{I}(\theta)$ averaged over
the azimuthal angle $\phi$, as a function of the polar angle $\theta$, at time
$\tau=800$. The dashed black line stands for the total intensity, the blue
line for the component $\ell^{\prime}=2$ and the red line for the component
$\ell^{\prime}=-3$: the latter is the only component present at
$\theta=\pi/2$. Simulations realized for $k_{0}\rho=5$, $\ell=2$ and
$\gamma=0.2$.

This bunching in phase leads to five-lobe radiation pattern in the azimuthal
direction, see panel (b). Because the Bessel function
$J_{2}(k_{0}\rho\sin\theta)$ vanishes for $k_{0}\rho=5$ and $\theta=\pi/2$, no
photons with $\ell^{\prime}=\ell=2$ are emitted in the transverse direction.
At this angle, only the component with $\ell^{\prime}=l+m=-3$ of the scattered
light is present. Thus, in this specific direction, light with a well-defined
OAM, which is different from that of the pump, can be collected.

In conclusion, we have demonstrated how orbital angular momentum can be
transferred coherently from a single Laguerre-Gaussian to a Bose-Einstein
condensate in a ring trap. The spontaneous formation of the azimuthal bunching
of ultracold atoms in the ring trap is at the core of the process. The process
is similar to the quantum regime of collective atomic recoil lasing [8, 13,
14] and to the superradiant Rayleigh scattering [9], but with transfer of
orbital angular momentum in units of $\hbar$ instead of linear momentum in
units of $\hbar k$. Differently from these linear configurations, the many
harmonics of the potential here allow to promote the atoms to one or more
higher OAM states. This mechanism is particularly promising for the creation
of states of matter with tailored OAM [15, 16, 17]: in the present scheme, the
choice of the ring radius, by tuning the trap frequencies, allows one to
manipulate the relative populations of different angular momentum states. The
emission of photons with angular momentum different from that of the pump and
with a specific spatial pattern offers the possibility to either use this
light to read the atomic state, or to use the atomic ring to produce light
with a new OAM. The atomic ring may also be used as a beam splitter, dividing
the light into two or more topological states, so robust phase qubits can be
generated. In this context, one may take advantage of the complete orthonormal
basis formed by the OAM states to design more elaborate quantum states for
quantum information processing [18].

The proposed setup correspond to ring traps with radii of a few optical
wavelengths, or subwavelength for the case $\rho<\lambda$ (yet this latter
case does not produce efficiently higher OAM light). Recent advances in the
control of quantum gases have seen the development of atom ring traps formed
from both magnetic and optical dipole potential, and more recent
implementations using hybrids of both, down to $\rho=9\mu m$ ring traps [19,
20]. Although much more challenging, nanometer-scale traps are also currently
under investigation [21].

## Acknowledgments

We thank André Cidrim for his help with editing figures. R.B. has received
support from the São Paulo Research Foundation (FAPESP) through Grants
No.2019/12842-2, 2019/13143-0 and 2018/15554-5, as well as from the National
Council for Scientific and Technological Development (CNPq) Grant Nos.
409946/2018-4 and 313886/2020-2. Part of this work was performed in the
framework of the European Training Network ColOpt, which is funded by the
European Union (EU) Horizon 2020 program under the Marie Sklodowska-Curie
action, grant agreement No. 721465.

## References

  * Beth [1935] R. A. Beth, Direct detection of the angular momentum of light, Phys. Rev. 48, 471 (1935).
  * Allen  _et al._ [1992] L. Allen, M. W. Beijersbergen, R. J. C. Spreeuw, and J. P. Woerdman, Orbital angular momentum of light and the transformation of laguerre-gaussian laser modes, Phys. Rev. A 45, 8185 (1992).
  * Krenn  _et al._ [2017] M. Krenn, M. Malik, M. Erhard, and A. Zeilinger, Orbital angular momentum of photons and the entanglement of laguerre–gaussian modes, Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 375, 20150442 (2017).
  * Andersen  _et al._ [2006] M. Andersen, C. Ryu, P. Cladé, V. Natarajan, A. Vaziri, K. Helmerson, and W. D. Phillips, Quantized rotation of atoms from photons with orbital angular momentum, Physical review letters 97, 170406 (2006).
  * Wright  _et al._ [2008] K. Wright, L. Leslie, and N. Bigelow, Optical control of the internal and external angular momentum of a bose-einstein condensate, Physical Review A 77, 041601 (2008).
  * Fickler  _et al._ [2016] R. Fickler, G. Campbell, B. Buchler, P. K. Lam, and A. Zeilinger, Quantum entanglement of angular momentum states with quantum numbers up to 10,010, Proceedings of the National Academy of Sciences 113, 13642 (2016).
  * Gisbert  _et al._ [2022] A. Gisbert, N. Piovella, G. Robb, and D. McLellan, Superradiant atomic recoil lasing with orbital-angular-momentum light, Physical Review A 105, 023526 (2022).
  * Bonifacio and Souza [1994] R. Bonifacio and L. D. S. Souza, Collective atomic recoil laser (carl) optical gain without inversion by collective atomic recoil and self-bunching of two-level atoms, Nuclear Instruments and Methods in Physics Research 341, 360 (1994).
  * 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 (1999).
  * [10] See supplemental material at [url will be inserted by publisher] for the derivation of eqs.(2), (5) and (8)., .
  * Ayllon  _et al._ [2019] R. Ayllon, J. T. Mendonça, A. T. Gisbert, N. Piovella, and G. R. M. Robb, Multimode collective scattering of light in free space by a cold atomic gas, Physical Review A 100, 023630 (2019).
  * Robb  _et al._ [2005] G. Robb, N. Piovella, and R. Bonifacio, The semiclassical and quantum regimes of super-radiant light scattering from a bose–einstein condensate, Journal of Optics B: Quantum and Semiclassical Optics 7, 93 (2005).
  * Piovella  _et al._ [2001] N. Piovella, M. Gatelli, and R. Bonifacio, Quantum effects in the collective light scattering by coherent atomic recoil in a bose–einstein condensate, Optics Communications 194, 167 (2001).
  * Slama  _et al._ [2007] S. Slama, S. Bux, G. Krenz, C. Zimmermann, and P. W. Courteille, Superradiant rayleigh scattering and collective atomic recoil lasing in a ring cavity, Physical Review Letters 98, 10.1103/physrevlett.98.053603 (2007).
  * Robb [2012] G. Robb, Superradiant exchange of orbital angular momentum between light and cold atoms, Physical Review A 85, 023426 (2012).
  * Das  _et al._ [2016] P. Das, M. E. Tasgin, and Özgür E Müstecaplıoğlu, Collectively induced many-vortices topology via rotatory dicke quantum phase transition, New Journal of Physics 18, 093022 (2016).
  * Kumar  _et al._ [2021] P. Kumar, T. Biswas, K. Feliz, R. Kanamoto, M.-S. Chang, A. K. Jha, and M. Bhattacharya, Cavity optomechanical sensing and manipulation of an atomic persistent current, Physical Review Letters 127, 113601 (2021).
  * Kapale and Dowling [2005] K. T. Kapale and J. P. Dowling, Vortex phase qubit: Generating arbitrary, counterrotating, coherent superpositions in bose-einstein condensates via optical angular momentum beams, Physical review letters 95, 173601 (2005).
  * Amico  _et al._ [2021] L. Amico, M. Boshier, G. Birkl, A. Minguzzi, C. Miniatura, L.-C. Kwek, D. Aghamalyan, V. Ahufinger, D. Anderson, N. Andrei, _et al._ , Roadmap on atomtronics: State of the art and perspective, AVS Quantum Science 3, 039201 (2021).
  * Eckel  _et al._ [2014] S. Eckel, J. G. Lee, F. Jendrzejewski, N. Murray, C. W. Clark, C. J. Lobb, W. D. Phillips, M. Edwards, and G. K. Campbell, Hysteresis in a quantized superfluid ‘atomtronic’circuit, Nature 506, 200 (2014).
  * Richberg  _et al._ [2021] R. Richberg, S. Szigeti, and A. Martin, Optical focusing of bose-einstein condensates, Physical Review A 103, 063304 (2021).


Standard Gaussian processes cannot be solved analytically for data with non-
Gaussian likelihoods, so we will use the approximate, variational inference
Gaussian process method from Hensman et al. (2014) and implemented in
GPyTorch333https://github.com/cornellius-gp/gpytorch (Gardner et al., 2021).

The Gaussian process assumption models the latent function in Eq. (26) as
(following the notation of Rasmussen & Williams, 2006)

$\displaystyle f(\mathbf{x})$
$\displaystyle\sim\mathcal{GP}(\mu,k(\mathbf{x},\mathbf{x}^{\prime}|\theta))\,,$
(28)

which is shorthand for

$\displaystyle\mathds{E}[f(\mathbf{x})]=\mu\,,$ (29)
$\displaystyle\mathrm{Var}[f(\mathbf{x}),f(\mathbf{x}^{\prime})]=k(\mathbf{x},\mathbf{x}^{\prime}|\theta)\,,$
(30)

and means that the values of $f$ are fully determined by the mean, $\mu$, and
the covariance function $k(\mathbf{x},\mathbf{x}^{\prime}|\theta)$ between
different inputs $\mathbf{x}$ and $\mathbf{x}^{\prime}$. We will be using the
radial basis function (or squared exponential) kernel for $k$:

$k(\mathbf{x},\mathbf{x}^{\prime}|\theta)=\sigma^{2}\prod_{i=0}^{d}\exp\left(-\frac{((\mathbf{x})_{i}-(\mathbf{x}^{\prime})_{i})^{2}}{2\ell_{i}^{2}}\right)\,,$
(31)

where $i$ runs over the $d=9$ dimensions of $\mathbf{x}$ and each dimension
has its own covariance lengthscale $\ell_{i}$, resulting in hyperparameter
$\theta=(\mu,\sigma^{2},\bm{\ell})$.

The power of Gaussian process regression stems from the conditioning property
of Gaussian distributions. In what follows, we assume $\mu=0$ for simplicity.
Given the assumed joint Gaussian distribution between function values at
$\mathsf{X}$ and $\mathsf{X}^{*}$, $p(\mathbf{f},\mathbf{f}^{*})$, which we
can write as

$p(\mathbf{f},\mathbf{f}^{*})=p\left(\begin{bmatrix}\mathbf{f}\\\
\mathbf{f}^{*}\end{bmatrix}\right)=\mathcal{N}\left(\begin{bmatrix}\mathbf{0}\\\
\mathbf{0}\end{bmatrix},\begin{bmatrix}\mathsf{K}_{\mathsf{X}\mathsf{X}}&\mathsf{K}_{\mathsf{X}\mathsf{X}^{*}}\\\
\mathsf{K}_{\mathsf{X}^{*}\mathsf{X}}&\mathsf{K}_{\mathsf{X}^{*}\mathsf{X}^{*}}\\\
\end{bmatrix}\right)\,,$ (32)

the conditional distribution $p(\mathbf{f}^{*}|\mathbf{f})$ is a new Gaussian
distribution given by

$p(\mathbf{f}^{*}|\mathbf{f},\theta)=\mathcal{N}\left(\mathsf{K}_{\mathsf{X}^{*}\mathsf{X}}\mathsf{K}_{\mathsf{X}\mathsf{X}}^{-1}\mathbf{f},\mathsf{K}_{\mathsf{X}^{*}\mathsf{X}^{*}}-\mathsf{K}_{\mathsf{X}^{*}\mathsf{X}}\mathsf{K}_{\mathsf{X}\mathsf{X}}^{-1}\mathsf{K}_{\mathsf{X}\mathsf{X}^{*}}\right)\,.$
(33)

Here we have introduced the $n\times n$ covariance matrix
$\mathsf{K}_{\mathsf{X}\mathsf{X}}$, with
$(\mathsf{K}_{\mathsf{X}\mathsf{X}})_{ij}=k(\mathbf{x}_{i},\mathbf{x}_{j})$
and $k$ given by Eq. (31), containing the covariance between different input
points in $\mathsf{X}$. Importantly, the probability distribution of
$\mathbf{f}^{*}$ for an arbitrary input location $\mathsf{X}^{*}$ depends
solely on the finite set of _measured_ inputs $\mathsf{X}$. Clearly, the
accuracy of the prediction $f(\mathbf{x}^{*})$ depends on the distance to the
nearest measured input $\mathbf{x}$ in $\mathsf{X}$ and the lengthscale
hyperparameter $\bm{\ell}$, with the function values $f^{*}$ regressing to the
mean $0$ and prior uncertainty $k(\mathbf{x}^{*},\mathbf{x}^{*})$ for
$\mathsf{K}_{\mathbf{x}^{*}\mathsf{X}}\to\mathbf{0}^{T}$. We can use
$p(\mathbf{f}^{*}|\mathbf{f},\theta)$ to predict $\mathbf{N}(\mathsf{X}^{*})$.

The optimal hyperparameters, $\theta$, for the simulated data, $\mathbf{N}$,
are found by maximizing

$p(\theta|\mathbf{N})=\frac{p(\theta)p(\mathbf{N}|\theta)}{p(\mathbf{N})}\,.$
(34)

where we introduced the marginal likelihood

$p(\mathbf{N}|\theta)=\int
p(\mathbf{N}|\mathbf{f})p(\mathbf{f}|\theta)\mathrm{d}\mathbf{f}\,,$ (35)

which cannot be solved analytically in the case of a negative binomial
likelihood.

This standard Gaussian process encounters two major difficulties. First, the
$\mathsf{K}_{\mathsf{X}\mathsf{X}}^{-1}$-term in Eq. (33) becomes
computationally expensive for datasets with large $n$. Second, non-Gaussian
likelihoods require approximations to optimize Eq. (34), since no closed-form
analytical solution exists. Both of these problems have been solved by the
sparse method using inducing variables and the variational free energy as
introduced by Titsias (2009) and applied to non-Gaussian likelihoods by
Hensman et al. (2014) and formalized by Matthews et al. (2016). We will
briefly introduce the necessary ingredients for this method.

The idea behind the method of Titsias (2009) is to introduce both an extra set
of $m\ll n$ inducing (or pseudo) inputs $\mathsf{Z}$ of the Gaussian process
such that $\mathbf{f}(\mathsf{Z})\equiv\mathbf{u}$ and an approximate
distribution over these function values, $q_{\psi}(\mathbf{u})$. The inducing
point locations $\mathsf{Z}$ and the parameters $\psi$ of the approximate
distribution family will be chosen in such a way that they optimally capture
the true posterior probability of the Gaussian process, i.e.
$q_{\psi}(\mathbf{f})\simeq p(\mathbf{f}|\mathbf{N})$. Assuming a Gaussian
distribution for $q(\mathbf{u})$ with

$q_{\psi}(\mathbf{u})=\mathcal{N}(\mathbf{m},\mathsf{S})\,,$ (36)

we get $\psi=(\mathbf{m},\mathsf{S})$ and we calculate the full approximate
distribution as

$q_{\psi}(\mathbf{f},\mathbf{u})=p(\mathbf{f}|\mathbf{u})q_{\psi}(\mathbf{u})\,.$
(37)

The optimization of $(\mathsf{Z},\mathbf{m},\mathsf{S})$ now needs to ensure
that

$\displaystyle p(\mathbf{f}|\mathbf{N})$ $\displaystyle\simeq\int
p(\mathbf{f}|\mathbf{u})q(\mathbf{u})\mathrm{d}\mathbf{u}$
$\displaystyle\Updownarrow
p(\mathbf{f}|\mathbf{u})=\mathcal{N}(\mathsf{K}_{\mathsf{X}\mathsf{Z}}\mathsf{K}_{\mathsf{Z}\mathsf{Z}}^{-1}\mathbf{u},\mathsf{D}_{\mathsf{X}\mathsf{X}})$
$\displaystyle=\mathcal{N}(\mathsf{K}_{\mathsf{X}\mathsf{Z}}\mathsf{K}_{\mathsf{Z}\mathsf{Z}}^{-1}\mathbf{m},\mathsf{D}_{\mathsf{X}\mathsf{X}}+\mathsf{K}_{\mathsf{X}\mathsf{Z}}\mathsf{K}_{\mathsf{Z}\mathsf{Z}}^{-1}\mathsf{S}\mathsf{K}_{\mathsf{Z}\mathsf{Z}}^{-1}\mathsf{K}_{\mathsf{Z}\mathsf{X}})\,,$
(38)

with
$\mathsf{D}_{\mathsf{X}\mathsf{X}}=\mathsf{K}_{\mathsf{X}\mathsf{X}}-\mathsf{K}_{\mathsf{X}\mathsf{Z}}\mathsf{K}_{\mathsf{Z}\mathsf{Z}}^{-1}\mathsf{K}_{\mathsf{Z}\mathsf{X}}$,
due to the conditioning property of Eq. (33) (see Chapter 4.3 of Matthews,
2017 for detailed explanations). Evaluating this expression only requires the
inverted $m\times m$ matrix $K_{\mathsf{Z}\mathsf{Z}}^{-1}$, significantly
reducing the computational cost of making model predictions.

To determine the optimal values $(\mathsf{Z},\mathbf{m},\mathsf{S})$, we
minimize the difference between the approximate distribution
$q_{\psi}(\mathbf{f},\mathbf{u})$ and the model posterior
$p(\mathbf{f},\mathbf{u}|\mathbf{N})$ through the Kullback-Leibler (KL)
divergence

$\displaystyle\mathcal{KL}[q(\mathbf{f},\mathbf{u})||p(\mathbf{f},\mathbf{u}|\mathbf{N})]=-\int
q(\mathbf{f},\mathbf{u})\log\left(\frac{p(\mathbf{f},\mathbf{u}|\mathbf{N})}{q(\mathbf{f},\mathbf{u})}\right)\mathrm{d}\mathbf{f}\mathrm{d}\mathbf{u}\,.$
(39)

Defining this equation as $\mathcal{K}$, we use Bayes’ theorem to rewrite
$p(\mathbf{f},\mathbf{u}|\mathbf{N})=p(\mathbf{N}|\mathbf{f})p(\mathbf{f}|\mathbf{u})p(\mathbf{u})/p(\mathbf{N|\theta})$,
making use of the fact that the observations are only conditionally dependent
on their corresponding function values $\mathbf{f}$. Also filling in Eq. (37),
we then find

$\displaystyle\mathcal{K}=$ $\displaystyle-\int
p(\mathbf{f}|\mathbf{u})q(\mathbf{u})\log\left(\frac{p(\mathbf{N}|\mathbf{f})p(\mathbf{u})}{p(\mathbf{N|\theta})q(\mathbf{u})}\right)\mathrm{d}\mathbf{f}\mathrm{d}\mathbf{u}$
$\displaystyle=$ $\displaystyle-\int q(\mathbf{f})\log
p(\mathbf{N}|\mathbf{f})\mathrm{d}\mathbf{f}+\int
p(\mathbf{f}|\mathbf{u})q(\mathbf{u})\log
p(\mathbf{N}|\theta)\mathrm{d}\mathbf{f}\mathrm{d}\mathbf{u}$
$\displaystyle+\int
p(\mathbf{f}|\mathbf{u})q(\mathbf{u})\log\left(\frac{q(\mathbf{u})}{p(\mathbf{u})}\right)\mathrm{d}\mathbf{f}\mathrm{d}\mathbf{u}$
$\displaystyle\mathcal{K}=$ $\displaystyle\log
p(\mathbf{N}|\theta)-\mathds{E}_{q(\mathbf{f})}\left[\log
p(\mathbf{N}|\mathbf{f})\right]+\mathcal{KL}[q(\mathbf{u})||p(\mathbf{u})]\,.$
(40)

We can rearrange terms in this expression and use the fact that the KL
divergence is strictly positive to arrive at the variational evidence lower
bound (ELBO), which provides a lower bound on the marginal likelihood—also
called the model evidence—as the name suggests

$\log
p(\mathbf{N}|\theta)\geq\mathcal{L}_{\mathrm{ELBO}}=\mathds{E}_{q(\mathbf{f})}\left[\log
p(\mathbf{N}|\mathbf{f})\right]-\mathcal{KL}[q(\mathbf{u})||p(\mathbf{u})]\,.$
(41)

Equality for this equation holds exactly when Eq. (39) equals zero, which is
the case when $q(\mathbf{f},\mathbf{u})=p(\mathbf{f},\mathbf{u}|\mathbf{N})$.
Assuming no covariance between $f_{i}$ and $N_{j\neq i}$, the likelihood
factors and we have

$\mathcal{L}_{\mathrm{ELBO}}=\sum_{i=1}^{n}\mathds{E}_{q(f_{i})}\left[\log
p(N_{i}|f_{i})\right]-\mathcal{KL}[q(\mathbf{u})||p(\mathbf{u})]\,,$ (42)

where the first term consists of a sum of one dimensional integrals which can
be computed easily using Gauss-Hermite quadrature, and the second term is the
KL divergence between two multivariate Gaussian distributions, since
$p(\mathbf{u})=\mathcal{N}(\mathbf{0},\mathsf{K}_{\mathsf{Z}\mathsf{Z}})$ due
to the Gaussian process assumption. Optimizing the ELBO is equivalent to
maximizing the marginal log-likelihood in Eq. (35).

We use the
ApproximateGP444https://docs.gpytorch.ai/en/latest/examples/04_Variational_and_Approximate_GPs/Non_Gaussian_Likelihoods.html
implementation of GPyTorch to model and optimize $f(\mathbf{x}_{i})$ with a
custom implementation of the negative binomial likelihood between
$N(\mathbf{x}_{i})$ and the measurements number counts from the simulations,
$N_{i}$.

##  Appendix C Emulator performance

Figure 15: _Top row:_ The residuals between the best-fit emulator and the
individual simulations in the hypercube for different apertures at fixed
redshift $z=0.43$ (_first three columns_) and at different redshifts for a
fixed aperture $(R_{1},R_{2},R_{\mathrm{m}})=(1,2,3)\,\mathrm{cMpc}$
(_rightmost column_). All halo samples were selected with
$m_{\mathrm{200c}}>10^{13.5}\,\mathrm{M_{\odot}}$. For the aperture
variations, we show the individual simulations and for the redshift variations
the 16th to 84th percentile scatter to avoid clutter. The thick, coloured line
indicates the median emulator residual and the shaded regions in the aperture
plots indicate the median of the bootstrapped simulation uncertainties. The
black, dash-dotted lines indicate the $\pm 2\,\mathrm{per\,cent}$ deviation.
The median emulator residual is within the simulation uncertainty, and
unbiased for low aperture masses and $R_{1}=0.5,1.5\,\mathrm{cMpc}$, while
biased low at the $\approx 1\,\mathrm{per\,cent}$ level for
$R_{1}=1.0\,\mathrm{cMpc}$. For abundant, low-aperture mass haloes the
emulator reaches an accuracy of $\approx 2\,\mathrm{per\,cent}$ for the bulk
of the simulations, rarely exceeding the $5\,\mathrm{per\,cent}$ level. The
most significant outlier, M046, which was run with a smaller box size, is
indicated with a red line in the different panels. _Bottom row:_ The Gaussian-
equivalent significance of the likelihood ratio between the simulated data and
the emulator. The median significance is mostly $\approx 1\,\sigma$. The
significance for the individual simulations can vary wildly, but, as the top
row shows, the fractional uncertainty remains reasonable.

The approximate Gaussian process does not sample the simulation inputs, but
instead optimizes the inducing point locations to accurately reproduce the
posterior of the full Gaussian process, i.e. Eq. (39). Hence, we will not
trivially reproduce the simulation aperture mass function. In our set-up, we
first normalize the input parameters, $\mathsf{X}$, so that all parameters lie
between $0$ and $1$. We use 500 inducing points in the ApproximateGP
variational distribution and minimize the marginal likelihood, approximated by
the gpytorch.mlls.VariationalELBO, with the Adam optimizer with a learning
rate of 0.01 and mini-batches of 512 observations each. These settings
resulted in the fastest loss function minimization in a coarse, manual search
for the optimal parameter settings. We resample the initial hyperparameters 5
times from their uniform priors to avoid local minima in the optimization. The
emulator parameters are specified by the inducing point locations,
$\mathbf{u}$, from Eq. (36) in Appendix B, the Gaussian process mean, $\mu$,
from Eq. (29), and the kernel lengthscales and normalization, $\bm{\ell}$ and
$\sigma$, respectively, from Eq. (31). We choose uniform priors
$\mathbf{u}\sim\mathcal{U}(0,1)$, $\mu\sim\mathcal{U}(-1,1)$,
$\bm{\ell}\sim\mathcal{U}(0.05,2.0)$, and
$\sigma^{2}\sim\mathcal{U}(0.05,2.0)$.

In the top row of Fig. 15, we show the resulting absolute deviation between
the emulated latent function, Eq. (26), and the normalized number density from
the simulation for different apertures and all cosmologies. The first three
columns correspond to the different aperture sizes at $z=0.43$, and the final
column shows the median and 16th to 84th percentile scatter for the emulator
at different redshifts. For low aperture masses, the emulator error rarely
exceeds the $5\,\mathrm{per\,cent}$ difference level, and the bulk of the
simulations have residuals within $\pm 2\,\mathrm{per\,cent}$ for the high
abundance aperture mass regime. The median deviation is biased slightly low
for $R_{1}=1.0\,\mathrm{cMpc}$, but it is within $\pm 2\,\mathrm{per\,cent}$
for all aperture sizes and all but the most massive haloes. The bulk of the
simulations lack haloes at the highest aperture masses resulting in the
noticeable downturn. While the fractional deviation becomes large, it is still
within the variance of the simulations, which is shown as the shaded region.

To quantify the quality of the fit in the high-mass tail, we show the
equivalent Gaussian significance of the deviation between the simulated data
and the emulator. We compute the significance by calculating the difference
between the log-likelihood of the measured number counts in the simulation
given the predicted aperture mass function of the emulator,
$\ln\mathcal{L}(N_{\mathrm{true}}|N_{\mathrm{pred}})$, and the log-likelihood
of the emulated number counts,
$\ln\mathcal{L}(N_{\mathrm{pred}}|N_{\mathrm{pred}})$, and converting this
probability ratio into the equivalent Gaussian confidence interval $n\sigma$
around the mean expectation, $\mu$, given by $\ln P(\mu+n\sigma)-\ln P(\mu)$.
We show the significance of the deviation between the emulator and the
simulation data in the bottom row of Fig. 15. Individual simulations behave
erratically for low aperture masses, rapidly oscillating between large and
small significance, but rarely exceeding $3\,\sigma$. The median significance
of all cosmologies, on the other hand, is consistently $\approx 1\,\sigma$.
For the high-aperture mass tail, with large fractional deviations between the
emulator and the simulations, the significance of the deviation $\lesssim
1\,\sigma$, indicating that the emulator captures the trends in the data to
within the shot-noise.

Figure 16: The residuals between the leave-one-out emulator and its held out
simulation for all simulations that are not at the edge of the parameter
space. The cosmological parameters are indicated in the top-left corner of
each panel. All emulators have apertures
$(R_{1},R_{2},R_{\mathrm{m}})=(1,2,3)\,\mathrm{cMpc}$ and have been calibrated
at $z=0.43$. All halo samples were selected with
$m_{\mathrm{200c}}>10^{13.5}\,\mathrm{M_{\odot}}$. Coloured lines indicate the
simulations that are within the first (blue shades) or last (red shades)
decile for the plotted cosmological parameter, with the saturation indicating
the ordering. The black line indicates the median result for all simulations
and the shaded region the 16th to 84th percentile scatter. The black, dash-
dotted lines indicate the $\pm 5\,\mathrm{per\,cent}$ region. Most emulators
are able to reproduce the held out simulation prediction within $\approx
5\,\mathrm{per\,cent}$. The most significant outliers, M012 and M033, are
indicated with red arrows and are simulations that are close to the edge of
the parameter space for several cosmological parameters.

Finally, we also perform a leave-one-out test on all simulations that are not
at the edge of the parameter space for any of the cosmological parameters. In
Fig. 16, we show how accurately the emulator predicts the aperture mass
function for all left-out simulations. We colour the lines for simulations
with cosmological parameters that are within the first or the last decile of
the hypercube with different shades of blue and red, respectively, with darker
shades indicating more significant outliers. The emulator can predict the
outcome of most simulations to within $\approx 5\,\mathrm{per\,cent}$ up to
the tail of the mass function. The most significant deviations are found for
simulations that are close to the edge of the cosmological parameter space in
one or more dimensions. The accuracy achieved by the emulator in the leave-
one-out test indicates that the emulator generalizes well beyond the trained
simulation inputs.