# Threshold studies for a hot beam superradiant laser including an atomic
guiding potential

Martin Fasser1 1Institut für Theoretische Physik, Universität Innsbruck,
Technikerstr. 21a, A-6020 Innsbruck, Austria  Christoph Hotter1 1Institut für
Theoretische Physik, Universität Innsbruck, Technikerstr. 21a, A-6020
Innsbruck, Austria  David Plankensteiner1 1Institut für Theoretische Physik,
Universität Innsbruck, Technikerstr. 21a, A-6020 Innsbruck, Austria  Helmut
Ritsch1 1Institut für Theoretische Physik, Universität Innsbruck,
Technikerstr. 21a, A-6020 Innsbruck, Austria

###### Abstract

Background:
Motivated by the outstanding short time stability and reliable continuous
operation properties of microwave clock masers, intense worldwide efforts
target the first implementation of their optical analogues based on narrow
optical clock transitions and using laser cooled dilute atomic gases. While as
a central line of research large efforts are devoted to create a suitably
dense continuous ultracold and optically inverted atom beam source, recent
theoretical predictions hint at an alternative implementation using a filtered
thermal beam at much higher density. Corresponding numerical studies give
encouraging results but the required very high densities are sensitive to beam
collimation errors and inhomogeneous shifts. Here we present extensive
numerical studies of threshold conditions and the predicted output power of
such a superradiant laser involving realistic particle numbers and velocities
along the cavity axis. Detailed studies target the threshold scaling as a
function of temperature as well as the influence of eliminating the hottest
part of the atomic distribution via velocity filtering and the benefits of
additional atomic beam guiding. Using a cumulant expansion approach allows us
to quantify the significance of atom-atom and atom-field correlations in such
configurations.
Methods:
To enable studies with realistic high particle number we implement a
simulation framework based on a first order as well as a second-order cumulant
expansion of the coupled atom-field dynamics with the atomic center of mass
motion represented by classical trajectories. We assume a velocity filtered
initial thermal distribution and guiding forces from prescribed optical
potentials along the cavity axis, while the transverse atomic motion is not
explicitly included as dynamical variable. Effectively, atomic motion along
the cavity axis induces a time varying atom-cavity coupling determined by the
cavity mode structure and the optical guiding potentials. Our model includes a
continuous effective incoherent pump mechanism. Using such simulations we
study how the intra-cavity photon number depends on the atom number, initial
velocity distribution and atomic guiding. Comparing the results in different
expansion orders including customized mixed order models allows to extract the
relative importance of atom-atom and atom-field correlations in different
operating regimes.
Results:
We predict necessary conditions to achieve a certain threshold photon number
depending on the atomic temperature and density. In particular, we show that
the temperature threshold can be significantly increased by using more atoms.
Interestingly, a velocity filter removing very fast atoms has only almost
negligible influence despite their phase perturbing properties. On the
positive side an additional conservative optical guiding towards cavity mode
antinodes leads to significantly lower threshold and higher average photon
number. Interestingly we see that higher order atom-field and direct atom-atom
quantum correlations play only a minor role in the laser dynamics, which is a
bit surprising in the superradiant regime.
Conclusions:
A hot atomic beam laser operated in the superradiant regime can achieve useful
power levels, if the inhomogeneous atomic broadening from the velocity
distribution along the cavity axis can be compensated by using sufficiently
more atoms. Velocity filtering to remove very hot atoms negatively influences
the average photon number but can potentially reduce power fluctuations as
well as a reduced linewidth and cavity pulling effects. Adding a confining
optical lattice potential, however, creates significant modifications of the
dynamics, especially for higher temperatures where our model simulations
predict significantly higher laser output power.

Corresponding author: Martin Fasser<EMAIL_ADDRESS>

## Introduction

Ever since the demonstration of the superior operation of an optical atomic
clock [1, 2, 3] with respect to microwave implementations people started
wondering about a possible active clock implementation in form of a
superradiant laser [4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19,
20]. In analogy to the very successful hydrogen micro-maser technology [21]
such an active device in principle promises at least similar stability and
accuracy at a reduced technical cost and fragility. Here an implementation
based on a dilute gas of cold clock atoms within a high-Q optical cavity was
predicted to be a very promising path to go [5]. As a particular example,
model setups based on a beam of inverted ultracold atoms traversing an optical
resonator [22, 23, 24] were predicted to yield unprecedented stability and
accuracy. However, as the operating optical wavelength is several orders of
magnitude shorter than the typical maser lines, the resulting technological
challenges are similarly larger.

Already the construction of a suitable beam source is virtually equivalent to
operating a continuous atom laser [25, 26] and thus poses a way more
challenging task than its microwave analog. As the cavity mode volume is also
way smaller the atoms have to be tightly confined before sending them through
the resonator. Yet another problem is creating atomic inversion with minimal
perturbation [27]. Despite intense efforts which recently lead to amazing
technical advancements on guided atom beams [28, 29, 30, 31] and pulsed
superradiant lasing [32, 33, 34, 35, 36, 37, 38, 39], a fully continuous
operating device [9, 19] has not been realized.

In recent calculations it was pointed out that a high precision could already
be achieved using a much hotter atomic ensemble with a sufficiently high
intra-cavity density [40, 41, 42, 24]. Here the sheer number of contributing
atoms allows to cross the laser threshold inducing an efficient collective
synchronization process towards a very narrow effective linewidth. This
approach at least in principle requires significantly less components and
would allow for a rather compact setup involving only a beam oven, an optical
resonator and a pump laser to create inversion. The precise conditions for a
concrete setup, however, cannot be reliably analytically predicted and thus
require large scale numerical simulations [24]. Here we set up a framework for
such simulations which can deal with realistically large atom numbers and
velocities. This allows to identify the minimal operating conditions and study
the influence of extra elements as velocity filters or guiding potentials.

## Model

