%~Mouliné par MaN_auto v.0.40.4 (2aa04cce) 2026-03-25 11:13:30 \documentclass[CRPHYS,Unicode,biblatex,published]{cedram} \addbibresource{CRPHYS_Delehaye_20250725.bib} \usepackage{subcaption} \usepackage{overpic} \usepackage{siunitx} \newcommand{\dif}{\mathop{}\!{\operatorfont{d}}} \newcommand{\bbR}{\mathbb{R}} \newcommand{\crit}{\mathrm{crit}} %\newcommand{\Ytterbium}{\prescript{171}{}{\smash{\mathrm{Yb}}\vphantom{\mathrm{Sr}}}} \newcommand{\Ytterbium}{{}^{171}\mathrm{Yb}} \newcommand{\StrontiumS}{{}^{87}\mathrm{Sr}} \newcommand{\StrontiumE}{{}^{88}\mathrm{Sr}} \newcommand{\groundstate}{{}^{1}\mathrm{S}_0} \newcommand{\excitedstate}{{}^{3}\mathrm{P}_0} \newcommand{\steady}{\mathrm{st}} \makeatletter \renewcommand\p@subfigure{} \makeatother %%%--------------------------------------------------------------------------- %%% DELIMITERS \usepackage{mathtools} \DeclarePairedDelimiter{\parens}{\lparen}{\rparen} \DeclarePairedDelimiter{\braces}{\{}{\}} \DeclarePairedDelimiter{\abs}{\lvert}{\rvert} \DeclarePairedDelimiter{\bracks}{[}{]} \DeclarePairedDelimiter{\angles}{\langle}{\rangle} \let\braket\angles \DeclarePairedDelimiter{\bra}{\langle}{\rvert} \DeclarePairedDelimiter{\ket}{\lvert}{\rangle} %%%--------------------------------------------------------------------------- \graphicspath{{./figures/}} \newcommand*{\mk}{\mkern -1mu} \newcommand*{\Mk}{\mkern -2mu} \newcommand*{\mK}{\mkern 1mu} \newcommand*{\MK}{\mkern 2mu} \hypersetup{urlcolor=purple, linkcolor=blue, citecolor=red} \let\oldtilde\tilde \renewcommand*{\tilde}[1]{\mathchoice{\widetilde{#1}}{\widetilde{#1}}{\oldtilde{#1}}{\oldtilde{#1}}} \title{A framework for continuous superradiant laser operation via sequential transport of atoms} \alttitle{Un cadre pour le fonctionnement continu d’un laser superradiant par transport séquentiel d’atomes} \author{\firstname{Jana} \lastname{El Badawi}\IsCorresp} \address{Université Marie et Louis Pasteur, SUPMICROTECH, CNRS, Institut FEMTO-ST (UMR 6174), F-25000 Besançon, France} \address{Université Marie et Louis Pasteur, CNRS, Institut UTINAM (UMR 6213), Observatoire des Sciences de l’Univers THETA, 41 bis avenue de l'Observatoire, F-25010 Besançon, France} \email{jana.elbadawi@femto-st.fr} \author{\firstname{Marion} \lastname{Delehaye}} \address[1]{Université Marie et Louis Pasteur, SUPMICROTECH, CNRS, Institut FEMTO-ST (UMR 6174), F-25000 Besançon, France} \email{marion.delehaye@femto-st.fr} \author{\firstname{Bruno} \lastname{Bellomo}} \address[2]{Université Marie et Louis Pasteur, CNRS, Institut UTINAM (UMR 6213), Observatoire des Sciences de l’Univers THETA, 41 bis avenue de l'Observatoire, F-25010 Besançon, France} \email{bruno.bellomo@univ-fcomte.fr} \thanks{This project has been supported by the EIPHI Graduate school (contract ANR-17-EURE-0002) and by the Bourgogne-Franche-Comté Region, by the ANR (project CONSULA, ANR-21-CE47- 0006-02), by First-TF Labex (contract ANR-10-LABX-48-01), and by Oscillator-IMP Equipex (contract ANR-11-EQPX-0033)} \CDRGrant[ANR]{ANR-17-EURE-0002} \CDRGrant[ANR]{ANR-21-CE47- 0006-02} \CDRGrant[ANR]{ANR-10-LABX-48-01} \CDRGrant[ANR]{ANR-11-EQPX-0033} \keywords{\kwd{Superradiance} \kwd{open quantum systems} \kwd{time and frequency metrology}} \altkeywords{\kwd{Superradiance} \kwd{systèmes quantiques ouverts} \kwd{métrologie temps-fréquence}} \editornote{Article submitted by invitation} \alteditornote{Article soumis sur invitation} \begin{abstract} We perform a theoretical study of a continuous superradiant laser supporting its experimental realization at FEMTO-ST using two sequentially-emitting ensembles of {}\textsuperscript{171}Yb atoms coupled to the same Fabry--Perot cavity. Using an open quantum system approach, we identify for the simplest case the parameter space where the laser reaches tens of picowatts of power with a sub-millihertz linewidth. Studying the impact of inhomogeneous frequency broadening and variations in atom-cavity coupling on the superradiant emission, we find the laser properties robust with respect to such perturbations, also thanks to the occurrence of synchronization of the atomic dipoles. We then consider a two-site configuration, in which atoms in each site are equally coupled to the cavity and have equal detunings, with different values for the two ensembles. We find for balanced and imbalanced atom numbers that synchronization leads in a certain parameter space to a single narrow spectral line whose central frequency follows the weighted average frequency. This result indicates that sequential loading can enable continuous superradiant emission for metrological applications, provided that the relative frequencies of the two ensembles are controlled to the level required by the target stability. \end{abstract} \begin{altabstract} Nous réalisons une étude théorique d’un laser superradiant continu en vue de son implantation expérimentale à FEMTO-ST utilisant deux ensembles d’atomes de {}\textsuperscript{171}Yb émettant séquentiellement dans une même cavité Fabry--Perot. En adoptant l'approche des systèmes quantiques ouverts, nous identifions dans le cas le plus simple l’espace des paramètres pour lequel le laser atteint des puissances de l’ordre de dizaines de picowatts avec une largeur de raie inférieure au millihertz. En étudiant l’impact de l’élargissement inhomogène de la fréquence et des variations du couplage atome-cavité sur l'émission superradiante, nous montrons que les propriétés du laser restent robustes face à de telles perturbations, notamment grâce à l’apparition d’une synchronisation des dipôles atomique. Nous considérons ensuite une configuration à deux sites, dans laquelle les atomes au sein de chaque site sont également couplés à la cavité et présentent le même désaccord en fréquence, avec des valeurs différentes pour les deux ensembles. Nous montrons pour des nombres d'atomes égaux et inégaux que la synchronisation peut conduire à une raie spectrale unique et étroite dont la fréquence centrale suit la moyenne pondérée des fréquences. Ce résultat indique que le chargement séquentiel peut permettre une émission superradiante continue pour des applications métrologiques, à condition que les fréquences relatives des deux ensembles soient contrôlées au niveau requis. \end{altabstract} \COI{The authors do not work for, advise, own shares in, or receive funds from any organization that could benefit from this article, and have declared no affiliations other than their research organizations.} \dateposted{2026-04-15} \begin{document} \maketitle \section{Introduction} The field of time and frequency metrology has seen remarkable advancements over recent decades. These developments have paved the way for a wide range of applications, from investigating variations in fundamental constants~\cite{Safronova_2019} to seismology~\cite{Tanaka2021}. Today's most precise timekeepers consist of the combination of a laser pre-stabilized on an ultra-stable, monolithic Fabry--Perot cavity, providing short-term stabilization down to $4 \times 10^{-17}$ at one second~\cite{Matei2017}, with an atomic ensemble that ensures long-term stabilization down to $4.8 \times 10^{-17} (\tau/\mathrm{s})^{-1/2}$~\cite{Oelker2019, Bothwell2022} and accuracy that can reach the $10^{-19}$ range~\cite{Aeppli2024, Hausser_2025}. Achieving even lower instabilities in optical clocks could unlock new applications, such as detecting gravitational waves~\cite{Kolkowitz2016} or searching for dark matter~\cite{Derevianko_2016, Filzinger2023}. However, these systems are currently limited by the fundamental thermal Brownian noise of the Fabry--Perot resonators. To overcome this limitation, several projects have emerged worldwide, including spectral hole burning~\cite{Thorpe2011}, Ramsey--Bordé interferometers~\cite{Olson2019}, and superradiant lasers~\cite{Meiser2009}. This article focuses on the project of a superradiant laser. The concept of a superradiant laser is based on superradiance, first introduced by Dicke in 1953~\cite{Dicke1953}, and further explored in the context of lasers in~\cite{Gross1982, Haake1993}. About 15 years ago, the idea of using such a laser as a new frequency reference emerged, with the potential to outperform current optical clocks by an order of magnitude~\cite{Meiser2009}. This idea has generated significant excitement, leading to several superradiant laser experiments aimed at metrological applications worldwide~\cite{Bohnet2012, Norcia2016, Norcia2016a, Norcia2017, Laske2019, Kristensen2023, Bohr_2024}. To date, all reported experiments operate in a pulsed manner, with the emission duration of one pulse below 1 second. Promising stabilities in the $10^{-16}$ range have been reported~\cite{Norcia2017}, and sub-natural linewidth emission has been observed~\cite{Kristensen2023}. However, the current pulsed nature of the emission remains a significant limitation for its use as a frequency reference. Intense theoretical efforts have been devoted to exploring possibilities for achieving continuous emission~\cite{steady_Holland, Kazakov2013, Hotter2022, Dubey2024}, which predict an appealing ultimate frequency instability at the $10^{-18}$ level at one second integration time~\cite{Kazakov2022}. One currently investigated possibility involves producing a beam of cold atoms loaded into a conveyor belt inside a bow-tie cavity where superradiance occurs~\cite{Dubey2024}. Other experiments focus on hot beams of atoms crossing a Fabry--Perot cavity to exploit the sub-natural linewidth emission achievable with such a system~\cite{Liu2020, Tang2022, Laburthe-Tolra2023, Oh2024}. Here, we theoretically investigate a system in which continuous superradiant emission could be achieved using two sequentially-emitting ensembles of $\Ytterbium$ atoms coupled to the same Fabry--Perot cavity, in order to realize predictions for the experiment under development at FEMTO-ST described in Section~\ref{sec:femto}. We use the framework of open quantum systems and the second-order cumulant expansion to describe the system, and we use the quantum regression theorem to obtain the emission spectrum in Section~\ref{sec:model}. In the first description of the superradiant laser system, a fixed number of identical $^{87}$Sr atoms was considered~\cite{Meiser2009}. This description captured the fundamental physics of the collective atomic dynamics and showed that it can produce a mHz linewidth laser with useful power. In Section~\ref{1class} we adopt a similar approach for our system of $^{171}$Yb atoms and show the characteristics of such a laser. In a recent study, this framework was extended to more realistic experimental conditions, in which the atoms have varying transition frequencies and are considered to be coupled to the cavity with varying strengths~\cite{Bychek_2021}. In Section~\ref{Xclass} we follow this approach using experimental parameters relevant to our system, and demonstrate the robustness of the emission with respect to such perturbations. Following this, we explore in Section~\ref{sec:2sites} the effects of a key feature of FEMTO-ST experiment, the loading of two sites labeled $A$ and $B$ inside the cavity, each loaded by a group of atoms with equal parameters. First, in Section~\ref{section:NA=NB}, we revisit the situation explored theoretically in~\cite{Minghui_NA} where the total number of atoms in each site is equal. In~\cite{Minghui_NA} this is explored by adiabatically eliminating the cavity field while here we keep the full description of the system. Second, in Section~\ref{section:NAdiffNB}, we consider the realistic experimental implementation for our system where the atoms are going to be transported to each site sequentially and consider that the total number of atoms in each site can vary. The characteristics of the laser is shown with the impact of atom number imbalance between the two sites. \section{Experimental concept}\label{sec:femto} At FEMTO-ST, we are establishing a superradiant experiment to explore the spectral capabilities of this new ultra-stable frequency reference. Our approach involves an original scheme to achieve sustained superradiant emission on the $\Ytterbium$ $\groundstate \rightarrow \excitedstate$ clock transition at 578 nm, which has a natural linewidth of 7 mHz. Our experimental apparatus consists of two chambers. In the first chamber, we prepare a magneto-optical trap of up to $10^8$ $\Ytterbium$ atoms, and we expect to transport up to $10^6$ atoms into one zone of a 56 mm-long Fabry--Perot cavity. This defines the practical operating range considered in this work. The cavity has a finesse of approximately $7000$ at the clock transition wavelength and the atoms are coupled to a single transverse cavity mode with an estimated mode diameter of about $145~\unit{\micro\meter}$. \begin{figure}[!h] \begin{subfigure}{0.31\textwidth} \begin{overpic}[width=\textwidth, trim={15cm 1cm 15cm 5cm}, clip]{rendered images/2025-07-24-MD-20-1-ss_fleche.png} \put(31,70){$A$} \put(62,70){$B$} \end{overpic} \caption{}\label{fig:femto_exp_a} \end{subfigure} \hfill \begin{subfigure}{0.31\textwidth} \begin{overpic}[width=\textwidth, trim={15cm 1cm 15cm 5cm}, clip]{rendered images/2025-07-24-MD-20-2-ss_fleche.png} \end{overpic} \caption{}\label{fig:femto_exp_b} \end{subfigure} \hfill \begin{subfigure}{0.31\textwidth} \begin{overpic}[width=\textwidth, trim={15cm 1cm 15cm 5cm}, clip]{rendered images/2025-07-24-MD-20-3-ss_fleche.png} \end{overpic} \caption{}\label{fig:femto_exp_c} \end{subfigure} \caption{Scheme of the experiment at FEMTO-ST. Black dots: atoms. Yellow lattice: cavity mode. Red lattices: trapping lattice within the cavity and optical conveyor belts. \eqref{fig:femto_exp_a}~Atoms are trapped in site $A$ of the cavity. \eqref{fig:femto_exp_b}~The number of atoms in site $A$ has decayed, hence a new ensemble of atoms is loaded into the site $B$ of the cavity by an optical conveyor belt. \eqref{fig:femto_exp_c}~Atoms are trapped in site $B$, and a new ensemble of atoms is prepared to be loaded in site $A$. }\label{fig:femto_exp} \end{figure} \begin{figure}[b] \includegraphics[width=0.3\textwidth]{all images/Yb_scheme.pdf} \caption{Proposed repumping scheme for the continuous superradiant laser based on $\Ytterbium$ atoms. The superradiant lasing occurs on the clock transition $\groundstate\leftrightarrow \excitedstate$. The intermediate states $\prescript{3}{}{\mathrm{P}_1}$ and $\prescript{3}{}{\mathrm{P}_2}$ participate in the repumping cycle which involves also higher auxiliary states, as shown schematically in the figure. }\label{fig:repumping-scheme} \end{figure} The experiment is based on two fixed zones in the cavity, referred to as ``site $A$'' and ``site $B$'', separated by approximately 3~mm. The two sites are sequentially populated using an optical conveyor belt in order to maintain the superradiant emission, as shown in Figure~\ref{fig:femto_exp}. In a typical experimental sequence, an atomic ensemble is first loaded into site $A$. This ensemble will be continuously repumped via a multi-step process which involves several auxiliary levels, as shown qualitatively in Figure~\ref{fig:repumping-scheme}, similarly to what is described in~\cite{Hotter2022}. This mechanism gives rise to an effective repumping rate used in the model described in Section~\ref{sec:model}. As the atom number decays over time due to residual heating and collisions with the background gas, we will prepare another ensemble to be loaded into site $B$ before the emission from site $A$ ends with a targeted refresh rate of $1$--$2$~Hz. In this paper, we investigate theoretically the spectral properties of the superradiant light emitted by such a system. $\Ytterbium$ has a nuclear spin $\mathrm{I}=1/2$, so both the ground state~$\groundstate$ and the excited clock state~$\excitedstate$ have two Zeeman sublevels. An external magnetic field will be applied orthogonally to the cavity axis, which will shift the sublevels symmetrically around the bare atomic frequency. The average of the emissions from the Zeeman sublevels allows for the reconstruction of a signal that is insensitive to the magnetic field~\cite{Norcia2017}. For typical magnetic fields of a few gauss, the resulting frequency shifts are on the order of a few hundreds Hz, which is larger than the targeted repumping rate, as found later in Section~\ref{1class}. In this paper, we consider a single pair of lasing states. \section{Description of the model}\label{sec:model} As model for the experimental setup described in the previous section and shown schematically in Figure~\ref{fig:femto_exp}, we consider an ensemble of $N$ cold $\Ytterbium$ atoms confined in an optical lattice and treated as two-level systems. Each atom in the ensemble is labeled by an index $i$, and is coupled to the cavity field with a position-dependent coupling strength $g_i$. The atomic transition frequency of the $i$-th atom is denoted by $\omega_i$ and the resonance frequency of the cavity by $\omega_c$. In particular, in the absence of individual atomic shifts, all the atoms are characterized by the unperturbed clock transition frequency $\omega_a$, and the cavity is assumed to be tuned in resonance with this transition. The unitary dynamics of the atoms and the cavity is then governed by the Hamiltonian (if not specified otherwise, we set $\hbar = 1$) \begin{equation} H=\sum_{i=1}^N\omega_i \sigma_i^{+} \sigma_i^{-}+ \omega_c a^{\dagger}a +\sum_{i=1}^N g_i\parens{a \sigma_i^{+}+a^{\dagger} \sigma_i^{-}}, \end{equation} where $\sigma_i^{+} = \ket{e}_i\bra{g}_i$ and $\sigma_i^{-} = \parens{\sigma_i^{+} }^{\dagger}$ are the raising and lowering operators for the $i $-th atom, $\ket{g}$ and $\ket{e}$ correspond to the atomic ground and excited states, $\groundstate$ and $\excitedstate$, respectively, and $a^{\dagger}$ and~$a$ are the photon creation and annihilation operators of the cavity mode. The atoms and cavity mode are subject to various dissipative processes. We consider a zone of the parameters in which the global dynamics for the atoms-cavity density matrix $\rho$ is governed by a master equation in the Lindblad form, \begin{equation} \frac{\dif}{\dif t}\rho=-i[H, \rho]+\mathcal{L}[\rho], \end{equation} where $\mathcal{L}[\rho]$ contains all the dissipative processes acting on the system via the operator $O$ at a rate $x$, described by the Lindblad superoperator $x \mathcal{D}[O] \rho= x\parens{2 O \rho O^{\dagger}-O^{\dagger} O \rho-\rho O^{\dagger} O}/2$. The dissipation includes the photons losses from the cavity at a rate $\kappa$, the spontaneous atomic decay at a single-atom rate $\gamma$, the individual incoherent repumping with a rate $R$, and inhomogeneous dephasing of each atom at a rate $1/T_2$. A cavity dephasing term denoted by the parameter $\xi$ which accounts for the effective noise induced by thermal fluctuations of the cavity mirrors is also included as in~\cite{Bychek_2021, Kazakov2022}. The single-atom rates $\gamma$, $R$, and $1/T_2$ are assumed to be identical for all atoms. These various dissipative processes contribute to the following total Lindblad dissipator, \begin{equation} \mathcal{L}[\rho] = \kappa\mathcal{D}[a] \rho + \xi \mathcal{D}[a^{\dagger} a] \rho + \sum_{i=1}^N \bracks[\Bigg]{\gamma \mathcal{D}\bracks{\sigma_i^{-}} \rho + R \mathcal{D}\bracks{\sigma_i^{+}} \rho + \frac{1}{2 T_2} \mathcal{D}\bracks{\sigma_i^{+} \sigma_i^{-}} \rho}. \end{equation} The zone of the parameters we work in justifies the use of a local approach to derive the total dissipator. In this approach, $\mathcal{L}[\rho]$ is obtained by the sum of the individual contributions from each dissipative source, as if it was the only source acting on the system~\cite{breuer2002theory, Rivas_2010, Bellomo2020}. Superradiant lasers operate in the so-called ``bad-cavity limit''~\cite{Meiser2009, BonifacioQST1}, in which the cavity decay rate $\kappa$ is much larger than the natural atomic decay rate $\gamma$, i.e., $\kappa\gg\gamma$. For our system the decay rate of the cavity is $\kappa = 2\pi \times 400~\unit{kHz} \approx 2.5 \times 10^6~\unit{rad.s^{-1}}$ which is many orders of magnitude more than the atomic decay rate $\gamma = 2\pi \times 7~\unit{mHz} \approx 4.4 \times 10 ^{-2}~\unit{rad.s^{-1}}$, corresponding to the ultra-narrow linewidth of the clock transition in the fermionic isotope of ytterbium $\Ytterbium$. The rate of the atom-cavity coupling for our system is $g = 2\pi \times 3.8~\unit{Hz} \approx 24~\unit{rad.s^{-1}}$, so we have $\kappa\gg g$ in our setup. The $1/T_2$ value of our experiment is not measured yet, so in the following we consider $1/T_2 \approx 1~\unit{rad.s^{-1}}$ to be a good approximation~\cite{Boyd_T2}. The parameters described above are fixed for our system, as they are determined by the experimental setup and physical characteristics of the components. A parameter that can be experimentally scanned over some range is that of the incoherent repumping rate $R$. This parameter describes the average rate at which atoms are repumped from the ground state $\groundstate$ to the excited state $\excitedstate$ of the clock transition. It provides the population inversion required for continuous superradiance. $R$ enters the model as simple rate parameter, however it is an effective description of a complex multi-level repumping scheme within the full atomic structure of $\Ytterbium$. The value of $R$ depends on the specific repumping configuration and can be calculated using the numerical method presented in~\cite{Hotter2022}, which models continuous multi-step repumping for the clock transition in bosonic $\StrontiumE$. In~\cite{Hotter2022}, it is stated that this method can extend well to the fermionic $\StrontiumS$ and to other alkaline-earth atoms, and that includes $\Ytterbium$ which is used in our system. As a result, the value of $R$ can have a tunable range based on the configuration of the multi-level system. The rate $\xi$ is associated in the model to a dephasing process for the cavity field, an additional dissipation term for the cavity other than the one governed by the photon loss rate $\kappa$. We follow the approach used in~\cite{Bychek_2021} and treat $\xi$ as a tunable parameter to check its impact on the characteristics of the laser. In all of this study, $\kappa$ stays the largest decay rate of the system. \subsection{Second-order cumulant expansion}\label{subsec:2ndcumulant} As the Hilbert space describing such a system grows exponentially with the atom number $N$, solving the Lindblad master equation becomes computationally complex for large $N$. In such cases, performing appropriate approximations is essential to simplify the problem. From the various approximation methods, the second-order cumulant expansion is commonly used~\cite{kubo_Ryogo}. As the $N$ atoms interact collectively with the cavity field, high-order atom correlations develop during the dynamics. Accurately considering such correlations is important for describing cooperative phenomena such as superradiant emission. However, accounting for all high-order correlations in a system of many atoms quickly becomes intractable, both computationally and analytically. The second-order cumulant expansion is a powerful approximation method that gives higher-order correlations in terms of lower-order ones. The key idea behind the second-order cumulant expansion is to truncate the correlations at second order by assuming that third- and higher-order cumulants are negligible. This assumption is valid in many physical systems where three-body, or higher, correlations are weak or can be captured by one- and two-body terms. The cumulants of arbitrary operators $X_1$, $X_2$, $X_3$ up to third order are defined as: \begin{equation} \begin{aligned} \braket{X_1}_c & = \braket{X_1}, \\ \braket{X_1 X_2}_c & = \braket{X_1 X_2} - \braket{X_1} \braket{X_2}, \\ \braket{X_1 X_2 X_3}_c & = \braket{X_1 X_2 X_3} - \braket{X_1 X_2} \braket{X_3} - \braket{X_1 X_3} \braket{X_2} - \braket{X_1} \braket{X_2 X_3} + 2 \braket{X_1} \braket{X_2} \braket{X_3}. \end{aligned} \end{equation} In this framework, the approximation $\braket{X_1 X_2 X_3}_c = 0$ allows us to express third-order correlations approximately as: \begin{equation} \braket{X_1 X_2 X_3}_c = 0 \quad \Longrightarrow \quad \braket{X_1 X_2 X_3} \approx \braket{X_1 X_2}\braket{X_3} + \braket{X_1 X_3}\braket{X_2} + \braket{X_1}\braket{X_2 X_3} - 2\braket{X_1}\braket{X_2}\braket{X_3}. \end{equation} This second-order cumulant expansion greatly reduces the computational cost while conserving the two-body quantum correlations. In the context of the linewidth of a superradiant laser, it was shown that this method, when compared with the full quantum solution, provides more accurate results than the mean-field and Langevin-based methods~\cite{PhD_SwadheenDubey}. The total phase invariance of the system is conserved under the second-order cumulant expansion and, as a result, the phase-dependent terms such as $\braket{a}$, $\braket{a^{\dagger}}$, and $\braket{\sigma^{\pm}}$ are zero~\cite{Meiser2009}. Starting from the master equation with the Hamiltonian of the system in the rotating frame of the cavity $H_{r} = -\sum_{i=1}^N\Delta_i\sigma_i^+\sigma_i^- +\sum_{i=1}^Ng_i(a\sigma_i^+ + a^{\dagger}\sigma_i^-)$, where $\Delta_i = \omega_c - \omega_i$ is the detuning between the $i$-th atom and the cavity, and applying the second-order cumulant expansion, we obtain the following closed set of equations that describe the system’s dynamics: \begin{equation}\label{eqs:Steady_state} \begin{aligned} \frac{\dif}{\dif t} \braket{a^{\dagger} a} & = -\kappa \braket{a^{\dagger} a} + \mathrm{i} \sum_{i=1}^N g_i \parens[\big]{\braket{a \sigma_i^{+}} - \braket{a^{\dagger} \sigma_i^{-}}}, \\ \frac{\dif}{\dif t} \braket{a \sigma_i^{+}} & = - \frac{1}{2}\parens*{\kappa + \gamma + R + \xi + \frac{1}{T_2} + \mathrm{i} 2\Delta_i } \braket{a \sigma_i^{+}} + \mathrm{i} g_i \braket{a^{\dagger} a} - 2\mathrm{i} g_i \braket{a^{\dagger} a} \braket{\sigma_i^{+} \sigma_i^{-}} - \mathrm{i} g_i \braket{\sigma_i^{+} \sigma_i^{-}} \\ & \quad - \mathrm{i} \sum_{j,j\neq i}^N g_j \braket{\sigma_{i}^{+} \sigma_{j}^{-}}, \\ \frac{\dif}{\dif t} \braket{\sigma_i^{+} \sigma_i^{-}} & = \mathrm{i} g_i \braket{a^{\dagger} \sigma_i^{-}} - \mathrm{i} g_i \braket{a \sigma_i^{+}} - (\gamma + R) \braket{\sigma_i^{+} \sigma_i^{-}} + R, \\ \frac{\dif}{\dif t} \braket{\sigma_{i}^{+} \sigma_{j}^{-}} & = -\mathrm{i}(\Delta_i - \Delta_j)\braket{\sigma_{i}^{+} \sigma_{j}^{-}} + \mathrm{i} g_i \braket{a^{\dagger} \sigma_j^{-}} - \mathrm{i} g_j \braket{a \sigma_i^{+}} - 2\mathrm{i} g_i \braket{a^{\dagger} \sigma_j^{-}} \braket{\sigma_i^{+} \sigma_i^{-}} \\ & \quad + 2\mathrm{i} g_j \braket{a \sigma_i^{+}} \braket{\sigma_j^{+} \sigma_j^{-}}- \parens[\Bigg]{\gamma + R + \frac{1}{T_2}} \braket{\sigma_{i}^{+} \sigma_{j}^{-}}. \end{aligned} \end{equation} In the following sections, we handle specific cases of this system and we only focus on its steady-state solutions. First in Section~\ref{1class}, we consider identical atoms equally coupled to the cavity mode, and this reduces Eqs.~\eqref{eqs:Steady_state} to a set of four equations. Note that in this case there is symmetry in the terms with respect to exchange of atoms. For instance, $\braket{a^{\dagger} \sigma_j^{-}} =\braket{a^{\dagger} \sigma_1^{-}}$ for all $j$, and $\braket{\sigma_i^{+}\sigma_j^-}=\braket{\sigma_1^{+} \sigma_2^{-}}$ for all $i \neq j$. The analytical steady-state solutions for this case are presented in Appendix~\ref{app:ss}, with further simplifications in the large $\kappa$ and large $R$ limit in Appendix~\ref{app:simple}, and large $N$ limit in Appendix~\ref{app:largeN}. Next in Section~\ref{Xclass}, we consider the case of non-identical atoms, where detunings can differ, $\Delta_i \neq \Delta_j$. To manage with the increasing number of equations, we adopt an approach commonly used in the literature~\cite{Bychek_2021}, where the atoms are grouped into $M$ frequency classes. Each class contains atoms that share the same detuned frequency, and the number of atoms are distributed according to a Gaussian frequency distribution. Similarly, in the case of atoms with different coupling strengths to the cavity mode, the atoms can be grouped in $K$ coupling strength classes. In the presence of both inhomogeneous broadening and varying coupling strengths, the system is described by $M \times K$ classes. The atoms are divided equally among $K$ coupling classes, and within each coupling class they are distributed into $M$ frequency classes. Last, in Section~\ref{sec:2sites} we consider the experimentally relevant scenario with two groups of atoms with their own detuning $\Delta_A$ and $\Delta_B$, and consider the case of equal number of atoms and the case of an imbalance in the number of atoms. \subsection{Spectrum of the cavity field}\label{sec:spectrum} Having formulated the dynamics within the second-order cumulant framework, we now focus on the characteristics of the spectrum of the cavity field. To calculate this spectrum, we make use of the Wiener--Khinchin theorem, according to which the spectrum $S(\omega)$ is given by the Fourier transform of the first-order correlation function $g^{(1)}(\tau) = \braket[\big]{a^{\dagger}(\tau+t) a(t)}$~\cite{Puri2001}, \begin{equation} S(\omega)=2 \, \mathbb{R} \braces*{ \int_0^{\infty} \dif \tau e^{-\mathrm{i} \omega \tau} g^{(1)}(\tau) }. \end{equation} To compute this correlation function, we apply the quantum regression theorem, which states that the time derivative of $\braket[\big]{a^{\dagger}(\tau+t) a(t)}$ follows the time derivative of $\braket[\big]{a^{\dagger}(t)}$, and we obtain using the second-order cumulant expansion: \begin{equation}\label{eq:time_evo} \begin{aligned} \frac{\dif}{\dif \tau} \braket[\big]{a^{\dagger}(\tau+t) a(t)} & = -\frac{\kappa+\xi}{2} \braket[\big]{a^{\dagger}(\tau+t) a(t)} + \mathrm{i} \sum_{i=1}^N g_i \braket[\big]{\sigma_i^{+}(\tau +t) a(t)}, \\ \frac{\dif}{\dif \tau } \braket[\big]{\sigma_i^{+}(\tau +t) a(t)} & = -\mathrm{i} g_i \braket[\big]{\sigma_i^z(\tau + t)} \braket[\big]{a^{\dagger}(\tau+t) a(t)} - \frac{1}{2}\parens[\Bigg]{\gamma+R+ \frac{1}{T_2} + 2 \mathrm{i}\Delta_i } \braket[\big]{\sigma_i^{+}(\tau +t) a(t)}. \end{aligned} \end{equation} To solve the above equations, the starting point $t=0$ is taken when the system is in steady-state. The term $\braket[\big]{\sigma_i^z(\tau+t)}$ is then considered a constant equal to the steady-state value $\braket{\sigma_i^z}^{\steady}$, which can be obtained from Eqs.~\eqref{eqs:Steady_state}. Analytical expressions for the spectrum are provided in Appendix~\ref{app:spectrum} for the one-class case. The cavity spectrum provides access to key observables of the system, such as the spectral linewidth and the power and the output of the cavity. The spectral linewidth corresponds to the full-width at half-maximum (FWHM) of the spectrum. The output power, on the other hand, is determined by the mean intracavity photon number $\braket{a^{\dagger}a}$ and is given by $P = \eta\hbar\omega_0\kappa\braket{a^{\dagger}a}$, where $\eta$ is the relative transmission of the output mirror of the cavity and $\omega_0$ is the central angular frequency of the spectrum~\cite{Kazakov2022}. In what follows, all power calculations are performed with $\eta=1$, that corresponds to ideal transmission. When all the atoms are on resonance with the cavity field, the spectrum is expected to be centered at $\omega_0=\omega_c$. \section{One-class case}\label{1class} In this section we consider the case of identical atoms ($\Delta_i = \Delta_j=\Delta)$ equally coupled to the cavity mode ($g_i = g_j = g$). In particular, we focus on the case when the atoms are on resonance with the cavity mode ($\Delta = 0)$ since the non-resonant case gives similar results if $\Delta$ is small enough. This follows directly from Eqs.~\eqref{Steady state relations}, which imply that for $\Delta\ll \parens[\big]{\gamma +\kappa + (1/T_2) +\xi +R}/2$, all the variables of Eqs.~\eqref{eqs:Steady_state} do not depend appreciably on $\Delta$. We also notice that in this regime $\mathbb{R}\bracks*{\braket{a\sigma^{+}}^\steady}$, which is zero for $\Delta=0$, is much smaller than $\mathbb{I}\bracks*{\braket{a\sigma^{+}}^\steady}$. The parameters are set equal to the experimental values provided in Section~\ref{sec:model}, $\gamma = 2\pi \times 7~\unit{mHz} $, $\kappa = 2\pi \times 400~\unit{kHz} $, $g = 2\pi\times 3.8~\unit{Hz}$, and $1/T_2 = 1~\unit{rad.s^{-1}}$, with $\xi = 0$ except when specified otherwise. In what follows, we study the steady-state solutions of the system described by Eqs.~\eqref{eqs:Steady_state}. The solutions are found numerically in the main text, though for this section they can be computed using the analytical derivations in Appendix~\ref{app:ss}. We also observe that according to the derivation of the spectrum for $\Delta=\xi=0$ in Eq.~\eqref{app:Spectrum special case}, when $\braket{\sigma^z}^\steady>0$ the spectrum has only one peak and is symmetric with respect to the cavity frequency. For a laser operating in the bad-cavity limit, the superradiant emission in the system is characterized by both a lower and upper threshold with respect to $R$~\cite{Meiser2009}. The dependence of the various variables with respect to $R$ is presented in Figure~\ref{2x2obs} for $N = 10^6$ atoms. The lower threshold corresponds to the point at which the repumping rate overcomes the atomic decay rate $R_{\min } \ge \gamma$, which is necessary to generate the population inversion. From Eqs.~\eqref{Steady state relations} it follows that $\braket{\sigma^z_1}^{\steady} < 0$ for $R<\gamma$, as it can be observed in Figure~\ref{2x2obs}(a). For such values of $R$ where the population inversion is not reached, the intracavity photon number is observed in Figure~\ref{2x2obs}(b) to be very small. Above this threshold $\braket{\sigma^z_1}^{\steady}$ crosses zero and $\braket{a^\dagger a}^\steady$ shows a sharp increase, indicating the onset of the superradiant emission. The atom-atom correlations $\braket{\sigma^+_1\sigma^-_2}^{\steady}$ in Figure~\ref{2x2obs}(c) and the imaginary part of the atom-field correlations $\braket{a\sigma^+_1}^{\steady}$ (real part is equal to zero for $\Delta = 0$) in Figure~\ref{2x2obs}(d), are close to zero for such small values of $R$. \begin{figure}[t] \includegraphics[width=0.7\textwidth]{all images/atomic_N1e6_colormix_2x2.pdf} \caption{Steady-state solutions of Eqs.~\eqref{eqs:Steady_state} as a function of $R$ for $N=10^6$ atoms represented by solid lines. The dotted lines refer to the large $\kappa$ and $R$ limit described in Appendix~\ref{app:simple}, while dashed lines correspond to the large $N$ limit reported in Appendix~\ref{app:largeN}. Red (resp.\ black) vertical dotted line: $R_{\min}$ (resp.\ $R_{\max}$). }\label{2x2obs} \end{figure} As the value of $R$ increases, all the above variables increase in modulus. In particular, the photon number increases linearly with $R$, as found in Eqs.~\eqref{ada large N} in the large $N$ limit for $R>\gamma$. Based on the approximate Eqs.~\eqref{Steady state relations simple large N} in the large $\kappa$ and $R$ limit, we observe that the atom-atom correlations and the photon number both reach a maximum at $R= 2g^2N/\kappa$, with $\braket{\sigma_i^+ \sigma_j^-}^{\steady} \approx 1/8$ and $\braket{a^{\dagger}a}^{\steady} \approx g^2N^2/(2\kappa^2)$. This value of $R$ is then defined as the optimal rate, $R_{\mathrm{opt}} \equiv 2g^2N/\kappa$, where superradiance manifests the most. At $R_\mathrm{opt}$, it holds $\braket{\sigma^z_1}^{\steady} \approx 1/2$, and the atom-field correlations reach a minimum, $\braket{a \sigma^+_1}^{\steady} \approx -gN/(4\kappa)$. For values of $R$ above this optimal rate, $\braket{\sigma^z}^\steady$ approaches 1 for $R =2 R_{\mathrm{opt}}= 4g^2N/\kappa$, and maintains this value for larger values of $R$, as shown in Eqs.~\eqref{Steady state relations small N} in the large $\kappa$ and $R$ limit. At $2 R_{\mathrm{opt}}$, the photon number, the atom-atom correlations, and the atom-field correlations are all close to zero. This value of $R$ is then the upper threshold of the superradiant emission and it is defined as $R_{\max} \equiv 2 R_{\mathrm{opt}}=4g^2N/\kappa = NC\gamma$, where the cooperativity parameter $C$ is defined as $C = 4g^2/ (\kappa\gamma)$~\cite{Kazakov2022}. The term $C\gamma$ is the single-atom decay rate defined in the atomic master equation following the adiabatic elimination of the cavity mode in the bad-cavity limit~\cite{intensity_Holland}. It follows that $NC\gamma$ characterizes the $N$-atom dynamics. The condition $R \ge R_{\max }$ then corresponds to the case when the repumping rate is equal or larger than this collective decay rate. In Figure~\ref{2x2obs} we also compare the solutions of Eqs.~\eqref{eqs:Steady_state} with their approximate expressions for the large $\kappa$ and $R$ limit and the large $N$ limit. The large $\kappa$ and $R$ limit is accurate for the atom-atom and the atom-field correlations for all values of $R$. This is not the case for the atomic inversion and for the photon number. For the atomic inversion, it holds for $R\ge R_{\min}$, and for the photon number for $R^* \le R\le R_{\max}$, where $R^*$ is a certain value above $R_{\min}$. The large $N$ approximation is well satisfied for all quantities up to some value of $R$, depending on the value of $N$ we are considering. For the photon number at $R\le R_{\min}$, the solution becomes unstable for values close to $R_{\min}$, and for this reason the corresponding region is not shown in the figure. \begin{figure}[h!] \includegraphics[width=0.7\textwidth]{all images/ata_analytical.pdf} \caption{Number of intracavity photons at the steady-state $\braket{a^{\dagger}a}^{\steady}$ as a function of the repumping rate $R$ for various values of $N$. The solid lines are the steady-state solutions of Eqs.~\eqref{eqs:Steady_state}, while the dotted lines refer to the large $\kappa$ and $R$ limit reported in Appendix~\ref{app:simple} and the dashed lines correspond to the large $N$ limit reported in Appendix~\ref{app:largeN}.}\label{fig:steady variables R} \end{figure} In order to observe the behavior of the photon number with the change in the atom number, we compare in Figure~\ref{fig:steady variables R} $\braket{a^{\dagger}a}^{\steady}$ for several values of $N$, with its approximate expressions for the large $\kappa$ and $R$ limit, and the large $N$ limit. As observed above, the two limits can be well satisfied for certain values of $R$. This is the case for $N=10^5$ and $10^6$, while the agreement degrades more and more for smaller values of $N$. We also notice that the maximum photon number when we move from $N=10^5$ to $10^6$, is increased by around two orders of magnitude, as expected from its quadratic scaling with $N$. The study of the steady-state solutions of the atomic and cavity variables pointed out that around $R_{\mathrm{opt}}$ the system could be exploited as a superradiant laser~\cite{Meiser2009}. In Figure~\ref{fig:power_and_linewidth_1class}, we analyze the steady-state laser characteristics for a large range of $R$ and $N$. Figure~\ref{fig:power_and_linewidth_1class}(a) shows the photon number and the output power, while Figure~\ref{fig:power_and_linewidth_1class}(b) displays the linewidth of the cavity field spectrum, all as functions of $R$ and $N$. To show the importance of the common cavity mode in the generation of atom-atom correlations, we compare the collective configuration with an alternative one consisting of $N$ independent single-atom systems, each coupled to its own cavity with parameters identical to that of the $N$-atom case. The intracavity photons from all single-atom systems are summed, and the corresponding total output power is given in Figure~\ref{fig:power_and_linewidth_1class}(c). In this configuration, the linewidth is equal for all single-atom systems and it corresponds to the linewidth of the total system, which is given in Figure~\ref{fig:power_and_linewidth_1class}(d). In the presence of strong collective effects in the $N$-atom system, significant differences between the two configurations are expected. \begin{figure}[t] \includegraphics[width=0.75\textwidth]{all images/NR_withcontours_bottop.pdf} \caption{(a),~(c)~Intracavity photon number $\braket{a^\dagger a}^{\steady}$ and output power $P$, and (b),~(d)~linewidth $\Delta\nu$ both as a function of atom number $N$ and repumping rate $R$. (a)--(b)~Collective dynamics with $N$ atoms coupled to a single cavity mode. Red dashed line: minimum repumping rate $R_{\min} = \gamma$ for superradiant emission onset; red dot-dashed line: maximum repumping rate $R_{\max} = N C \gamma$ above which superradiance is suppressed; light blue dotted line: critical atom number $N_{\crit} = 2/(C \gamma T_2)$ for the onset of collective behavior~\cite{Meiser2009}. (c)--(d) Non-collective dynamics of $N$ single-atom systems each coupled to an individual cavity with identical parameters as in (a)--(b). Various contour lines (linking points with the same value) are depicted by white dot-dashed lines.}\label{fig:power_and_linewidth_1class} \end{figure} In the collective dynamics with $RR_{\min}$, this power shows negligible variation with an increase in $R$. The linewidth for this system is broad at small $R$ and it increases linearly with $R$. This behavior is in contrast to that observed in the collective dynamics for large~$N$, where the power scales with $N^2$ and the linewidth becomes very narrow. This clearly shows the importance of collective coupling, enabled by the common cavity mode, in order to have a superradiant emission process that enhances the laser performance. With respect to the atom number impact on the dynamics, it was found in~\cite{Meiser2009} that in the limit of $T_2 \gamma \ll 1$, there is a critical number of atoms $N_{\crit}=2 /\parens{C \gamma T_2} \approx 2.2 \times 10^3$ required to reach the collective dynamics of the superradiant laser. This critical number is marked in Figure~\ref{fig:power_and_linewidth_1class}(a)--(b) by a light blue line, and for $N\gamma$. Here we focus on the case $R\gg \{\gamma, T_2^{-1}, \xi \}$ and note that $\Delta$ is also expected to be much smaller than $\kappa$. Since the number of atoms is expected to be very large, we can use the approximations $N-1 \approx N$ and $N-2 \approx N$. With such conditions, the equation governing the steady-state value $\braket{a^{\dagger} a}^{\steady}$ is given by the positive solution of a second order equation, $\mathrm{\tilde{a}} \braket{a^{\dagger} a}^{\steady~2}+\mathrm{\tilde{b}}\braket{a^{\dagger} a}^{\steady} +\mathrm{\tilde{c}}=0 $, where (using $N-2\approx N$) \begin{equation} \mathrm{\tilde{a}} = \frac{2 g \kappa (\kappa +R)}{N R^2}, \qquad \mathrm{\tilde{b}}=(\kappa +R)\parens*{\frac{\kappa}{4 g N}-\frac{ g}{R}}, \qquad \mathrm{\tilde{c}}=-g. \end{equation} In this limit we then get the following steady-state values of the variables, considering the condition $N\ge R \kappa/(4 g^2)$ in which we first obtain $\braket{a^{\dagger} a}^{\steady}$ and using $N\gg8$: \begin{equation}\label{Steady state relations simple large N} \begin{aligned} \braket{a^{\dagger} a}^{\steady} & \approx \frac{N R}{2 \kappa}-\frac{R^2}{8 g^2}, \\ \braket{\sigma^z}^{\steady} & \approx 1 -\frac{2 \braket{a^{\dagger} a}^{\steady} \kappa }{R N}\approx\frac{R \kappa}{4 g^ 2 N}, \\ \mathbb{R}\bracks[\big]{\braket{a\sigma^{+}}^{\steady}}& \approx -\frac{\braket{a^{\dagger} a}^{\steady} \Delta \kappa }{g N (\kappa +R)}\approx \frac{\Delta R \parens{\kappa R-4 g^2 N}}{8 g^3 N (\kappa +R)}, \\ \mathbb{I}\bracks[\big]{\braket{a\sigma^{+}}^{\steady}} & = -\frac{\kappa\braket{a^{\dagger} a}^{\steady} }{2 g N}\approx \frac{R \parens{\kappa R-4 g^2 N}}{16 g^3 N}, \\ \braket{\sigma_{i}^{+} \sigma_{j}^{-}}^{\steady}& \approx\frac{\braket{a^{\dagger} a}^{\steady} }{N R^2}\bracks*{\kappa R-\frac{2 \braket{a^{\dagger} a}^{\steady} \kappa ^2}{N}} \approx \frac{\kappa R \parens{4 g^2 N-\kappa R}}{32 g^4 N^2}. \end{aligned} \end{equation} Several considerations can be derived from the above equations. For the atom-atom correlations $\braket{\sigma_{i}^{+} \sigma_{j}^{-}}^{\steady}$, one finds that they reach a maximum at $N=\kappa R/(2 g^2)$ (equivalently, for fixed $N$, at $R=2 g^2 N/\kappa$), where the value of the maximum is $1/8$. At this particular value of $N$, we have $\braket{a^{\dagger} a}^{\steady}\approx g^2 N^2/(2 \kappa^2)$ (this is a maximum with respect to $R$ at $R=2 g^2 N/\kappa$), $\braket{\sigma^z}^{\steady}\approx1/2$, $\mathbb{R}\bracks[\big]{\braket{a\sigma^{+}}^{\steady}}\approx-\frac{\Delta g N}{2 (\kappa^2+2 g^2 N)}$, and $\mathbb{I}\bracks[\big]{\braket{a\sigma^{+}}^{\steady}}\approx-g N/(4 \kappa)$ (this is a minimum with respect to $R$ at $R=2 g^2 N/\kappa$). As for the threshold $R_{\min}=\gamma$, Eq.~\eqref{Steady state relations} shows that for $R<\gamma$ we have $\braket{\sigma_{i}^{+} \sigma_{j}^{-}}^{\steady}<0$. In the above equations for the limit considered here the steady-state correlations are positive, $\braket{\sigma_{i}^{+} \sigma_{j}^{-}}^{\steady}>0$. Using the equations derived for the large $N$ limit in the next subsection (Eqs.~\eqref{SpmijN}), it follows that for a large enough value of $N$ one has $\braket{\sigma_{i}^{+} \sigma_{j}^{-}}^{\steady}>0$ for $R>\gamma$ provided that $\kappa R>4 g^2$ (provided the system is in the same limit considered in this subsection). When $N \leq R\kappa /(4 g^2)$ and at the same time $N \gg 8$, we have \begin{equation}\label{Steady state relations small N} \braket{a^{\dagger} a}^{\steady} \approx 0, \qquad \braket{\sigma^z}^{\steady} \approx 1, \qquad \mathbb{R}\bracks[\big]{\braket{a\sigma^{+}}^{\steady}} \approx 0, \qquad \mathbb{I}\bracks[\big]{\braket{a\sigma^{+}}^{\steady}} \approx 0, \qquad \braket{\sigma_{i}^{+} \sigma_{j}^{-}}^{\steady}\approx 0. \end{equation} \begingroup \allowdisplaybreaks \section{Simplified expressions: large \texorpdfstring{$N$}{N} limit}\label{app:largeN} The steady solutions in Eqs.~\eqref{Steady state relations} can be used to compute the expression of $\braket{a^{\dagger} a}^{\steady}$ and the other variables in the limit of a very large number of atoms $N$: %\begin{align} %\braket{a^{\dagger} a}^{\steady}_{N\rightarrow \infty} & = %\begin{dcases}\label{ada large N} % \begin{aligned}[b] % & \textstyle \kappa R (\gamma +T_2^{-1} +R) % \\ & \textstyle \quad \times \frac{\bracks[\Big]{4 T_2^{-1} (\gamma +R)+16 \gamma R -\frac{(R-\gamma) (\gamma +R) (\gamma +T_2^{-1} +R) \parens*{4 \Delta^2+(\gamma +\kappa +T_2^{-1} +\xi +R)^2}}{g^2 (\gamma +\kappa +T_2^{-1} +\xi +R)}}}{4 N (R-\gamma)^3 (\gamma +\kappa +T_2^{-1} +R)^2} % \\ & \textstyle \quad + \frac{R (\gamma +T_2^{-1} +R)}{(\gamma -R) (\gamma +\kappa +T_2^{-1} +R)}, % \end{aligned} %& \text{for $R < \gamma$}, %\\ \textstyle \sqrt{\frac{ N \gamma (2 \gamma +T_2^{-1})}{2 \kappa (2 \gamma +\kappa +T_2^{-1})}}, %& \text{for $R = \gamma$}, %\\ \textstyle \frac{N (R-\gamma)}{2 \kappa }, %& \text{for $R > \gamma$}, %\end{dcases} %\\ %\braket{\sigma^z}^\steady_{N\rightarrow \infty} & = %\begin{dcases}\label{sigmaz large N} % \textstyle \frac{R-\gamma }{\gamma +R} + \frac{2 \kappa R (\gamma +T_2^{-1} +R)}{N (R-\gamma) (\gamma +R) (\gamma +\kappa +T_2^{-1} +R)}, %& \text{for $R < \gamma$}, %\\ \textstyle \sqrt{\frac{\kappa (2 \gamma +T_2^{-1})}{2 \gamma N (2 \gamma +\kappa +T_2^{-1})}}, %& \text{for $R = \gamma$}, %\\ \begin{aligned}[b] % & \textstyle \kappa (\gamma +T_2^{-1} +R) % \\ & \textstyle \quad \times \frac{\bracks[\Big]{(\gamma +\kappa +T_2^{-1} +\xi +R)^2+4 \Delta^2 -\frac{4 g^2 (\gamma +\kappa +T_2^{-1} +\xi +R)}{R-\gamma }}}{4 g^2 N (\gamma +\kappa +T_2^{-1} +R) (\gamma +\kappa +T_2^{-1} +\xi +R)}, % \end{aligned} %& \text{for $R > \gamma$}, %\end{dcases} %\\ %\mathbb{R}\bracks[\big]{\braket{a\sigma^{+}}^{\steady}_{ N\rightarrow \infty}} & = %\begin{dcases} % \textstyle \frac{\Delta \kappa R (\gamma +T_2^{-1} +R)}{g N (R-\gamma) (\gamma +\kappa +T_2^{-1} +R) (\gamma +\kappa +T_2^{-1} +\xi +R)}, %& \text{for $R < \gamma$}, %\\ \textstyle -\frac{(\Delta \kappa) }{g (2 \gamma +\kappa +T_2^{-1} +\xi)} \sqrt{\frac{\gamma (2 \gamma +T_2^{-1})}{2 \kappa N (2 \gamma +\kappa +T_2^{-1})}}, %& \text{for $R = \gamma$}, %\\ \begin{aligned}[b] % & \textstyle \Delta \kappa (\gamma +R) (\gamma +T_2^{-1} +R) % \\ & \textstyle \quad \times \frac{(R-\gamma) \parens*{4 \Delta^2+(\gamma +\kappa +T_2^{-1} +\xi +R)^2}-4 g^2 (\gamma +\kappa +T_2^{-1} +\xi +R)}{8 g^3 N (R-\gamma) (\gamma +\kappa +T_2^{-1} +R) (\gamma +\kappa +T_2^{-1} +\xi +R)^2} % \\ & \textstyle \quad + \frac{\Delta (R-\gamma)}{2 g (\gamma +\kappa +T_2^{-1} +\xi +R)}, % \end{aligned} %& \text{for $R > \gamma$}, %\end{dcases} %\\ %\mathbb{I}\bracks*{\braket{a\sigma^{+}}^\steady_{ N\rightarrow \infty} } & = %\begin{dcases}\label{SpmijN} % \textstyle \frac{\kappa R (\gamma +T_2^{-1} +R)}{2 g N (R-\gamma) (\gamma +\kappa +T_2^{-1} +R)}, %& \text{for $R < \gamma$}, %\\ \textstyle - \frac{1}{2 g}\sqrt{\frac{\gamma \kappa (2 \gamma +T_2^{-1})}{2 N (2 \gamma +\kappa +T_2^{-1})}}, %& \text{for $R = \gamma$}, %\\ \begin{aligned}[b] % & \textstyle - \kappa (\gamma +R) (\gamma +T_2^{-1} +R) % \\ & \textstyle \quad \times \frac{\frac{4 g^2 (\gamma +\kappa +T_2^{-1} +\xi +R)}{R-\gamma }-4 \Delta^2-(\gamma +\kappa +T_2^{-1} +\xi +R)^2}{16 g^3 N (\gamma +\kappa +T_2^{-1} +R) (\gamma +\kappa +T_2^{-1} +\xi +R)}-\frac{R-\gamma }{4 g}, % \end{aligned} %& \text{for $R > \gamma$}, %\end{dcases} %\\ %\braket{\sigma_{i}^{+} \sigma_{j}^{-}}^\steady_{ N\rightarrow \infty} & = %\begin{dcases} % \textstyle - \frac{\kappa R}{N (\gamma +R) (\gamma +\kappa +T_2^{-1} +R)}, %& \text{for $R < \gamma$}, %\\ \textstyle - \frac{\kappa }{2 N (2 \gamma +\kappa +T_2^{-1})}, %& \text{for $R = \gamma$}, %\\ \begin{aligned}[b] % & \textstyle \frac{\kappa (R-\gamma) \parens*{4 \Delta^2+(\gamma +\kappa +T_2^{-1} +\xi +R)^2}}{8 g^2 N (\gamma +\kappa +T_2^{-1} +R) (\gamma +\kappa +T_2^{-1} +\xi +R)} % \\ & \textstyle \quad -\frac{\kappa (4 g^2) (\gamma +\kappa +T_2^{-1} +\xi +R)}{8 g^2 N (\gamma +\kappa +T_2^{-1} +R) (\gamma +\kappa +T_2^{-1} +\xi +R)}, % \end{aligned} %& \text{for $R > \gamma$}. %\end{dcases} %\end{align} \begin{alignat}{2} \braket{a^{\dagger} a}^{\steady}_{N\rightarrow \infty} & = \begin{dcases}\label{ada large N} \begin{aligned}[b] & \textstyle \kappa R (\gamma +T_2^{-1} +R) \\ & \textstyle \quad \times \frac{\bracks[\Big]{4 T_2^{-1} (\gamma +R)+16 \gamma R -\frac{(R-\gamma) (\gamma +R) (\gamma +T_2^{-1} +R) \parens*{4 \Delta^2+(\gamma +\kappa +T_2^{-1} +\xi +R)^2}}{g^2 (\gamma +\kappa +T_2^{-1} +\xi +R)}}}{4 N (R-\gamma)^3 (\gamma +\kappa +T_2^{-1} +R)^2} \\ & \textstyle \quad + \frac{R (\gamma +T_2^{-1} +R)}{(\gamma -R) (\gamma +\kappa +T_2^{-1} +R)}, \end{aligned} \\ \textstyle \sqrt{\frac{ N \gamma (2 \gamma +T_2^{-1})}{2 \kappa (2 \gamma +\kappa +T_2^{-1})}}, \\ \textstyle \frac{N (R-\gamma)}{2 \kappa }, \end{dcases} & \quad \begin{aligned} & \vphantom{\textstyle T_2^{-1}} \\ & \vphantom{\textstyle \frac{\bracks[\Big]{\frac{ \parens*{T_2^{-1}}}{T_2^{-1}}}}{T_2^{-1}}} \\ & \vphantom{\textstyle \frac{T_2^{-1}}{T_2^{-1}}} \text{for $R < \gamma$}, \\ & \vphantom{\textstyle \sqrt{\frac{T_2^{-1}}{T_2^{-1}}}} \text{for $R = \gamma$}, \\ & \vphantom{\textstyle \frac{N (R-\gamma)}{2 \kappa }} \text{for $R > \gamma$}, \end{aligned} \\ \braket{\sigma^z}^\steady_{N\rightarrow \infty} & = \begin{dcases}\label{sigmaz large N} \textstyle \frac{R-\gamma }{\gamma +R} + \frac{2 \kappa R (\gamma +T_2^{-1} +R)}{N (R-\gamma) (\gamma +R) (\gamma +\kappa +T_2^{-1} +R)}, \\ \textstyle \sqrt{\frac{\kappa (2 \gamma +T_2^{-1})}{2 \gamma N (2 \gamma +\kappa +T_2^{-1})}}, \\ \begin{aligned}[b] & \textstyle \kappa (\gamma +T_2^{-1} +R) \\ & \textstyle \quad \times \frac{\bracks[\Big]{(\gamma +\kappa +T_2^{-1} +\xi +R)^2+4 \Delta^2 -\frac{4 g^2 (\gamma +\kappa +T_2^{-1} +\xi +R)}{R-\gamma }}}{4 g^2 N (\gamma +\kappa +T_2^{-1} +R) (\gamma +\kappa +T_2^{-1} +\xi +R)}, \end{aligned} \end{dcases} & \quad \begin{aligned} & \vphantom{\textstyle \frac{R-\gamma }{\gamma +R} + \frac{2 \kappa R (\gamma +T_2^{-1} +R)}{N (R-\gamma) (\gamma +R) (\gamma +\kappa +T_2^{-1} +R)},} \text{for $R < \gamma$}, \\ & \vphantom{\textstyle \sqrt{\frac{\kappa (2 \gamma +T_2^{-1})}{2 \gamma N (2 \gamma +\kappa +T_2^{-1})}},} \text{for $R = \gamma$}, \\ & \vphantom{\textstyle \kappa (\gamma +T_2^{-1} +R)} \\ & \vphantom{\textstyle \quad \times \frac{\bracks[\Big]{(\gamma +\kappa +T_2^{-1} +\xi +R)^2+4 \Delta^2 -\frac{4 g^2 (\gamma +\kappa +T_2^{-1} +\xi +R)}{R-\gamma }}}{4 g^2 N (\gamma +\kappa +T_2^{-1} +R) (\gamma +\kappa +T_2^{-1} +\xi +R)},} \text{for $R > \gamma$}, \end{aligned} \\ \mathbb{R}\bracks[\big]{\braket{a\sigma^{+}}^{\steady}_{ N\rightarrow \infty}} & = \begin{dcases} \textstyle \frac{\Delta \kappa R (\gamma +T_2^{-1} +R)}{g N (R-\gamma) (\gamma +\kappa +T_2^{-1} +R) (\gamma +\kappa +T_2^{-1} +\xi +R)}, \\ \textstyle -\frac{(\Delta \kappa) }{g (2 \gamma +\kappa +T_2^{-1} +\xi)} \sqrt{\frac{\gamma (2 \gamma +T_2^{-1})}{2 \kappa N (2 \gamma +\kappa +T_2^{-1})}}, \\ \begin{aligned}[b] & \textstyle \Delta \kappa (\gamma +R) (\gamma +T_2^{-1} +R) \\ & \textstyle \quad \times \frac{(R-\gamma) \parens*{4 \Delta^2+(\gamma +\kappa +T_2^{-1} +\xi +R)^2}-4 g^2 (\gamma +\kappa +T_2^{-1} +\xi +R)}{8 g^3 N (R-\gamma) (\gamma +\kappa +T_2^{-1} +R) (\gamma +\kappa +T_2^{-1} +\xi +R)^2} \\ & \textstyle \quad + \frac{\Delta (R-\gamma)}{2 g (\gamma +\kappa +T_2^{-1} +\xi +R)}, \end{aligned} \end{dcases} & \quad \begin{aligned} & \vphantom{\textstyle \frac{\Delta \kappa R (\gamma +T_2^{-1} +R)}{g N (R-\gamma) (\gamma +\kappa +T_2^{-1} +R) (\gamma +\kappa +T_2^{-1} +\xi +R)},} \text{for $R < \gamma$}, \\ & \vphantom{\textstyle -\frac{(\Delta \kappa) }{g (2 \gamma +\kappa +T_2^{-1} +\xi)} \sqrt{\frac{\gamma (2 \gamma +T_2^{-1})}{2 \kappa N (2 \gamma +\kappa +T_2^{-1})}},} \text{for $R = \gamma$}, \\ & \vphantom{\textstyle \Delta \kappa (\gamma +R) (\gamma +T_2^{-1} +R)} \\ & \vphantom{\textstyle \quad \times \frac{(R-\gamma) \parens*{4 \Delta^2+(\gamma +\kappa +T_2^{-1} +\xi +R)^2}-4 g^2 (\gamma +\kappa +T_2^{-1} +\xi +R)}{8 g^3 N (R-\gamma) (\gamma +\kappa +T_2^{-1} +R) (\gamma +\kappa +T_2^{-1} +\xi +R)^2}} \\ & \vphantom{\textstyle \quad + \frac{\Delta (R-\gamma)}{2 g (\gamma +\kappa +T_2^{-1} +\xi +R)},} \text{for $R > \gamma$}, \end{aligned} \\ \mathbb{I}\bracks*{\braket{a\sigma^{+}}^\steady_{ N\rightarrow \infty} } & = \begin{dcases}\label{SpmijN} \textstyle \frac{\kappa R (\gamma +T_2^{-1} +R)}{2 g N (R-\gamma) (\gamma +\kappa +T_2^{-1} +R)}, \\ \textstyle - \frac{1}{2 g}\sqrt{\frac{\gamma \kappa (2 \gamma +T_2^{-1})}{2 N (2 \gamma +\kappa +T_2^{-1})}}, \\ \begin{aligned}[b] & \textstyle - \kappa (\gamma +R) (\gamma +T_2^{-1} +R) \\ & \textstyle \quad \times \frac{\frac{4 g^2 (\gamma +\kappa +T_2^{-1} +\xi +R)}{R-\gamma }-4 \Delta^2-(\gamma +\kappa +T_2^{-1} +\xi +R)^2}{16 g^3 N (\gamma +\kappa +T_2^{-1} +R) (\gamma +\kappa +T_2^{-1} +\xi +R)}-\frac{R-\gamma }{4 g}, \end{aligned} \end{dcases} & \quad \begin{aligned} & \vphantom{\textstyle \frac{\kappa R (\gamma +T_2^{-1} +R)}{2 g N (R-\gamma) (\gamma +\kappa +T_2^{-1} +R)},} \text{for $R < \gamma$}, \\ & \vphantom{\textstyle - \frac{1}{2 g}\sqrt{\frac{\gamma \kappa (2 \gamma +T_2^{-1})}{2 N (2 \gamma +\kappa +T_2^{-1})}},} \text{for $R = \gamma$}, \\ & \vphantom{\textstyle - \kappa (\gamma +R) (\gamma +T_2^{-1} +R)} \\ & \vphantom{\textstyle \quad \times \frac{\frac{4 g^2 (\gamma +\kappa +T_2^{-1} +\xi +R)}{R-\gamma }-4 \Delta^2-(\gamma +\kappa +T_2^{-1} +\xi +R)^2}{16 g^3 N (\gamma +\kappa +T_2^{-1} +R) (\gamma +\kappa +T_2^{-1} +\xi +R)}-\frac{R-\gamma }{4 g},} \text{for $R > \gamma$}, \end{aligned} \\ \braket{\sigma_{i}^{+} \sigma_{j}^{-}}^\steady_{ N\rightarrow \infty} & = \begin{dcases} \textstyle - \frac{\kappa R}{N (\gamma +R) (\gamma +\kappa +T_2^{-1} +R)}, \\ \textstyle - \frac{\kappa }{2 N (2 \gamma +\kappa +T_2^{-1})}, \\ \begin{aligned}[b] & \textstyle \frac{\kappa (R-\gamma) \parens*{4 \Delta^2+(\gamma +\kappa +T_2^{-1} +\xi +R)^2}}{8 g^2 N (\gamma +\kappa +T_2^{-1} +R) (\gamma +\kappa +T_2^{-1} +\xi +R)} \\ & \textstyle \quad -\frac{\kappa (4 g^2) (\gamma +\kappa +T_2^{-1} +\xi +R)}{8 g^2 N (\gamma +\kappa +T_2^{-1} +R) (\gamma +\kappa +T_2^{-1} +\xi +R)}, \end{aligned} \end{dcases} & \quad \begin{aligned} & \vphantom{\textstyle - \frac{\kappa R}{N (\gamma +R) (\gamma +\kappa +T_2^{-1} +R)},} \text{for $R < \gamma$}, \\ & \vphantom{\textstyle - \frac{\kappa }{2 N (2 \gamma +\kappa +T_2^{-1})}, } \text{for $R = \gamma$}, \\ & \vphantom{\textstyle \frac{\kappa (R-\gamma) \parens*{4 \Delta^2+(\gamma +\kappa +T_2^{-1} +\xi +R)^2}}{8 g^2 N (\gamma +\kappa +T_2^{-1} +R) (\gamma +\kappa +T_2^{-1} +\xi +R)}} \\ & \vphantom{\textstyle \quad -\frac{\kappa (4 g^2) (\gamma +\kappa +T_2^{-1} +\xi +R)}{8 g^2 N (\gamma +\kappa +T_2^{-1} +R) (\gamma +\kappa +T_2^{-1} +\xi +R)},} \text{for $R > \gamma$}. \end{aligned} \end{alignat} In the limit of large $N$, the above expressions provide useful predictions for experimental setups which operate in this regime. For example, since the output power of the laser is proportional to $\braket{a^{\dagger} a}^{\steady}$, it scales linearly with the number of atoms $N$ when $R>\gamma$. In the same range of $R>\gamma$, we have that $\braket{\sigma_z}^{\steady}$ is expected to go to zero (corresponds to equally populated ground and excited atomic states), the atom-field correlations $\braket{a\sigma^{+}}^{\steady}$ go to a constant value, and the atom-atom correlations $\braket{\sigma_{i}^{+} \sigma_{j}^{-}}^{\steady}$ go to zero. \endgroup \begin{figure}[!h] \includegraphics[width=0.7\textwidth]{all images/ata_analytical_2.pdf} \caption{Number of intracavity photons at the steady-state $\braket{a^{\dagger}a}^{\steady}$ as a function of atom number $N$ for various values of $R$, for $\gamma = 2\pi \times 7~\unit{mHz} $, $\kappa = 2\pi \times 400~\unit{kHz} $, $T_2^{-1} = 1~\unit{rad.s^{-1}}$, $g = 2\pi\times 3.8~\unit{Hz}$, $\Delta = 0$, and $\xi = 0$. }\label{fig:steady variables N} \end{figure} In Figure~\ref{fig:steady variables N}, we show $\braket{a^{\dagger}a}^{\steady}$ as a function of $N$ for three values of $R$, smaller than, equal to, and larger than $\gamma$. Each case is compared with the corresponding large $N$ limit expressions given in Eqs.~\eqref{ada large N} (for $R<\gamma$ we prefer to consider only the large $N$ constant value $\frac{R (\gamma +T_2^{-1} +R)}{(\gamma -R) (\gamma +\kappa +T_2^{-1} +R)}$), and in the case $R=10^4 \gamma$ we also consider the large $\kappa$ and $R$ expression, given in Eqs.~\eqref{Steady state relations simple large N}. This is valid for $N\ge R \kappa/(4 g^2)\approx 4.85\times 10^5$. For smaller values of $N$, in this limit we find $\braket{a^{\dagger}a}^{\steady}\approx 0$ (Eqs.~\eqref{Steady state relations small N}), which cannot be reported in the log-log scale used in the figure. The results show that for large $N$, $\braket{a^{\dagger}a}^{\steady}$ approaches a constant value for $R<\gamma$, increases as $\sqrt{N}$ for $R=\gamma$, and increases linearly in $N$ for $R>\gamma$. The last case $R>\gamma$, is the one in which the superradiant laser operates and then, for any fixed value of $R$, the power of the laser is expected to become proportional to $N$ for values of this variable large enough to have access to this limit. \begin{figure}[t] \includegraphics[width=0.7\textwidth]{all images/atomic_N1e6_colormix2.pdf} \caption{Atomic variables as a function of $N$ for $R=10^4 \gamma$, for $\gamma = 2\pi \times 7~\unit{mHz} $, $\kappa = 2\pi \times 400~\unit{kHz} $, $T_2^{-1} = 1~\unit{rad.s^{-1}}$, $g = 2\pi\times 3.8~\unit{Hz}$, $\Delta = 0$, and $\xi = 0$. }\label{fig:atomic steady variables N} \end{figure} \newpage In Figure~\ref{fig:atomic steady variables N}, we show for the case $\Delta=\xi=0$, the steady-state variables which involve the atoms, $\braket{\sigma^{z}}^{\steady}$, $\braket{\sigma_{i}^{+} \sigma_{j}^{-}}^{\steady}$, and $\mathbb{I}\bracks*{\braket{a \sigma^{+}}^{\steady}}$ (the real part is zero for $\Delta=0$) as a function of $N$ for $R=10^4 \gamma$. This is compared with both the large $\kappa$ and $R$ expressions given in Eqs.~\eqref{Steady state relations small N} for $N< R\kappa/(4 g^2)\equiv N^* \approx 4.85 \times 10^5$ and in Eqs.~\eqref{Steady state relations simple large N} for larger values of $N$, and the large $N$ limit expressions given in Eqs.~\eqref{ada large N}. We can observe how the large $\kappa$ and $R$ expressions work well for any $N$ while the large $N$ expressions are valid above specific values of $N$. The variable $\braket{\sigma^{z}}^{\steady}$ here starts from 1 for $N< N^*$ and decreases to zero for larger values of $N$. Both $\braket{\sigma_{i}^{+} \sigma_{j}^{-}}^{\steady}$ and $\mathbb{I}\bracks*{\braket{a \sigma^{+}}^{\steady}}$ start from zero for $N< N^*$. For larger values of $N$, $\braket{\sigma_{i}^{+} \sigma_{j}^{-}}^{\steady}$ increases up to a maximum of $1/8$ for $N =2 N^*$ and then decreases back to zero. For this specific value of $N$, we have $\braket{\sigma^{z}}^{\steady}\approx 1/2$, $\braket{\sigma_{i}^{+} \sigma_{j}^{-}}^{\steady}\approx 1/8$, and $\mathbb{I}\bracks*{\braket{a \sigma^{+}}^{\steady}}\approx -R/(8 g) \approx -2.30$. For values of $N>N^*$, $\mathbb{I}\bracks*{\braket{a \sigma^{+}}^{\steady}}$ decreases towards a constant minimum value approximately equal to $-R/(4 g) \approx -4.61$. This analysis is based on the approximated expressions that are valid in the large $\kappa$ and $R$ limit. \section{Cavity field spectrum}\label{app:spectrum} In the case when all the atoms have the same parameters, the $N+1$ equations in Eqs.~\eqref{eq:time_evo} reduce to two equations: \begin{equation}\label{app:2eqs:matrix} \frac{\dif}{\dif \tau} \begin{pmatrix} \braket[\big]{a^\dagger (\tau+t) a(t)} \\ \braket[\big]{\sigma^\dagger(\tau+t) a(t)} \end{pmatrix} = \begin{pmatrix} c_1 & \mathrm{i}c_2 \\ \mathrm{i}c_3 & c_4 \end{pmatrix} \begin{pmatrix} \braket[\big]{a^\dagger (\tau+t) a(t)} \\ \braket[\big]{\sigma^\dagger(\tau+t) a(t)} \end{pmatrix}, \end{equation} where the coefficients are defined as $c_1 = -(\kappa+\xi)/2$, $c_2 = gN$, $c_3 = -g\braket{\sigma^z}^\steady$, and $c_4 = c_{4\mathbb{R}}+\mathrm{i} c_{4\mathbb{I}} = -(\gamma + R + 1/T_2)/2 - \mathrm{i}\Delta$. The system can be solved by applying the Laplace transform, using the standard operational rule $\mathcal{L}\bracks[\big]{\frac{\dif}{\dif \tau} f(\tau)}=s F(s)- f(0)$, with the initial conditions $\braket[\big]{a^{\dagger}(0) a(0)} = \braket{a^{\dagger}a}^{\steady}$ and $\braket[\big]{\sigma^{+}(0) a(0)} = \braket{\sigma^{+}a}^{\steady}$. Solving the resulting set of equations gives \begin{equation}\label{Spectrum1} \braket[\big]{a^{\dagger}(s) a(0)} = \frac{(s - c_4)\braket{a^{\dagger} a}^{\steady} + c_2\braket{\sigma^{+}a}^{\steady}}{(s-c_1)(s-c_4) - c_2c_3}. \end{equation} Using Eqs.~\eqref{Steady state relations}, the above expression leads to \begin{equation}\label{app:ata:spectrum:N} \braket[\big]{a^{\dagger}(\omega) a(0)} =\braket{a^{\dagger} a}^{\steady}\frac{\frac{1}{2} (\gamma +\kappa +T_2^{-1} +R)+i \parens[\Big]{\frac{\Delta (\gamma +T_2^{-1} +\xi +R)}{\gamma +\kappa +T_2^{-1} +\xi +R}+\omega } }{\parens[\big]{\frac{\kappa +\xi }{2}+i \omega } \bracks*{\frac{1}{2} (\gamma +T_2^{-1} +R)+i (\Delta +\omega) }-g^2 N \braket{\sigma^z}^{\steady}}, \end{equation} where we used $s=\mathrm{i} \omega$ to correctly compute $S(\omega)$, which is given by $2 \mathbb{R}\braces[\big]{ \braket[\big]{a^{\dagger}(\omega) a(0)} }$. The above expression shows that in the large $N$ limit, where $\braket{a^{\dagger} a}^{\steady}$ increases linearly with $N$, as shown in Eqs.~\eqref{ada large N}, and $\braket{\sigma^z}^{\steady}$ scales inversely proportional to $N$, as shown in Eqs.~\eqref{sigmaz large N}, the spectrum takes the form of a function of $\omega$ multiplied by $N$. It follows that in this limit the width of the spectrum becomes independent on $N$. This behavior is observed in Figure~\ref{fig:power_and_linewidth_1class}(b) for large values of $N$ and for values of $R$ larger than $\gamma$ up to some value (for $R$ very large we can expect that the limit of large $N$ occurs for values of $N$ larger than the ones considered in this figure). The independence of the linewidth from $N$ is reflected by the vertical behavior of the white contour lines at constant values of $R$. We note that the linewidth stays independent of $N$ also in the case considered in Appendix~\ref{app:simple} for $N\ge R \kappa/(4 g^2)$, since $\braket{\sigma^z}^{\steady}$ is inversely proportional to $N$, as shown in Eqs.~\eqref{Steady state relations simple large N}. From Eq.~\eqref{app:ata:spectrum:N}, the spectrum can be put under the final form \begin{equation}\label{eq:Spectrum} \begin{split} S(\omega) & = \braket{a^{\dagger} a}^{\steady}\left[(\gamma +\kappa +T_2^{-1} +R)\left(\frac{(\kappa +\xi) \parens*{4 \Delta ^2+(\gamma +T_2^{-1} +R) (\gamma +\kappa +T_2^{-1} +\xi +R)}}{4 (\gamma +\kappa +T_2^{-1} +\xi +R)} \right.\right. \\ & \quad \left.\left. - g^2 N \braket{\sigma^z}^{\steady} \vphantom{\frac{(\kappa +\xi) \parens*{4 \Delta ^2+(\gamma +T_2^{-1} +R) (\gamma +\kappa +T_2^{-1} +\xi +R)}}{4 (\gamma +\kappa +T_2^{-1} +\xi +R)}} \right) + 2 \Delta \xi \omega +\xi \omega^2\right]\Bigg/\bracks[\Bigg]{\frac{(\kappa +\xi)^2}{16} \parens*{4 \Delta ^2+(\gamma +T_2^{-1} +R)^2} \\ & \quad -\frac{1}{2} g^2 (\kappa +\xi) (\gamma +T_2^{-1} +R)N \braket{\sigma^z}^{\steady} + g^4 N^2 \braket{\sigma^z}^{\steady~2} +\Delta \parens*{\frac{(\kappa +\xi)^2}{2} +2 g^2 N \braket{\sigma^z}^{\steady}}\omega \\ & \quad + \parens*{\frac{1}{4} \parens*{4 \Delta ^2+(\kappa +\xi)^2+(\gamma +T_2^{-1} +R)^2}+2 g^2 N \braket{\sigma^z}^{\steady}} \omega ^2 +2 \Delta \omega ^3 +\omega ^4}. \end{split} \end{equation} An alternative way to manipulate Eq.~\eqref{Spectrum1} consists in transforming it into \begin{equation} \braket[\big]{a^{\dagger}(s) a(0)} = \frac{(s - c_4)\braket{a^{\dagger} a}^{\steady} + c_2\braket{\sigma^{+}a}^{\steady}}{(s - \lambda_+)(s - \lambda_-)}, \end{equation} where \begin{equation}\label{app:lambda_pm} \lambda_{\pm}=\frac{1}{2} \parens*{(c_1+c_4) \pm \sqrt{\parens*{c_1+c_4 }^2-4\parens*{c_1 c_4 -c_2 c_3}}} = a_{\pm} + \mathrm{i} b_{\pm}. \end{equation} By using the method of partial fractions, $\braket[\big]{a^{\dagger}(s) a(0)}$ can be then expressed as \begin{equation}\label{Partial fractions} \braket[\big]{a^{\dagger}(s) a(0)} = \frac{X}{(s - \lambda_+)} + \frac{Y}{(s - \lambda_-)}, \end{equation} where \begin{equation} X = \frac{c_4 \braket{a^{\dagger} a}^{\steady} - c_2 \braket{\sigma_A^{+} a}^{\steady} - \lambda_{+} \braket{a^{\dagger} a}^{\steady}}{\lambda_{-} - \lambda_{+}}, \quad Y = \frac{-c_4 \braket{a^{\dagger} a}^{\steady} + c_2 \braket{\sigma_A^{+} a}^{\steady} + \lambda_{-} \braket{a^{\dagger} a}^{\steady}}{\lambda_{-} - \lambda_{+}}. \end{equation} From the above equation we can find that $X + Y = \braket{a^{\dagger} a}^{\steady}$, with $\braket{a^{\dagger} a}^{\steady}$ a real value. It then follows that $X = X_{\mathbb{R}} + \mathrm{i} X_{\mathbb{I}}$ and $Y = Y_{\mathbb{R}} - \mathrm{i} X_{\mathbb{I}}$. The inverse Laplace transform of Eq.~\eqref{Partial fractions} then gives \begin{equation}\label{inverseLT} \braket[\big]{a^{\dagger}(\tau) a(0)} = Xe^{\lambda_+\tau} + Ye^{\lambda_-\tau}. \end{equation} The structure of $\lambda_{\pm}$ determines that the second term in the above Eq.~\eqref{inverseLT} decays to zero faster than the first~\cite{PhD_DavidTieri, PhD_SwadheenDubey}. To obtain the final expression for the spectrum we have to compute $S(\omega) \equiv 2 \mathbb{R}\braces[\big]{ \braket[\big]{a^{\dagger}(\omega) a(0)} }$. From Eq.~\eqref{Partial fractions}, we have \begin{equation} \begin{split} \braket[\big]{a^{\dagger}(\omega) a(0)} & = \frac{X_{\mathbb{R}} + \mathrm{i} X_{\mathbb{I}}}{-a_+ +i(\omega - b_+)} + \frac{Y_{\mathbb{R}} - \mathrm{i} X_{\mathbb{I}}}{-a_- +i(\omega - b_-)} \\ & \begin{multlined}[t][.7\displaywidth] = \frac{-a_+X_{\mathbb{R}}+X_{\mathbb{I}}(\omega-b_+) - \mathrm{i}a_+X_{\mathbb{I}} -\mathrm{i} X_{\mathbb{R}}(\omega-b_+)}{a_+^2 +(\omega-b_+)^2} \\ + \frac{-a_-Y_{\mathbb{R}}-X_{\mathbb{I}}(\omega-b_-) + \mathrm{i}a_-X_{\mathbb{I}} -\mathrm{i} Y_{\mathbb{R}}(\omega-b_-)}{a_-^2 +(\omega-b_-)^2}, \end{multlined} \end{split} \end{equation} leading to \begin{equation} S(\omega) = 2\frac{-a_+X_\mathbb{R} + X_{\mathbb{I}}(\omega - b_+)}{a_+^2 + (\omega - b_+)^2} + 2\frac{-a_- Y_\mathbb{R} - X_{\mathbb{I}}(\omega - b_-)}{a_-^2 + (\omega - b_-)^2}. \end{equation} For $X_{\mathbb{I}} =0$, we have two Lorentzian functions centered at $b_+$ and $b_-$. In the general case $X_{\mathbb{I}} \neq 0$, each Lorentzian is modified by a linear function at the numerator. \subsection*{Cavity field spectrum for \texorpdfstring{$\Delta=\xi=0$}{Delta=xi = 0}} In the case $\Delta=\xi=0$, Eq.~\eqref{eq:Spectrum} simplifies to \begin{equation}\label{app:Spectrum special case} \begin{split} S(\omega) & = \braket{a^{\dagger} a}^{\steady} \bracks[\Bigg]{(\gamma +\kappa +T_2^{-1} +R)\parens[\Bigg]{\frac{\kappa (\gamma +T_2^{-1} +R) (\gamma +\kappa +T_2^{-1} +R)}{4 (\gamma +\kappa +T_2^{-1} +R)}-g^2 N \braket{\sigma^z}^{\steady} }}\Bigg/\bracks[\Bigg]{\frac{\kappa^2}{16} \\ & \quad \times (\gamma +T_2^{-1} +R)^2 -\frac{1}{2} g^2 \kappa (\gamma +T_2^{-1} +R)N \braket{\sigma^z}^{\steady}+g^4 N^2 \braket{\sigma^z}^{\steady~2} \\ & \quad + \parens*{\frac{1}{4} \parens*{\kappa^2+(\gamma +T_2^{-1} +R)^2}+2 g^2 N \braket{\sigma^z}^{\steady}}\omega ^2+\omega ^4}, \end{split} \end{equation} where $\braket{a^{\dagger} a}^{\steady}$ and $\braket{\sigma^z}^{\steady}$ are the steady-state values for $\Delta=\xi=0$. The dependence of the above spectrum on $\omega$ has the form \begin{equation}\label{app:Spectrum special case form} S(\omega)=\frac{X^\prime}{Y^\prime+Z\omega ^2+\omega ^4}, \end{equation} where the form of $X^\prime$, $Y^\prime$, and $Z$ can be obtained from Eq.~\eqref{app:Spectrum special case}. The spectrum is then symmetric with respect to $\omega=0$ and for $Z>0$, which is valid for $\braket{\sigma^z}^{\steady}>0$, condition satisfied in the regime treated in Appendix~\ref{app:simple}, $S(\omega)$ has a maximum equal to $X^\prime/Y^\prime$ for $\omega=0$. The central frequency of the spectrum is then equal to zero in the rotating frame of the cavity meaning that it is equal to the cavity frequency. The linewidth $\Delta\nu$ which corresponds to the FWHM of the spectrum, is determined by the equation $S(\omega)=X^\prime/(2Y^\prime)$ leading to \begin{equation}\label{app:FWHM special case form} \Delta\nu=\sqrt{2} \sqrt{\sqrt{4 Y^\prime+Z^2}-Z}. \end{equation} We then find that a possible approach to reduce $\Delta\nu$ is to operate where $Y^\prime\approx 0$ or $Y^\prime\ll Z$. In the large $\kappa$ and $R$ limit in Appendix~\ref{app:simple}, $Y^\prime$ takes the form $\frac{\kappa^2R^2 }{16} -\frac{1}{2} g^2 \kappa R N \braket{\sigma^z}^{\steady}+g^4 N^2 \braket{\sigma^z}^{\steady~2} $. For $N\le R \kappa/(4 g^2)$, we have for the atomic variable $\braket{\sigma^z}^{\steady}\approx 1$, which gives $Y^\prime\approx 0$ at $N= R \kappa/(4 g^2)$, and for $N\ge R \kappa/(4 g^2)$, we have $\braket{\sigma^z}^{\steady}\approx R \kappa/(4 g^ 2 N)$ and $Y^\prime\approx 0$. In the large $N$ limit of Appendix~\ref{app:largeN}, with $R>\gamma$, we consider the large $\kappa$ and $R$ limit, such that the large $N$ approximation stays valid. For the case $\kappa \gg R$, we find that $Y^\prime=g^4$ and $Z=\kappa^2/4$. This implies that $Y^\prime\ll Z^2$, which allows us to develop Eq.~\eqref{app:FWHM special case form} to obtain \begin{equation}\label{app:FWHM special case form large N} \Delta\nu\approx 2 \sqrt{\frac{Y^\prime}{Z}}=\frac{4 g^2 }{\kappa}. \end{equation} This is equal to the one reported in~\cite{Meiser2009, Kazakov2022}, where the minimal linewidth of a superradiant laser is estimated by $\frac{4 g^2}{\kappa}$. \printbibliography \end{document}