Figure 1: _Schematic of the system._ Atoms are incoherently pumped with rate
and couple to an optical cavity. The clusters move according to a Gaussian
velocity distribution with standard deviation v, which relates to a certain
temperature $T$. This motion is included along the cavity axis with a position
dependent coupling strength $g$. The atoms decay with the rate and the photons
leak out of the cavity with a rate .

We consider an ensemble of $N$ two-level atoms with a transition frequency of
${}_{\text{a}}$ coupled to a single mode cavity with resonance frequency
${}_{\text{c}}$, see Figure 1. The coherent dynamics of this system is
governed by the Tavis Cummings Hamiltonian, which in the rotating frame of the
atomic transition frequency ${}_{\text{a}}$ reads

$H=\Delta
a^{\dagger}a+\sum^{N}_{j=1}g(x_{j})(a^{\dagger}{}^{\mathrm{ge}}_{j}+{}^{\mathrm{eg}}_{j}a),$
(1)

with the atomic transition operator of the $j$-th atom
${}^{\mathrm{ge}}_{j}=|g\rangle_{j}\langle e|$, the cavity photon creation
(annihilation) operator $a^{\dagger}$ ($a$), and the cavity-atom detuning
$\Delta={}_{\text{c}}-{}_{\text{a}}$. Atoms at a position $x$ along the cavity
axis couple to the cavity according to the mode function $g\cdot
f(x)=g\cdot\cos(2\pi x/\lambda)$, where $g$ is the maximal cavity coupling
strength and the wavelength of the cavity photons. The motion of the atoms is
treated classically by an initial atomic velocity drawn from a Gaussian
distribution with width v corresponding to a temperature $T$. In the beginning
we treat freely moving atoms, later on we investigate the influence of an
external optical potential. In both cases the atomic motion is fully
determined by the initial velocity and position.

Additional to the coherent dynamics we need to account for the dissipative
processes, which are described by Liouvillian term $\mathcal{L}[\rho]$ in the
master equation for the system density matrix

$\dot{\rho}=i[\rho,H]+\mathcal{L}[\rho].$ (2)

In the Born-Markov approximation the Liouvillian super-operator for a jump
operator $J_{i}$ with a corresponding rate i can be written in Lindblad form
as

$\mathcal{L}[\rho]=\sum_{i}\frac{{}_{i}}{2}\left(2J_{i}\rho
J_{i}^{\dagger}-J_{i}^{\dagger}J_{i}\rho-\rho J_{i}^{\dagger}J_{i}\right).$
(3)

The considered dissipative processes are listed in table 1, including cavity
photon losses, individual spontaneous decay as well as incoherent pumping of
the two-level atoms.

Table 1: _Dissipative Processes._ The system features a damped cavity mode, atomic decay and an incoherent pump. $i$ | $c_{i}$ | i | Description
---|---|---|---
1 | $a$ |  | cavity photon losses
2 | ${}_{j}^{\mathrm{ge}}$ |  | decay from $|2\rangle_{j}$ to $|1\rangle_{j}$
3 | ${}_{j}^{\mathrm{eg}}$ |  | excitation of the $j$-th atom

As we want to simulate systems with up to one million atoms, we cannot solve
the full master equation due to exponential scaling of the size of the Hilbert
space with the particle number. To this end we employ a mean field approach
for operator averages. For a not explicitly time-dependent operator
$\mathcal{O}$ we can calculate the time-evolution of its expectation value as

$\frac{\mathrm{d}}{\mathrm{d}t}\langle\mathcal{O}\rangle=\mathrm{Tr}\\{\mathcal{O}\dot{\rho}\\}=i\langle[H,\mathcal{O}]\rangle+\langle\mathcal{\bar{L}}[\mathcal{O}]\rangle$
(4)

with

$\mathcal{\bar{L}}[\mathcal{O}]=\sum_{i}\frac{{}_{i}}{2}\left(2J_{i}^{\dagger}\mathcal{O}J_{i}-J_{i}^{\dagger}J_{i}\mathcal{O}-\mathcal{O}J_{i}^{\dagger}J_{i}\right),$
(5)

where we inserted the master equation (2) into Eq. (4) and used the cyclic
property of the trace. To truncate the set of equations we approximate higher
order quantum correlations with the cumulant expansion method [43, 44].

The first order approximation (mean field) is sufficient to accurately
describe our system (see Appendix), the corresponding set of equations is

$\displaystyle\dot{\langle a\rangle}$
$\displaystyle=-\Bigl{(}\frac{\kappa}{2}+i\Delta\Bigl{)}\langle a\rangle-
ig\sum_{{m}}f(x_{{m}})\langle{}^{\mathrm{ge}}_{{m}}\rangle$ (6a)
$\displaystyle\dot{\langle{}^{\mathrm{ge}}_{{m}}\rangle}$
$\displaystyle=-\frac{\gamma+\nu}{2}\langle{}^{\mathrm{ge}}_{{m}}\rangle+2igf(x_{{m}})(\langle{}^{\mathrm{ee}}_{{m}}\rangle-\frac{1}{2})\langle
a\rangle$ (6b) $\displaystyle\dot{\langle{}^{\mathrm{ee}}_{{m}}\rangle}$
$\displaystyle=-(\gamma+\nu)\langle{}^{\mathrm{ee}}_{{m}}\rangle+igf(x_{{m}})(\langle
a\rangle^{*}\langle{}^{\mathrm{ge}}_{{m}}\rangle-\langle{}^{\mathrm{eg}}_{{m}}\rangle^{*}\langle
a\rangle)+\nu,$ (6c)

where the last two equations exist for each atom from $m=1$ to $N$. The number
of equations in this approximation scales linearly with the atom number $N$.
Still, to simulate a system with several hundred thousand atoms we need
further approximations. To this end, we assume that groups of $K$ atoms behave
in a sufficiently identical way, such that we can describe them with the same
equations. This means, we group atoms into clusters. For the set of equations
in 6, this means that the index $m$ now runs only from $m=1$ to the number of
clusters $N_{\mathrm{Cl}}$, but the source term
$-ig\sum_{{m}}f(x_{{m}})\langle{}^{\mathrm{ge}}_{{m}}\rangle$ in Eq. (6a) gets
enhanced by the clustersize $K$. In the following, we use a clusternumber of
$N_{\mathrm{Cl}}=400$, which is a good compromise between simulating the real
setup and still having a feasible computation time (see Appendix).

## Photon number dependence on temperature

Figure 2: _Photon number dependence on temperature and atom number._ (a)
steady state photon number (mean field, see Appendix for comparison to higher
orders) as a function of the standard deviation of the velocity distribution v
for different atom numbers. (b) Photon number normalized to
${}_{\mathrm{v}}=0$ (maximum). The red dots indicate the value at which the
curves have dropped to half of their initial value. After normalizing, the
curves do not longer align, hinting to a nonlinear relation between atom and
photon number. (c) Dependence of the threshold value of the velocity spread
${}_{\mathrm{v}}^{2}\sim T$ on the number of atoms. Larger atom numbers allow
for larger possible values of v, while maintaining a lot of photons inside the
cavity. The parameters are $\nu=0.5$, $\Delta=0$, $g=0.00136$ and $\kappa=1$.
We neglect since it is much smaller than , see Eq. (6).

First, we investigate the influence of the atom number $N$ and the temperature
$T$ on the steady state photon number. The results in mean field (for a
comparison to higher orders see Appendix) are depicted in Figure 2.
Unsurprisingly, we always end up with the highest number of photons for the
lowest temperature. The drop-off of the photon number for higher temperature
can be explained by the Doppler effect, as higher velocities of the individual
atoms lead to them being shifted out of resonance. Comparing each curve to the
others, we can investigate the significant influence of the atom number on the
photon number: Higher atom numbers yield higher photon numbers. While this is
not surprising, the crucial aspect of this result becomes apparent, when we
normalize the curves to their photon number value for zero temperature (Figure
2(b)): The lacking coincidence of the curves hints to a nonlinear relation
between atoms and photons. However, for low atom numbers we see the curves
dropping to virtually 0 photons as early as for
${}_{\mathrm{v}}/(\lambda\kappa)\approx 0.05$, whereas for higher atom numbers
such a temperature has hardly any effect on the photon number. For example,
for $10^{6}$ atoms there are still approximately $95\%$ of the photons at zero
temperature in the cavity. The red dots indicate the points, at which the
photon number has dropped to half their maximal value at zero temperature for
each individual curve. These values are depicted in Figure 2(c), they
initially follow a quadratic relation, while for higher atom numbers the
relation seems linear. From that graph we can clearly see, that using higher
atom numbers allows a higher temperature while still having a decent amount of
photons inside the cavity. This could be beneficial to experimental setups,
where it is generally difficult to cool the atoms down to low temperature.
However, it remains to be seen what this means to the linewidth of the
outcoming photons, as higher temperature usually leads to Doppler broadening
of the linewidth.

## Adding a velocity filter

Figure 3: _Photon number for different velocity cutoffs._ (a) Mean steady
state photon number. (b) Normalized to maximum photon number. Cutting off the
the high velocity part of the distribution always yields a smaller steady
state photon number, but it remains at a high level for increasing v. We only
depict the results in mean field, as the results in mixed order almost
perfectly align. The parameters are the same as in Figure 2.

An interesting idea is to consider a cutoff in the velocity distribution for
the atoms. Experimentally, this may be achieved by using velocity filter for
the atomic beam to get rid of the atoms with very high velocity, creating a
smaller but colder ensemble. This could be more advantageous than trying to
narrow the velocity distribution by cooling, as cooling the atoms proposes an
experimental challenge. To simulate such a filter, we start with the same
number of atoms as before, but exclude all atoms from the dynamics with
velocities higher than the cutoff. We depict the results in Figure 3 for a
cutoff at ${}_{\text{v}}$ ($\approx 68\%$ slowest atoms) and at
${}_{\text{v}}/2$ ($\approx 38\%$ slowest atoms). To compare, we also depict
the previous results without cutoff. In Figure 3(a) the steady state photon
number is shown for the different cutoffs. By cutting off the high velocity
part of the atomic distribution, we lose only fast atoms, but it always leads
to a decrease in photon number. A velocity filter therefore only reduces the
photon number. It might still be useful with respect to the linewidth, which
is not calculated at this point.

However, the photon number seems to be more stable with increasing
temperature, which can be seen in Figure 3(b), where we normalize the graphs
to their highest photon number (at zero temperature). While the photon number
significantly drops already for low temperature (e.g. ${}_{\mathrm{v}}\approx
0.1$) for the uncut case, it stays relatively stable for the other cases,
where we exclude the fast atoms.

In summary, including a velocity filter will negatively affect the photon
number, but allow it to stay at a relatively stable value for low
temperatures. However, it remains to be seen what effect such a velocity
filter has on the linewidth of the atoms. In general, one would expect a
velocity filter to be beneficial for the linewidth reduction, as the excluded,
fast atoms increase the linewidth by means of the Doppler effect.

## Dynamics including an optical lattice potential

So far, we assumed that the atoms move freely without any potential. After
some time at non-zero temperature, the atoms will be randomly distributed over
the mode function. To compare that situation to the zero-temperature case with
non-moving atoms, we always place the atoms equidistantly distributed over the
mode function. Therefore, the differences in photon number stem from the
movement of the atoms itself, not from the somewhat artificial location of
atoms on the mode function. The coupling between atoms and cavity field can be
controlled (atleast for low temperature) by including an optical lattice
potential, as the atoms are trapped at the potential minimums. Aligning these
potential minimums with the antinodes of the mode function leads to enhanced
coupling compared to the random sampling of the mode function in the case
where atoms move freely. The potential $V(x)$ has the form $V(x)=\Delta
E\cdot\cos^{2}(\frac{4\pi}{{}_{\mathrm{latt}}}x)$ with spatial periodicity
latt and energy barrier $\Delta E$. By using the Hamiltonian formalism one
ends up with two linear differential equations of first order for each
cluster. These additional equations can simply be added to the equations in
mean field (6) for simulation purposes. To keep the computation time low, we
now simulated only $N_{\mathrm{Cl}}=100$ clusters instead of 400 and only
averaged over 10 trajectories instead of 50.
For the optimal case, we assume an optical potential with exactly half the
wavelength of the mode function, that way the minima of the optical potential
are situated at the maximum values of $|f(x_{m})|$. The results of such a
simulation are depicted in Figure 4(a). To compare to the previous results, we
include the graph that we already calculated, which is the result for an
energy barrier of $\Delta E=0$. The energy is expressed in terms of the recoil
energy $E_{\mathrm{rec}}=\frac{\mathord{\mathchar 0\relax\symmathdesignA
07Eh}^{2}k^{2}}{2m}$, where $k=2\pi/\lambda$, $v_{\mathrm{rec}}$ is the
associated recoil velocity and $m$ is the mass of the atoms. To get exemplary
values for these energy scales we use the ${}^{40}\mathrm{Ca}$
intercombination line between ${}^{1}\mathrm{S}_{0}$ and
${}^{3}\mathrm{P}_{1}$ with $\lambda=657\,\mathrm{nm}$ and $\Gamma=2\pi\cdot
375\,\mathrm{Hz}$, the cavity linewidth is $\kappa=2\pi\cdot
120\,\mathrm{kHz}$. To compare the two energy scales of the energy barrier and
the velocity of the atoms, we color the data points according to the fraction
of trapped atoms at each particular energy barrier and velocity. For high
energy barriers and low velocities we have a lot of trapped atoms, these data
points are colored blue, whereas for low energy barriers and high velocities
almost no atoms are trapped, these points are colored yellow. As we can see
from Figure 4, the results change significantly by including a barrier.

Figure 4: _Quasi-steady state photon number including a optical lattice
potential with wavelength ${}_{\mathrm{latt}}=0.5\,\lambda$ (a) and
${}_{\mathrm{latt}}=0.99\,\lambda$ (b) with varying energy barriers._ The
color of the data points indicates the fraction of trapped atoms, where at
blue most of the atoms are trapped, while at yellow most of the atoms can
escape their trap. (a) ${}_{\mathrm{latt}}=0.5\,\lambda$: Due to the choice of
latt, for low temperatures most of the atoms are trapped at the antinode of
the mode-function, leading to a higher atom-cavity coupling and higher photon
number. For higher temperatures, most of the atoms escape their trap, but the
photon number stays significantly above the barrierless case. This is due to
the energy barriers, which slow the atoms down, yielding lower Doppler effect
and therefore higher photon number. (b) ${}_{\mathrm{latt}}=0.99\,\lambda$:
For low temperatures the results coincide, regardless of the chosen energy
barrier. This means, that the atom-cavity coupling is the same as in the
barrierless case, where we distributed the atoms equidistantly over a mode
function. For higher temperatures, most of the atoms escape their trap, but
the photon number behaves similarly to the case with
${}_{\mathrm{latt}}=0.5\,\lambda$. Once again, the atoms get slowed by the
barrier, eventually leading to significantly higher photon number compared to
the barrierless case. The parameters are the same as in Figure 2.

Let us first look at the case of low velocities/low temperatures: When the
barrier is sufficiently high, then almost all the atoms are trapped at the
antinodes of the mode function, leading to high coupling between atoms and
cavities and eventually to higher photon number. For higher
velocities/temperature, the number of photons stays significantly above the
comparison graph for no energy barrier. On the one hand, this is because there
are still atoms, that are trapped near the antinode of the mode function. On
the other hand, there is a different effect at play, of which we can get a
better insight when choosing the wavelength of the lattice slightly different
from the wavelength of the mode function. Such a system may be of interest, as
there usually is a mismatch between the effective optical lattice periodicity
and the optical wavelength.
The results can be seen in Figure 4(b), where we chose
${}_{\mathrm{latt}}=0.99\,\lambda$. Note, how the results converge for very
low temperature, regardless of the value of the energy barrier $\Delta E$.
This means, that the coupling is the same as in the barrierless case, where we
distributed the atoms equidistantly over the mode function. So by choosing the
wavelength of the optical lattice only slightly different to the wavelength of
the mode function, we once again sample an average over the mode function, as
long as the atoms sit on sufficiently many sites. Therefore the whole
advantage of keeping atoms near the antinode of the mode function disappears
and the coupling between the atoms and the cavity is the same compared to the
coupling in the barrierless case. However, increasing the temperature, we
still see a significant difference: Even though the average coupling is the
same, the photon number is much higher. The atoms are decelerated, whenever
they go over a potential barrier, yielding a reduced average velocity, which
in turn decreases the Doppler effect and increases the photon number.

In summary, for low velocities/temperatures, the change in photon number
depends mostly on the change of coupling of the atoms, which again depends on
the location of the minima of the optical potential compared to the antinodes
of the mode function. For higher velocities/temperatures, the other effect of
having a lower average velocity due to the barriers comes also into play. The
lower average velocity leads to a decreased Doppler effect and therefore a
significant increase of photon number.

## Conclusions

We studied incoherently pumped two-level-atoms inside a single mode cavity. In
particular, we investigated the influence of temperature on the photon number.
For high temperatures, where the atoms are fast, the photon number is
suppressed because of the Doppler effect. However, we saw that increasing the
atom number, to get sufficient collective atom-cavity-coupling, can compensate
this suppression. Including a velocity filter leads to a lower amount of atoms
participating in the dynamics and therefore, as the number of atoms is a
significant parameter, to a lower photon number. We plan to calculate the
spectrum to see if excluding the fast atoms has a decreasing effect on the
linewidth.

With an optical lattice potential the coupling between the atoms and the
cavity field can be controlled. For low temperatures (compared to the energy
barrier), most atoms are trapped and the atom-cavity coupling (and in
consequence the photon number) can be engineered. However, for high
temperatures, even though most of the atoms are not trapped, the photon number
significantly increases compared to the situation without a potential. In some
cases, where virtually no photons are generated anymore for a certain v
without the potential, we still have a macroscopic photon occupation number,
even if v is 4 times as high, by including the optical potential. Therefore
including an optical potential becomes very beneficial for higher
temperatures.

## Appendix. Equations in higher orders

As mentioned in the main text, for our parameter regime the physics is
accurately described by mean field. To verify this, we derive the equations
for higher orders in the cumulant expansion and compare the numerical results
with the mean field. For second order, one derives equations for two-operator
averages. These will depend on other averages of two-operator and three-
operator products. In order to close the set of differential equations one
needs to derive equations for the two-operator averages and approximate the
three-operator averages according to [43]

$\langle abc\rangle\approx\langle ab\rangle\langle c\rangle+\langle
ac\rangle\langle b\rangle+\langle cb\rangle\langle a\rangle-2\langle
a\rangle\langle b\rangle\langle c\rangle.$ (7)

Using the phase invariance of the system [5, 44] one ends up with the set of
equations in second order of the cumulant expansion:

$\displaystyle\dot{\langle a^{\dagger}a\rangle}=$
$\displaystyle-\Bigl{(}\kappa-2i\Delta\Bigl{)}\langle
a^{\dagger}a\rangle+igK\sum_{{m}}f(x_{{m}})(\langle
a{}^{\mathrm{eg}}_{{m}}\rangle-\langle a^{\dagger}{}^{\mathrm{ge}}_{{m}})$
(8a) $\displaystyle\dot{\langle a^{\dagger}{}^{\mathrm{ge}}_{{m}}\rangle}=$
$\displaystyle-\Bigl{(}\frac{\gamma+\kappa+\nu}{2}-i\Delta\Bigl{)}\langle
a{}^{\mathrm{eg}}_{{m}}\rangle+igf(x_{{m}})(\langle{}^{\mathrm{ee}}_{\mathrm{m}}\rangle-\langle
a^{\dagger}a\rangle)+ig\sum_{{j=1}}^{N}f(x_{{j}})\langle{}^{\mathrm{eg}}_{{j}}{}^{\mathrm{ge}}_{{m}}\rangle$
(8b) $\displaystyle+2igf(x_{{m}})\langle
a^{\dagger}a\rangle\langle{}^{\mathrm{ee}}_{{m}}\rangle$
$\displaystyle\dot{\langle{}^{\mathrm{ee}}_{{m}}\rangle}=$ $\displaystyle
igf(x_{{m}})(\langle a^{\dagger}{}^{\mathrm{ge}}_{{m}}\rangle-\langle
a{}^{\mathrm{eg}}_{{m}}\rangle)-(\gamma+\nu)\langle{}^{\mathrm{ee}}_{{m}}\rangle+\nu$
(8c)
$\displaystyle\dot{\langle{}^{\mathrm{ge}}_{k}{}^{\mathrm{eg}}_{l}\rangle}=$
$\displaystyle-(\nu+\gamma)\langle{}^{\mathrm{ge}}_{k}{}^{\mathrm{eg}}_{l}\rangle-
igf(x_{{k}})\langle a{}^{\mathrm{eg}}_{l}\rangle+igf(x_{{l}})\langle
a{}^{\mathrm{eg}}_{k}\rangle-\nu\langle{}^{\mathrm{ge}}_{k}{}^{\mathrm{eg}}_{l}\rangle$
(8d) $\displaystyle+2igf(x_{{k}})\langle{}^{\mathrm{ee}}_{k}\rangle\langle
a{}^{\mathrm{eg}}_{j}\rangle-2igf(x_{{l}})\langle{}^{\mathrm{ee}}_{j}\rangle\langle
a{}^{\mathrm{eg}}_{k}\rangle.$

The advantage of this second order is, that we can include atom-field and
atom-atom correlations, however, we pay the price in form of a higher
computation time, as the number of equations now scales quadratically with the
cluster number instead of linear as above. The terms involving two atom
operators $\langle{}_{i}{}_{j}\rangle$ are responsible for the quadratic
scaling, as both $i$ and $j$ are between 1 and $N$. In an attempt to take the
advantage of the upsides of both first order and second order, we approximate
two-operator products if they both are atom operators by splitting them, but
keep the correlations between atom and field, we call this particular
truncation "mixed order":

$\displaystyle\dot{\langle a\rangle}=$
$\displaystyle-\Bigl{(}\frac{\kappa}{2}-i\Delta\Bigl{)}\langle a\rangle-K\cdot
ig\sum_{{m}}f(x_{{m}})\langle{}^{\mathrm{ge}}_{{m}}\rangle$ (9a)
$\displaystyle\dot{\langle a^{\dagger}a\rangle}=$
$\displaystyle-\Bigl{(}\kappa-2i\Delta\Bigl{)}\langle
a^{\dagger}a\rangle+igK\sum_{{m}}f(x_{{m}})(\langle
a{}^{\mathrm{eg}}_{{m}}\rangle-\langle a^{\dagger}{}^{\mathrm{ge}}_{{m}})$
(9b) $\displaystyle\dot{\langle aa\rangle}=$
$\displaystyle-\Bigl{(}\kappa-2i\Delta\Bigl{)}\langle
aa\rangle-2igK\sum_{{m}}f(x_{{m}})\langle a{}^{\mathrm{ge}}_{{m}}\rangle$ (9c)
$\displaystyle\dot{\langle a{}^{\mathrm{eg}}_{{m}}\rangle}=$
$\displaystyle-\Bigl{(}\frac{\gamma+\kappa+\nu}{2}-i\Delta\Bigl{)}\langle
a{}^{\mathrm{eg}}_{{m}}\rangle+igf(x_{{m}})\langle a^{\dagger}a\rangle-
ig\sum_{{j=1}}^{N_{\mathrm{Cl}}}\sum_{i=1}^{K}f(x_{{j}})\langle{}^{\mathrm{eg}}_{{m}}{}^{\mathrm{ge}}_{{j,i}}\rangle$
(9d) $\displaystyle-2igf(x_{{m}})\langle
a^{\dagger}a{}^{\mathrm{ee}}_{{m}}\rangle$ $\displaystyle\dot{\langle
a{}^{\mathrm{ge}}_{{m}}\rangle}=$
$\displaystyle-\Bigl{(}\frac{\gamma+\kappa+\nu}{2}-i\Delta\Bigl{)}\langle
a{}^{\mathrm{ge}}_{{m}}\rangle-igf(x_{{m}})\langle aa\rangle-
ig\sum_{{j=1}}^{N_{\mathrm{Cl}}}\sum_{i=1}^{K}f(x_{{j}})\langle{}^{\mathrm{ge}}_{{m}}{}^{\mathrm{ge}}_{{j,i}}\rangle$
(9e) $\displaystyle+2igf(x_{{m}})\langle aa{}^{\mathrm{ee}}_{{m}}\rangle$
$\displaystyle\dot{\langle{}^{\mathrm{ee}}_{{m}}\rangle}=$ $\displaystyle
igf(x_{{m}})(\langle a^{\dagger}{}^{\mathrm{ge}}_{{m}}\rangle-\langle
a{}^{\mathrm{eg}}_{{m}}\rangle)-(\gamma+\nu)\langle{}^{\mathrm{ee}}_{{m}}\rangle+\nu$
(9f) $\displaystyle\dot{\langle{}^{\mathrm{ge}}_{{m}}\rangle}=$
$\displaystyle-\frac{\gamma+\nu}{2}\langle{}^{\mathrm{ge}}_{{m}}\rangle-
igf(x_{{m}})\langle a\rangle+2igf(x_{{m}})\langle
a{}^{\mathrm{ee}}_{{m}}\rangle$ (9g) $\displaystyle\dot{\langle
a{}^{\mathrm{ee}}_{{m}}\rangle}=$ $\displaystyle-(\gamma+\nu)\langle
a{}^{\mathrm{ee}}_{{m}}\rangle-\frac{\kappa}{2}\langle
a{}^{\mathrm{ee}}_{{m}}\rangle+\nu\langle a\rangle-igf(x_{{m}})\langle
aa{}^{\mathrm{eg}}_{{m}}\rangle$ (9h) $\displaystyle-
ig\sum_{{j=1}}^{N_{\mathrm{Cl}}}\sum_{i=1}^{K}f(x_{m})\langle{}^{\mathrm{ee}}_{{m}}{}^{\mathrm{ge}}_{{j,i}}\rangle+igf(x_{{j}})\langle
a^{\dagger}a{}^{\mathrm{ge}}_{{m}}\rangle.$

For the sake of better readability we refrain from carrying out all the
expansions of three-operator expectation values like $\langle
aa{}^{\mathrm{ee}}_{\mathrm{m}}\rangle$. In order to end up with a closed set
of equations, one needs to approximate the three-operator expectation values
by two-operator expectation values according to Eq. (7). By neglecting the
atom-atom correlations we get rid of the quadratic scaling, but we are still
able to keep higher order correlations, namely the ones between atom and
field.

Grouping the atoms into clusters is more involved in the second order and
mixed order as opposed to first order, one has to carefully take into account
for example the correlations within a cluster [17, 45].

## Determining the steady state photon number

Figure 5: _Determining the steady state photon number._ Mean photon number
time evolution for different v with (a) 200 atom clusters and (b) 1000 atom
clusters. After a certain time, the system reaches its quasi steady state. The
fluctuations become smaller if we divide the velocity distribution into a
larger number of atom clusters. The other parameters are the same as in Figure
2 for $N=10^{6}$.

In order to study the photon number values depending on different parameters
we initialize the set of differential equations with a certain parameter set
and let the system evolve until it reaches its quasi-steady-state. The
velocities are drawn from a Gaussian distribution with width ${}_{\text{v}}$
corresponding to a temperature. We average over the photon number from the
timepoint, at which the quasi-steady-state has been reached, see Figure 5.
Using more clusters leads to much smaller fluctuations, however, this also
requires more computation time. Finding a compromise between high
clusternumbers and reasonable computation time is discussed in the next
section. Moreover, one can imagine that certain samples of the distribution
exhibit special properties, in order to average them out we use 50 (unless
indicated otherwise) different samples and also average over them.

## Comparison expansion orders and number of clusters

Limited computation time forces us to make compromises, as both higher
expansion orders and higher number of clusters describe the physical reality
better, but also require more computation time.

In the following we compare the results for different expansion orders and
number of cluster. In Figure 6(a) the solid line depicts the first order, the
dotted line mixed order and the dashed line represents the data in second
order. We see that all three lines almost perfectly align. Therefore, using
higher expansion orders does not lead to significant differences and the first
order approximation captures the essential physics in our parameter regime.

In Figure 6(b) we depict the results in first order for different
clusternumbers ranging from $N_{\mathrm{Cl}}=50$ to $N_{\mathrm{Cl}}=1000$.
Here we see, that varying the number of clusters $N_{\mathrm{Cl}}$ does have
an influence on the results. We chose $N_{\mathrm{Cl}}=400$, as the difference
between choosing $N_{\mathrm{Cl}}=400$ and as high as $N_{\mathrm{Cl}}=1000$
is only marginal, while the computation time significantly increases.

Figure 6: _Expansion orders and number of clusters._ (a) The solid line
represents mean field, the dotted line is mixed order and the dashed line
depicts the result in second order, we use $N_{\mathrm{Cl}}=25$. The results
for all orders agree very well. (b) The chosen number of clusters
$N_{\mathrm{Cl}}$ makes a difference, especially for medium temperatures.
However, increasing the number of clusters from $N_{\mathrm{Cl}}\approx 400$
has only a marginal effect on the results. The other parameters are the same
as in Figure 2 for $N=10^{6}$.

## Acknowledgements

We acknowledge funding from the Austrian Science Fund (FWF) doctoral college
No. DK-ALM W1259-N27 (M.F.) and the European Union’s Horizon 2020 research and
innovation program under the Grant Agreement No. 820404 iqClock (C. H., D. P.,
H. R.).

## Data availability

## Underlying data

Figshare: Superradiant$\\_$laser$\\_$Figures.
https://doi.org/10.6084/m9.figshare.c.6781920.v1 [46].
This project contains the following underlying data:

  * •

Data used in Figures 2-6. All data are stored using DelimitedFiles.jl package
in Julia.

  * •

Data are stored together with the source code (see below), specific files to
readout and plot the data are provided.

Data are available under the terms of the Creative Commons Zero "No rights
reserved" data waiver (CC0 1.0 Public domain dedication).

## Software availability

  * •

Source code available from: https://github.com/martinf97/SRL_paper

  * •

Archived source code at time of publication:
https://doi.org/10.5281/zenodo.8232295 [47]

  * •

License: MIT License

## References

  * [1] M. Takamoto, F.-L. Hong, R. Higashi, and H. Katori, “An optical lattice clock,” Nature, vol. 435, no. 7040, pp. 321–324, 2005.
  * [2] B. Bloom, T. Nicholson, J. Williams, S. Campbell, M. Bishof, X. Zhang, W. Zhang, S. Bromley, and J. Ye, “An optical lattice clock with accuracy and stability at the 10- 18 level,” Nature, vol. 506, no. 7486, pp. 71–75, 2014.
  * [3] A. D. Ludlow, M. M. Boyd, J. Ye, E. Peik, and P. O. Schmidt, “Optical atomic clocks,” Reviews of Modern Physics, vol. 87, no. 2, p. 637, 2015.
  * [4] F. Haake, M. I. Kolobov, C. Fabre, E. Giacobino, and S. Reynaud, “Superradiant laser,” Phys. Rev. Lett., vol. 71, pp. 995–998, Aug 1993.
  * [5] D. Meiser, J. Ye, D. Carlson, and M. Holland, “Prospects for a millihertz-linewidth laser,” Physical review letters, vol. 102, no. 16, p. 163601, 2009.
  * [6] J. Chen, “Active optical clock,” Chinese Science Bulletin, vol. 54, no. 3, pp. 348–352, 2009.
  * [7] D. Meiser and M. Holland, “Steady-state superradiance with alkaline-earth-metal atoms,” Physical Review A, vol. 81, no. 3, p. 033847, 2010.
  * [8] D. Meiser and M. Holland, “Intensity fluctuations in steady-state superradiance,” Physical Review A, vol. 81, no. 6, p. 063827, 2010.
  * [9] J. G. Bohnet, Z. Chen, J. M. Weiner, D. Meiser, M. J. Holland, and J. K. Thompson, “A steady-state superradiant laser with less than one intracavity photon,” Nature, vol. 484, no. 7392, pp. 78–81, 2012.
  * [10] V. Vuletic, “An almost lightless laser,” Nature, vol. 484, no. 7392, pp. 43–44, 2012.
  * [11] T. Maier, S. Kraemer, L. Ostermann, and H. Ritsch, “A superradiant clock laser on a magic wavelength optical lattice,” Optics express, vol. 22, no. 11, pp. 13269–13279, 2014.
  * [12] Y. Zhang, Y.-X. Zhang, and K. Mølmer, “Monte-carlo simulations of superradiant lasing,” New Journal of Physics, vol. 20, p. 112001, nov 2018\.
  * [13] K. Debnath, Y. Zhang, and K. Mølmer, “Lasing in the superradiant crossover regime,” Physical Review A, vol. 98, no. 6, p. 063837, 2018.
  * [14] C. Hotter, D. Plankensteiner, L. Ostermann, and H. Ritsch, “Superradiant cooling, trapping, and lasing of dipole-interacting clock atoms,” Optics express, vol. 27, no. 22, pp. 31193–31206, 2019.
  * [15] A. Gogyan, G. Kazakov, M. Bober, and M. Zawada, “Characterisation and feasibility study for superradiant lasing in 40ca atoms,” Opt. Express, vol. 28, pp. 6881–6892, Mar 2020.
  * [16] Y. Zhang, C. Shan, and K. Mølmer, “Ultranarrow superradiant lasing by dark atom-photon dressed states,” Physical Review Letters, vol. 126, no. 12, p. 123602, 2021.
  * [17] A. Bychek, C. Hotter, D. Plankensteiner, and H. Ritsch, “Superradiant lasing in inhomogeneously broadened ensembles with spatially varying coupling [version 2; peer review: 2 approved],” Open Research Europe, vol. 1, no. 73, 2021.
  * [18] G. A. Kazakov, S. Dubey, A. Bychek, U. Sterr, M. Bober, and M. Zawada, “Ultimate stability of active optical frequency standards,” Phys. Rev. A, vol. 106, p. 053114, Nov 2022.
  * [19] S. L. Kristensen, E. Bohr, J. Robinson-Tait, T. Zelevinsky, J. W. Thomsen, and J. H. Müller, “Subnatural linewidth superradiant lasing with cold ${}^{88}\mathrm{Sr}$ atoms,” Phys. Rev. Lett., vol. 130, p. 223402, May 2023\.
  * [20] J. Zhang, T. Shi, J. Miao, and J. Chen, “The development of active optical clock,” AAPPS Bulletin, vol. 33, no. 1, p. 10, 2023.
  * [21] H. M. Goldenberg, D. Kleppner, and N. F. Ramsey, “Atomic hydrogen maser,” Phys. Rev. Lett., vol. 5, pp. 361–362, Oct 1960.
  * [22] G. A. Kazakov and T. Schumm, “Active optical frequency standard using sequential coupling of atomic ensembles,” Physical Review A, vol. 87, no. 1, p. 013821, 2013.
  * [23] G. A. Kazakov and T. Schumm, “Active optical frequency standards using cold atoms: Perspectives and challenges,” in 2014 European Frequency and Time Forum (EFTF), pp. 411–414, 2014.
  * [24] M. Tang, S. A. Schäffer, and J. H. Müller, “Prospects of a superradiant laser based on a thermal or guided beam of ${}^{88}\mathrm{Sr}$,” Phys. Rev. A, vol. 106, p. 063704, Dec 2022.
  * [25] T. Salzburger and H. Ritsch, “Atom-photon pair laser,” Phys. Rev. A, vol. 75, p. 061601, Jun 2007.
  * [26] M. Holland, K. Burnett, C. Gardiner, J. I. Cirac, and P. Zoller, “Theory of an atom laser,” Phys. Rev. A, vol. 54, pp. R1757–R1760, Sep 1996.
  * [27] C. Hotter, D. Plankensteiner, G. Kazakov, and H. Ritsch, “Continuous multi-step pumping of the optical clock transition in alkaline-earth atoms with minimal perturbation,” Optics Express, vol. 30, no. 4, pp. 5553–5568, 2022.
  * [28] S. Bennetts, C.-C. Chen, B. Pasquiou, and F. Schreck, “Steady-state magneto-optical trap with 100-fold improved phase-space density,” Phys. Rev. Lett., vol. 119, p. 223202, Dec 2017.
  * [29] C.-C. Chen, S. Bennetts, R. G. Escudero, B. Pasquiou, and F. Schreck, “Continuous guided strontium beam with high phase-space density,” Phys. Rev. Applied, vol. 12, p. 044014, Oct 2019.
  * [30] R. G. Escudero, C.-C. Chen, S. Bennetts, B. Pasquiou, and F. Schreck, “Steady-state magneto-optical trap of fermionic strontium on a narrow-line transition,” Phys. Rev. Res., vol. 3, p. 033159, Aug 2021.
  * [31] C.-C. Chen, R. González Escudero, J. Minář, B. Pasquiou, S. Bennetts, and F. Schreck, “Continuous bose–einstein condensation,” Nature, vol. 606, no. 7915, pp. 683–687, 2022.
  * [32] M. A. Norcia, M. N. Winchester, J. R. Cline, and J. K. Thompson, “Superradiance on the millihertz linewidth strontium clock transition,” Science advances, vol. 2, no. 10, p. e1601231, 2016.
  * [33] M. A. Norcia and J. K. Thompson, “Cold-strontium laser in the superradiant crossover regime,” Physical Review X, vol. 6, no. 1, p. 011025, 2016.
  * [34] M. A. Norcia, J. R. Cline, J. A. Muniz, J. M. Robinson, R. B. Hutson, A. Goban, G. E. Marti, J. Ye, and J. K. Thompson, “Frequency measurements of superradiance from the strontium clock transition,” Physical Review X, vol. 8, no. 2, p. 021036, 2018.
  * [35] T. Laske, H. Winter, and A. Hemmerich, “Pulse delay time statistics in a superradiant laser with calcium atoms,” Physical Review Letters, vol. 123, no. 10, p. 103601, 2019.
  * [36] S. A. Schäffer, M. Tang, M. R. Henriksen, A. A. Jørgensen, B. T. Christensen, and J. W. Thomsen, “Lasing on a narrow transition in a cold thermal strontium ensemble,” Physical Review a, vol. 101, no. 1, p. 013819, 2020.
  * [37] M. Tang, S. A. Schäffer, A. A. Jørgensen, M. R. Henriksen, B. T. Christensen, J. H. Müller, and J. W. Thomsen, “Cavity-immune spectral features in the pulsed superradiant crossover regime,” Physical Review Research, vol. 3, no. 3, p. 033258, 2021.
  * [38] C. Hotter, L. Ostermann, and H. Ritsch, “Cavity sub- and superradiance for transversely driven atomic ensembles,” Phys. Rev. Res., vol. 5, p. 013056, Jan 2023.
  * [39] E. Bohr, S. L. Kristensen, C. Hotter, S. A. Schäffer, J. Robinson-Tait, J. W. Thomsen, T. Zelevinsky, H. Ritsch, and J. H. Müller, “Collectively enhanced ramsey readout by cavity sub-to superradiant transition,” arXiv preprint arXiv:2306.12544, 2023.
  * [40] H. Liu, S. B. Jäger, X. Yu, S. Touzard, A. Shankar, M. J. Holland, and T. L. Nicholson, “Rugged mhz-linewidth superradiant laser driven by a hot atomic beam,” Physical Review Letters, vol. 125, no. 25, p. 253602, 2020\.
  * [41] S. B. Jäger, H. Liu, A. Shankar, J. Cooper, and M. J. Holland, “Regular and bistable steady-state superradiant phases of an atomic beam traversing an optical cavity,” Phys. Rev. A, vol. 103, p. 013720, Jan 2021.
  * [42] S. B. Jäger, H. Liu, J. Cooper, T. L. Nicholson, and M. J. Holland, “Superradiant emission of a thermal atomic beam into an optical cavity,” Phys. Rev. A, vol. 104, p. 033711, Sep 2021.
  * [43] R. Kubo, “Generalized cumulant expansion method,” Journal of the Physical Society of Japan, vol. 17, no. 7, pp. 1100–1120, 1962.
  * [44] D. Plankensteiner, C. Hotter, and H. Ritsch, “Quantumcumulants.jl: A julia framework for generalized mean-field equations in open quantum systems,” 2021\.
  * [45] M. Fasser, “Simulating a superradiant laser based on a thermal atomic beam,” Master’s thesis, University of Innsbruck, Innsbruck, Austria, 3 2022.
  * [46] M. Fasser, “Superradiant_laser_figures,” Figshare, https://doi.org/10.6084/m9.figshare.c.6781920.v1, 2023\.
  * [47] M. Fasser, “Srl_paper,” Zenodo, https://doi.org/10.5281/zenodo.8232295, 2023\.